Belytschko T. - Introduction (779635), страница 42
Текст из файла (страница 42)
Belytschko, Lagrangian Meshes, December 16, 1998NI ( ξ ) =1(1+ ξ I ξ )(1 + η I η )4(E4.2.2)where (ξ I ,η I ) , I= 1 to 4, are the nodal coordinates of the parent element shown in Fig. 4.4.They are given byξ −1 1 1 −1−1 1 1 [ξiI ] = ηI = −1 I(E4.2.3)Since (E4.2.1) also holds for t=0, we can write X (ξ) XI = N I (ξ ) Y(ξ ) YI (E4.2.4)where X I ,YI are the coordinates in the undeformed configuration. The nodal velocities are givenbyv x (ξ,t ) v xI (t ) = N I (ξ )vy (ξ ,t ) vyI ( t) (E4.2.5)which is the material time derivative of the expression for the displacement.Rate-of-Deformation and Internal Nodal Forces. The map (E4.2.1) is not invertible for theshape functions given by (E4.2.2).
Therefore it is impossible to write explicit expressions for theelement coordinates in terms of x and y, and the derivatives of the shape functions are evaluated byusing implicit differentiation. Referring to (4.4.47) we have[NIT, x = NI , x[]NI , y = N TI ,ξ x,−1ξ = NI , ξξ, x ξ, y NI ,η η, x η, y ](E4.2.6)The Jacobian of the current configuration with respect to the element coordinates is given by x,ξx,ξ = y,ξx, η x I = [ x iI ] ∂NI ∂ξ j = N I ,ξy, η y I [[] x I N I ,ξNI , η = y I N I ,ξ]x I NI ,η yI N I ,η (E4.2.7a)For the 4-node quadrilateral the above is x ( t)ξI (1+ η Iη) x I (t )η I (1 +ξI ξ)x,ξ = ∑ II =1 y I ( t) ξI (1+η Iη) yI (t )η I (1+ ξIξ ) 4(E4.2.7b)In the above, the summation has been indicated explicitly because the index I appears three times.As can be seen from the RHS, the Jacobian matrix is a function of time. The inverse of Fξ isgiven by4-39T.
Belytschko, Lagrangian Meshes, December 16, 1998x ,−1ξ =1Jξ y,η −y,ξ− x,η , J = x, ξ y,η − x,η y,ξx, ξ ξ(E4.2.7c)The gradients of the shape functions for the 4-node quadrilateral with respect to theelement coordinates are given by∂N1∂N2TN, ξ = [∂NI ∂ξ i ] =∂N3 ∂N4∂ξ∂ξ∂ξ∂ξξ1( 1+η1η) η1(1+ ξ1ξ) ∂η∂η 1 ξ2 ( 1+η2η) η 2 (1+ ξ2 ξ )=∂η 4 ξ3 (1+ η3η) η3(1+ ξ3ξ ) ∂ηξ4 ( 1+η4η ) η 4 (1+ ξ4ξ )∂N1∂N2∂N3∂N4The gradients of the shape functions with respect to the spatial coordinates can then be computedbyB I = NIT, x ξ1(1+ η1η) η1 (1+ ξ1ξ) y,ηT −1 ξ2 ( 1+η2η ) η 2(1+ ξ2ξ ) 1= NI , ξ x ,ξ =ξ3 ( 1+η3η) η3 (1+ ξ3ξ ) Jξ −x,ηξ4 (1+ η4η ) η 4(1+ ξ4ξ )−y, ξ x,ξ (E4.2.8a)and the velocity gradient is given by Eq.
(4.5.3)L = v IB IT = v I NIT, x(E4.2.8b)For a 4-node quadrilateral which is not rectangular, the velocity gradient, and hence the rate-ofdeformation, is a rational function because Jξ = det x, ξ appears in the denominator of x , ξ and( )hence in L. The determinant Jξ is a linear function in (ξ, η) .The nodal internal forces are obtained by (4.5.6), which gives( ) [f intIT= f xIf yI]int= ∫ B TI σdΩ =Ω∫[N I, xΩσ xx σ xyN I, y dΩσ xy σ yy ](E4.2.9)The integration is performed over the parent domain.
For this purpose, we usedΩ = Jξ adξdη(E4.2.10)where a is the thickness. The internal forces are then given by (4.4.11), which when written outfor two dimensions gives:(f intI ) = [ f xITf yI]int[= ∫ NI , x∆σ xx σ xy NI , y aJξ d∆σ xy σ yy ]4-40(E4.2.11)T. Belytschko, Lagrangian Meshes, December 16, 1998where NI ,i is given in Eq. (E4.2.8a). Equation (E4.2.14) applies to any isoparametric element intwo dimensions. The integrand is a rational function of the element coordinates, since Jξ appearsin the denominator (see Eq.
(4.2.8a)), so analytic quadrature of the above is not feasible.Therefore numerical quadrature is generally used. For the 4-node quadrilateral, 2x2 Gaussquadrature is full quadrature. However, for full quadrature, as discussed in Chapter 8, the elementlocks for incompressible and nearly incompressible materials in plane strain problems. Therefore,selective-reduced quadrature as described in Section 4.5.4, in which the volumetric stress isunderintegrated, must be used for the four-node quadrilateral for plane strain problems when thematerial response is nearly incompressible, as in elastic-plastic materials.The displacement for a 4-node quadrilateral is linear along each edge.
Therefore, theexternal nodal forces are identical to those for the 3-node triangle, see Eqs. (E4.1.29-E4.1.33).Mass Matrix. The consistent mass matrix is obtained by using (4.4.52), which gives N1 N 2˜M=[N N3 1Ω0N 4∫N2N3N4 ]ρ0 dΩ0(E4.2.12)We usedΩ 0 = Jξ0 (ξ ,η )a0dξdη(E4.2.13)where Jξ0 (ξ,η ) is the determinant of the Jacobian of the transformation of the parent element to the˜ wheninitial configuration a is the thickness of the undeformed element.
The expression for M0evaluated in the parent domain is given byN12N1 N2+1+ 1N22˜ =M∫ ∫ symmetric−1 −1N1 N3N2 N3N32N1 N4 N2 N4 ρ0a0 J ξ0 (ξ,η) dξdηN 3N4N 24 (E4.2.14)The matrix is evaluated by numerical quadrature. This mass matrix can be expanded to an 8x8matrix using the same procedure described for the triangle in the previous example.A lumped, diagonal mass matrix can be obtained by using Lobatto quadrature with thequadrature points coincident with the nodes. If we denote the integrand of Eq.
(E4.2.14) bym(ξI ,η I ) , then Lobatto quadrature gives4M=∑ m(ξI ,η I )(E4.2.15)I =14-41T. Belytschko, Lagrangian Meshes, December 16, 1998Alternatively, the lumped mass matrix can be obtained by apportioning the total mass of theelement equally among the four nodes. The total mass is ρ 0 A0 a0 when a0 is constant, so dividingit among the four nodes givesM=1ρ Aa I4 0 0 0 4(E4.2.16)where I 4 is the unit matrix of order 4.Example 4.3. Three Dimensional Isoparametric Element. Develop the expressions forthe rate-of-deformation, the nodal forces and the mass matrix for three dimensional isoparametricelements.
An example of this class of elements, the eight-node hexahedron, is shown in Fig. 4.5.8z5ζ7φ (x(ξ,0), t)564638117242yηξ3xFig. 4.5. Parent element and current configuration for an 8-node hexahedral element.Motion and Strain Measures. The motion of the element is given byxx ( t) I y = N I (ξ )y I (t ) z z (t ) I ξ = (ξ , η, ζ )(E4.3.1)where the shape functions for particular elements are given in Appendix C.
Equation (E4.3.1) alsoholds at time t=0, so4-42T. Belytschko, Lagrangian Meshes, December 16, 1998 XX I Y = N I (ξ ) YI Z Z I(E4.3.2)The velocity field is given byv v x xI vy = NI (ξ)v yI v v z zI (E4.3.3)The velocity gradient is obtained from Eq. (4.5.3), giving[BIT = N I ,xNI ,yN I ,zv xI L = v IB = v yI N I ,xv zI [TIvxI NI , x= vyI NI , x vzI N I , x](E4.3.4)NI ,yNI ,z](E4.3.5)v xI N I , y v xI N I , z vyI N I , y v yI N I , z v zI NI , y vzI N I ,z (E4.3.6)The derivatives with respect to spatial coordinates are obtained in terms of derivatives with respectto the element coordinates by Eq.
(4.4.37).NIT, x = N IT,ξ x −1,ξx,ξ ≡ Fξ x,ξ =x I N IT,ξ(E4.3.7)x I= yI N I ,ξz I[N I, ηN I ,ζ](E4.3.8)The deformation gradient can be computed by Eqs. (3.2.10), (E4.3.1) and (E4.3.7):F=( )∂xT0= x I N I, X = x I N TI ,ξX ,−1ξ ≡ x I NI , ξ Fξ∂X−1(E4.3.9)whereX, ξ ≡ Fξ0 = XI NIT, ξ(E4.3.10)The Green strain is then computed by Eq. (3.3.5); a more accurate procedure is described inExample 4.12.4-43T. Belytschko, Lagrangian Meshes, December 16, 1998Internal Nodal Forces.
The internal nodal forces are obtained by Eq. (4.5.6):(f ) = [ fint TIxI, f yI , fzI] =∫BintTIΩ[σdΩ = ∫ N I, x∆NI, yNI , zσ xx σ xy σ xz σ xy σ yy σ yz Jξ d∆σ xz σ yz σ zz ](E4.3.11)The integral is evaluated by numerical quadrature, using the quadrature formula (4.5.26).External Nodal Forces. We consider first the nodal forces due to the body force.(4.4.13), we havefiIext = N I ρbidΩ = N I ( ξ)ρ(ξ)bi ( ξ)J ξd∆∫∫ΩBy Eq.(E4.3.12a)∆b (ξ) f xI ext 1 1 1 x fyI = ∫ ∫ ∫ N I ( ξ)ρ(ξ )by (ξ) Jξ dξdηdζf b ( ξ)− 1−1−1 zI z (E4.3.12b)where we have transformed the integral to the parent domain. The integral over the parent domainis evaluated by numerical quadrature.To obtain the external nodal forces due to an applied pressure t =− pn , we consider asurface of the element.
For example, consider the external surface corresponding with the parentelement surface ζ =− 1; see Fig. 4.6. The nodal forces for any other surface are constructedsimilarly.On any surface, any dependent variable can be expressed as a function of two parentcoordinates, in this case they are ξ and η . The vectors x,ξ and x,η are tangent to the surface. Thevector x,ξ × x,η is in the direction of the normal n and as shown in any advanced calculus text, itsmagnitude is the surface Jacobian, so we can write()pndΓ = p x,ξ × x,η dξdη(E4.3.13)For a pressure load, only the normal component of the traction is nonzero. The nodal externalforces are then given by1 1fiIext= ∫ ti NI dΓ = − ∫ pniN I dΓ= − ∫ΓΓ∫ p eijk x j ,ξ xk ,η N I dξdη(E4.3.14)−1−1where we have used (E4.3.13) in indicial form in the last step. In matrix form the above isf Iext =− ∫ pNI x ,ξ × x ,η dΓ(E4.3.16)Γ4-44T.
Belytschko, Lagrangian Meshes, December 16, 1998We have used the convention that the pressure is positive in compression. We can expand theabove by using Eq. (4.4.1) to express the tangent vectors in terms of the shape functions andwriting the cross product in determinant form, giving ex= f xIe x + f yIe y + f zI e z =− ∫ ∫ pN I det xJ NJ ,ξ−1 −1 xK NK ,η1 1fextIeyy J NJ, ξyK N K, ηez z J N J ,ξ dξdηzK N K ,η (E4.3.17)This integral can readily be evaluated by numerical quadrature over the loaded surfaces of theparent element.Example 4.4. Axisymmetric Quadrilateral. The expressions for the rate-of-deformationand the nodal forces for the axisymmetric quadrilateral element are developed.
The element isshown in Fig. 4.7. The domain of the element is the volume swept out by rotating thequadrilateral 2π radians about the axis of symmetry, the z-axis. The expressions in indicialnotation, Eqs. (4.5.3) and (4.5.6), are not directly applicable since they do not apply to curvilinearcoordinates.z43Ω12θxyrFig.