Belytschko T. - Introduction (779635), страница 97
Текст из файла (страница 97)
In addition, we define the standard B matrixby57neNDv = ∑D(N Iv I) ≡ Bv(1.5.21)I=1Substituting Eqs. (20) into (18) and using the orthogonality condition for t as before gives0 = δΠ = δv T∫Ω B Ts dΩ − δv Tf ext(1.5.22)Exploiting the arbitrariness of δv we obtain the discrete equilibrium equationsf int - f ext = 0(1.5.23)f int = ∫Ω B t s dΩwhere the stress is given by some nonlinear constitutive equation(1.5.24)s = (e° ,s ,...) = (Bv,s ,...)(1.5.25)The above formulation is applicable to problems with both material and geometricnonlinearities.
In applying the assumed strain stabilization procedure, it is convenient touse a corotational formulation as discussed in Section 1.4.7, where the Cauchy stressesand velocity strains are expressed in terms of a coordinate system (x, y) which rotates withthe element. As with mixed method stabilization of Section 1.4, a corotational coordinatesystem also assure that the element is frame invariant.The internal forces in a corotational formulation are given by∫fint = Ω B t s dΩ(1.5.26)where the superposed tildes indicate quantities expressed in terms of the corotationalcoordinates. The counterpart of (20) is°e = Bvand the rate form of the constitutive equation can be written(1.5.27)°s = Ce(1.5.28)where C is a matrix which depends on the stress and other state variables; for anincrementally isotropic hypoelastic material, C is given by (1.4.40b).The above form of a stress-strain law is objective (frame-invariant).
The spin is thengiven by∂v∂vω=1 y − x(1.5.29)2 ∂x∂yIn developing the hourglass resistance based on physical parameters, two assumptionsmust be made:1. the spin is constant within the element2. the material response tensor C is constant within the element.The velocity for the 4 node quadrilateral is given by a form identical to (1.4.9b)vi = (s t + xbtx + ybty + hgt )v i(1.5.30)58The spin (29) is then given byω = 12 (btx v y − bty v x + qy h,x − qx h,y ))t(1.5.31)tqx = g v xqy = g v y(1.5.32)Because of the orthogonality property (1.4.16), the average spin is given by1ω o = A∫Ω eωdΩ = 12(btx v y − bty v x )(1.5.33)This corresponds to the spin at the center of the element.
It can be seen from (31) that thestronger the hourglass mode, the more assumption 1 is violated.To illustrate the remainder of the development, the special case, e2 =-e1 , e3 =0 isconsidered. The corotational components of the velocity strain are then given by thecounterpart of (8).°εxttbx +e 1 h,x g°t=εy- e 1 h,x gtby°2εxyt- e 1 h,y gttby +e 1 h,y gtbxvx= Bv(1.5.34)vyFor an anisotropic material, the stress rate is then given bys=so+s1°o= Ce(C11 -C12 )( qxh,x-qyh,y)+ e1 (C -C )( q h, -q h, ) 22 21 0 y y x x (1.5.35)It can be seen that the corotational stress rate always has the same distribution within theelement, so the stress also has the same form; at any point in time, s o is the constant part ofthe element stress field evaluated at the quadrature point, and s 1 is the nonconstant part.Taking advantage of this form of the stress field, and inserting (34) and (35) into(26), and taking advantage of the orthogonality properties of hx and hy (1.4.16) and thefact that C is constant in the element, givesotfint = AB s o + fstabwhere fstab are the hourglass (stabilization) nodal forces, which are given byfstab =Qx gQy gwhere(1.5.36)(1.5.37)59Qx= e 21 C11 -C12 -C21 +C22Hxx qx - Hxy qy(1.5.38)Hyy qy - Hxy qxQyoand B is the constant part of B which is given by btxoB = 0 bty0btybtx(1.5.39a)11btx = 2A[y24 , y31 , y42 , y13 ]bty = 2A[x42 , x13 , x24 , x31 ]The nodal force vector is arranged by components:f t = [f tx , f ty ]f tx = [fx1 , fx2 , fx3 , fx4 ] f ty = [fy1 , fy2 , fy3 , fy4 ](1.5.39b)(1.5.40)For an isotropic material, Qx and Qy can be written in terms of the constants given inTable 3.Qx= 2µc1 Hxx +c2 Hyy qx + c3 Hxy qy(1.5.41)cH+cHq+cHqQy1 yy 2 xx y3 xy xAs with mixed method stabilization, the shear modulus in a nonlinear isotropic process isgiven by (1.4.9a).Table 4 is a flowchart outlining the procedure to evaluate nodal forces in an explicitprogram with time step ∆t.
Implementation in a static program simply requires thereplacement of the products of rates and ∆t by an increment in the corresponding integral;for example, ∆s replaces s ∆t.Table 4. Element nodal force calculation1. update corotational coordinate system2. transform nodal velocities v and coordinates x to corotational coordinate systemoo3. compute strain-rate at quadrature the point by e=B v (Eq. (39) gives B )4.
compute stress-rate by constitutive law and update stress (note: s =∆s /∆t)5. compute generalized hourglass strain rates by Eq. (32)6. compute the generalized hourglass stresses rates by (38) and update thegeneralized hourglass stresses7. compute fint by Eqs. (36) and (37)8. transform fint to global system and assemble60Remark 3.1 The stress rate in (36) corresponds to the Green-Naghdi rate if the corotationalcoordinate system is rotated by ω∆t in each time step.Remark 3.2 If the Jaumann rate is used in conjunction with a fixed coordinate system, thestress field loses the form of (35) and other approximations are needed.Remark 3.3 Because of the assumption of a constant spin and material response in theelement, deviations from this assumption are directly proportional to the strength of thehourglass modes (see for example (31-33)); thus in h-adaptive methods, it isadvantageous to refine by fission those elements which exhibit substantial hourglassenergy, as advocated in Belytschko, Wong, and Plaskacz (1989).Remark 3.4 If a Jaumann rate is used in a fixed coordinate system, the stress field does notmaintain the distribution (35) This is one reason that the corotational form is preferred.1.5.6 Assumed Strain with Multiple Integration Points.
In the development above,stabilization forces are obtained for a reduced one-point integration element. One-pointintegration was chosen because it is usually advantageous to keep the number of stressevaluations to a minimum; however, there is a correlation between the number ofintegration points needed in a mesh and the nonlinearity of the stress field. An example ofthis is the dynamic cantilever beam of Section 1.6.3.
For elastic material, a very accuratesolution can be obtained with only one element through the depth of the beam, because theaxial stress varies linearly through the depth. For elastic-plastic material, many elementsare need through the depth to obtain a reasonably accurate solution, because the axial stressvaries nonlinearly through the depth. The number of integration points can be increasedby refining the mesh, or by increasing the number of integration points in each element.The latter method has the advantage of being able to increase the number of quadraturepoints without reducing the stable time step of an explicit method.The assumed strain fields developed above can be used with any number ofintegration points without encountering locking since the strain fields have zero dilatationalstrain throughout the element domain for incompressible material.
The element forcevector for multi-point integration using an assumed strain field is analogous to (1.2.31b)and is given byintfe =nQ∑α=1twα J xα B xα s xα(1.5.42)where B xα and s xα are the corotational counterparts of (8) and (25) evaluated at aquadrature point, xα. Stabilization forces, may or may not be necessary with (42)depending on the location of the integration points.The g terms in (8) assure rank sufficiency, as long as h, x and h, y are not too small.If we consider the rectangular element in Fig. 12 with a corotational coordinate system, thereferential axes are parallel to the corotational axes, soξ, x = 1η, y = 1ab11h, x = a η h, y = b ξη, x = ξ, y = 0(1.5.43b)(1.5.43a)61From (43b), it is apparent that h, x = 0 along the η axis and h, y = 0 along the ξ axis.Therefore if the integration points are all located on one of the referential axes, stabilizationforces will be needed in either the x or y directions to maintain rank sufficiency.ηyξ2bx2aFigure 12.
A rectangular element in the corotational coordinate systemFull 2x2 integration using Eq. (42) is rank sufficient, but nearly the same results areobtained with two integration points using a modified form of Eq. (42) given byintfe = 2 J 02∑α=1tB xα s xα(1.5.44)In Eq. (44), J(0) is the Jacobian evaluated at the origin of the referential coordinatesystem, and the two integration points are either x 1 = (-1 3, -1 3), x 2 = (+1 3,+1 3), or else x 1 = (-1 3, +1 3), x 2 = (+1 3, -1 3). The choice of the pair ofintegration points makes little difference in the solution.
This 2-point integration scheme issimilar to the IPS2 element reported in Liu et al. (1988). The formulation here differs byusing an assumed strain field is used to improve accuracy. Using, the QBI strain field, aflexural-superconvergent 2-point element is obtained.In Section 1.6.3, we observe that the ASQBI element with 1-point integrationprovides an accurate coarse mesh solution with elastic material; however, with elasticplastic material, the coarse mesh solution is poor. We can therefore attribute the error inthe elastic-plastic solution to an insufficient number of integration points.
This large erroris not surprising if we consider the nature of the solution. The plastic deformation of abeam in bending initiates at the top and bottom surfaces of a beam where the axial stress isgreatest. With 1-point integration, the only stress evaluation is at the center of the element,so while the stress state at the integration point remains within the yield surface, the stressstate may be outside the yield surface at other points in the element domain. For coarsemesh bending, the error is large, resulting in too little plastic deformation.The 2-point integration scheme of Eq.
(44), and 2x2 integration by Eq. (42) improveon 1-point integration by placing integration points nearer the edge of the element. InSection 1.6.3, the effect of multiple stress evaluations is demonstrated by the solution ofan elastic-plastic cantilever beam. Results for the 2 and 4-point integration schemes aregiven in Tables 11a through 11d. Both use the QBI strain field, so the 2 point scheme iscalled ASQBI(2 pt), and 2x2 integration is called ASQBI(2x2).
Both of these elementshave flexural-superconvergence as does the 1-point element with ASQBI stabilization, sothe elastic part of the solution is solved very accurately. Therefore, the difference in thesolutions of these three elements with elastic-plastic material can be attributed to the effectof multiple stress evaluations on the nonlinear part of the solution.621.6 Numerical ResultsThe numerical examples reported here include linear and nonlinear problems.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
