Belytschko T. - Introduction (779635), страница 80
Текст из файла (страница 80)
For explicitcalculations, the following constraints on the parameters are used:α =0(7.16.33a)β=0(7.16.33b)W.K.Liu, Chapter 74312γ ≥(7.16.33c)The flowchart of the computational procedure for the class of pressure-insensitive materialsis as follows:Step 1. Initialization. Set n=0, input initial conditions.Step 2.
Time stepping loop, t ∈[0,t max ].Step 3. Integrate the mesh velocity to obtain the mesh displacement and spatial coordinates.Step 4. Calculate incremental hydrostatic pressure by integration Eq.(7.16.7e) with s and zreplaced by p and u respectively:(a) the rate of pressure due to convection,(b) the rate of pressure due to deformation.Step 5. Calculate incremental deviatoric stresses, yield stresses, and back stresses byintegration Eq.(7.16.7e) which stress update procedures have been discussed in detail inlast section:(a) the rate of stresses due to convection,(b) the rate of stresses due to rotation,(c) the rate of stresses due to deformation.Step 6. Compute the internal force vector.Step 7. Compute the acceleration by the equations of motion, Eq.(7.16.32a).Step 8. Compute the density by the equation of mass conservation, Eq.(7.16.31a).Step.7.
Integrate the acceleration to obtain the velocity.Step 10. If (n + 1)∆t > t max , stop; otherwise, replace n by n+1 and go to Step 2(2) Central difference methodThe central difference can also be applied for the time integration. Same as the standardcentral difference procedure, we can obtain displacement and acceleration vectors at each i1time step, while get velocity vector at each i + 2 time step. The momentum equation will beintegrated as an example to illustrate this explicit scheme.1After obtaining the velocity vector at (n − 2 )∆t and the displacement vector at n∆t , theacceleration vector at n∆t is :intv,t[χ ]n = (Mn ) −1 (f extn − f n − Ln− 1 v n− 1 )22With central difference scheme, the velocity and displacement of next time step are:v n+ 1 = v n− 1 +∆tv,t[ χ]n22anddn+1 = dn +∆ tv n+ 12For operator splitting update, the Lagrangian part for strains is calculated with a usualcentral difference scheme as:W.K.Liu, Chapter 744( L)ε(t ∗ ) = ε n +∆ tε,t[x]n+12whereε (,t[L)x]n+12= (v,x )n+ 12After considering the Godunov scheme for the convect part, the full upwind integration for1D 2-node element can be written as:∆t 1 + sign(c1 )1 − sign(c 2 )ε n+1 = ε(t∗ ) − [c1(ε e − ε e −1 ) + c2(εe +1 − εe )]h22For elastic-plastic problem, the radial return method can be used to determine the correctstates of stress and strain.7.10 ALE Governing Equation7.10.1 Slightly Compressible Viscous Flow with Moving Boundary ProblemIn this section, we develop the governing equations for a slightly compressibleNewtonian fluid.
In a generalized Newtonian fluid, the stress is a function of the rate-ofdeformation. Therefore, the stress is independent of the history of deformation. The mostwell known case is a linear Newtonian fluid, where the stress is a linear function of therate-of-deformation.This class of materials simplifies the implementation since the constitutive equation isindependent of the strain history.
Therefore, the constitutive equation can be used in itsstrong form.The formulation is restricted to isothermal processes; therefore energy equation is notneeded. The continuity equation (7.7.15), momentum equation (7.7.25) and the meshupdate equation (7.2.19) in ALE description are:∂ρ∂ρ∂v+ ci+ρ i = 0∂t χ∂xi∂x iρ∂vi∂tχ+ ρc j∂vi ∂σ ij=+ ρbi∂x j ∂x j∂x iw = ci∂χ j j(7.10.1a)(7.10.1b)(7.10.1c)In the constitutive equation, the hydrostatic stress is independent of deformation while theshear stress, which depend on the rate of deformation.σ ij = − pδij + 2µDij(7.10.1d)W.K.Liu, Chapter 745whereDij =1 ∂vi ∂v j and µ = µ (Dij )+2 ∂x j ∂xi (7.10.1e)where µ is the dynamic viscosity which is shear rate dependent.
The functions µ forgeneralized Newtonian models are presented in Table 7.2.W.K.Liu, Chapter 746Model1-D Viscosity3-D GeneralizationNewtonianµ 0 = constants = 2µ0 DPower Lawµ = mDn −1TruncatedPower LawCarreau-ABingham2tr(D)2if D ≤ D0µ = µ0if2tr(D)2 ≤ D0 :s = 2µ0 Dif D ≥ D0if2tr(D)2 ≥ D0 :µ = µ0 (Carreau[s = 2m 2tr(D)2 s = 2µ0 γ˙ 0D n−1)D0[µ − µ∞= 1+ (λ D)2µ0 − µ∞]n −1]Dn−1Ds = 2µ∞ D(n−1)/2[+2(µ 0 − µ ∞ ) 1 + 2λ 2 tr(D)2[µ∞ = 0s = 2µ0 1 + 2λ 2 tr(D)2]( n−1)/2]( n−1)/2D1tr(s)2 ≤ τ 022s =01if tr(s)2 ≥ τ 022τ0s = 2µ p 1 +D2tr(D)2 ifµ = ∞ τ ≤ τ0µ = µp +τ0Dτ ≥ τ0where s and D are the deviatoric part of the stress and stretch tensors, respectivelyTable 7.2DW.K.Liu, Chapter 747We will first consider a viscous, barotropic fluid, so that the pressure depends onlyon the density, p = p(ρ) . The material time derivative of the pressure p gives:p˙ =∂p ˙ρ∂ρ(7.10.2a)Now we can define the bulk modules B by:B ∂p=ρ ∂ρ(7.10.2b)Using this definition in Eq.
(7.10.2a) yields:p˙ =B˙ρ˙ p˙ρ or=ρρ B(7.10.2c)Substituting Eq. (7.10.2c) into Eq. (7.10.1a), the continuity equation may be rewritten as(Liu and Ma, 1982):1 Dp ∂vi+=0B Dt ∂xi(7.10.3a)or, by introducing Eq. (7.2.20b) in (7.10.3a), as:1 ∂p1 ∂p ∂vi+ ci+=0B ∂t χ B ∂xi ∂xi(7.10.3b)The objective here is to find the density, the material velocity and the mesh velocityby solving Eqs. (7.10.1a-e).
Prior to presenting the strong from of the governingequations and boundary conditions, the finite element mesh updating procedure will first bediscussed in the next section.In each Eq (7.10.1) a convective term is present; thus, one of the drawbacks of anEulerian formulation is retained in the ALE methods. The major advantage of an ALEapproach is that it simplifies the treatment moving boundaries and interfaces.7.11 Mesh Update Equations7.11.1 IntroductionThe option of arbitrarily moving the mesh in the ALE description offers interestingpossibilities.
By means of ALE, moving boundaries (which are material surfaces) can betracked with the accuracy characteristic of Lagrangian methods and the interior mesh can bemoved so as to avoid excessive element distortion and entanglement. However, thisrequires that an effective algorithm for updating the mesh, i.e. the mesh velocities vˆ , mustW.K.Liu, Chapter 748be prescribed. The mesh should be prescribed so that mesh distortion is avoided and sothat boundaries and interfaces remain at least partially Lagrangian.In this section, we will describe several procedures for updating the mesh. Thematerial and mesh velocities are related by Eq.
(7.2.19); hence, once one of them isdetermined, the other is automatically fixed. It is important to note that, if vˆ is given, dˆand aˆ can be computed and there is no need to evaluate w. On the other hand, if vˆ isconsidered the unknown but w is given, Eq. (7.2.19) must be solved to evaluate vˆ beforeupdating the mesh. Finally, mixed reference velocities can be given (i.e. a component of vˆcan be prescribed and w in the other(s)). Finding the best choice for these velocities and analgorithm for updating the mesh constitutes one of the major hurdles in developing aneffective implementation the ALE description.
Depending on which velocity ( vˆ or w ormixed) is prescribed, three different cases may be studied.7.11.2 Mesh Motion Prescribed a PrioriThe case where the mesh motion vˆ is given corresponds to an analysis where thedomain boundaries are known at every instant. When the boundaries of the fluid domainhave a known motion, the mesh movement along this boundary can be prescribed a priori.7.11.3 Lagrange-Euler Matrix MethodThe case where the relative velocity w is arbitrarily defined is a format for apropos byHughes et al (1981). Let w be:wi =∂χ i∂t()= δij − α ij v jX(7.11.1)where δ ij is the Kroneker delta and α ij is the Lagrange-Euler parameter matrix such thatα ij = 0 if i ≠ j and α ii is real (underlined indices meaning no sum on them).
In general,the α' s can vary in space and be time-dependent; however α ij is usually taken as timeindependent. According to Eq. (7.11.1) the relative velocity w becomes a linear functionof the material velocity and it was chosen because, if α ij = δ ij , w = 0 and the Lagrangiandescription is obtained, whereas, if α ij = 0 , w = v and the Eulerian formulation is used.The Lagrange-Euler matrix needs to be given once and for all at each grid point.Since w is defined by Eq. (7.11.1), the other velocities are determined by Eq.(7.2.19), which become, respectively.ci =∂xiδ − α jk vk∂χ j jk()()(7.11.2a)andvˆi = vi − δ jk − α jk vk∂xi∂χ j(7.11.2b)W.K.Liu, Chapter 749The latter equations must be satisfied in the referential domain along its boundaries.Substituting Eq.
(7.2.8a) into (7.11.2b) yields a basic equation for mesh rezoning:∂xi∂tχ()+ δ jk − α jk vk∂xi− vi = 0∂χ j(7.11.3)The explicit form of Eq. (7.11.3) in 1D, 2D and 3D are listed:1D Form∂x∂x+ (1 − α ) v−v =0∂t χ∂χ(7.11.4)2D Form∂x1∂x∂x+ (1 − α11 ) v1 1 + (1 − α 22 ) v2 1 − v1 = 0∂t χ∂χ1∂χ 2(7.11.5a)∂x2∂t(7.11.5b)χ+ (1− α 11) v1∂x2∂x+ (1 − α 22 ) v2 2 − v2 = 0∂χ1∂χ 23D Form∂x1∂x∂x∂x+ (1 − α11 ) v1 1 + (1 − α 22 ) v2 1 + (1 − α 33 ) v3 1 − v1 = 0∂t χ∂χ1∂χ 2∂χ 3(7.11.6a)∂x2∂t∂x2∂x∂x+ (1 − α 22 ) v2 2 + (1− α 33 ) v3 2 − v2 = 0∂χ1∂χ 2∂χ 3(7.11.6b)∂x3∂x∂x∂x+ (1 − α11 ) v1 3 + (1− α 22 ) v2 3 + (1 − α 33 ) v3 3 − v3 = 0∂t χ∂χ1∂χ 2∂χ 3(7.11.6c)χ+ (1− α 11) v1Remark 1:Equation (7.11.3) differs only in its last term from the one proposed by Hughes et al(1981). This difference is not noticeable if the Lagrange-Euler parameters α ij are chosenequal to zero or one.
Moreover, Eq. (7.11.3) includes the Jacobian matrix (i.e. ∂xi / ∂χ j )that is missing in the Liu and Ma (1982) formulation. Finally, Eq. (7.11.3) is a transportequation without any diffusion so the classic numerical difficulties associated with transportequations are expected.W.K.Liu, Chapter 750The ALE technique with a mesh update based on the Lagrange-Euler parameters isvery useful in surface wave problems. We assume that the free surface is oriented relativeto the global coordinates so that it can be written as x3s = x3s (x1 , x2 ,t).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
