Belytschko T. - Introduction (779635), страница 70
Текст из файла (страница 70)
The results of () and () are in factidentical: for a system which has no velocity dependent terms, the perturbationanalysis is identical to taking a derivative of the forces, and the stability of theresult of the perturbation analysis simply depends on the sign of the derivative ofthe forces, which is the second derivative of the work potential.The linearized test for stability used in Example 8.?.1 is not a foolproof test forstability. For example, if y0 = −0.99b , the test indicates that the equilibrium pointis unstable. However, a numerical dynamics solution when started at tat point,will only oscillate with an amplitude of 0.002, which to most engineers would notbe an instability.
The test as posed in Example 1 checks whether any perturbationwill grow at all and conform to the criterion () based on a linearized analysis,which need not conform to any physical notion of instability. In contrast to theexponential instabilities seen in stability analysis of numerical methods, theinstabilities in physical systems will not exhibit continuing growth. What it isdoes predict accurately is that when the dynamics is added to the system, thesystem will not oscillate about the unstable equilibrium point in response to aperturbation but move to oscillating about a nearby point on a stable equilibriumpath.Example 6.5. Consider a linear stability analysis of the beam element shownin Fig.
6.5E. Node 2 is clamped, node 1 is free to rotate and move in the xdirection. Find the equilibrium equation and the equilibrium branches of thesystem.6-65T. Belytschko & B. Moran, Solution Methods, December 16, 199821PlBAµ1CPsolution 1solution 2u1 = 15 IlAθ1Figure 6.5E. Beam model used for stability analysis and equilibrium paths.The displacement boundary conditions imply thatux1 = uy1 = θ1 = uy2 = 0(E6.5.1)Therefore, the only nonzero degrees-of-freedom are uy1 ≡ u1 and θ1 .equations of equilibrium can be deduced from Example ??? to beEA2EA 2u1 −θ =Fl15 1−The(E6.5.2)2EA3EAP 4EI 2EAθ = 0θ1u1 + −u1 + l151535 1(E6.5.3)The above system of two nonlinear algebraic equations in two unknownspossesses two solutions:6-66T. Belytschko & B.
Moran, Solution Methods, December 16, 1998PlEA(E6.5.4)Solution 2: u1 =2l 2 Plθ +15 1 EA(E6.5.5)u1 =Pl 2 15Iθ +28 1 Al(E6.5.6)Solution 1: θ1 = 0, u1 =These two curves are plotted in Figure 6. It can be seen a pitchfork bifurcationoccurs atu1 =15IAl(E6.5.7)This is the critical point for this beam. The corresponding load can be found bysubstituting (???) and θ1 = 0 into Eq. (E6.5.2), which givesFcrit =15EIl2(E6.5.8)The linearized stability of any of the equilibrium paths can be examined byconsidering the linearized equations of motion about a point on the path:()M∆˙d˙ + Kmat + K geo ∆d = 0(E6.5.9)where ∆d here is the displacement from the path.
The equations can be writtenout by using the mass matrix given in Eq. (9.3.18) and the material and tangentstiffnesses given in Eqs. (???) and (???). The resulting equations are AE +˙˙2100∆uρ0 l 0 A0 l ∆u1 12 ˙ ˙ + 4EI ∆θ1 = 0420 0 αl ∆ θ1l (E6.5.10)We will examine the stability of two of the paths for u1 > 15I Al ; the path PA andthe path PC.The problem parameters are Young’s modulus E, the moment of the crosssection I, and the original length of the beam l o . The beam is modeled by asingle element with a linear axial displacement field and a cubic transversedisplacement field.
This is a standard beam element described in Chapter 9. Theunknowns are dT = ux uy θ , where θ is the rotation of the node; nodalsubscripts have been dropped because they all refer to node 1.[]NUMERICAL STABILITY6-67T. Belytschko & B. Moran, Solution Methods, December 16, 1998At this point it is worthwhile to comment on the differences betweenphysical stability and numerical stability. Physical stability pertains to thestability of an solution of a model, whereas numerical stability pertains to thestability of the numerical solution. Numerical instabilities arise from thediscretization of the model equations, whereas physical instabilities areinstabilities in the solutions of the model equations independent of the numericaldiscretization.
Numerical stability is usually only examined for processes whichare physically stable. Very little is known how a “stable” numerical proceduresbehave in physically unstable processes. This shortcoming has importantpractical ramifications, because many computations today simulate physicalinstabilities, and if we cannot guarantee that our methods track these instabilitiesaccurately, then these simulations may be suspect.Numerical stability of a time integration procedure is defined inanalogously to stability of solutions, Eq.
( 6.5.1-2). A numerical procedure isstable if small perturbations of initial data result in small changes in the numericalresponse. More formally, the numerical procedure is stable ifu nA − u nB ≤ Cε ∀n > 0(6.5.29)u 0A − u 0B ≤ ε(6.5.30)whenLATERIt is of interest to note that numerical stability of a process that isphysically unstable cannot be examined by this definition, i.e. we cannot sayanything about the stability of a numerical procedure when applied an equationthat exhibits unstable response. The reason can be seen as follows.
If a system isunstable, then the solution to the system will not satisfy (). Therefore, even if thenumerical solution procedure is stable, it will not satisfy ().General results for numerical stability of time integrators are largely basedon the analysis of linear systems. These results are extrapolated to nonlinearsystems by applying them to the linearized equations. Therefore, we will firstdescribe the stability theory which is used to obtain critical time steps for linearsystems. Next we described the procedures for applying these results to nonlinearsystems.
In conclusion, we will describe some results on stability of timeintegrators which apply directly to nonlinear systems. However, we stress that atthe present time there is no stability theory which encompasses the nonlinearproblems which are routinely solved by nonlinear finite element methods, andmost of our insight into stability stems from the analysis of linear models.Numerical Stability of Linear Systems. Most of the theory of stabilityof numerical methods is concerned with linear systems.
The idea is that if anumerical method is unstable for linear systems, it will of course be unstable fornonlinear systems also, since linear systems are a subset of nonlinear systems.Luckily, the converse has also turned out to be true: numerical methods which arestable for linear systems in almost all cases turn out to be stable for nonlinearsystems. Therefore, the stability of numerical procedures for linear systemsprovides a useful guide to their behavior in both linear and nonlinear systems.6-68T. Belytschko & B. Moran, Solution Methods, December 16, 1998To begin our exploration of stability of numerical procedures, and inparticular the stability of time integrators, we first consider the equations of heatconduction:Mu˙ + Ku= f(6.5.31)where M is the capacitance matrix, K is the conductance matrix, f is the forcingterm and u is a matrix of nodal temperatures. This system is chosen as a startingpoint because it is a first order system of ordinary differential equations, while theequations of motion are second order in time.To apply the definition of stability, we consider two solutions for the same systemwith the same discrete forcing function but slightly different initial data.
The twosolutions satisfy the same equation with the same f, soMu˙ A + KuA = fMu˙ B + KuB = f(6.5.32)If we now take the difference of the two equations, we obtainMd˙ + Kd= 0d = u A − uB(6.5.33)We now consider a two-step family of time integrators:dn +1 = dn + (1−α )∆t˙dn +α∆t ˙dn +1(6.5.34)Since (6.5.33) holds at time steps n and n +1, we can multiply them respectively( )by 1 −α ∆t and α∆t , respectively(1−α )∆tM˙dn + (1−α )∆tKdn = 0,α∆tMd˙ n+1 +α∆tKdn+ 1 = 0(6.5.35)Adding the above two equations and using (6.5.34) to eliminate the derivatives,we obtain(M +α∆tK)dn +1 = (M + (1− α )∆t )Kdn(6.5.36)This equation is in general amplification matrix form: it gives the numericalsolution at times step n +1 in terms of the solution at time step n . Anamplification matrix A is a matrix which gives the solution at time step n +1 in ofthe solution at time step n bydn +1 = Adn(6.5.37)The generalized amplification matrix form isBdn +1 = Adn(6.5.38)We shall now show that the time integrator is stable if the eigenvalues of thegeneralized amplification matrix form lie within the unit circle in the complexplane.6-69T.
Belytschko & B. Moran, Solution Methods, December 16, 1998For this purpose, we need to recall the eigenvalue problem associated with(6.5.33):Kyi = λiMyi(6.5.39)where λi are the eigenvalues and y i the eigenvectors of the system. We recall thatthe matrix M is positive definite and symmetric, whereas the matrix K is positivesemidefinite and symmetric. Because of the symmetry of the matrices, theeigenvectors of (6.5.39) are orthogonal with respect to M and K , which can bewritten asy j Myi = δij ,y jKyi = λiδ ij ( nosumon i)(6.5.40)and from the positiveness of the matrices the eigenvalues are nonnegative.
Thegeneralized amplification equation is associated the generalized eigenvalueproblemAzi =µi Bzi(6.5.41)The eigenvalues of the above system will be shown to control the stability of thetime integrator. In general, these eigenvalues may be complex. Stability thenrequires that the moduli of all of the eigenvalues be less or equal to 1. Otherwiseat least one component of the solution grows exponenetially like zn , so thesolution is unstable. In other words, if we consider the complex plane as shownin Fig, X, then the eigenvalues must lie within or on the unit circle for thenumerical method to be stable.The eigenvectors span the space Rn , so any vector d ∈Rn D can be written as alinear combination of the eigenvalues, see XXX,.
The eigenvectors of (6.5.41)and are identical to the eigenvectors of the (6.5.39) and the eigenvalues are relatedby the following:if A = a1M + a2 K and B = b1M + b2 K then µ =a1 +a2 λib1 +b2 λi(6.5.42)This is shown as follows. Since the eigenvectors y i span the space, we canexpand the eigenvectors z i in terms of y i byz i = ci yi(6.5.43)Substituting the above into (6.5.41), premultiplying by y j and using theorthogonality relations (6.5.40) givesa1 + a2λi = µi (b1 +b2 λi )(6.5.44)from which the last equation in (6.5.42) follows immediately.6-70T.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
