Belytschko T. - Introduction (779635), страница 66
Текст из файла (страница 66)
It can easily be shown via(???) that the above is identical to (E6.2.12).Geometric Stiffness Matrix in Total Lagrangian Form. Thegeometric stiffness is developed from (6.4.15):K IJ = HIJ IH = ∫ B0T SB0 dΩ 0(E6.2.28)Ω0where the B0 matrix is given in (4.6.3), soH=∫Ω01 −11S11] [−1 +1]dΩ 0[l0 +1l0(E6.2.29)Assuming that the stress is constant givesˆ = S11 A0 +1 −1Hl0 −1 +1(E6.2.30)Expanding the above, we obtain the geometric stiffnessK geo +1 0 −1 0 A0S11 0 +1 0 −1=l0 −1 0 +1 0 0 −1 0 +1(E6.2.31)The total tangent stiffness is then given by the sum of the material and geometricstiffnesses, (E6.2.17).6.3.7.Constraints.
Three types of methods are frequently used for treatingthe constraint Eq. (6.3.10). They are:1. penalty methods2. Lagrangian multiple methods3. augmented Lagrangian methodsThese methods originate in constrained optimization theory. As will beseen, they can readily be adapted to the solution of the nonlinear algebraicequations that correpond to the momentum or equilibrium equations, Eq. (6.3.10).To motivate these methods, we begin with a description of how they are appliedto the nonlinear minimization problem, Eq. (6.3.27) (note that while Eq. (6.3.27)6-48T. Belytschko & B.
Moran, Solution Methods, December 16, 1998[...] the stationary prints are determined in the problem is often called aminimization problem because often only the stable equilibrium solutions are ofinterest).The problem then is to solve:r( d) = 0 subject to gI ( d) = 0, I =1 to ncwhere r ∈R n, d ∈Rn . The following notation is usedW,a =∂W∂dara, b = AabGIa =or∂gI= gI , a∂daG Ia =∂gI= gI, a∂da∂r= A = M˙d˙ + f int −f ext∂d[]We will use the conventions GI = GI1 , GI 2, ..., GInc and H I ∈R nc × Rnc , asbefore. Recall thatW,a = ra = faint − faext(6.3.40d)in a conservation problem and that W,ab , the JacobianW,ab = Aab(6.3.40e)matrix of the system.We will also examine the less general problem of finding stationary thepoints of()()W d = 0 subject to g I d = 0(6.3.41)Lagrange Multiplier Method. In the Lagrange method, the constraints areappended to the objective function with the Lagrangian multipliers.
Theminimization Eq. (6.3.41) becomes equivalent to finding the stationary points ofW + λ I gI ≡ W + λT g(6.3.42)Note that at a minimum of W, the augmented function given above has a saddlepoint.The stationary points are found by setting to zero the derivatives of theabove with respect to da and λI :W,a +λ I gI, a ≡ ra +λ I gI, a = 0(6.3.43a)6-49T.
Belytschko & B. Moran, Solution Methods, December 16, 1998gI = 0(6.3.43b)The above is the system of equations for the constrained problem. The constraintintroduces extra forces λI gI, a , which are linear combinations of the Lagrangemultipliers.The linear model of (6.3.43) is the first two terms of a Taylor expansion of(6.3.43), givingra + λ I gI , a + ra, b ∆db + λ I gI, a + λ I gI , ab ∆db = 0(6.3.44a)g I + gI , a∆da = 0(6.3.44b)Using matrix notation we can write the above asA + λI H I GGT ∆d −r − λT G = 0 ∆λ −g (6.3.45)So the linear model has nc additional equations due to the constraint.
Even whenthe A is positive definite, the augmented system of equations will not be positivedefinite because of the zeroes on the diagonal in the lower right hand corner of thematrix. For a linear statics problem with a linear constraints Gd= g , the abovebecomesK GT d f ext = G 0 λ g since(6.3.46)1. A = K for linear statics;2. H I = 0 for linear constraints, see Eq. (6.3.40c);3. the starting value is zero and ∆d = d , ∆λ = λ , and the constraint isGd = 0 .For the general problem with nonconservative materials, dynamics, etc., theLagrangian multiplier method is formulated as follows.
The stationary conditionEq. (6.3.43) can be written( )δW + δ λI gI = 0(6.3.47)From Eq. (B4.6.1) and Eq. (6.3.1)6-50T. Belytschko & B. Moran, Solution Methods, December 16, 1998δW =δW ext −δW int + δW inert(= δd T −f ext + f int + s DM˙d˙)(6.3.48)= δdT r = δdaraSubstituting Eq. (6.3.48) into Eq. (6.3.47) and writing out the differentials on thesecond term givesδda ra +δλ TI gI + λTI gI , aδda = 0∀δdaδλ I(6.3.49)Using the arbitrariness of the differentials in the above implies Eq. (4.3.44-45).Thus the same structure is obtained for a nonconservative dynamic problem. Thelinearization procedure leads to the same equations, Eq.
(6.3.45). While thedevelopment has been given for the virtual work δW , it applies identically tovirtual power.Penalty Method. Again, we first consider conservative problems where thesolution is determined by minimization. In the penalty method, the constraint isenforced by adding the square of the constraints to the poential, so we minimizeas modified potential1W (d) = W (d) + βgI (d) gI (d)2(6.3.50)where β is a penalty parameter.
The penalty parameter is generally orders ofmagnitude greater than other parameters of the problem. The idea is that if β islarge enough, the minimum of W (d) cannot be attained without satisfying theconstraints.The stationary (or minimum) conditions giver + βg T G = 0(6.3.51)(ra, b + βgI , b gI , a +g I gI, ab )db = 1(ra + βgI gIa )(6.3.52)W ,a = W,a + βgI gI, a = 0orThe linear model isor in matrix form()A ∆d = A + βG T G + gI H I ∆d b = −r + βgT G(6.3.53)The size of this system is not increased over the unconstrained system.
For linearconstraints, if A > 0 , A > 0, i.e. the augmented system if positive definite if theoriginal Jacobian matrix is positive definite. The major drawback of penaltymethods is that they impair the conditioning of the equations.6-51T. Belytschko & B. Moran, Solution Methods, December 16, 1998The discrete equations for nonconservative systems are obtained by thesame procedure as for the Lagrange multipliers.
Write the stationary conditions indifferential form10 =δW =δW + βδ (gI gI )2(6.3.54)Now apply Eq. (6.3.48) to replace δW in the above. The discrete equations andlinear model are then given by Eqs. (6.3.51) and (6.3.52), respectively.6.5. STABILITY: CONTINUATION AND ARCLENGTHMETHODSStability.In nonlinear problems, stability of solutions is of considerableinterest.
There are many definitions of stability: stability is a concept thatdepends on the beholder and his objectives. However, some general definitionsare widely accepted. We will here describe a theory of stability that originatesfrom Liapunov(??) and is widely used throughout mathematical analysis, seeSaybol(??) for a very lucid account of its computtional application to a variety ofproblems. We will focus on its application to finite element methods.We will first give a definition of stability and explore its implications.Consider a process that is governed by an evolution equation such as the equationof motion or the heat conduction equations.
Let the solution for the initial()()conditions d A 0 = d 0A be denoted by d A t . Now consider additional solutions()for initial conditions d B 0 = d 0B , where d 0B are small perturbation of d 0A . Thismeans that d 0B is close to d 0A in some norm, i.e.dA0 − d 0B ≤ ε(6.5.1)A solution is stable when for all initial conditions that satisfy (6.5.1), the solutionsatisfies() ()dA t − dB t ≤ Cε∀t > 0(6.5.2)To explain this definition, we specify the norm to be the l 2 norm.
Note that allinitial conditions which satisfy (6.5.1) lie in a neighborhood of d 0A . It is oftensaid that the initial conditions lie in a ball around d 0A . The definition then states()that for all initial conditions which lie in the ball around d 0A , the solutions d B t( ) for all time.must lie in a ball around the solution d A tThis definition isillustrated for a system with two dependent variables in Fig. 6.7. The right-hand6-52T. Belytschko & B. Moran, Solution Methods, December 16, 1998side shows the behavior of a stable system. Here we have only shown twosolutions resulting from perturbations of the initial data, since it is impossible toshow an infinite number of solutions.
The leftr hand side shows an unstablesystem. It suffices for a single solution starting in the ball about d 0A to diverge toindicate an unstable solution.The applicability of this definition to processes we intuitively considerstable and unstable can be seen by the following examples. Consider a beam abeam loaded axially by a vertical load as shown in Fig. 6.8. We first consider thenumerical response when the beam is perfectly straight. The lateral response in()this case is the path is denoted by d A t , and as can be seen, the lateraldisplacement is zero even though the load eventually exceeds the buckling load.If you don't believe this, try it.
The beam will usually not buckle in anincremental solution or a dynamic solution with explicit or implicit integration.Only if roundoff error introduces a "numerical imperfection" will the perfectlystraight beam buckle. We then plot the lateral displacement of the beam as weperturb the location of node 2, which can be considered an initial condition on thedisplacement of that node. The resulting paths are also shown in Fig. 6.8. It canbe seen that when the load is below the buckling load, the paths for different()initial conditions remain in a ball about the d A t . However, when the loadexceeds the buckling load, the solutions for different initial conditions in thelocation of point A diverge drastically.