Belytschko T. - Introduction (779635), страница 93
Текст из файла (страница 93)
Adding these two rows corresponds to adding 2 generalized strains.The matrices for the one-point quadrature QUAD4 are thenB=b Tx0C11 C12 C13000bTyC11 C22 C2300bTybTxC13 C23 C3300gT0000CQ00gT0000CQC=(1.3.1a)s T = σx , σy , σxy , Qx , Qy(1.3.1b)eT = ε x , ε y , 2ε xy , qx , qy(1.3.1c)32where C is the constitutive matrix augmented by two rows and columns in terms of aconstant to be determined (CQ); s and e are the stress and strain matrices augmented by thegeneralized stresses and strains (Qx, Qy) and (qx, q y), respectively.To maintain linear consistency for the element, these additional generalized strainsshould vanish when the nodal displacements (or velocities) emanate from linear fields.Consistency and stability (rank sufficiency) are essential requirements for a soundnumerical method.The requirement that qx = qy = 0 for linear fields impliesgT uLin = gT α 0 s + α 1 x + α 2 y = 0∀α i(1.3.2)The above must be satisfied for both ux and uy so we have not specified the component;uLin is taken from Eq.
(1.2.15b). The above can be interpreted as an orthogonalitycondition: g must be orthogonal to all linear fields.1.3.1 The Gamma Vector. We first define a set of four vectors, b*.b* ≡ bx , by , s , hTo obtain g, two properties of b* are used:1.
the vectors bi are biorthogonal to x j2. the vectors b*i are linearly independent.The biorthogonality property, given bybti x j = δij (i, j) = 1 to 2is an identity which holds for all isoparametric elements:(1.3.3)(1.3.4)∂Nx j = δij(1.3.5)∂xiThe demonstration of this identity is based on the isoparametric mapping, Eq.(1.2.3c), which when combined with (5) gives∂x∂Nx j = j = δ ij∂xi∂xi(1.3.6)where the last equality expresses the fact that in two dimensions, for example,∂x/∂x=∂y/∂y=1, ∂x/∂y=∂y/∂x=0. Eq. (5) holds for the derivatives of the shape functionsat any point.
In particular, it also holds for the point ξ = η = 0 in QUAD4. Additionalorthogonality conditions,bti s = 0bti h = 0ht s = 0i=1 to 2(1.3.7)can easily be verified by arithmetic using the definitions of these vectors.The linear independence of the 4 b*i vectors is demonstrated as follows. Assume b*iare linearly dependent.
Then it follows that there exists α i ≠ 0 such that33α 1bx + α 2by + α 3s + α 4h = 0(1.3.8)Premultiplying the above by s t , and using (7) yields α 3 = 0. Similarly premultiplying byht yields α 4 = 0. Then premultiplying x t yieldsα1 + α3x ts + α4x th = 0(1.3.9)and since it has just been determined that α 3 = α 4 = 0, it follows that α 1 = 0. Similarly,premultiplying by y t shows that α 2 = 0.
Thus α i = 0, for i = 1 to 4, and it follows that thevectors b*i are linearly independent.The preceding developments are now used as tools for the construction of g, via theconsistency requirement (2). Since the vectors b*i are linearly independent, they span R4.It follows that any vector in R4, including g can be expressed as a linear combination ofb*i :g = β1 bx + β2 by + β3 h + β4 s(1.3.10)where β i are constants to be determined by the linear consistency requirement (2).Substituting (10) into (2) and collecting the coefficients of α i yieldsα o β 1 btx s + β 2 bty s + β 3 ht s + β 4 s t s+ α 1 β 1 btx x + β 2 bty x + β 3 ht x + β 4 s t x(1.3.11)+ α 2 β 1 btx y + β 2 bty y + β 3 ht y + β 4 s t y = 0Since the above must vanish for all α i, each coefficient of α i must vanish. Taking thecoefficient of α o and simplifying by means of Eqs.
(4) and (7) givesβ 4 s t s = 4β 4 = 0Using (12) and (6) to simplify the coefficient of α 1 in (11) gives(1.3.12)β1 + β3 ht x = 0and a similar procedure for the coefficient of α 2 gives(1.3.13)β2 + β3 ht y = 0(1.3.14)Using Eqs. (13) and (14) in (10) to express β 1 and β 2 in terms of β 3 and using (12) yieldsg = β3 h - ht x bx - ht y by(1.3.15)The constant β 3 remains undetermined, since for any value of β 3 the vector g is orthogonalto all linear fields. It will be convenient later to have g t h=1, so we set β 3 = 1/4. Thevalue of β 3 = 1 was used in Flanagan and Belytschko (1981) because the reference elementwas a unit square; this changes some of the subsequent constants but not the substance ofthe development.
In this description we choose β 3 = 1/4, which givesg = 1 h - ht x bx - ht y by4The above expression can be written in indicial notation as(1.3.16a)34org = 1 h - ht x i bi4(1.3.16b)γ I = 1 hI - hJxiJ biI4Using (4), (7), and (16a) the following are easily verified by:g tx = g ty = g ts = 0g th = 1(1.3.16c)(1.3.17)1.3.2 Stabilization Forces and Stiffness Matrix.
Replacing B(0) and s (0) in Eq. (1.2.39)by the augmented matrices B and s gives the nodal forcesf int = A bx00bybybxσxσyσxy+AQx gQy g(1.3.18)= AB T(0)s (0) + f stab(1.3.19)The first term in the internal force is obtained by one-point quadrature. The generalizedstresses and are strains are obtained by the stress-strain law, s = Ce, and the straindisplacement equation e = Bd:Qx = CQqxg tuQy = CQqy(1.3.20)g tuqx =x qy =yThe stiffness matrix is obtained by substituting replacing B(0) and C in (1.2.42) bythe augmented matrices B and C which givesstabK e = K 1pte + Ke(1.3.21a)tK 1pte = AB (0)CB(0)(1.3.21b)where gg t 0 QK stab=AC(1.3.21c)e 0 gg t stabK 1ptis obtained from thee is the stiffness matrix obtained from one-point quadrature.
K eadditional generalized strains which were introduced to stabilize the element and is oftencalled a stabilization matrix. The stabilization matrix is of rank 2. Combined with the onepoint quadrature stiffness, it yields a matrix of rank 5, which is the correct rank for theQUAD4.This form of the linearly consistent generalized strains occurs in many stabilizationprocedures for underintegrated elements, and will be designated as the g vector.
Note thatthe vector is not completely determined by linear consistency: an unspecified constant β 3remains. This vector is orthogonal to the nodal displacements which emanate from a linear35field for arbitrarily shaped quadrilaterals. When the element is rectangular, ht x = ht y = 0and g = β 3h.1.3.3 Scaling the Stabilization Forces. Since the constants CQ in Eq.
(1) are not truematerial constants, it is important to provide formulas for these constants which provideapproximately the same degree of stabilization regardless of the geometry and materialproperties of the element. The basic objective is to obtain a scaling which perturbs theelement sufficiently to insure the correct rank but not to overwhelm the one-pointquadrature stiffness.One procedure for selecting CQ is to scale the maximum eigenvalue of the stabilizationstiffness to the maximum eigenvalue of the underintegrated stiffness. In fact, it would bedesirable to shift eigenvalues associated with hourglass modes out of the spectrum that is ofinterest in the response.
The hourglass modes in a fully integrated element are smaller thanthe 1-point quadrature element eigenvalues. To avoid locking, the stabilization should be asmall fraction of the one point quadrature eigenvalue.According to Flanagan and Belytschko (1981), the maximum eigenvalue for the 1point quadrature version of QUAD4 for an isotropic material is bounded by1 222Ac b ≤ λ max ≤ Ac b2b = ∑ bTi bi(1.3.22a)c2 = λ + 2µ(1.3.22b)i=1The eigenvalues of K e are given by the eigenvalue problemKx=λx(1.3.23)The eigenvalue associated with the stabilization can be estimated by letting x = dHx in theRayleigh quotient, which with (20) and the orthogonality properties (4) and (7) givestλ = x Kx =x txACQht ggt h(1.3.24)ht hwhere the second equality follows because K 1pte h = 0.
Using Eq. (17), it can be seen thatλ = ACQ 4(1.3.25)Using Eqs. (22) and (25) it follows that the eigenvalue associated with the stabilization isscaled to the lower bound on the maximum eigenvalue of the element ifCQ = 2α s c2 b = 2α s (λ+2µ)2∑ btibii=1(1.3.26)where α s is a scaling parameter.In Flanagan and Belytschko (1981), the hourglass control parameter was scaled bythe dynamic eigenvalue, i.e., the frequency, of the element. However, since the hourglasscontrol is strictly an element-stiffness related stress, this seems inconsistent and a purestiffness scaling is more appropriate.361.4 Mixed Method Hourglass StabilizationThe mixed variational principles are another vehicle for developing four-nodequadrilaterals which do not lock.
Furthermore, they can be used to develop rankcompensation procedures for underintegrated elements (stabilization matrices) which donot involve any arbitrary perturbation parameters and are based on the material propertiesand geometry of the element.1.4.1 Displacement Field of QUAD4 in Terms of Biorthogonal Basis. Before developingelements based on a mixed method, the displacement field in the element is expressed in aform in which the parts which cause locking can easily be identified and corrected.Furthermore, when this expression for the displacement field is used, the stiffness matrixcan be obtained in closed form without any numerical integration.
This is useful forunderstanding its properties and for implementation in vector method programs.The procedure described here is based on Belytschko and Bachrach (1986). As apreliminary to developing this expression (which will be called the Belytschko-Bachrachform), the basis vectors b** and x * are defined so they are biorthogonal in R4:x *β t b**α = δαβα,β = 1 to 4(1.4.1)b** = bx , by , g, s *(1.4.2a)x * = x, y, h, s(1.4.2b)wheres * = 1 s - s t x bx - s t y by(1.4.3)4The vector s * is obtained by orthogonalizing s to x and y.
The arbitrary constant whichemerges is chosen to be 1/4 so that s t s * = 1. Most of the identities involved in (1) havealready been proven; see Eqs. (1.3.4-7); the remaining ones are easily verified using (3)and (1.3.7).To develop the Belytschko-Bachrach form, we start by expressing the displacementfield asu x,y = α o + α 1x + α 2y + α 3ξη(1.4.4)Only a single component is considered since the procedure for both components isidentical. Evaluating the above at the 4 nodes givesuI = u xI,yI = α o + α 1 xI + α 2 yI + α 3 ξIη I(1.4.5)It is easily shown that ξIη I = hI, so writing the above in matrix form givesu = αos + α1x + α2y + α3h(1.4.6)*which is a linear combination of the linearly independent x vectors. Linear independenceof x * follows from the biorthogonality of x * and b** .37We now exploit this biorthogonality to evaluate α i.
Premultiplying (6) by s * t andinvoking the orthogonality conditions, s * t x = s * t y = s * t h = 0 yieldsα o = s *t u(1.4.7)Similarly premultiplying respectively by btx , bty and g t and using the biorthogonality(4) yieldsα 1 = btx u(1.4.8a)α 2 = bty u(1.4.8b)α3 = g tu(1.4.8c)Substituting the above into (4) yieldsu x,y = s *t + xbtx + ybty + ξηgt u(1.4.9a)The two components of the displacement field can be expressed in the same formux x,y = s t + xbtx + ybty + hgt ux(1.4.9b)uy x,y = s t + xbtx + ybty + hgt uy(1.4.9c)h=ξη(1.4.9d)This is the same interpolation as the standard isoparametric form (1.2.1), however, thisexpression will more clearly reveal what causes locking and how to eliminate it.1.4.2 Volumetric Locking.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
