Belytschko T. - Introduction (779635), страница 95
Текст из файла (страница 95)
Substituting (28c),(37), and (40a) into (26b), we obtainQxQyq= 4µ qxy(1.4.41)It can be seen already that because of the way the assumed strain field was designed, thenonconstant part depends only on the shear modulus µ and is independent of the bulkmodulus.Evaluating the stiffness by (26c),Ke = B t 0 CB 0 + KstabeKstab=ec1 ggtc2 ggtc2 ggtc3 ggt(1.4.42a)(1.4.42b)B(0) is given by (1.2.39). Constants for OI stabilization as well as those for QBI aregiven in Table 1.
QBI stabilization is derived in the same way as OI with E given by (32).Table 1. Constants for the mixed method stabilization matrixStabilizationOIQBIc14µHxx H*2µ(1+ν)Hxx H**c24µHxy H*2µν(1+ν)Hxy H**c34µHyy H*2µ(1+ν)Hyy H**Note. ν ≡ ν for plane stress and ν ≡ ν/(1-ν) for plane strain;Hxx HyyHxx HyyH* ≡H** ≡Hxx Hyy - H2xyHxx Hyy - νH2xy46Explicit approximate expressions for Hij can be obtained by integrating (28) inclosed form assuming the Jacobian to be constant within the element domain.Hxx =Hyy =t2t2t2t2L 1yL 1xtHxywhere+ L 2y3A+ L 2x3At(1.4.43)tt- L 1x L 1y - L 2x L 2y=3AtL 1 = -1, +1, +1, -1(1.4.44)tL2= -1, -1, +1, +1.The quantities used to evaluate Hij can also be used to evaluate the element area byttttA = 1 (L 1 x)( L 2 y) - ( L 1 y)( L 2 x)(1.4.45)4The first part of the stiffness corresponds to the one-point quadrature stiffness.
Thesecond part is the stabilization or rank compensating stiffness, which is of rank 2 and thusincreases the rank of the total stiffness from 3 (rank of one-point quadrature stiffness) to 5.This is the correct rank of QUAD4 according to Eq. (1.2.26a).The form given in (35) can be considered a canonical form for the stiffness matrix ofQUAD4 if the constants ci are arbitrary. For any ci, this element stiffness will satisfy thepatch test. The constants ci can be varied to improve the performance of the element forspecific problem classes, but as shown by numerical studies in Belytschko and Bachrach(1986), the rate of convergence will always be the same, provided the element does notlock.
When the stabilization matrix is independent of the bulk modulus λ + 23 µ, theelement will not lock volumetrically.This development also provides guidance about the design of stabilization proceduresin nonlinear problems which are based on material properties.
If the current linearizedvalue of µ can be estimated, then (34) provides a stress-strain relation between Qi and qi.1.4.6 Frame Invariance. The elimination of the nonconstant part of the shear strain as isdone with the OI and QBI elements improves their performance in bending problems. Thecantilever test problems of the next section demonstrates the excellent coarse mesh bendingaccuracy of QBI; however, elimination of the nonconstant shear strain also causes theelement to lose frame invariance.
For most problems, the effect is negligible, but forcoarse mesh bending, the effect can be significant.Elements based on the OI or QBI assumed strain field can be made frame invariant byevaluating the stabilization matrix using an orthogonal, local coordinate system that isaligned with the element. The local coordinate system, called the (x, y) system, is relatedto the global coordinate system by a rotation matrix R.47The nodal coordinates, x and y, evaluated in the local (x, y) coordinate system, arerenamed x and y.
They are obtained byx R (1.4.46)y where R consists of the standard two dimensional rotation matrix arranged in an 8x8 matrixto transform each pair of nodal coordinates, xI and yI.x =y cosθI4x4R=sinθI4x4(1.4.47)-sinθI4x4 cosθI4x4where I4x4 is a rank 4 identity matrix.The angle between the global and local coordinate system can be defined byttanθ = L 1 y L t1 x(1.4.48)This definition aligns the x axis with the referential ξ axis of the element as shown in Fig.8. This definition may not be appropriate for anisotropic material. This point is discussedfurther in the explicit formulation that follows.3ηyy2xθξ4x1Figure 8. Local (x, y) coordinate system aligned with ξ axis of an elementWe evaluate bx , by , and g by substituting x and y for x and y in (1.2.40) and(1.3.15).tbx = 1 y24 , y31 , y42 , y132Atby = 1 x42 , x13 , x24 , x312AyIJ = yI-yJ xIJ = xI-xJg = 1 h+ ht x bx + ht y by4(1.4.49a)(1.4.49b)(1.4.49c)(1.4.50)48Hats are added to terms to indicate that local coordinates are used in their evaluation.
Hxx ,Hyy , and Hxy , are evaluated by Eq. (44) with x and y substituted for x and y. Likewise,the constants in Table 1 are evaluated in terms of Hxx , Hyy , and Hxy and are renamed c1 ,c2 , and c3 . The stabilization matrix is analogous to (42b) and is given bystabKetc1 gg=tc2 gg(1.4.51)ttc2 ggc3 ggIn order to add the element stabilization matrix to the global stiffness matrix, it must betransformed back to the global coordinate system byK stab= Rt Kstabee RIt is simple enough to evaluateKstabe=(1.4.52)K stabec*1 ggtc*2 ggtc*3 ggc*2 ggin closed form astt(1.4.53)wherec*1 = c1 cos 2 θ + c3 sin 2 θ - 2c2 sinθcosθc*2 = c2 cos 2 θ - sin 2 θ + c1 - c3 sinθcosθc*3 = c3 cos 2 θ + c1 sin 2 θ + 2c2 sinθcosθ(1.4.54)1.4.7 Hourglass Control Procedure.
We seek to evaluate internal forces directly by thefirst equality in Eq. (26c) which in corotational coordinates isfintt=Bs(1.4.55)For the OI and QBI strain and stress fields, Eq. (55) can be shown to take the form of aone-point element plus stabilization forces:fintt= AB 0 s 0 + fstab(1.4.56)A procedure for large deformation, nonlinear problems based on the previous analysis ofthe mixed element is described. The mixed field OI given in Section 1.4.5 will be used forthis purpose. Implementations based on other assumed strain fields can be developedsimilarly.The development hinges on the fact that the linear theory developed in Section 1.4.5is identical to the nonlinear theory if all variables are interpreted as rates.
Thuscomponents of the generalized hourglass strain rates can be obtained by Eq. (26a) writtenin the form49ttqx = 1 Hxx Hyy g v x + Hxy Hyy g v yHt(1.4.57)tqy = 1 Hxy Hxx g v x + Hxx Hyy g v yHwhere the superposed carets indicate that the quantities are evaluated using a corotationalcoordinate system. This formulation results in a frame invariant element. The corotationalcoordinate system is equivalent to the local coordinate system presented in the previoussection, however it is embedded in the element and rotates with the element as the elementdeforms.
The corotational coordinate system can be embedded by various techniques, andfor anisotropic materials, it is advantageous to embed the coordinate system so that itcoincides with axes of orthotropy or other directional features of the material.The corotational coordinate system can also be used to evaluate the rate-ofdeformation and update the stress at the quadrature point. An advantage of the corotationalsystem is that a frame invariant stress rate is not needed for large deformation problems.The stress-strain law for the generalized hourglass strain rates and stress rates isgiven byQi = 4µqi(1.4.58)This relation involves the shear modulus µ and assumes an isotropic material response.The shear modulus, µ is obtained by taking the ratio of the effective deviatoric stress ratesand strain rates.2µ =s ijs ij12(1.4.59a)e ije ijs ij = σij - 1 pδij3(1.4.59b)e ij = εij - 1 εkk δij3In two-dimensional plane stress problems222s ijs ij = σx - σx σy + σx + 3σxy222e ije ij = εx - εx εy + εx + 3εxyThe stabilization stresses are then updated by(1.4.59c)(1.4.60a)(1.4.60b)50t n+1n+1Qi=nQi +Qi dt(1.4.61a)tnwhich for the central difference method givesn+1Qinn+1/2= Qi + ∆t Qi(1.4.61b)The stabilization internal forces, evaluated by (26b), are then given byfstab= (Hxx Qx -Hxy Qy )gx(1.4.62a)fstab= (Hyy Qy -Hxy Qx )gy(1.4.62b)Not that the stress and strain which is used to evaluate the shear modulus is marked withhats to indicate that these are corotational quantities.
This is not necessary since the shearmodulus is an invariant quantity for isotropic material.The assumptions made in this development is that the material response is uniformover the element and the deviatoric response is isotropic. The second assumption can beavoided by using a C matrix based on a fully anisotropic C. However, this entailsavailability of C in the computational process, and in procedures such as radial return forelastoplasticity, C is not available. The first assumption is more troublesome; as theelastic-plastic front passes across an element, one-point quadrature is not as effective inresolving the behavior along the boundary.
This effect has been noted in Liu et al. (1988).Usually, however, the substantially reduced cost of one-point quadrature elements allowsmore elements to compensate for this effect. Adaptive schemes with automatic meshrefinement in zones of rapidly varying material behavior are also effective.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
