Belytschko T. - Introduction (779635), страница 57
Текст из файла (страница 57)
Initial conditions must be given for thevelocitites, the Cauchy stresses, and all state variables of the materials in themodel. The displacements are initially considered to vanish.Flowchart inorrect, half missing on time steps, not n orderBox 6.1Flowchart for Explicit Time Integration1. Initial conditions and initialization:set v 0 , σ0 , and other material state variables;d0 = 0, n = 0, t = 0; compute M(2. getf f n , ∆tcrit)3.
compute accelerations a n = M −1f n4. compute kinetic energy and check energy balance, see Section ??5. update nodal velocities: v n + 2 = vn + 12 ∆tn a n16-4T. Belytschko & B. Moran, Solution Methods, December 16, 19986. enforce velocity boundary conditions:n+1if node I on Γvi : viI 2 = vi x I , tn + 1 2 n + 21 n + 217. update nodal displacements: dn +1 = dn +∆ t8. update counter and time: n ← n +1, t ←t +∆t9. update nodal velocities: v n +1 = vn+ 12v+ 1 ∆tn a n210. output, if simulation not complete, go to 2(Subroutine getf f n , ∆tcrit)0. initialization: f n = 0, ∆tcrit = ∞1.
compute external nodal forces f ext ,n which are global2. loop over elements ei. GATHER element nodal displacements and velocitiesii. feint,n = 0iii. loop over quadrature points ξQ1. if n=0, go to 8n − 212. compute measures of deformation: D( )(ξ Q ), Fn (ξQ ),E n( ξQ )3.
compute stress σ n ξQ by constitutive equation4. feint,n ← feint,n +B T σ n wQ JξQEND quadrature point loopiv. compute external nodal forces on element, feext , nv. fen = feext ,n − feint, neeevi. compute ∆tcrit, if ∆tcrit<∆ tcrit then ∆t crit = ∆tcritvii. SCATTER fen to global f n3.
END loop over elementsIn this algorithm, the accelerations are first integrated to obtain thevelocities. The integration of the velocities is broken into two half-steps so thatthe velocities are available at an integer step in the computation of the energybalance. The displacements are computed in each time step by integrating thevelocities.The main part of the procedure is the calculation of the nodal forces fromthe nodal displacements at a given time step, which is performed in getf. In thissubroutine, the equations governing a continuum are used along with thegather/scatter procedures:6-5T. Belytschko & B.
Moran, Solution Methods, December 16, 19981. the nodal displacements of the element are extracted from the globalmatrix of nodal displacements by the “gather “ operation;2. the strain measures are computed at each quadrature point of theelement;3. the stresses are computed by the constitutive equation at eachquadrature point;4.
the internal nodal forces are computed by integrating the product of theB matrix and the stresses over the domain of the element with theCauchy stress;5. the nodal forces of the element are scattered into the global array.In the first time step, the strain measures and the stress are not computed. Instead,as shown in the flowchart, the initial stresses are used to obtain the internal nodalforces.The flowchart shows the algorithm with the matrix form of the internalforce computation, in which the stress tensor is stored as a square matrix and theB matrix is used.
The change to the Voigt form only requires the use of a columnmatrix for the stresses and the B matrix, (4.5.14). Similarly, the internal forcecomputation can be changed to the total Lagrangian format by replacing thediscrete values of the integrand in step 10 by the integrands of (B4.8.2).Most essential boundary conditions are easily handled in explicit methods.For example, if the velocities or displacements are prescribed as functions of timealong any boundary, then the velocity/displacement boundary conditions can beenforced by setting the nodal velocities according to the data:(viIn = vi x I , tn)(6.2.7)If the data is not available on the nodes, the least square procedure given inSection 2.4.5 can be used to fit the nodal values.The velocity boundary conditions can also be enforced in local coordinatesystems as shown in the Box 6.1.
In that case, the equations of motion at thesenodes must be expressed in the local coordinate system, so the nodal forcecomponents must be expressed in the local coordinate systems before assemblyand time integration. The boundary condition is also enforced in the localcoordinate system. The orientation of the local coordinate system may vary withtime but the time integration formulas must then be modified to account for theadditional terms in the equations of motion.When essential boundary conditions are given as linear or nonlinearalgebraic equations relating the displacements, the implementation is morecomplicated.
One approach is to use a linearization of the constraint. Considerfor example the nonlinear constraintg(d(t )) = 0(6.2.8)where g(d(t )) is a linear or nonlinear algebraic function of the nodaldisplacements. If the constraint involves integral or differential relationships,such as a dependence on the velocities, it can be put in the above form by using6-6T. Belytschko & B. Moran, Solution Methods, December 16, 1998difference equations or a numerical approximation of the integral.
The above canbe linearized as follows:( )() vn +1/ 2 = 0 ∂G dn ∂G dn +1 ∂da + ∂da(6.2.9)aAfter a large number of time steps, linearizations such as the above combinedwith a central difference update of the displacements may substantially violate theconstraint. This drift in the enforcement of the constraint can be avoided bycorrecting the linearized update so that the constraint is enforced exactly at thenext time step, n+1. When accurate treatment of the constraints is important,techniques for differential-algebraic equations should be used, Petzold (??).As can be seen from the flowchart, an explicit method is easilyimplemented. Furthermore, explicit time integration is very robust, by which wemean that the explicit procedure seldom aborts due to failure of the numericalalgorithm.
The salient disadvantage of explicit integration, the price you pay forthe simplicity of the method and its avoidance of the solution of equations, is theconditional stability of explicit methods. If the time step exceeds a critical value∆tcrit , the solution may grow unboundedly and will in any case be erroneous.The critical time step is also called the stable time step.
The critical timestep for a model depends on the mesh and the material properties. For low orderelements, we will show in Section X that the critical time step for linear responseis given by∆tcrit = minlece(6.2.10)where l e is a characteristic length of element e and ce the wavespeed of elemente. Thus the critical time step decreases with mesh refinement and increasingstiffness of the material.
The cost of an explicit simulation is independent of thefrequency content which is of interest and depends only on the size of the modeland the time of the simulation relative to the critical time step given by (6.2.10).The time step is calculated in the flowchart on an element basis.
For eachelement, a critical time step is calculated, and if it is smaller than that calculatedfor all previous elements in that time step, it is reset. The theoretical justificationfor setting the critical time step on an element basis and other approaches aredescribed in Section 6.??.6.3EQUILIBRIUMINTEGRATION.SOLUTIONSANDIMPLICITTIME6.3.1.Equilibrium and Transient Problems. We will combine thedescription of the solution of the equilibrium equations with time integration byimplicit methods because they share many common features. To begin, we write6-7T. Belytschko & B.
Moran, Solution Methods, December 16, 1998the discrete momentum equation at time step n+1 in a form applicable to bothequilibrium and dynamic problems:()()()0 = r dn +1 , t n+1 = sD M˙d˙ n +1 − f n+1 = s DMan+1 − f ext dn +1,t n +1 + f int d n +1 (6.3.1)where sD is a switch which is set by:0 for a static(equilibrium)problemsD = 1 for a dynamic(transient) problem(6.3.2)The column matrix r( dn +1 , t n+1 ) is called a residual. When sD = 0 , the above arethe equilibrium equations at the next step. In addition, the displacement boundaryconditions must be met; these can be written as a set of nc nonlinear algebaricequations()G i dn +1 = gi , i = 1 to nc(6.3.2b)Differential and integral constraints are put in discrete form by usingdiscretizations of the derivatives and integrals, respectively.
In most cases thedisplacement boundary conditions are linear algebraic equations, but we havewritten the general form (6.3.2b) because complex boundary conditions are oftenneeded in nonlinear problems.When the accelerations vanish or are negligible, a system is in equilibriumand the solution of the resulting equations is called an equilibrium solution. Theequilibrium equations are given by (6.3.1) with sD = 0 :()()(0 = r dn +1 , t n+1 = fint dn +1, tn +1 − f ext d n+1 , t n+1)(6.3.3)In equilibrium problems, the residuals correspond to the out-of-balance forces;problems in which the accelerations can be neglected are called static problems.The governing equations for both the implicit update of the equations ofmotion and the equilibrium equations are a set of nonlinear algebraic equations inthe nodal displacements, dn+1 .
In equilibrium problems with rate-independentmaterials, t need not be the real time. Instead it can be any monotonicallyincreasing parameter which describes the changing load. If the constitutiveequation is a differential or integral equation, it must also be discretized in time toobtain a set of algebraic equations for the system.6.3.2a.Newmark β-equations.
We will now show that the discreteequations obtained with an implicit time integrator are nonlinear algebraicequations in the unknowns dn+1 . For this purpose we consider a popular class oftime integrators called the Newmark β-method. In this time integration formula,the updated displacements and velocities are given bydn +1 = ˜dn + β∆t 2a n +1(6.3.4)6-8T. Belytschko & B. Moran, Solution Methods, December 16, 1998˜dn = dn +∆tvn + ∆t 2 (1− 2β )a n2(6.3.5)v n+1 = ˜v n + γ∆tan+1(6.3.6)˜v n = v n + (1− γ )∆ta n(6.3.7)Here β and γ are parameters whose useful values are summarized in Box 6.2. Inwriting the time integration formulas, we have segregated the historical values ofthe nodal variables, i.e. those pertaining to time step n, in ˜v n and ˜dn .
Theresulting formulas correspond to the predictor-corrector form given by Hughesand Liu( ). This segregation of the historical terms is convenient for the algebraicoperations which follow and for the construction of explicit-implicit timeintegration procedures.Equation (6.3.4) can be solved for the updated accelerations for β > 0 ,giving(n +1 ˜ n+11a n+1 = β ∆t−d2 d)(6.3.8)Substituting (6.3.8 ) into (6.3.1) gives()(s0 = r = β∆tD 2 M( dn+1 − ˜dn ) − f ext dn+1 , t n+1 + f int dn+1 , t n +1)(6.3.9)which is a set of nonlinear algebraic equations in the nodal displacements dn+1 .Eq.(6.3.9) applies to both the static and dynamic problems. Therefore, in bothcases we consider the discrete problem to be()()find dn +1 so that r dn +1 = 0 subject to g dn +1 = 0((6.3.10))where r dn +1 is given by Eq.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
