Belytschko T. - Introduction (779635), страница 94
Текст из файла (страница 94)
The four node quadrilateral locks in plane strain forincompressible materials when it is fully integrated. The cause of volumetric locking canbe explained by considering a mesh of elements in plane strain with fixed boundaries ontwo sides as shown in Fig. 6. Consider the element in the lower left-hand corner, element1. The nodal displacements of the element for an incompressible material must preservethe total volume of the element (or to be specific, the area in plane strain, since constantvolume implies that the area be constant). If we consider small displacements, the onlydisplacements of node 3 which maintain constant area areux3 = −αa(1.4.10)uy3 = + αbwhere α is an arbitrary parameter; all other nodal displacements of element 1 are zero dueto the boundary conditions.
Differentiating (9b) and (9c), we obtain the dilatationthroughout the element.∂h∂h∆ = ux,x + uy , y = b tx ux + bty uy + g t ux + g t uy(1.4.11a)∂x∂ySubstituting Eq. (10) into (11a), the constant part of the dilatation drops out leaving38∂h∂h1∆ = 4α(b − a )(1.4.11b)∂y∂xwhich only vanishes everywhere except along a line that passes through the origin.For rectangular elements as in Fig.
7, Eq. (11b) simplifies toα∆ = ab(bx−ay)(1.4.12)where (x, y) is a local coordinate system. The volumetric strain is non zero everywhereexcept along the line y = (b/a)x; therefore, for an incompressible material the volumetricenergy will be infinite if the strain energy is evaluated exactly as is the case in a fullyintegrated element. Thus node 3 will not be able to move; nodes 2 and 3 then provide arigid boundary for the left hand side of element 2, and it can similarly be shown that byusing these arguments for element 2 that node 6 cannot move. This argument can then berepeated for all nodes of the mesh to show that deformation of the mesh is impossible.This argument also applies to meshes of skewed elements as in Fig.
6.Fixed boundary:ux = uy = 04b3112aFigure 6. Mesh of quadrilateral elements with fixed boundaries on two sides7b9834y43112x625aFigure 7. Partial mesh of rectangular elements fixed on two sidesAnother way to examine this behavior is consider an arbitrary element deformation asexpressed by Eq. (4).ux = α 0x + α 1x x + α 2x y + α 3x h39(1.4.13)uy = α 0y + α 1y x + α 2y y + α 3y hThe dilatation can be evaluated by differentiating.∆ = ux,x + uy,y = α 1x + α 2y + α 3x∂h+ α 3y∂h(1.4.14)∂x∂yWe can evaluate the change in area of the element by integrating the dilatation over theelement domain.⌠∫Ω e∆dA = (α 1x + α 2y + α 3x∂h+ α 3y∂h)dA(1.4.15)∂x∂y⌡Ω eWe can show algebraically using (1.2.7b) and (1.2.8a-d) an important property of h(ξ,η):⌠ ∂h dΩ = ⌠ ∂h dΩ = 0∂x⌡Ωe⌡Ωe∂yTherefore (15) is trivial to integrate and the change in element area is∫Ω e∆dA = (α1x + α2y )A(1.4.16)(1.4.17)which is zero only for α 2y = −α 1x .
If we now consider volume preserving elementdeformation, i.e. α 2y = −α 1x , the dilatation is∆ = α 3x∂h+ α 3y∂h(1.4.18)∂x∂yThis dilatation will be non zero everywhere within the element except along the curveα 3x ∂h ∂x = -α 3y ∂h ∂y, even though the overall element deformation is volumepreserving.
Thus it becomes apparent that locking arises from the inability of the elementto represent the isochoric field associated with the hourglass mode, as reflected by α 3x andα 3y in the above equations. To eliminate locking, the strain field must be designed, orprojected, so that in the hourglass mode the dilatation in the projected strain field vanishesthroughout the element. In more general terms this may be stated as follows: to avoidlocking, the strain field must be isochoric throughout the element for anydisplacement field which preserves the total volume of the element.
In particular, inthe quadrilateral, the dilatation must vanish in the entire element for the hourglass mode,because this displacement mode is equivoluminal.1.4.3 Variational principle The weak form corresponding to the Hu-Washizu variationalprinciple is given for a single element domain by0 = δπ(u,e,s ) = ∫Ωeδe t CedΩ +δ ∫Ω s t (e- D u)dΩ - δWexte40= ∫Ω [δe t (Ce-s ) - δs t (e- D u) + δ( D u)t s ]dΩ - δdt f exte(1.4.19)where δ denotes a variation, e is the interpolated strain, and s the interpolated stress. Duis the symmetric displacement gradient which would be equivalent to the strain in adisplacement method. In mixed elements that are derived from the Hu-Washizu variationalprinciple, the displacement gradient is projected on a smaller space to avoid locking.
Theterm δWext designates the external work, f ext the external nodal forces. The domainchosen for (19) is a single element, but an arbitrary domain can also be assumed ifconnectivity is introduced into the subsequent development.The isoparametric shape functions are used to interpolate the displacement field,which when integrated, gives the symmetric displacement gradient asDu = Bd(1.4.20)We introduce additional interpolants for the strains and stresses.e = Ee(1.4.21)s = Ss(1.4.22)where the interpolation matrices, E and S , and the augmented strains and stresses, e and swill be defined subsequently.
Substituting (20), (21), and (22) into (19), we obtain0 = ∫Ω [δe t Et (CEe-S s ) - δs t S t (Ee-Bd) + δdt B t S s ]dΩ - δdt f ext (1.4.23)eBy invoking the stationary condition on (19), we obtaintCe = E sEe = BdtB s = f ext(1.4.24a)(1.4.24b)(1.4.24c)C≡E t CE dΩ(1.4.25a)S t E dΩ(1.4.25b)S t B dΩ(1.4.25c)whereΩeE≡ΩeB≡ΩeWe obtain expressions for e, s , and the stiffness matrix from (24a-c).-1e = E Bd-1 ts = E Ce(1.4.26a)(1.4.26b)41t-1 tt-1f ext = B s = B E CE Bd ≡ K ed(1.4.26c)1.4.4 Strain Interpolations to Avoid Locking. The strain-field associated with thedisplacement field (9) can be obtained by straightforward differentiation which gives:∂uxbtx +∂x∂uyDu ==∂y∂ux∂y+∂uyεox + qx=εoy + qy2εoxy + qx∂h∂x0bty +∂x∂h∂h∂ygt0bty +gt btx +∂h∂y∂h∂xgtux= Bduy(1.4.27a)gt∂h∂x∂h(1.4.27b)∂y+ qy∂h∂y∂xwhere the naughts indicate the constant part of the strain field.In Section 1.4.2, it was demonstrated that QUAD4 elements of incompressiblematerial can lock when the dilatational energy at any point other than the origin is includedin the stiffness matrix.
It was also shown that this is caused by the dilatational fieldassociated with the hourglass modes, dHx and dHy, which in a fully integrated elementalways leads to non-vanishing dilatation. From Eq. (27b), it can be seen that thehourglass modes generate the nonconstant part of the volumetric field.In constructing a strain interpolant which will not lock volumetrically, we then havetwo alternatives:1. the nonconstant terms of the first two rows of Eq. (27b) can be dropped2.
the first two rows can be modified so that no volumetric strains occur in thehourglass modes.The first alternative leads to the strain fieldεoxe=εoy2εoxy + qx ∂h ∂y + qy ∂h ∂x(1.4.28a)42This can be written in the form of Eq. (17) by lettingE=et =1000001000001∂h ∂y∂h ∂xεox ,εoy , 2εoxy , qx ,qy(1.4.28b)(1.4.28c)The second alternative is to define the strain field byεox + qx ∂h ∂x - qy ∂h ∂ye=εoy + qy ∂h ∂y - qx ∂h ∂x(1.4.29)2εoxy + qx ∂h ∂y + qy ∂h ∂xIn Eq. (29), the dilatation ε x + ε y still vanishes in the hourglass mode, since regardless ofthe value of qx and qy, Eq. (29) yieldsεx + εy = εox + εoy(1.4.30)The question then arises as to which of the two alternatives, (28) or (29), ispreferable.
The field in (29) is frame invariant whereas (28) is not. However, thecomputations associated with (28) are simpler. However, neither of these are particularlyattractive for most problems from the viewpoints of accuracy and efficiency.For elements which involve beam bending, the performance of the element can beimproved strikingly by omitting the nonconstant part of the shear field.
This shear strainfield cannot be combined with the extensional strains in (28) because the strain field wouldthen only contain three independent functions, and the element would be rank deficient.Therefore this shear strain field is combined with the extensional strains in (29), whichgivesεox + qx ∂h ∂x - qy ∂h ∂ye=εoy + qy ∂h ∂y - qx ∂h ∂x(1.4.31a)εoxyFor this elementE=100∂h ∂x -∂h ∂y010-∂h ∂x ∂h ∂y00100(1.4.31b)43with e given in (28c). The strain field (31) leads to the "Optimal Incompressible" or OIelement in Belytschko and Bachrach.
This element performs well in beam bendingproblems when one set of element sides are parallel to the axis of the beam and theelements are not too distorted.The performance of QUAD4 in bending can be enhanced even further for isotropic,elastic problems by using a strain field which depends on Poisson's ratio as follows:εox + qxe=εoy + qy∂h∂x∂h∂y- νqy- νqx∂h∂y∂h(1.4.32)∂x2εoxyν for plane stresswhere ν ≡ ν/(1-ν) for plane strainThis is the field called "Quintessential Bending and Incompressible" or QBI in Belytschkoand Bachrach (1986), which has great accuracy in bending with linear elastic material.
Fora rectangle, as shown by Fröier et al. (1974), this element corresponds to the incompatiblequadrilateral of Wilson et al. (1973).1.4.5 Stiffness Matrix for OI Element. In order to gain more insight into these mixedelements and to see how they are used to construct stabilization (rank-compensating)matrices which do not involve arbitrary parameters, the stiffness matrix for the OI element,which is based on the strain field (31b) will be developed.The stress field is chosen to bes=st =100∂h ∂x00100∂h ∂y00100σox ,σoy , σoxy , Qxs ≡SsQy(1.4.33.a)(1.4.33b)Using (16), and integrating (25b) and (25c), we obtain E and B as follows:AI3x30 3x20 2x3Hxx -Hxy-Hxy HyyE=(1.4.34)44where A is the area of the element, and∂h ∂xi ∂h ∂xj dΩHij =(1.4.35)ΩeThe B matrix is given byB=Abtx00AbtyAbtyAbtxHxx gt00Hyy gt(1.4.36)The inverse of E is given by-1E =1 I3x3A0 2x30 3x2(1.4.37)1 Hyy HxyH Hxy HxxwhereH = H xx Hyy - H2xyUsing (26a), (36) and (37), we can evaluate e to be(1.4.38)εox = btx uxεoy = bty uy2εoy = btx uy +bty ux(1.4.39)qx = 1 Hxx Hyy gt ux + H xy Hyy gt uyHqy = 1 Hxy Hxx gt ux + H xx Hyy gt uyHFrom (39) we can see that in the mixed element, the constant parts of the strain fieldcorresponds exactly to the constant part which emanates from the displacement field D ugiven in (27).
The nonconstant part depends strictly on the hourglass mode (anycomponent of ux or uy which is not orthogonal to g). The effect of the projection is tomodify this part of the strain field so that the volumetric strains vanish.45To complete the evaluation of the element stiffness, we obtain C by integrating(25a).AC3x30 3x20 2x34µHxx -4µHxy-4µHxy 4µHyyC=(1.4.40a)For a linear isotropic material, C is given by λ+2µ λ 0C=λλ+2µ00µ0(1.4.40b)For plane strain, λ=λ=2µν/(1-2ν); for plane stress λ=2µν/(1-ν).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
