Hutton - Fundamentals of Finite Element Analysis (523155), страница 64
Текст из файла (страница 64)
Hence, 36 integrationsare required to obtain the complete stiffness matrix. The integrations are straightforward but algebraic tedious. Here, we develop only a single term of the stiffness matrix in detail, then discuss the more-efficient numerical methods used infinite element software packages.If we carry out the matrix multiplications just indicated, the first diagonalterm of the stiffness matrix is found (after a bit of algebra) to be(e)k 11Etb=16a(1 + 2)1 1Eta(s − 1) dr ds +32b(1 + )1 1(r − 1) 2 dr ds2−1 −1−1 −1(9.65)and this term evaluates to(e)k 1111Etb2(s − 1) 3 Eta2(r − 1) 3 = + 32b(1 + )16a(1 + 2)33−1−1 Etb16Eta16=+16a(1 + 2) 332b(1 + ) 3(9.66)Hutton: Fundamentals ofFinite Element Analysis9.
Applications in SolidMechanicsText© The McGraw−HillCompanies, 20049.4 Isoparametric Formulation of the Plane Quadrilateral ElementNote that the integrands are quadratic functions of the natural coordinates. Infact, analysis of Equation 9.64 reveals that every term of the element stiffnessmatrix requires integration of quadratic functions of the natural coordinates.From the earlier discussion of Gaussian integration (Chapter 6), we know that aquadratic polynomial can be integrated exactly using only two integration (orevaluation) points. As here we deal with integration in two dimensions, we mustevaluate the integrand at the Gauss points√√33ri = ±sj = ±33with weighting factors W i = W j = 1 .
If we apply the numerical integration(e), we obtain, as expected, the result identical to thattechnique to evaluation of k 11given by Equation 9.66. More important, the Gauss integration procedure can beapplied directly to Equation 9.64 to obtain the entire element stiffness matrix ask (e) = t ab22W i W j [B(r i , s j )] T [D][B(r i , s j )](9.67)i=1 j=1where the matrix triple product is evaluated four times, in accordance with thenumber of integration points required.
The summations and matrix multiplications required in Equation 9.67 are easily programmed and ideally suited todigital computer implementation.While written specifically for the four-node rectangular element, Equation 9.67 is applicable to higher-order elements as well. Recall that, as the polynomial order increases, exact integration via Gaussian quadrature requires increasein both number and change in value of the integration points and weighting factors. By providing a “look-up” table of values fashioned after Table 6.1, computerimplementation of Equation 9.67 can be readily adapted to higher-order elements.We use the triangular element to illustrate plane stress and the rectangularelement to illustrate plane strain.
If the developments are followed clearly, it isapparent that either element can be used for either state of stress. The only difference is in the stress-strain relations exhibited by the [D] matrix. This situationis true of any element shape and order (in terms of number of nodes and order ofpolynomial interpolation functions). Our use of the examples of triangular andrectangular elements are not meant to be restrictive in any way.9.4 ISOPARAMETRIC FORMULATION OFTHE PLANE QUADRILATERAL ELEMENTWhile useful for analysis of plane problems in solid mechanics, the triangularand rectangular elements just discussed exhibit shortcomings. Geometrically, thetriangular element is quite useful in modeling irregular shapes having curvedboundaries.
However, since element strains are constant, a large number of smallelements are required to obtain reasonable accuracy, particular in areas of high347Hutton: Fundamentals ofFinite Element Analysis3489. Applications in SolidMechanicsCHAPTER 9Text© The McGraw−HillCompanies, 2004Applications in Solid Mechanicssv34yu31(1, 1)3u4rv2v1x(⫺1, 1)43v4u121(⫺1, ⫺1)u2(a)2(1, ⫺1)(b)Figure 9.6(a) A four-node, two-dimensional isoparametric element. (b) The parent element innatural coordinates.stress gradients, such as near geometric discontinuities. In comparison, the rectangular element provides the more-reasonable linear variation of strain components but is not amenable to irregular shapes. An element having the desirablecharacteristic of strain variation in the element as well as the ability to closely approximate curves is the four-node quadrilateral element.
We now develop thequadrilateral element using an isoparametric formulation adaptable to eitherplane stress or plane strain.A general quadrilateral element is shown in Figure 9.6a, having elementnode numbers and nodal displacements as indicated. The coordinates of node iare (xi, yi) and refer to a global coordinate system. The element is formed bymapping the parent element shown in Figure 9.6b, using the procedures developed in Section 6.8. Recalling that, in the isoparametric approach, the geometricmapping functions are identical to the interpolation functions used to discretizethe displacements, the geometric mapping is defined by4x =N i (r, s)x ii=1(9.68)4y=N i (r, s) yii=1and the interpolation functions are as given in Equation 9.60, so that the displacements are described as4u(x , y) =N i (r, s)u ii=1(9.69)4v(x , y) =N i (r, s)vii=1Hutton: Fundamentals ofFinite Element Analysis9.
Applications in SolidMechanicsText© The McGraw−HillCompanies, 20049.4 Isoparametric Formulation of the Plane Quadrilateral ElementNow, the mathematical complications arise in computing the strain componentsas given by Equation 9.55 and rewritten here as ∂∂u0 ∂x ∂x εx ∂ ∂v 0 u={ε} = ε y=(9.70)∂y v ∂y␥x y ∂∂ ∂u + ∂v ∂y∂x∂y ∂xUsing Equation 6.83 with = u, we have∂u∂u ∂x∂u ∂y=+∂r∂ x ∂r∂ y ∂r(9.71)∂u ∂x∂u ∂y∂u=+∂s∂x ∂s∂y ∂swith similar expressions for the partial derivative of the v displacement. WritingEquation 9.71 in matrix form ∂u ∂u ∂x ∂y ∂u ∂r = ∂r ∂r ∂ x = [J ] ∂ x(9.72) ∂x ∂y ∂u ∂u ∂u ∂y∂y∂s∂s ∂sand the Jacobian matrix is identified as[J ] =J11J21J12J22∂x ∂r= ∂x∂s∂y∂r ∂y ∂s(9.73)as in Equation 6.83. Note that, per the geometric mapping of Equation 9.68, thecomponents of [ J ] are known as functions of the partial derivatives of the interpolation functions and the nodal coordinates in the x y plane.
For example,J11 =∂x=∂r4i=1∂ Ni1x i = [(s − 1)x 1 + (1 − s)x 2 + (1 + s)x 3 − (1 + s)x 4 ]∂r4(9.74)a first-order polynomial in the natural (mapping) coordinate s. The other termsare similarly first-order polynomials.Formally, Equation 9.72 can be solved for the partial derivatives of displacement component u with respect to x and y by multiplying by the inverse ofthe Jacobian matrix. As noted in Chapter 6, finding the inverse of the Jacobianmatrix in algebraic form is not an enviable task. Instead, numerical methods areused, again based on Gaussian quadrature, and the remainder of the derivationhere is toward that end.
Rather than invert the Jacobian matrix, Equation 9.72349Hutton: Fundamentals ofFinite Element Analysis3509. Applications in SolidMechanicsCHAPTER 9Text© The McGraw−HillCompanies, 2004Applications in Solid Mechanicscan be solved via Cramer’s rule. Application of Cramer’s rule results in ∂u ∂r J12 ∂u∂u J22 ∂u1∂s∂r==[J22 −J12 ]∂u |J ||J |∂x∂s J11 ∂u ∂r ∂u J ∂u 21∂u1∂s==[−J21 +J11 ] ∂r|J ||J |∂y ∂u ∂sor, in a more compact form,∂u ∂x ∂u ∂y=1J22−J|J |21∂u∂r−J12∂uJ11 ∂s(9.75)(9.76)The determinant of the Jacobian matrix | J | is commonly called simply theJacobian.Since the interpolation functions are the same for both displacement components, an identical procedure results in∂v ∂v ∂x 1∂rJ22 −J12=(9.77)∂v ∂v J11 |J | −J21∂y∂sfor the partial derivatives of the v displacement component with respect to globalcoordinates.Let us return to the problem of computing the strain components perEquation 9.70.
Utilizing Equations 9.76 and 9.77, the strain components areexpressed as∂u ∂u ∂u∂r ∂r ∂x∂u∂u εxJ−J00∂v121 22∂s∂s{ε} = ε y ===[G](9.78)00−JJ2111∂y |J | ∂v ∂v ␥x y−JJJ−J21112212∂u∂v ∂r ∂r +∂v∂v∂y∂x∂s∂sHutton: Fundamentals ofFinite Element Analysis9.
Applications in SolidMechanicsText© The McGraw−HillCompanies, 20049.4 Isoparametric Formulation of the Plane Quadrilateral Elementwith what we will call the geometric mapping matrix, defined asJ−J0022121 [G] =(9.79)J11 00−J21|J |−J21J11J22 −J12We must expand the column matrix on the extreme right-hand side of Equation 9.78 in terms of the discretized approximation to the displacements. ViaEquation 9.69, we have ∂u ∂ N1 ∂ N2 ∂ N3 ∂ N4 u1 0000 ∂ru ∂r∂r∂r∂r2∂u ∂ N1 ∂ N2 ∂ N3 ∂ N4u3 0000 ∂s = u4∂s∂s∂s ∂s∂ N1 ∂ N2 ∂ N3 ∂ N4 ∂v v1 0000v ∂r∂r∂r∂r2∂r ∂ N1 ∂ N2 ∂ N3 ∂ N4 v3 ∂v 0000v4∂s∂s∂s∂s∂s(9.80)where we reemphasize that the indicated partial derivatives are known functionsof the natural coordinates of the parent element.
For shorthand notation, Equation 9.80 is rewritten as∂u ∂r∂u∂s =[P] {␦}(9.81)∂v ∂r ∂v ∂sin which [ P ] is the matrix of partial derivatives and {␦} is the column matrix ofnodal displacement components.Combining Equations 9.78 and 9.81, we obtain the sought-after relation forthe strain components in terms of nodal displacement components as{ε} = [G ][ P ]{␦}(9.82)and, by analogy with previous developments, matrix [B] = [G ][ P ] has beendetermined such that{ε} = [B]{␦}(9.83)and the element stiffness matrix is defined by (e) k= t [B] T [D][ B] d AA(9.84)351Hutton: Fundamentals ofFinite Element Analysis3529. Applications in SolidMechanicsCHAPTER 9Text© The McGraw−HillCompanies, 2004Applications in Solid Mechanicswith t representing the constant element thickness, and the integration is performed over the area of the element (in the physical xy plane).
In Equation 9.84,the stiffness may represent a plane stress element or a plane strain element, depending on whether the material property matrix [D] is defined by Equation 9.6or 9.54, respectively. (Also note that, for plane strain, it is customary to take theelement thickness as unity.)The integration indicated by Equation 9.84 are in the x-y global space, butthe [B] matrix is defined in terms of the natural coordinates in the parent elementspace. Therefore, a bit more analysis is required to obtain a final form.
In thephysical space, we have d A = dx dy , but we wish to integrate using the naturalcoordinates over their respective ranges of −1 to +1. In the case of the four-noderectangular element, the conversion is straightforward, as x is related only to rand y is related only to s, as indicated in Equation 9.61.
In the isoparametric caseat hand, the situation is not quite so simple. The derivation is not repeated here,but it is shown in many calculus texts [1] thatd A = dx dy = | J | dr dshence, Equation 9.84 becomes1 1 (e) T[B] T [D][ B] | J | dr dsk= t [B] [D][ B] | J | dr ds = t(9.85)(9.86)−1 −1AAs noted, the terms of the [B] matrix are known functions of the naturalcoordinates, as is the Jacobian | J |. The terms in the stiffness matrix representedby Equation 9.86, in fact, are integrals of ratios of polynomials and the integrations are very difficult, usually impossible, to perform exactly.
Instead, Gaussianquadrature is used and the integrations are replaced with sums of the integrandevaluated at specified Gauss points as defined in Chapter 6. For p integrationpoints in the variable r and q integration points in the variable s, the stiffnessmatrix is approximated byk (e) = tpqW i W J [B(r i , s j )] T [D][ B(r i , s j )]| J (r i , s j )|dr ds(9.87)i=1 j=1Since [B] includes the determinant of the Jacobian matrix in the denominator, thenumerical integration does not necessarily result in an exact solution, since theratio of polynomials is not necessarily a polynomial. Nevertheless, the Gaussianprocedure is used for this element, as if the integrand is a quadratic in both r ands, with good results. In such case, we use two Gauss points for each variable, asis illustrated in the following example.EXAMPLE 9.3Evaluate the stiffness matrix for the isoparametric quadrilateral element shown in Figure 9.7 for plane stress with E = 30(10) 6 psi, = 0.3, t = 1 in. Note that the propertiesare those of steel.Hutton: Fundamentals ofFinite Element Analysis9.