Static Solution by Implicit Time-Integration

The static behavior of many structures can be characterized by a load-deflection or force-displacement response. If the response plot is nonlinear, the structure behavior is nonlinear. From computational point of view the resolution of a nonlinear problem is much more complex with respect to the linear case. However, the use of relatively recent resolution methods based on sparse iterative techniques allows saving substantially in memory.

Linear Static Solver

A linear structure is a mathematical model characterized by a linear fundamental equilibrium path for all possible choices of load and deflection variables.

This implies that:

The response to different load systems can be obtained by superposition,
Removing all loads returns the structure to the reference position.

The requirements for such a model to be applicable are:

Perfect linear elasticity for any deformation,
Infinitesimal deformation,
Infinite strength.

Despite of obvious physically unrealistic limitations, the linear model can be a good approximation of portions of nonlinear response. As the computational methods for linear problems are efficient and low cost, Radioss linear solvers can be used to find equilibrium of quasi-linear systems. The Preconditioned Conjugate Gradient method is the iterative linear solver available in Radioss. The algorithm enables saving a lot of memory for usual application of Radioss as a sparse storage method is used. This means that only the non-zero terms of the global stiffness matrix are saved. In addition, the symmetry property of both stiffness and preconditioning matrices is worthwhile to save memory.

The performance of conjugate gradient method depends highly to the preconditioning method. Several options are available in Radioss using the card /IMPL/SOLV/1. The simplest method is a so-called Jacobi method in which only the diagonal terms are taken into account. This choice allows saving considerable memory space; however, the performance may be poor. The incomplete Choleski is one of the best known effective preconditioning methods. However, it can result in negative pivots in some special cases even if the stiffness matrix is definite positive. This results a low convergence of PCG algorithm. The problem can be resolved by using a stabilization method. ¹ Finally, the Factored Approximate Inverse method may be the best choice which is used by default in Radioss.

Nonlinear Static Solver

As explained in the beginning of this chapter, a nonlinear behavior is characterized by a nonlinear load-deflection diagram called path. The tangent to an equilibrium path may be formally viewed as the limit of the ratio force increment on displacement increment. This is the definition of a stiffness or more precisely the tangent stiffness related to a given equilibrium state. The reciprocal ratio is called flexibility. The sign of the tangent stiffness is closely associated with the stability of an equilibrium state. A negative stiffness is necessary associated with unstable equilibrium. A positive stiffness is necessary but not sufficient for stability.

The problem of nonlinear analysis can be viewed as that of minimising the total potential energy

Π

which is a function of the total displacement

X

. A truncated Taylor series then leads to:(1)

Π_{n} (X + δ X) = Π_{0} (X) + \frac{\partial Π}{\partial X} δ X + \frac{1}{2} δ X^{T} \frac{\partial^{2} Π}{\partial X^{2}} δ X + ...

Where the subscripts

n

and

0

denote respectively final and initial configurations. The term

\frac{\partial Π}{\partial X}

can be identified as the out-of-balance forces or gradient

F

, of the total potential energy which is the difference between the internal force vector

F_{int}

and the external force vector

F_{ext}

. The term

\frac{\partial^{2} Π}{\partial X^{2}}

describes the tangent stiffness matrix

K_{T}

. The principle of minimum energy and the equilibrium of stable state give:(2)

δ Π = Π_{n} (X + δ X) - Π_{0} (X) = 0

Which is implied in 式 1:(3)

K_{T} δ X = F_{ext} {(X)-F}_{int} (X)

The tangent matrix

K_{T}

should be defined positive at the equilibrium point for stable case:(4)

δ X^{T} K_{T} δ X > 0

式 3 can handle the solution of the nonlinear problem when an incremental method is used. The solution methods are generally based on continuous incremental and corrective phases. The most important class of corrective methods concerns the Newton-Raphson method and its numerous variants as modified, modified-delayed, damped, quasi and so forth. All of these Newton-like methods require access to the past solution. In the following section the conventional and modified Newton methods under general increment control are studied.

Newton and Modified New Methods

As you will often prefer to trace the complete load/deflection response or in other words, the equilibrium path, it is useful to combine the incremental and iterative solution procedures. You can recall that the purpose is to solve Nonlinear Static Solver, 式 3 which can be written in residual form:(5)

R (X, λ) = F_{ext} (X, λ {)-F}_{int} (X) = 0

with

F_{ext} (X, λ) = F_{C} + λ F_{ext} (X)

. This equation represents a system of n algebraic nonlinear equations depending on only one loading parameter

λ

. If the loading depends to only one loading variable independent to the state of deflection, you have:(6)

F_{ext} (X, λ) = λ F_{ext} (X)

Several techniques are available to resolve 式 5. In some situations, the parameter $λ$ is fixed, and the equations are resolved to determine n components of $X$ in order to verify 式 5. In this case, the technique is called load control method. Another technique called displacement control consists in fixing a component of $X$ and searching for $λ$ and 'n-1' other components of displacement vector $X$ . A generalization of displacement control technique will enable to imply several components of displacement vector by using an Euclidian norm. The method is called arc-length control and intended to enable solution algorithms to pass limit points (that is, maximum and minimum loads). The techniques making possible to obtain the load-deflection curve by finding point by point the solution are called piloting techniques.

When the piloting technique is chosen for a given step, the associated solution is obtained by an iterative resolution of so-called Newton-Raphson methods. At iteration

i

, the residual vector

R^{i}

is:(7)

R^{i} = F_{ext} {(X}^{i}, λ^{i} {),-F}_{int} {(X}^{i})

A correction

Δ X

and

Δ λ

can be considered with:(8)

R^{i + 1} = R^{i} + {[\frac{\partial R}{\partial X}]}^{i} ΔX + {[\frac{\partial R}{\partial λ}]}^{i} Δ λ

Combining 式 8 with 式 7, you obtain:(9)

K_{T}^{i} Δ {X-F}_{ext}^{i} Δ λ = R^{i}

as $R^{i + 1} = 0$ and:

$X^{i + 1} = X^{i} + Δ X$

(10)

λ^{i + 1} = λ^{i} + Δ λ

The tangent matrix

K_{T}^{i}

is obtained by assembling the elementary matrices

k_{T}^{i}

. It corresponds to:(11)

K_{T}^{i} = \frac{\partial F_{int}}{\partial X} - \frac{\partial F_{ext}}{\partial X}

Using load control technique, the standard Newton-Raphson method resolves 式 9 to 式 11 by applying a known load increment

Δ λ

as illustrated in 図 1. The tangent matrix is updated and triangulized at each iteration. This insures a quadratic convergence to exact solution.

However, it is possible to save computation time which depends on the size of the problem and on the degree of the nonlinearity of the problem. The method is called modified Newton-Raphson which is based on the conservation of the tangent matrix for all iterations (図 2). This method can also be combined with the acceleration techniques as line-search explained in Line Search Method to Optimize the Resolution.

The convergence criteria may be based on Euclidian norm of residual forces, residual displacements or energy where an allowable tolerance is defined.

Line Search Method to Optimize the Resolution

The Newton-Raphson resolution of Newton and Modified New Methods, 式 9 implies updating the variables at each iteration with Newton and Modified New Methods, 式 10. The new estimation of

X^{i + 1}

does not satisfy Newton and Modified New Methods, 式 9 only if

R^{i + 1} = 0

. In order to reduce the number of iterations the line-search method is used. The line-search technique is an important feature of most numerical techniques used in optimization problems. ² The method consists in introducing a parameter

α

, such as:(12)

X^{i + 1} = X^{i} + α Δ X

Where, $α$ is obtained to minimize the total potential energy or to satisfy the principle of virtual works. The techniques to determine $α$ use often a Raleigh-Ritz procedure with only one unknown parameter.

The principle of virtual work can be written in the general form:(13)

$W (X, δ X) = δ X • R (X) = 0$ For all kinematical acceptable $δ X$

Considering 式 12, write:(14)

δ X = α Δ X

and:(15)

$δ W = δ α Δ X • R (X^{i} + α Δ X) = 0$ for all $α$

Then,

α

is determined from:(16)

Δ X • R (X^{i} + α Δ X) = 0

which leads to a three-order polynomial equation in

α

for elastic materials:(17)

C {}_{1}+ C {}_{2}α + C {}_{3}α^{2} + C_{4} α^{3} = 0

The coefficients $C {}_{1}$ , $C {}_{2}$ , $C {}_{3}$ and $C {}_{4}$ can be expressed in terms of displacements $X^{i}$ and the increment of displacements $Δ X$ .

Arc Length Method

To obtain the load-deflection behavior of a structure, the load or the displacement of a given point of the structure must be parameterized. Up to now, you have parameterized the load by the time

t

. However, a single parameter is not always sufficient to control in an optimum way the time step. On the other hand, it is not possible to pass limit points with "snap-through" and "snap-back" when using load-controlled or displacement-controlled techniques. This is due to the fact that the increase in load or in a given displacement component may result a dynamic response losing a part of load-deflection curve as shown in 表 1 (a) and (b).

表 1. Various Load-deflection Curves and Step-by-Step Solution by Arc-length Method
(a) Snap-through	(b) Snap-backs
(c) Arc-length method: intersection of the equilibrium branch with the circle about the last solution	(d) Buckling with or without imperfections

The tracing of equilibrium branches are quite difficult. In arc-length method, instead of incrementing the load parameter, a measure of the arc length in the displacement-load parameter space is incremented. This is accomplished by adding a controlling parameter to the equilibrium equations.

The arc-length method was originally introduced by Riks ³ and Wempner. ⁴ Considering a function

f

implying several components of the displacement vector

X

, the arc-length method consists in determining in each step the Euclidian norm of the increase in

X

:(18)

f = 〈 X^{i + 1} - X^{n} 〉 {X^{i + 1} - X^{n}} - {(Δ S_{n})}^{2} = 0

This leads to:(19)

〈 X^{i + 1} - X^{n} 〉 {X^{i + 1} - X^{n}} = {(Δ S_{n})}^{2}

And (20)

a {(Δ λ)}^{2} + b (Δ λ) + c = 0

With:

$a = 〈 Δ X_{F} 〉 {Δ X_{F}}$ ; ${Δ X_{F}} = {[K_{T}^{n}]}^{- 1} {F_{e x t}}$

$b = 2 〈 Δ X_{F} 〉 {Y}$ ; ${Y} = {Δ X_{F}} + {X^{i} - X^{n}}$

$c = 〈 Y 〉 {Y} - {(Δ S_{n})}^{2}$

In each of the Newton-Raphson iterations, 式 20 must be resolved to select a real root. If there is no root, $Δ S_{n}$ should be reduced. The most closed root to the last solution is retained in the case of two real roots.

表 1(c) illustrates the intersection of the equilibrium branch with the circle about the last solution.

¹ Ajiz M.A. and Jennings A., 「A robust incomplete Choleski-conjugate gradient algorithm」, Int. J. Method Eng., Vol. 20, pp. 949-966, 1984.

² Flectcher R., 「Practical methods of optimization」, 2nd edition, Wiley, 1987.

³ Riks E., 「An incremental approach to the solution of snapping and buckling problems」, Int. J. Solids & Structs, Vol. 15, pp. 529-551, 1979.

⁴ Wempner G.A., 「Discrete approximations related to nonlinear theories of solids」, Int. J. Solids & Structs, Vol. 7, pp. 1581-1599, 1971.