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

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

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

(7.16.16b), both the stress gradient, which will be denoted by γ , and its spatialderivatives are required. To obtain the gradient in Eq.(7.16.16b), a classical least-squaresproject is employed to obtain a smoothed field of stress gradient. Via divergence theorem,we have∫Ω NγγdΩ = − ∫ σ∇N γ dΩ + ∫ N γ σndΓΩ(7.16.17)Γwhere n is the outward unit normal in the current configuration. After the regularassembling procedure, we obtain the linear set of equations:Mγ γ = Φ(7.16.18)where M γ is a consistent pseudo-mass matrix defined as:M γ = [M γIJ ] = ∫ NIγ N Jγ dΩ(7.16.19a)Ωγ is the vector of nodal smoothed values of the stress gradient, and the vector Φ is definedas:[Φ = [Φ I ] = ∑ − ∫ σ∇N γI dΩ + ∫ NγI σndΓeΩΓ](7.16.19b)W.K.Liu, Chapter 738To make the algorithm explicit, the lumped mass matrix is preferred instead of consistentone.

After doing so, the solution of stress gradient of γ can be achieved straightforward.Then Eq.(7.16.16) can be solved. To obtain the stress value at Gauss points, thecollocation technique can be applied to handle the weak form of Eq.(7.16.16).(2) Godunov-like update:In this phase, the convection equation of :σ ,t[ χ] + ci σ ,i = 0(7.16.20)will be solved. With the help of the stress-velocity product Y = σc , Eq.(7.16.20) can berewritten as:σ ,t[ χ] + Yi,i = σci,i(7.16.21)To apply Godunov method, the weak form is presented∫Ω δσσ ,t[χ ]dΩ = ∫Ω δσσ ci,idΩ − ∫ΓtδσYi ni dΓ(7.16.22)where the test functions δσ are constant within any single element.

Since both σ and δσare constants within an element e , Eq.(7.16.22) results inσ ,t[ χ] = −1 Ns∑ f (σ c − σ )[1−sign( f s )]2Ω s =1 s s(7.16.23)with the upwind consideration, where the element e has volume Ω and Ns faces, σ sc is thestress component in the contiguous element across face s, and f s is the flux of convectivevelocity c across face s,f s = ∫ ci ni ds .s(7.16.24)To apply the Godunov update to the situation of multi-point quadrature, we can divideevery finite element into various subelements, each of them corresponding to the influencedomain of a Gauss point. In every subelement, σ is assumed to be constant, andrepresented by the Gauss-point value. Because of this , σ is a piece wise constant filedwith respect to the mesh of subelements, and Eq.(7.16.24) can be integrated to update thevalue of σ for each subelement.7.16.3 Stress Update Procedures in 1-D Case:In this section, we will compare the performance of different update procedure.

Below, wewill emphasize on the stress update. For illustrative purposes, we consider onedimensional(1-D) case. In a 1-D case, the shape functions and the corresponding testfunctions for density, velocity, energy and stress-velocity product may be chosen to be thepiecewise linear C 0 functions such that1N1 = 2 (1− ξ)1N2 = 2 (1+ ξ )(7.16.25a)(7.16.25b)W.K.Liu, Chapter 739where ξ ∈[−1,1], while the functions for deviatoric stress and pressure can be C −1 , or inparticular, constant in each element. The test and trial functions for all variables areidentical. The full upwind method can be applied for all the matrices involving convectioneffects.For a uniform mesh when the 1-D rod is divided into M segments of equal size of h, wherethe elements and the nodes are numbered sequentially from 1 to M and M+1, respectively.Let cm designate the convective velocity at node m and sm the stress in element m.

Forsimplicity, all the nodal convective velocities are considered to be positive.7.16.3.1 Application of SUPG in 1-D Case:The stress update is according to Eq.(7.16.7e). The matrices and vectors appeared inEq.(7.16.7e) are as the following.The generalized mass matrix is:1 Os1M = hO 1 M × M(7.16.26a)The transpose of the divergence operator matrix is −1 1O OT−1 1G =O O −1 1 M× ( M+1)(7.16.26b)The matrix M y is diagonal by locating the integration points at the nodes11diag(M y ) = h{2 ,1,L,1,L1, 2 }T( M+1)×1(7.16.26c)The transport of stress vector isLys = 16 h{(2c1 + c2 )s1, L,(c m−1 + 2c m )sm −1 + (2c m + cm+1 )sm ,L,(c M + 2cM +1 )sM }T( M+1)×1(7.16.26d)if exact integration is used;or1Lys = 2 h{0,L,(cm + c m+1 )sm ,L, (c M + c M+1 )sM }T( M+1)×1if full upwind is used, where 1 < m < M hereafter.(7.16.26e)W.K.Liu, Chapter 740The generalized diffusion vector isDs = {(−c1 + c2 )s1 ,L,(−cm + cm+1 )sm ,L,(−cM + c M+1 )sM }TM×1(7.16.26f)and the rate of change of stress due to material deformation ( the rotation of stress vanishesin 1-D case) isz = h{s˙1 ,L, s˙m ,L, s˙M }TM×1(7.16.26g)where s˙ = C1111v(1,1)By substituting Eqs.(7.16.12a-g) into Eq.(7.16.7e), the rate of change of stress in ALEdescription can be shown to be:(−3c1 )s1 + (2c2 + c3 )s2   s1,[x]  s1,[χ] M (− c1 + c2 )s1 −(c+2c)sm−1mm−1  M  M M −(c−c)s11mm+1 msst,[χ]=  sm,[χ] = h  (−cm + cm+1 )sm  − 6h +  m,[x] +(2cm+1 + cm+ 2 )sm +1MM  M M(−c + csM,[χ]  −(cM −1 + 2c M )sM −1  sM,[x ] MM +1 )s M M ×1 +3c sM +1 M(7.16.27a)if exact integration is used;or s1,[χ]  M st,[χ]=  sm,[χ]  M sM,[χ]  (−c1 + c2 )s1  s1,[x]  (− c1 + c2 )s1 M M M 1  (−cm−1 + c m )sm −1  1= h  (−cm + cm+1 )sm  − 2h  +(cm + cm +1 )sm  +  sm,[x]  (7.16.27b)MM  M −(c+c)sM−1MM−1(−c + c sM,[x] MM +1 )s M M ×1+2(c M + c M+1 )sM if full upwind is used to evaluate Lys .The second bracket on the right-hand side of Eq.(7.16.27a) shows the central differencing(or simple averaging) effects for the transport of stresses, while Eq.(7.16.27b) exhibits thedonor-cell differencing.

This can be further clarified by lettingc1 =L= cm =L= cM +1 = c(constant) ,thenW.K.Liu, Chapter 7 s1,[χ]  M st,[χ]=  sm,[χ]  M sM,[χ] 41 −s1 + s2   s1,[x] M  M  c= − 2h −sm−1 + sm+1  +  sm,[x] M  M  −s M −1 + s M  s M,[x] M ×1(7.16.28a)if exact integration is used; and s1,[χ]  M st,[χ]=  sm,[χ]  M sM,[χ] M ×1  s1,[x]   M   +  sm,[x]   M  M  s M,[x] s1Mc= − h  −sm−1 + smM−s+ 2sM −1(7.16.28b)if full upwind is used.Eq.(7.16.28a) shows that the transport of the stresses at odd and even elements tends to bedecoupled, therefore physically unrealistic oscillations would be expected when the simpleaveraging method is employed to evaluate the spatial derivatives of stresses.7.16.3.2.

Application of Operator Split in 1D Case:Also, for illustrative purposes, a 1D rod is considered which is divided into M segments ofequal size of h and assume:c1 =L= cm =Lc M+1 = c(constant) > 0 ,(1) Lax-Wendroff update:We can see that this procedure of update will not work for constant stress and linear shapefunction since the RHS of Eq.(7.16.19) will be zero.

In addition, we can see that for LaxWendroff update, the shape function of s , N s , must be in the same order as N γ .Assuming both of them are linear shape functions, we can obtain:11diag(M y ) = h{2 ,1,L,1,L1, 2 }T( M+1)×1Φ=(7.16.29a)1[−s + s ,−s + s ,− s + s ,L, −sm−1 + sm+1 ,L− sM −1 + sM +1 ](TM +1)×12 1 2 1 3 2 4(7.16.29b)(2) Godunov-like update:For constant stress and linear shape function, it is easy to obtain the update equation forelement nW.K.Liu, Chapter 7s(nt+∆t ) = s(nt) −42c∆t ( t)(t)(s − sn−1) +∆t s˙n(t +∆t /2)h n(7.16.30)7.16.4 Explicit Time Integration Algorithm:For simplicity, the coupled equations (7.16.7c-e) will be integrated by an explicit scheme.Lumped mass matrices are used to enhance the computational efficiency.

Both predictorcorrector method( Hughes & Liu, 1978) and standard central difference(Huerta & Casadei,1994) method can be applied here for explicit time integration. Below, ( )n and ( )n +1 willdenote the matrices at times tn = n∆t and tn +1 = (n + 1)∆t respectively, where ∆t is thetime increment.(1) Predictor-corrector methodThis kind of predictor-corrector method is similar to the Newmark algorithm. The majordifference is that the former algorithm is explicit, while the latter is implicit.Mass equations:ρ,t[χ ]n+1 = − (Mρn )−1 (Lρnρ˜ n+1 + Kρn ρ˜ n +1 )(7.16.31a)ρ˜ n+1 = ρ n + (1− α )∆tρ ,t[χ ]n(7.16.31b)ρ n+1 = ρ˜ n+1 + α∆tρ,t[ χ]n +1(7.16.31c)Momentum equations:int˜ n+1 )v,t[χ ]n+1 = (Mn ) −1 (f extn+1 − f n − Ln v(7.16.32a)v˜ n+1 = v n + (1− γ )∆tv,t[ χ]n(7.16.32b)v n+1 = v˜ n+1 + γ∆tv,t[ χ]n +1(7.16.32c)Eq.(7.16.32a) needs to be used in conjunction with1d˜ n+1 = dn +∆ tv n + ( 2 − β )∆t 2 v,t[ χ]n(7.16.32d)dn+1 = d˜ n+1 + β∆t 2 v,t[ χ]n +1(7.16.32e)to calculate the f intn .In the above equations, α , β and γ are the computational parameters.

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

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