Belytschko T. - Introduction (779635), страница 18
Текст из файла (страница 18)
Belytschko, Chapter 2, December 16, 1998forces for the Eulerian formulation are obtained in the following on an element level. Wedefine the inertial nodal forces by Eq. (2.7.8) which givesδ Peiner = δv Te feinert = δv Te ∫ ρN (N˙v +v, xv ) AdxT(2.11.14)ΩeFrom the above, it follows that the inertial nodal forces are given byfeiner = M ev˙ e + fetran(2.11.15)whereMe =∫TρN NAdx, fetran =Ωe∫ ρv, x vAdx(2.11.16)ΩeThe transport nodal forces have not been written in matrix form; they are quadratic in thenodal velocities. This term is needed in the Eulerian formulation because the nodes arefixed in space, so the time derivatives of the nodal velocities correspond to spatialderivatives.
The mass matrix differs from the mass matrix in the Lagrangian meshes in thatit is a function of time: as the density in the element changes, the mass matrix will changecorrespondingly.Example. Two-Node Eulerian Finite Element. The finite element equations aredeveloped for a one-dimensional, two node element with linear velocity, density and stressfields. The element, shown in Fig. 2.8, is of length l = x2 − x1 and unit cross-sectionalarea.
As can be seen, the spatial configuration does not change with time since it is anEulerian element. The map between element and spatial coordinates is given by x1 x (ξ ) = [1− ξ ξ ] ≡ N(ξ )x e x2 e(2.11.17)The density, velocity and stress are also interpolated by the same linear shape functionsρ(ξ ) = N(ξ)ρev (ξ ) = N (ξ )veσ (ξ ) = N (ξ)σe(2.11.18)Superscripts are not appended to the shape functions because all variables are interpolatedby the same shape functions.Density Equation.
The element matrices for the discrete continuity equation are given byM ρex21− ξ l 2= ∫ N N dx = ∫ [1−ξ ξ ]ldξ = 16 ξ x01T1g ρex2= ∫Nx11− ξ (ρv), x l dξ ξ 012(2.11.19)1T(ρv ), x dx = ∫ 2-59(2.11.20)T. Belytschko, Chapter 2, December 16, 1998The g ρe vector is usually evaluated by numerical quadrature. For linear interpolants it isgiven byg eρ =2v + v 1 2ρ + ρ 1(ρ2 − ρ1) 1 2 + ( v2 − v1 ) 1 2 6v1 + 2v2 6 ρ1 +2ρ 2 (2.11.21)The above matrix vanishes when the density and velocity are constant in the element.Stress Equation.
The element matrix for the stresses M σe = Meρ . The vector gσe is givenbyx211− ξ gσe = ∫ N T vσ , x dx = ∫ (v1 (1− ξ ) +v2 ξ )(σ2 −σ 1 ) dξξx10 2v + v 1= (σ 2 − σ1 ) 1 2 6v1 +2v2 (2.11.22)In summary, the finite element equations for the Eulerian formulation consists of threesets of discrete equations: the continuity equation, the constitutive equation, and themomentum equation, or equation of motion. The momentum equation is similar to theupdated Lagrangian form, except that the inertial term includes a transport term and varieswith time.
All nodal forces are defined over fixed intervals in space. The semidiscreteforms of the continuity and constitutive equations are first order ordinary differentialequations. We have only developed the discrete equations for the case where the endpointsare fixed.2 . 1 2 Solution MethodsWe have seen so far that the momentum equation can be discretized with a Lagrangianmesh in the form˙˙ = f ext – f int = fMu(2.12.1)These are ordinary differential equations in time.In order to enable some nonlinear problems to be solved at this point, we now describethe simplest solution method, explicit time integration of the equations of motion for aLagrangian mesh.
The most widely used explicit method is based on the central differenceformulas. Explicit integration can be simplified further by replacing M by a diagonal orlumped mass matrix.We start at time t=0 using time steps ∆t , so that at time step n, t = n∆t . The value ofa function at n∆t is denoted by a superscript n, i.e., un ≡ u( n∆t ) . In the central differencemethod, the velocities are approximated byn+1/ 2n –1/ 2u(t +∆t/ 2)−u(t− ∆t/ 2)u˙ n = v n+1/2 = u ∆– ut=∆t2-60(2.12.2)T. Belytschko, Chapter 2, December 16, 1998where the second equality is included to clarify the notation.
Half time step values are usedfor the velocities. The accelerations are given byn+1/2 –v n–1/2˙u˙ n ≡ a n = v∆t(2.12.3)In each case, the value of the derivative at the center of a time interval is obtained fromthe difference of the function values at the ends of the interval, hence the name centraldifference formulas. The flow chart for an explicit program is then given by the followingBox.Box 2.5 Flowchart for Explicit Time Integration of Lagrangian Mesh1. Initial conditions and initialization: set v0 , σ 0 ; n = 0,t = 0; compute M2.
get f n3. compute accelerations a n = M −1f n 1 if n = 0114. update nodal velocities: v n + 2 = vn + 2 −α + α∆tan : α = 21 if n > 05. enforce essential boundary conditions: if node I on Γv : vnI = v( x I ,tn )6. update nodal displacements: un +1 = un +∆tv n + 27. update counter and time: n ← n +1, t ←t +∆t8. output, if simulation not complete, go to 21Module: get f1. GATHER element nodal displacements une and velocities v ne +1/ 22.
compute measure of deformation3. compute stress by constitutive equation σ en4. compute internal nodal forces by equation in Box.5. compute external nodal forces on element and fe = feext − feint5. SCATTER element nodal displacements to global matricesUpdating for the displacements by Eq. (6) then does not require any solution ofalgebraic equations. Thus, in a sense, explicit integration is simpler than static linear stressanalysis. As can be seen from the flowchart, most of the explicit program is astraightforward interpretation of the governing equations and the time integration formulas.The program begins with the enforcement of the initial conditions; procedures for fittingdifferent initial conditions have already been described.
The first time step is somewhatdifferent from the others because only a half-step is taken. This enables the program tocorrectly account for the initial conditions on the stresses and velocities.Most of the programming and computation time is in computing the element nodalforces, particularly the internal nodal forces. The nodal forces are computed element-byelement. Prior to starting the element computations, the element nodal velocities anddisplacements are gathered from the global arrays. As can be seen from the flowchart, thecomputation of the internal nodal forces involves the application of the equations which areleft in strong form, the strain equation and the constitutive equation, followed by theevaluation of the internal nodal forces from the stress, which emanates from the weak form2-61T.
Belytschko, Chapter 2, December 16, 1998of the momentum equation. When the computation of the element nodal forces iscompleted, they are scattered to the global array according to their node numbers.The essential boundary conditions are enforced quite easily as shown. By setting thenodal velocities equal to the prescribed nodal velocities at all nodes on prescribed velocityboundaries, the correct displacements result, since the velocities are subsequently integratedin time.
The placement of this step in the flowchart insures that the correct velocities areavailable in the nodal force computation. The initial velocities must be compatible with theboundary conditions; this is not checked in this flowchart but would be checked in aproduction program. The reaction forces can be obtained by outputting the total nodalforces at the prescribed velocity nodes.It can be seen from the flowchart that the traction boundary conditions enter onlythrough the external nodal forces. Therefore, for a traction-free boundary, nothing need bedone: the homogeneous traction boundary condition is enforced naturally in a weak senseby the finite element solution.
However, the traction boundary condition is only satisfiedapproximately.Stability Criterion. The disadvantage of explicit integration is that the time step mustbe below a critical value or the solution "blows up" due to a numerical instability. This isdescribed in detail in Chapter 6. Here we limit ourselves to pointing out that the criticaltime step for the 2-node elements described in this Chapter is given by∆ t crit = cl(2.12.4)where l is the current length of the element and c is the wave speed given byc 2 = Eση / ρ(2.12.5)where Eση is the modulus in the elastic relation between Cauchy stress and rate-ofdeformation. For nonlinear materials, Eση is replaced by the current tangent modulustσ˙EσD=Dx(2.12.6)A computer program for the explicit integration of the one-dimensional updated andtotal Lagrangian formulation is given in Appendix B.Appendix A. Derivation of Conservation Equations in 1D.In this Appendix, simple derivations of the conservation equations in one dimensionwill be given.
These are "engineering" derivations which develop these equations insimplest terms, and they lack the mathematical rigor and generality associated with thederivations found in texts on continuum mechanics.We first derive the equation of conservation of mass, also called the continuityequation.
Consider a segment of the rod shown in Fig. 2.7, which in the initial, orundeformed, state is of length ∆X, cross-sectional area A 0 and density ρ 0. In the deformedstate this Lagrangian segment has length ∆x, area A and density ρ. The subdomain is amaterial, or Lagrangian subdomain, in that all material points remain in the subdomain and2-62T. Belytschko, Chapter 2, December 16, 1998the right and left end points are the same material points.
No flow of material occurs in orout of the segment. Therefore, by mass conservation, the mass in the undeformed segmentmust equal the mass of the deformed segment:ρA∆x = ρ0 A0 ∆X(A.1)tcurrentconfiguration( Aσ )x∆ xA ρb( Aσ )x+ ∆ x∆xreferenceconfigurationx, X( Ao P) X∆XAo ρo b( Ao P)X +∆ X∆XFig. 2.9. A segment of a rod in the reference (initial, undeformed) and current configurations showing allforces acting on the segment.Dividing by ∆X and taking the limit as ∆X → 0 givesρ0 A0 = ρA ∂x = ρAF or ρ0 ( X ) A0 ( X) = ρ( X,t ) A( X,t )F ( X,t )∂X(A.2)The above is one form of the equation of mass conservation.
On the right hand side theindependent variables have been indicated to stress that this equation only holds whenexpressed in terms of the material coordinates; the variables on the RHS are functions oftime, whereas the variables on the LHS are independent of time.To obtain another form of this equation, we note that the volumes of the segment arerelated by the Jacobian by dΩ = JdΩ 0 . Since dΩ = Adx and dΩ 0 = A0 dX , it follows thatJ = AA F0(A.3)Substituting the above into (A.1) gives the another form of the mass conservation equation2-63T.
Belytschko, Chapter 2, December 16, 1998ρ0 ( X ) = ρ( X,t ) J ( X, t )(A.4)The above equation also applies in multi-dimensional problems.Momentum Equation. The momentum equation is derived by considering the segmentof the rod shown in Fig. 2.9. The forces on the deformed segment are shown in Fig. 2.9and consist of the forces arising from the stress, which act on the right and left hand end ofthe segment, and the body force; the body force is distributed over the entire segment andits net resultant is placed at the center of the segment.
On the left hand end of the segmentthe force is ( Aσ ) x , i.e. the product of the stress and the current area at the point x . Theforce due to the stress on the right hand end is given by ( Aσ ) x + ∆x . The resultant force dueto the body force is obtained by multiplying b( x + ∆x2 ) by the mass of the segment ρA∆x .So if we write Newton's second law for the segment we have–Aσ x + ( Aσ ) x + ∆x + (ρAb) x + ∆x ∆x = ( ρA˙u˙ ) x + ∆x ∆x2(A.5)2where the LHS is the sum of the resultant forces from the stress and the body force and theRHS is the product of the mass of the segment and its acceleration.The forces due to the stresses are now expanded by a Taylor's series about themidpoint of the segment, with the product Aσ treated as a single function, which gives( Aσ ) x + ∆x = ( Aσ) x + ∆ x / 2 +( Aσ ) x = ( Aσ ) x +∆ x / 2 −∂ ( Aσ )∆x+O ∆x 2∂x x +∆ x / 2 2∂ (Aσ )∂x( )x +∆ x / 2(∆x+ O ∆x22(A.6a))(A.6b)The use of a Taylor series expansion of course presupposes that the function is smoothenough so that the first derivative exists; this is not the case wherever the stress or the areais discontinuous.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
