Belytschko T. - Introduction (779635), страница 40

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 40)

The most widely used procedure for numericalintegration in finite elements is Gauss quadrature. The Gauss quadrature formulas are, see forexample Dhatt and Touzot (1984, p.240), Hughes (1977, p. 137)nQ∫ −1 f (ξ )dξ = Q=∑1wQ f (ξQ )1(4.5.24)where the weights wQ and coordinates ξQ of the nQ quadrature points are available in tables; ashort table is given in Appendix C. Equation (4.5.24) integrates f (ξ ) exactly if it is a polynomialof order m ≤ 2nQ −1. Equation (4.5.24) is tailored to quadrature over the parent element domains,since it is over the interval [-1, 1].To integrate over a two dimensional element, the procedure is repeated over the seconddirection, yielding the followingnQ1 nQ21 1∫∆ f ( ξ)d∆ = ∫ ∫ f (ξ,η)dξdη = Q∑=1 Q∑=1wQ wQ f (ξ Q ,ηQ )1−1−1121(4.5.25)22In three dimensions, the Gauss quadrature formula isn Q1 n Q2 nQ31 1 1∫∆ f ( ξ)d∆ = ∫ ∫ ∫ f (ξ)dξdηdζ = Q∑=1Q∑=1Q∑=1wQ wQ wQ f ( ξQ ,ηQ1−1−1−1122312,ζ Q3)(4.5.26)3For example, in integrating the expression for the internal nodal forces over the biunit squareparent element in two dimensions, we usefint= ∫ B {σ}Jξ d∆ =T∆=nQ1 n Q2∑ ∑ wQ wQQ1 =1 Q2 =111 1∫ ∫BT−1−1T2({σ} Jξ dξdη){ (B ξQ1 , ηQ2 σ ξQ1 , ηQ2)} Jξ (ξQ , ηQ )1(4.5.27)2To simplify the notation for multi-dimensional quadrature, we often combine the weights intoa single weight and write the quadrature formula in any dimension as4-27T.

Belytschko, Lagrangian Meshes, December 16, 1998wQ f (ξ Q )∫∆ f (ξ)d∆ = ∑Q(4.5.28)where wQ is a products of the weights for one-dimensional quadrature wQ .The number of quadrature points used in nonlinear analysis is generally based on the samerules as for linear analysis; the number of quadrature points is chosen to exactly integrate the nodalinternal forces for a regular element.

A regular form of an element is one that can be obtained byonly stretching but not shearing the parent element; for example, a rectangle for two-dimensionalisoparametric elements. To choose the number of quadrature points for the internal nodal forcesfor a 4-node quadrilateral, we use the following argument.

The rate-of-deformation and the Bmatrix are linear in this element since the velocities are bilinear. If the stress is linearly related tothe rate-of-deformation, it will vary linearly within the quadrilateral element. The integrand for theinternal nodal forces is approximately quadratic, since it is a product of the B-matrix and thestresses. By the above rule for Gauss quadrature, two quadrature points are then needed in eachdirection for exact quadrature of a quadratic function, so 2x2 quadrature on a quadrilateralintegrates the internal nodal forces exactly on a regular element.

Quadrature formulas whichintegrate the nodal internal forces almost exactly for a linear constitutive equation are called fullquadrature.Gauss quadrature is very powerful for smooth functions which are polynomials or nearlypolynomials. In linear finite element analysis, the integrand in the expression for the stiffnessmatrix consists of polynomials for rectangular elements and is smooth and nearly a polynomial forisoparametric elements. In nonlinear analysis, the integrand is not always smooth. For example,for an elastic-plastic material, the stress may have a discontinuous derivative in space at the surfaceseparating elastic and plastic material. Even if the stress-strain law is smooth for an elastic-plasticmaterial, the derivative of the stress with respect to the strain changes drastically when the responsechanges from elastic to plastic, so the effect is the same.

Therefore, the errors in Gauss quadratureof an element that contains an elastic-plastic interface are likely to be large. However, higher orderquadrature is not recommended for circumventing these errors, since it often leads to stiff behavioror locking of elements.4.5.5.Selective-ReducedIntegration.

For incompressible or nearly incompressiblematerials, full quadrature of the nodal internal forces may cause an element to lock, i.e. thedisplacements are very small and do not converge or converge very slowly. The easiest way tocircumvent this difficulty is to use selective-reduced integration.In selective-reduced integration, the pressure is underintegrated, whereas the remainder of thestress matrix is fully integrated. For this purpose, the stress tensor is decomposed into thehydrostatic component or pressure p, which is the trace of the stress tensor1σ ijdil = pδ ij = σ kkδ ij3(4.5.29a)and the deviatoric components:σ ijdev = σij − pδ ij(4.5.29b)The rate-of-deformation is similarly split into dilatational and deviatoric components which aredefined by4-28T. Belytschko, Lagrangian Meshes, December 16, 19981Dijdev = Dij − Dkkδ ij31Dijdil = Dkk δij3(4.5.30)It is noted that the dilatational and deviatoric components are orthogonal to each other so that thetotal virtual internal power as defined in Eq.

(4.3.19) becomes:δ P int = ∫ δ Dijσ ij dΩ=13Ω∫ δD pdΩ + ∫δDΩσ ijdevdΩdevijii(4.5.31)ΩAfter expressing the rate-of-deformation in terms of the shape functions by (4.4.7b) and (4.5.30),the dilatational and deviatoric integrands become:δDii p = δviI NI , i p(4.5.32a)andδDijdevσ ijdev =11NI , jδviI + N I ,iδv jI )σ devN I ,kδvkIδ ij σijdev(ij −23(4.5.32b)Using the symmetry of σ ijdev and the fact that the trace of σ ijdev vanishes, the deviatoricintegrand simplifies to:δDijdevσ ijdev =δviI N I ,j σ ijdev(4.5.33)Selective-reduced integration consists of full integration on the deviatoric term and reducedintegration on the dilatational term on δPint .

Thus, for a four-node quadrilateral, selective-reducedintegration gives:41δ P int = δviI  Jξ (0) N I, i (0) p(0 ) + ∑ wQJξ ξ Q N I, j ξ Q σ jidev ξQ  3Q =1( )( )( )(4.3.34a)Hence the selective-reduced integration internal force expression is:( )fiIintT41= f Iiint =  Jξ (0) NI , i (0) p(0 ) + ∑ wQ Jξ ξQ NI , j ξ Q σ jidev ξQ  3Q =1( )( )( )(4.3.34b)where, as indicated in the above, the single quadrature point for the reduced quadrature is thecentroid of the parent element. The deviatoric part is integrated by full quadrature using two pointsin each direction, for a total of four quadrature points; this is called 2x2 quadrature. This scheme issimilar to the scheme used in linear analysis of incompressible materials.

Selective-reducedschemes for other elements can be developed by similarly modifying selective-reduced integrationschemes given linear finite element texts; see Hughes(1979) for selective-reduced integrationprocedures for linear problems.4.5.6.Element Force and Matrix Transformations. Often element nodal forces andelement matrices must be expressed which in terms of alternate degrees of freedom, i.e. for a4-29T.

Belytschko, Lagrangian Meshes, December 16, 1998different set of nodal displacements. In the following, transformations are developed for nodalforces and element matrices.ˆ . We wish toConsider an element or assemblage of elements with nodal displacements dˆ byexpress the nodal forces for the nodal displacements d which are related to dˆdddd=Tdtdtorδˆd = Tδd(4.5.35)The nodal forces associated with d are then given byf = TT ˆf(4.5.36)This transformation holds because the nodal foces and velocities are assumed to be conjugate inpower, see Section 2.4.2. It is proven as follows. An increment of work is given byˆ TˆfδW =δdT f =δ d∀δd(4.5.37)Either set of nodal forces and virtual displacements must give an increment in work, since work isa scalar independent of the coordinate system or choice of generalized displacements.

Substituting(4.5.35) into (4.5.37) givesδdT f =δd T TTˆf∀δd(4.5.38)Since (4.5.38) holds for all δd , Eq.(4.5.36) follows.The mass matrix can be transformed similarly. We first consider the case where T isindependent of time. then the mass matrix for the two set of degrees of freedom is related byˆ TM = TT M(4.5.39)This is shown as follows. By Eq.(4.5.36)f inert = TTˆf inert(4.5.40)and using (4.4.18) for the two sets of degrees of freedom, we haveˆ ˆ˙vM˙v = T T M(4.5.41)If T is independent of time, from (4.5.35), ˙vˆ = T˙v , and substituting this into the above and usingthe fact that this must hold for the arbitrary nodal accelerations, we obtain (4.5.39).

If the T matrix˙ v , sois time dependent, then ˙vˆ = T˙v + Tˆ T˙d˙ + TT MˆT˙ ˙df inert = TT M(4.5.42)˙ usually depends on the nodal velocities, so in this case terms which are not linear inThe matrix Tthe velocities occur in the equations of motion.4-30T. Belytschko, Lagrangian Meshes, December 16, 1998A transformation similar to (4.3.39) holds for the linear stiffness matrix and the tangentstiffness discussed in Chapter 6:ˆ T,K = T TKˆ tan TKtan = T T K(4.5.43)These transformations enable us to evaluate element matrices in coordinate systems which simplifythe procedure as in example 4.6.

They are also useful for treating slave nodes, see example 4.5Example 4.1. Triangular 3-node element. The triangular element will be developed usingtriangular coordinates (also called area coordinates and barycentric coordinates). The element isshown in Figure 4.2. It is a 3-node element with a linear displacement field; the thickness of theelement is a. The nodes are numbered counterclockwise in the parent element, and they must benumbered counterclockwise in the initial configuration, otherwise the determinant of the mapbetween the initial and parent domains will be negative.y23x ξi , t1x2XξYξ2( 0, 1)i3(1, 0)2ξ131parent element(triangular coordinates)1XFig.

4.2. Triangular element showing node numbers and the mappings of the initial and current configurations tothe parent element.The shape functions for the linear displacement triangle are the triangular coordinates, soNI = ξI . The spatial coordinates are expressed in terms of the triangular coordinates ξI by4-31T. Belytschko, Lagrangian Meshes, December 16, 1998 x  x x 2   1 y  =  y1 y 21   1 1x3 ξ1  y3 ξ2 1  ξ3(E4.1.2)where we have appended the condition that the sum of the triangular element coordinates is one.The inverse of Eq. (E4.1.2) is given byξ y 1  1  23ξ2  =y31ξ  2A  y 3 12x32x13x21x2 y3 − x3y 2  x  x3 y1 − x1 y3 y x1 y2 − x2 y1 1(E4.1.4a)where we have used the notationx IJ = x I − x JyIJ = yI − y J(E4.1.3)and2A = x32 y12 − x12 y32(E4.1.4b)where A is the current area of the element.

As can be seen from the above, in the triangular threenode element, the parent to current map (E4.1.2) can be inverted explicitly. This unusualcircumstance is due to the fact that the map for this element is linear. However, the parent tocurrent map is nonlinear for most other elements, so for most elements it cannot be inverted.The derivatives of the shape functions can be determined directly from (E4.1.4a) by inspection:[N I , j ] = [ξI , j ]ξ 1,x= ξ2 , xξ 3, xyξ1, y  1  23ξ2, y  =y2 A  31ξ3, y  y12x32 x13 x21 (E4.1.5)We can obtain the map between the parent element and the initial configuration by writing Eq.(E4.1.2) at time t = 0, which gives X  X   1 Y  =  Y1 1   1X2Y21X3 ξ1  Y3 ξ2 1 ξ3 (E4.1.6)The inverse of this relation is identical to (E4.1.4) except that it is in terms of the initial coordinatesξ  Y23 11 ξ2  = Y31ξ  2A0  Y 3 12X32X13X 21X2Y3 − X3 Y2   x  X3 Y1 − X1Y3   y X1Y2 − X2Y1 12A0 = X32Y12 − X12Y32(E4.1.7a)(E4.1.7b)4-32T.

Характеристики

Список файлов книги