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

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

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

Note that the spatialderivatives of the deviatoric stress are involved in the convection term.When C−1 functions are used to interpolate the element stresses, the ambiguity of the stressderivatives at the element interface renders the calculation of the spatial derivatives of stressa difficult task. As mentioned in Remark 7.16.1.3, this is remedied by replacingEq.(7.16.1d) by a set of coupled equations, Eqs.(7.16.1h) and (7.16.1i), and thecorresponding matrix equations have been given in Eqs.(7.16.5a) and (7.16.5b).

Thestress-velocity product y can be eliminated by inverting M y in Eq(7.16.5b) andsubstituting into Eq(7.16.5a):M s s,t[ χ ] + GT (M y )−1 Lys − Ds = z(7.16.7a)where G T (M y )−1 Ly s can be identified as the convective term, and the upwind techniquesmentioned in Remark 7.16.1.4 should be applied to evaluate Lys . When c = 0 , i.e.,χ = x , Eq.(7.16.7a) degenerates to the usual stress updating formula in the Lagrangiandescription,M s s,t[ χ ] = z(7.16.7b)REMARK 7.16.2.1. The conservation equations are listed here to show the similarityamong the equations:M ρρ,t[χ ] + Lρ ρ + Kρ ρ = 0(7.16.7c)W.K.Liu, Chapter 733Mv ,t[ χ ] + Lv + f int = f ext(7.16.7d)M s s,t[ χ ] + GT (M y )−1 Lys − Ds = z(7.16.7e)REMARK 7.16.2.2.

In the above, Eq.(7.16.5b) is eliminated because it can beincorporated in Eq.(7.16.7e). This procedure is analogous to that by Liu & Chang(1985)where a fluid-structure interaction algorithm is described.REMARK 7.16.2.3. All the path-dependent material properties, such as yield strains,effective plastic strains, yield stresses, and back stresses, should be convected viaEq.(7.16.7e) with s replaced by these properties, respectively, and with z appropriatelymodified.7.16.2.2 Stress Update Procedure for SUPG Method:In a nonlinear displacement finite element formulation, when applied to elastic-plasticmaterials with kinematic hardening, the velocities are stored at nodes while the stresshistories, back stresses, and yield radii are available only at quadrature points.

In order toestablish the nodal values for the stress-velocity product, a weak formulation is a logicalnecessity. In addition, based on the one-dimensional study(Liu et. al. 1986), in which theupwind procedure is used to define this intermediate variable, the artificial viscositytechnique(streamline upwind) is considered here as a generalization of this upwindprocedure to multi-dimensional cases.The relation for the stresses-velocity product of Eq.(7.16.1i), is modified to accommodatethe artificial viscosity tensor Aijkm byyijk = sij ck − Aijkm,m(7.16.8a)The ingredients of the artificial viscosity tensor consist of a tensorial coefficient multipliedby the stress:Aijkm = µkm sij(7.16.8b)where the tensorial coefficient is constructed to act only in the flow direction (streamlineupwind effect)µ km = µ ck cm cn cn(7.16.8c)and the scalar µ is given byµ =NSD∑ αicihiNSD(7.16.8d)i=1Here hi is the element length in the i-direction; NSD designates the number of spacedimensions; and α i is the artificial viscosity parameter given by1αi =  2 1− 2for ci > 0for ci < 0(7.16.8e)W.K.Liu, Chapter 734The weak form corresponding to Eq.

(7.16.8a) can be obtained by multiplying by testfunction for the stress-velocity product and integrating over the spatial domain Ω :∫Ω δyijk yijkdΩ = ∫Ω δyijksijck dΩ − ∫Ω δyijk Aijkm,m dΩ(7.16.9a)This equation may be written as∫Ω δyijk yijkdΩ = ∫Ω δyijksijck dΩ + ∫Ω δyijk,m AijkmdΩ(7.16.9b)by applying the divergence theorem and by assuming no traction associated with theartificial viscosity on the boundary. The expression for Aijkm , Eq. (7.16.8b and c), can besubstituted into this equation to yield∫Ω δyijk yijkdΩ = ∫Ω (δyijk + δyijk )sij ck dΩ(7.16.9c)δyijk = δyijk,m µ cm c nc n(7.16.9d)wherecan be viewed as a modification of the Galerkin finite element method because of thetransport nature of the stress-velocity product.The shape functions for the stress-velocity product can be chosen to be the standard C 0functions. The number and position of quadrature points for Eq.(7.16.9c and d) should beselected to be the Gauss quadrature points, since the stress histories in Eq.(7.16.9c) areonly available at these points.Following the procedures given above, the stress-velocity product can be defined at eachnodal point and it can be substituted into the constitutive equation of Eq.(7.16.5a) tocalculate the rate of change of stresses with the same procedure as that of without artificialviscosity.

Note that the interpolation functions for stresses are integrated over the spatialelement domain. The task of selecting the number of quadrature points for the displacementfinite element poses another important issue. For example, the locking phenomenon forfully integrated elements arises when the material becomes incompressible(Liu, et.al.,1985). While selective reduced integration can overcome this difficulty, it is just as costlyas full quadrature. To alleviate this computational hurdle, the use of one-point quadraturecombined with hourglass control is developed by Belytschko et.al. (1984). In addition, thenonlinear two-quadrature point element (Liu et.al., 1988) appears to be another candidatefor large-scale computations because it exhibits nearly the same accuracy as the selectivereduced integration element while with only one-third of the cost.

These two kind ofelements, as well as the other displacement elements, can be readily adopted in the ALEcomputations.Following the procedure described by Liu et.al(1986), the displacement element is dividedinto M sub-domain, where M denotes the number of quadrature points. Each sub-domainis designated by Ω I (I=1,...,M), which contains the quadrature points x I , and no twosub-domains overlap. Associated with Ω I , a stress interpolation function NIs is assignedand its value is prescribed to be unity only at the quadrature point x = x I such thatNIs (x = x I ) = 1(7.16.10a)W.K.Liu, Chapter 735The test function in Ω I is chosen to be the Dirac delta functionNIs = δ(x − x I )(7.16.10b)Substitution of these functions into the constitutive equation represents a mathematicalrequirement that the residual of the weal from vanishes at each collocative quadrature point.Because the collocation point is located right at the quadrature point, the algebraic equationsresulting from Eq.(7.16.5a) can be easily worked out without numerical integration andgiven as below:General mass matrix isMIJs = ∫ NIs N Js dΩ x = I M× M(7.16.11a)Ωwhere the subscripts I and J range from 1 to M, the number of stress quadrature points perelement.

The transpose of the divergence operator matrix reads,yGTIJ = ∫ N Is N J,xdΩΩFor 2D 4-node element, it will be: N1,x (ξ1 ) N1,x (ξ 2 )GT = M N (ξ ) 1,x MN3, x (ξ1 )N3,x (ξ2 )MN3,x (ξ M )N1,y (ξ1 )N1, y (ξ 2 )MN1,y (ξ M )N3,y (ξ1 )N3,y (ξ2 )MN3,y (ξ M )N2, x (ξ1 )N2, x (ξ2 )MN2, x (ξ M )N4,x (ξ1 )N4, x (ξ2 )MN4, x (ξ M )N2,y (ξ1 )N2, y (ξ 2 )MN2, y (ξ M )N4, x (ξ1 ) N4, y (ξ 2 ) MN4, y (ξ M )(7.16.11b)M ×8The generalized diffusion matrix for stress isDIJ =∫Ω NI ck,k NJ dΩssorc k,k (ξ1 )ck ,k (ξ 2 )D= Ock ,k (ξ M ) M × MThe generalized stress vector iszijI = ∫ NIs{CijklDkl + τ kj Wik + τ ki W jk }dΩΩor(7.16.11c)W.K.Liu, Chapter 7 (Cijkl Dkl + τ kj Wik + τ ki W jk ) ξ1  (Cijkl Dkl + τ kj Wik + τ ki W jk )ξ 2z=M (C D + τ W + τ W ) kj ikki jk ξ M  ijkl klM×136(7.16.11d)7.16.2.3.

Stress Update Procedure for Operator Split Method:In addition to the methods shown in the last section to solve the fully coupled equations,another approach referred to as an operator split is an alternative way to apply FEM tosolve this problem numerically(Benson, 1989). Conceptually, the approach is simple.“Splitting” stands for “decomposing” a set of PDE operators into several sets of simplePDE operators which will be solved sequentially. An operator split decouples the variousphysical phenomena in the governing equations to obtain simpler equations to be solvedmore easily.

In exchange of certain loss of accuracy, the operator split offers a genericadvantage: simpler equations leads to simpler and stable algorithms, specifically designedfor each decoupled equation according to the different physical characteristics.To illustrate the operator split concept, we consider the transport equation of onecomponent of the Cauchy stress by:σ ,t[ X] = σ ,t [ χ] +ciσ ,i = q(7.16.12)where the expression of q will be determined by the constitutive law.

An example of q isthe right hand of Eq.(7.16.1h).The operator split technique is to split Eq.(7.16.12) into 2 phases. The first phase of theoperator split is to solve a Lagrangian step without considering the convective effect asbelow:σ ,t[ χ] = q(7.16.13)In this Lagrangian phase, it is integrated in time to update stresses from σ (t) (stress at timet) to σ ( L) (denotation of the Lagrangian updated stress), neglecting the convective termswhich is equivalent to assuming that mesh points χ move with material particles X , that isχ = X. Thus this Lagrangian phase can proceed in the same way as the usual UpdatedLagrangian procedure. In addition, it is well known that the stresses are obtained at Gausspoints.The second phase is to deal with the convective term that has not been taken into accountduring the Lagrangian phase, where the governing PDE is :σ ,t[ χ] + ci σ ,i = 0(7.16.14)and during which phase, the stresses are updated from σ ( L) to σ (t+∆t) .According to the two phases strategy above, the constitutive equation is split into theparabolic equation of the Lagrangian phase and the hyperbolic equation of the convectionphase.Here, we follow [12] to apply explicit methods to integrate the convection equation.

Tocompute gradients of the stress fields on the surface of elements, two different approachesW.K.Liu, Chapter 737have been taken: i) use an explicit smoothing procedure (Lax-Wendroff update); or (ii) usean algorithm that circumvents the computation of the stress gradient (Godunov-liketechnique).(1) Lax-Wendroff update:The key point of the Lax-Wendroff method is to replace the time derivatives of dependingvariables with spatial derivatives using the governing equations. For the partial differentialequations of Eq.(7.16.14), we have:( L)(t +∆t /2) ( L )σ ,t[σ,iχ] = − ci(7.16.15a)and( L)(L)( t+∆t/2) ( L )( L)σ ,ttσ,i ),t = −ci(t +∆t /2) (σ ,t[[χ ] = (σ,t[χ ] ),t = −(ciχ] ),i( L)= −ci(t+∆t/2) (−c (tj +∆t /2)σ ,(jL) ),i = ci(t +∆t /2)c (tj +∆t /2)σ ,ij(7.16.15b)After substituting the above two equations into the Taylor series expansion of σ (t+∆t) withrespect to the time:1(L)( L)2σ (t+∆t) = σ t + σ,t[χ]∆t + 2 σ ,tt[χ ]∆t(7.16.16a)we obtain the update equation of the Lax-Wendroff method as:1σ (t+∆t) = σ t − c(i t+∆t /2) σi( L )∆t + 2 c(i t+∆t/2) c(j t+∆t/2) σ (,ijL )∆t 2(7.16.16b)where ci(t +∆t /2) is the convective velocity evaluated at the mid-step.In Eq.

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

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