Solution of Post-Buckling&Limit Load Problems,Without Inverting the Tangent Stiffness Matrix&Without Using Arc-Length Methods

T.A.Elgohary,L.Dong,J.L.Junkinsand S.N.Atluri

1 Solving NAEs without Inverting the Jacobian Matrix

A large number of problems in engineering and applied sciences,such as large deformation solid mechanics, fluid dynamics,post-buckling of structural frames,plates,and shells,etc,as characterized by nonlinear differential equations,will lead to a system of nonlinear algebraic equations(NAEs)after discretization:

Fi(x1,...,xn)=0,i=1,...,n,or in their vector-form:

where,a represents the initial guess of the solution,B is the Jacobian(tangent stiffness)matrix given byandkdenotes the number of iteration.

In computational solid mechanics,the trend over the past 30–40 years has been to directly derive the tangent stiffness matrix,B,(rather than forming the nonlinear equations,F(x)=0)through incremental finite element methods[Cris field(1983);Kondoh and Atluri(1986);Riks(1972);Wempner(1971)].Recently,however,[Dai,Yue,and Atluri(2014)]directly derived the system of equations,F(x)=0,for a von Kármán plate theory using the Galerkin method.

As motivated by Eq.2,[Hirsch and Smale(1979)]also introduced the continuous Newton’s method:

Figure 1:Newton’s Method Fails at Limit Points

However,in order to find the solution,both the iterative as well as the continuous Newton’s methods require the the inversion of the Jacobian matrix.On the one hand,inverting the Jacobian matrix in each iteration is computationally very expensive.On the other hand,for complex problems where the Jacobian matrix may be singular,such as near the limit-load points for geometrically nonlinear frames or in elastic-plastic solids,the iterative as well as the continuous Newtons’s methods become problematic,as shown in Fig.1.Various variants of the arc-length methods have been widely used for marching through the limit-points,in post-buckling analyses,such as those presented by[Wempner(1971);Riks(1972);Cris field(1983);Lam and Morley(1992)].These methods generally involve complex procedures by appending various constraints,and monitoring the eigenvalues of the Jacobian matrices.It will be advantageous to have a method to find the solutions for post-buckling problems of structures without inverting the Jacobian matrix,without using the arc-length method,and without worrying about initial guesses for the Newton’s methods.

In order to eliminate the need for inverting the Jacobian matrix in the Newton’s iteration procedure,[Liu and Atluri(2008a)]proposed an alternate first-order system of nonlinear ODEs,for the solution of the NAEs,F(x)=0,by postulating an evolutionary equation for x,thus:

whereνis a nonzero constant andq(t)may in general be a monotonically increasing function oft.In their approach,the termν/q(t)plays the major role of being a stabilizing controller to help obtain the solution even for a bad initial guess,and speeds up the convergence.This Fictitious Time Integration Method was successfully applied to the solution of various engineering problems in[Liu and Atluri(2008a,b,c)].In spite of its success,the Fictitious Time Integration Method was postulated only based on an engineering intuition,does not involve the Jacobian matrix at all,and was shown to only have local convergence.

The homotopy method,as firstly introduced by[Davidenko(1953)],represents one of the best methods to enhance the convergence from a local convergence to a global convergence.Previously,all the homotopy methods were based on the construction of a vector homotopy function H(x,t),which serves the objective of continuously transforming a function G(x)into F(x)by introducing a homotopy parametert(0≤t≤1).The homotopy parametertcan be treated as a time-like fictitious variable,and the homotopy function can be any continuous function such that

Two of the most popular vector homotopy functions are the Fixed-point Homotopy Function:

and the Newton Homotopy Function:

By using the vector homotopy method,the solution of the NAEs can be obtained by numerically integrating:

As can be seen in Eq.7,the implementation of the Vector Homotopy Method necessitates the inversion of the matrixat each iteration.In order to remedy the shortcoming of the Vector Homotopy Method,[Liu,Yeih,Kuo,and Atluri(2009)]proposed to solve the system of NAEs by constructing a Scalar Homotopy Function h(x,t),such that h(x,0)=0?‖G(x)‖=0 and h(x,1)=0?‖F(xiàn)(x)‖=0.As an example,the following Scalar Fixed-point Homotopy Function was introduced in[Liu,Yeih,Kuo,and Atluri(2009)]:

However,it may be more convenient to to define homotopy functions witht∈[0,∞]instead oft∈[0,1],and use a positive and monotonically increasing functionQ(t)to enhance the convergence speed.In this paper,we consider the following two scalar homotopy functions,as denoted by the Scalar Fixed-point Homotopy Function:

and the Scalar Newton Homotopy function:

respectively.

By selecting a driving vector u so that the evolution of˙x is parallel to u,the system of NAEs can be solved by numerically integrating:

With different scalar homotopy functionsh(x,t),differentQ(t),and different driving vectors,u,Eq.11 leads to different variants of scalar homotopy methods,see[Liu,Yeih,Kuo,and Atluri(2009);Ku,Yieh,and Liu(2010);Dai,Yue,and Atluri(2014)].In this paper,we select u such thatThus,ifhfis to be used,Eq.11 leads to:

In this manuscript,Q(t)=etis used for simplicity,while various possible choices can be found in[Dai,Yue,and Atluri(2014)].

For both of these two methods,the inversion of the Jacobian matrix is not involved.However,if a scalar equation is to be considered,Eq.13 will fail at whereB=0.Therefore,Eq.13 should only be used for non-scalar equations,while Eq.12 can be used for general scalar and vector NAEs.

2 Solving the Williams’Toggle with a Scalar Homotopy Method

2.1 The Classical Williams’Solution

The classical toggle problem,as introduced by[Williams(1964)],comprising of two rigidly jointed equal members of lengthland angleβwith respect to the horizontal axis,and subjected to an externally applied vertical loadWat the apex,is shown in Fig.2.

The structure deforms in a symmetrical mode as shown in Fig.3 with the deflected position of the neutral axis of memberrsdenoted byr′s′,

Figure 2:The classical Williams’toggle

Figure 3:Symmetrical deformation of toggle

Following the same assumptions and nomenclature in[Williams(1964)],the externally applied loadWcan be expressed in terms of the deformation at the apex,δ,through the following series of equations,

where,Fis the component of the reaction force perpendicular to the undeflected position of the neutral axis denoted here byrsandPis the component of the force at the end of the member parallel tors.Pis expressed in terms ofδas:

where,AEis the extensional rigidity of the member.Fis then expressed in terms ofδusing nonlinear elastic stability theory as:

Combining Eq.14 through Eq.17,for a given loadW,the vertical displacement,δ,of the Williams’toggle can be found by solving the following scalar nonlinear algebraic equation:

In this study,the set of parametersl,EI,AEare considered to be the same as those presented in[Williams(1964)],and are given in Tab.1.By changing the height of the apex of the toggle,three cases of interest are generated,as shown in Fig.4.The first case takeslsinβ=0.32,which represents the original plot in[Williams(1964)].The second and third cases,lsinβ=0.38 andlsinβ=0.44,respectively,show the effect of raising the apex on the load-deflection curve of the toggle,and introduce limit points in the scalar NAE at which the Jacobian is singular.

2.2 Solving Williams’Equation with a Scalar Homotopy Method

In order to characterize the deflectionδresulting from a specific external load,the scalar NAE(Eq.18)must be solved.To better understand the limitations on solving the NAE utilizing the classical Newton’s method,the behavior of the Jacobian derived analytically from Eq.18 is shown in Fig.5 for the three values oflsinβintroduced earlier.

Table 1:Parameters set in[Williams(1964)]

Figure 4:Three Cases of Load-Deflection Curves for Williams’Toggle

Limit-points are those where the Jacobian becomes close to zero and thus classical Newton’s method will fail.For this end,the previously introduced Scalar Homotopy Methods are used to avoid inverting the Jacobian for the solution of the NAE.As discussed in Section 1,the Scalar Newton Homotopy Method(Eq.13)only works for a system of NAEs.Thus,the Scalar Fixed-point Homotopy Method(Eq.12)is adopted here for the solution of the scalar NAE of the Williams’toggle.Setting the tolerance to 10?6for the original Williams’toggle withlsinβ=0.32,the Williams’equation can be solved for an arbitrary load,here chosen asW=27.11 lb.A comparison between Newton’s Method and the Scalar Fixed-point Homotopy Method is shown in Tab.2

Figure 5:Scalar Jacobian Evaluation of Williams’Scalar NAE

Table 2:Solution of Original Williams’Equation with No Limit Points

The comparison shows the fast convergence speed of Newton’s method achieving the required accuracy in just 9 iterations whereas it took the Scalar Fixed-point Homotopy Method 30 iterations to achieve similar accuracy results.This case oflsinβ=0.32 as shown in Fig.5 has no singularities in the Jacobian,thus the superior performance of Newton’s method is expected.Fig.6 show the evolution of the solution and the fast convergence of Newton’s Method as compared to the Scalar Fixed-point Homotopy Method.

A second case is considered forlsinβ=0.44 where the Williams’equation is solved for an externally applied load selected near the limit point,W=43.79 lb.Both Newton’s Method and the Scalar Fixed-point Homotopy Method are utilized,and the results are shown in Tab.3

After 1000 iterations Newton’s method did not converge to the solution whereas the Scalar Fixed-point Homotopy Method achieved the required accuracy in 345 iterations.Fig.7 shows a comparison between the evolution of the solution for the two methods.It is shown that Newton’s method will keep fluctuating about the solution and not converge to achieve the required accuracy,whereas the Scalar Fixed-point Homotopy Method converges to the solution with the required high accuracy.

Figure 6:Vertical Deflection vs.No.Iterations,lsinβ=0.32

Table 3:Solution of Williams’Equation for Loading near Limit Point

Figure 7:Vertical Deflection vs.No.Iterations,lsinβ=0.44

3 Application to Finite Element Analyses

3.1 A Generalized Finite Element Model for Frame Structures

The currently adapted Scalar Homotopy Methods can be easily combined with general purpose nonlinear finite element programs,by taking the directly derived tangent stiffness matrix at each iteration as the Jacobian matrix,and taking the difference between generalized internal force vector and the external force vector(the residual)as F(x).In this manuscript,explicitly derived tangent stiffness matrices and nodal forces of large-deformation beam-column members are adopted following the work of[Kondoh and Atluri(1986)].The basic derivations of[Kondoh and Atluri(1986)]are briefly reviewed here.

First,the nomenclature and the sign convention used in the derivation for a general beam column member are shown in Fig.8.The functionsw(z)andu(z)describe the displacement at the centroidal axis of the element along thezand thexaxes,respectively.The anglesare the angles between the tangent to the deformed centroidal axis and the line joining the two nodes of the deformed element at nodes 1 and 2,respectively.M1andM2are the bending moments at the two nodes andNis the axial force in the beam member.The total rotation of the beam member is then given by,

The total stretch/deformation of the beam member is then expressed in terms of the displacements at the two nodes as,

Figure 8:Kinematics&Nomenclature for a Beam Member

The axial force and bending moment are non-dimensionalized through,

The non-rigid rotation and the non-dimensional bending moment are decomposed into symmetric and anti-symmetric parts given by,

The relation between the generalized displacements and forces at the nodes of the beam member is given by,

where,haandhsare given by,

The kinematics variables can then be expressed in a vector form for a beam membermas,

The increment of the internal energy of a beam member is then expressed in terms of the increment of kinematics variables vector,the tangent stiffness matrixand the internal force vectoras,

The tangent stiffness matrix,and the internal force vector,for the membermare given by,

where the elements constructing Eq.29 and Eq.30 are given by,

The load-deflection curve generated using the finite element model in Eq.28 is compared against the original Williams’problem withlsinβ=0.32 in Fig.9.Other cases withlsinβ=0.38 andlsinβ=0.44 are shown in Fig.10 and Fig.11,respectively.The Scalar Fixed Point Homotopy Method,Eq.9,is used to generate the load-deflection curves for the finite element model for all three cases.As shown,the finite element model accurately describe the load-deflection characteristics of the Williams’toggle as it agrees well with the solutions of the scalar NAE presented in[Williams(1964)]and summarized in Eqs.14–18.The Scalar Fixed point Homotopy Method method successfully solved the FEM equations capturing the load-deflection relation around the limit points,at which the Newton’s method fails to find the solution,as will be shown in this next subsection.

3.2 Solution of the Finite Element Model Using Scalar Homotopy Methods

Figure 9:Load-Deflection,Williams’Equation&Finite Element,lsinβ =0.32

Figure 10:Load-Deflection,Williams’Equation&Finite Element,lsinβ =0.38

Figure 11:Load-Deflection,Williams’Equation&Finite Element,lsinβ =0.44

Table 4:Solution of Finite Element Model for Loading near Limit Point

The Scalar Fixed-point Homotopy Method,Eq.12,and the Scalar Newton Homotopy Method,Eq.13,are both applied to the finite element model to solve for the deflection given a specific load.As done in the previous section the case oflsinβ=0.44 is examined with the same value of the load applied near the limit point.Setting the tolerance for the relative residual error to be 10?6,the two methods are compared with Newton’s method and the results are shown in Tab.4.Both Scalar Homotopy Methods proved superior to the Newton’s method,as both converged to the solution with the required accuracy whereas the Newton’s method failed to find the solution after 1000 iterations.A zoomed in plot is shown in Fig.12 to illustrate the oscillating behavior of Newton’s method and its failure to find the solution.The Scalar Fixed-point Homotopy Method converged in 160 iterations(Fig.13),which is about one third the number of iterations required by the Scalar Newton Homotopy Method,(Fig.14).This makes the Scalar Fixed-point Homotopy Method more suitable for solving the problem of Williams’toggle,whereas the Scalar Newton Homotopy Method provides a valid alternative to obtain the solution.The Scalar Homotopy Methods developed in this work and in previous works are suitable to solve general nonlinear finite element models with very high accuracy,without inverting the tangent stiffness matrix,and without having to use the computationally expensive arc-length methods.

Figure 12:Residual Error in Newton’s Method

Figure 13:Residual Error in Scalar Fixed-point Homotopy Method

Figure 14:Residual Error in Scalar Newton Homotopy Method

Table 5:Solution of Finite Element Model for Loading near Limit Point

In order to illustrate the efficiency of the scalar homotopy methods the tolerance for the relative residual error is relaxed to match existing finite element solvers(0.1%).For this case an external load of 44 lb.is applied and the results are shown in Tab.5.The two methods achieved the required accuracy within 14 iterations,which demonstrates the power of the scalar homotopy methods in solving engineering problems and the fast convergence that can be achieved when addressing such problems.The Newton’s method failed after 1000 iterations,with the same oscillatory non-convergent behavior shown in Fig.12.Figs.15-16 show the path to convergence of the Scalar Fixed-point Homotopy Method and the Scalar Newton Homotopy Method,respectively.

Figure 15:Residual Error in Scalar Fixed-point Homotopy Method

Figure 16:Residual Error in Scalar Newton Homotopy Method

4 Conclusion

The Scalar Homotopy Method is applied to the solution of post-buckling and limit load problems of plane frames considering geometrical nonlinearities.Explicitly derived tangent stiffness matrices and nodal forces of large-deformation beam-column members are adopted following the work of[Kondoh and Atluri(1986)].By using the Scalar Homotopy Method,nodal displacements of the equilibrium state are iteratively solved for,without inverting the Jacobian(tangent stiffness)ma-trix and without using complex arc-Length methods.This simple method thus saves computational time and avoids the problematic behavior of the Newton’s method when the Jacobian matrix is singular.While the simple Williams’toggle is considered in this paper,extension to general finite element analyses of space frames,plates,shells and elastic-plastic solids will be considered in forthcoming studies.

Acknowledgement:This work is supported by the Texas A&M Institute for Advanced Study(TIAS).It was initiated while S.N.Atluri visited TIAS briefly in January,2014.The support of various U.S.government agencies during 1970-2014 is also thankfully acknowledged.Messrs Le and Riddick of ARL are thanked for their encouragement.

Cris field,M.A.(1983):An arc-length method including line searches and accelerations.International journal for numerical methods in engineering,vol.19,no.9,pp.1269–1289.

Dai,H.;Yue,X.;Atluri,S.N.(2014):Solutions of the von kármán plate equations by a galerkin method,without inverting the tangent stiffness matrix.Journal of Mechanics of Materials and Structures,vol.9,pp.195–226.

Davidenko,D.F.(1953):On a new method of numerical solution of systems of nonliear equations.Dokl.Akad.Nauk SSSR,vol.88,pp.601–602.

Hirsch,M.W.;Smale,S.(1979):On algorithms for solving f(x)=0.Communications on Pure and Applied Mathematics,vol.32,no.3,pp.281–312.

Kondoh,K.;Atluri,S.N.(1986): A simplified finite element method for large deformation,post-buckling analyses of large frame structures,using explicitly derived tangent stiffness matrices.International journal for numerical methods in engineering,vol.23,no.1,pp.69–90.

Ku,C.-Y.;Yieh,W.;Liu,C.-S.(2010): Solving non-linear algebraic equations by a scalar newton-homotopy continuation method.International Journal of Nonlinear Sciences and Numerical Simulation,vol.11,no.6,pp.435–450.

Lam,W.;Morley,C.(1992): Arc-length method for passing limit points in structural calculation.Journal of structural engineering,vol.118,no.1,pp.169–185.

Liu,C.-S.;Atluri,S.N.(2008a): A novel time integration method for solving a large system of non-linear algebraic equations.CMES:Computer Modeling in Engineering&Sciences,vol.31,no.2,pp.71–83.

Liu,C.-S.;Atluri,S.N.(2008b):A novel fictitious time integration method for solving the discretized inverse sturm-liouville problems,for specified eigenvalues.CMES:Computer Modeling in Engineering&Sciences,vol.36,no.3,pp.261–286.

Liu,C.-S.;Atluri,S.N.(2008c):A fictitious time integration method(ftim)for solving mixed complementarity problems with applications to non-linear optimization.CMES:Computer Modeling in Engineering&Sciences,vol.34,no.2,pp.155–178.

Liu,C.S.;Yeih,W.;Kuo,C.L.;Atluri,S.N.(2009): A scalar homotopy method for solving an over/under determined system of non-linear algebraic equations.CMES:Computer Modeling in Engineering and Sciences,vol.53,no.1,pp.47–71.

Riks,E.(1972): The application of newton’s method to the problem of elastic stability.Journal of Applied Mechanics,vol.39,no.4,pp.1060–1065.

Wempner,G.A.(1971): Discrete approximations related to nonlinear theories of solids.International Journal of Solids and Structures,vol.7,no.11,pp.1581–1599.

Williams,F.W.(1964):An approach to the non-linear behaviour of the members of a rigid jointed plane framework with finite deflections.The Quarterly Journal of Mechanics and Applied Mathematics,vol.17,no.4,pp.451–469.