ANALYSIS REFERENCE Chapter 5.

Algorithm

Section 1 Simultaneous Equation Solver

The simultaneous equation solver is a method for finding the solution u of a linear matrix equation in the form
shown in equation (5.1.1):

Ku = p (5.1.1)

Simultaneous equation solvers are not only used for linear static structural analysis, but also for all types of
analysis such as eigenvalue analysis, dynamic analysis, nonlinear analysis, etc. General solvers include the
Gauss elimination method, direct solver based on the decomposition method, and the iterative solver, which
converges to a solution that minimizes iterative calculations. Direct solver is generally used for structural
analysis because it is not affected by the numerical properties of matrices and can find the solution safely.
However, when the size of the problem increases, the memory capacity and computation amount tend to
increase rapidly. Hence, the iterative solver is recommended for large problems because it requires relatively
less memory capacity. However for structural analysis, the iterative solver may not provide the wanted
solution due to the numerical properties of matrices, and the number of iterative calculations needed to obtain
the converging solution may be large. FEA NX provides a function that automatically determines the direct
solver or iterative solver, depending on the size of the problem in question.

In direct solver, the simultaneous equation solution is found in two steps. The first step is matrix
decomposition and the second step is the forward-backward substitution (FBS) process. The LU solver,
generally used for asymmetric matrices, can be applied to matrix decomposition in the following form for the
symmetric stiffness matrix K obtained in finite element analysis.

LLT u  p or LDLT u  p (5.1.2)

L : Lower triangular matrix

D : Diagonal matrix

Generally, the matrix decomposition method that includes D is needed when the stiffness matrix is not
T
definitely positive. FEA NX uses the LL form matrix decomposition method (Cholesky decomposition
method) for linear static structural analysis. For eigenvalue analysis or nonlinear analysis, the positive definite
T
condition cannot be guaranteed and so, the LDL form matrix decomposition method is used.

When applying the direct solver, the sparse matrix needs to be applied appropriately. Generally, the stiffness
matrix K generated in finite element analysis is a sparse matrix with multiple '0', and the required memory

306 | Section 1. Simultaneous Equation Solver

 ANALYSIS REFERENCE Chapter 5. Algorithm

capacity and computation amount differs greatly depending on how this sparsity is used. Hence, FEA NX not
only provides the direct solver for general dense matrices that do not use the sparse matrix, but also provides
the multi-frontal solver, which appropriately uses the sparse matrix to greatly reduce the memory capacity
and computation amount.
The multi-frontal solver requires reordering of DOFs to minimize the memory capacity and computation
amount using the sparse matrix, and matrix decomposition is performed by separating the matrix into
multiple fronts according to this reordered information. Figure 5.1.1 displays the effective computation order
of a rectangular mesh generated by DOF reordering. The algorithm used to implement DOF reordering is a
recursive bisection, and forward substitution is done in the same order as the matrix decomposition whilst
backward substitution is done in the opposite order.

Figure 5.1.1 Matrix decomposition

1 2 1 4 1 2 1
order of multi-frontal solver

3 3

1 2 1 1 2 1

The multi-frontal solver used in FEA NX does not assemble and save the stiffness matrix of the entire region
individually and hence, requires less memory capacity than the general multi-frontal solver. The out-of-core
analysis function is supported to provide additional hard disk memory automatically during memory shortage
when solving large problems.

Also, implementation of the multi-frontal solver uses the computation ability of the Graphics Processing Unit
(GPU) to process calculations. The recent demand for complex problems highlights the importance of the
simultaneous equation solver performance, which is the core of finite element analysis. The GPU consists of
multiple computation units (cores) and provide a much higher computational performance than the CPU. The
GPU is applied to real matrix decomposition, which takes the longest computation time, to provide an overall
improved computational performance.

Section 1. Simultaneous Equation Solver | 307

 ANALYSIS REFERENCE Chapter 5. Algorithm

The iterative solver is a method of reducing the error of the approximate solution through iterative
calculations and so, it is very important to reduce the convergence error using only a small number of
calculations. Generally, the number of iterative calculations is determined by the preconditioning method.
FEA NX uses the SA(smoothed aggregation) AMG(algebraic multi-grid)1 methods which are preconditioning
methods that are known to be stable, regardless of the element shape. The number of calculations of the AMG
method is not greatly affected by the number of DOF because it uses a multi-grid, and this method displays
stable convergence when used on elements that have displacement and rotation nodal DOFs such as shell
elements. The multi-grid is composed automatically for the iterative solver using the AMG method, and this is
created by the representative DOFs of the adjacent node set and each node set.

Figure 5.1.2 Example of node set

for multi-grid composition

As explained above, the performance of the direct solver and iterative solver differ on the size of the problem
and FEA NX provides an automatic selection function to determine the solver. When using the automatic
section function, the direct solver using the dense matrix is selected for small size problems, the multi-frontal
solver is selected for medium size problems and the AMG iterative solver is selected for large size problems.

The automatic selection criterion is determined by considering the following points.

► When empirical condition is known : Determined with reference to the number of user input nodes or
elements
► When empirical condition is unknown : Determined within the program with reference to the number of
model DOF and system memory size

308 | Section 1. Simultaneous Equation Solver

ANALYSIS REFERENCE Chapter 5. Algorithm

Section 2 Eigenvalue Extraction

Eigenvalue extraction is a fundamental algorithm of normal mode analysis, and the eigenvalue extraction
problems in normal mode analysis have the following form.

Ki  i Bi  0 (no summation) (5.2.1)

K : Stiffness matrix
B : Mass matrix ( M ) when performing normal mode analysis

The eigenvalue extraction method in FEA NX is coupled and changes with the simultaneous equation solver.
The Lanczos resampling is used for the multi-frontal solver (the default value of the simultaneous equation
solver), and eigenvalue extraction or direct solver using the dense matrix is used for the dense matrix solver.
Each method has the following characteristics.

• Lanczos resampling
► Appropriate for large sized problems.
► Because the eigenvalue can be omitted, use of the Sturm sequence check option is recommended.

• Direct solver using the dense matrix

► The performance can suddenly decrease when the number of DOF is around 103 and so, it is appropriate
for small scale test models.
► The eigenvalue is not omitted.

Lanczos resampling is a method of finding the approximate eigenvalue using the tridiagonal matrix that arises
1
when generating the Krylov subspace span(V1, V2 ,..., Vk ) . For effective eigenvalue calculation, the block
2
tridiagonal matrix can be used, and because the tridiagonal matrix size is maintined similar to the number of
eigenvalues, the computation speed is very fast and it is appropriate for large scale problems. However,
eigenvalue omission can occur and so, it is useful to use the checking option.
The direct solver using the dense matrix goes through the stiffness matrix decomposition, tridiagonal matrix
generation and eigenvalue calculation processes. Tridiagonal matrix generation and eigenvalue calculation is
done for the entire matrix and eigenvalue omission does not occur. However, it is inappropriate for solving
large size problems.

1
Hughes, T.J.R., The Finite Element Method, Prentice-Hall International, Inc., New Jersey, 1987
2
Cullum, J. and Donath, W., “A Block Lanczos algorithms for computing the q algebraically largest eigenvalues and a
corresponding eigenspace of large real symmetric matrices,” Proc. 1974 IEEE Conference on Decision and Control, IEEE
Computer Society, 1974

Section 2. Eigenvalue Extraction | 309

 ANALYSIS REFERENCE Chapter 5. Algorithm

Eigenvalue computation range

For normal mode analysis, the number of eigenvalues and its range considers the modal participation factor
(equation 5.3.1) or modal effective mass (equation 5.3.3), or can be determined with reference to the frequency
region of interest. If the number and range of eigenvalues is determined, it can be set using the following
inputs.

Table 5.2.1 Setting the number Number of

and range of eigenvalues Variable setting ( v1 , v2 , N input or not input) Eigenvalue range
eigenvalues
v1 , v2 , N v1  v  v2 Maximum N
v1 , not input, N v1  v Maximum N

not input, v2 , N v  v2 Maximum N

not input, not input, N   v   Maximum N

v1 , v2 , not input v1  v  v2 All eigenvalues

v1 , not input, not input v1  v All eigenvalues
not input, v2 , not input v  v2 All eigenvalues
not input, not input, not input   v   All eigenvalues

The v1 , v2 inputs above are the frequency (Hz) in normal mode analysis.

Eigenvalue computation results

Eigenvectors, which are the results of the eigenvalue problem, satisfy equation (5.2.2), even if its size changes.

a(Ki  i Bi )  Ki  i Bi  0 (5.2.2)

i  ai

Hence, a method is needed to express the size of the calculated eigenvalue consistently. FEA NX applies the
eigenvector normalization process such that the following equation is satisfied, depending on the analysis
type.

iT Mi  1 (5.2.3)

The eigenvalue calculation algorithm is only an approximate solution, even when the direct solver for the
dense matrix is used, and its accuracy cannot be guaranteed. Therefore, FEA NX selects the following values as
the eigenvalue calculation results to check the accuracy of the calculated eigenvalue and eigenvector.

310 | Section 2. Eigenvalue Extraction

 ANALYSIS REFERENCE Chapter 5. Algorithm

Table 5.2.2 Calculation results Result article Calculation

excluding eigenvalue and
Generalized mass bi  iT Bi
eigenvector
Generalized stiffness ki  iT Ki

iT1Ki iT1Bi
Orthogonality loss  i  max( , )
ki bi
Ki  i Bi
Error measure ei 
Ki

Section 2. Eigenvalue Extraction | 311

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 3 Effective Mass and Mode

Superposition
After calculating the natural frequency, natural period and mode shape using mode analysis, these can be used
3.1 to compute useful information such as modal effective mass or modal participation factor. The i th modal
Effective Mass participation factor is expressed as i and can be calculated as follows:

1 T
i  i MT ,   1, 2,3, 4,5, 6 (no summation)
mi (5.3.1)
mi   Mi
i
T
(generalized mass)

 : DOF direction (1~3 : displacement, 4~6 : rotation)

Here, T is the matrix that represents the size of the directional stiffness behavior and it is defined for each
node to have the following property:

1 0 0 0 z  z0 y0  y   e1 
0  
 1 0 z0  z 0 x  x0  e2 
0 0 1 y  y0 x0  x 0  
e3 

    , e    (5.3.2)
0 0 0 1 0 0  e4 
0 0 0 0 1 0  e5 
  
0 0 0 0 0 1   e6 


x0 , y0 , z0 represent the center of rotation. FEA NX sets it as an arbitrary node or the center of mass for the
entire model.

The modal effective mass is also defined for each direction and can be simply calculated using the modal
participation factor as follows:

  ( i ) mi
mieff (5.3.3)
2

Adding the effective mass for all modes is the same as the mass of the entire model, excluding the nodes that
have assigned constraint conditions.

312 | Section 3. Effective Mass and Mode Superposition

ANALYSIS REFERENCE Chapter 5. Algorithm

The mode superposition method can be applied to dynamic response analysis. Mode superposition uses the
3.2 eigenmode found from eigenvalue analysis (instead of directly solving the linear dynamic equilibrium equation)
Mode Superposition to solve the size reduced mode equilibrium equation as shown below:

Mu(t )  Cu(t )  K u(t )  f (t ) (5.3.4)

The spatial coordinate system displacement u  t  can be expressed as a combination of the modal
displacement ξ  t  using the eigenmode shape Φ as follows:

u t   Φξ t  , Φ  1 2 ... N  (5.3.5)

Using this, the dynamic equilibrium equation (5.3.4) can be expressed in the modal coordinate system as
follows:

[ΦT MΦ]ξ(t )  [ΦT CΦ]ξ(t )  [ΦT KΦ]ξ(t )  ΦT f (t ) (5.3.6)

Generally when mode superposition is applied, the high order modes are excluded and only partial low order
modes are used to compose the eigenmode shape Φ and so, equation (5.3.6) is an approximation of equation
(5.3.4). Hence, if an insufficient number of eigenmodes is included in the calculation for expressing the actual
physical displacement, the accuracy of the calculated results can fall greatly.

The mode equilibrium equation (5.3.7) is expressed independently for each mode when the modal damping
matrix ΦT CΦ is '0', as shown below:

mii (t )  kii (t )  pi (t ) (5.3.7)

mi : i th modal mass pi : i th modal load
ki : i th modal stiffness i : i th modal displacement

Using the mode superposition method above, the equilibrium equation can be reduced to have the same
number of variables as the number of calculated eigenmodes, and analysis can be performed effectively when
the mode equilibrium equation is fully separated between modes.

Damping term treatment

If the modal damping matrix is diagonalized and the coupling removed for the modal equilibrium equation
(5.3.7) that is reduced using the eigenmode, it can be expressed as a separated form for each mode like
equation (5.3.8).

mii (t )  bii (t )  kii (t )  pi (t ) (5.3.8)

Section 3. Effective Mass and Mode Superposition | 313

 ANALYSIS REFERENCE Chapter 5. Algorithm

bi : i th modal damping

Or, it can be expressed as follows.

1 (5.3.9)
i (t )  2 iii (t )   2i (t )  pi (t )
i
mi

 i  bi  2mii  : Modal damping ratio

  ki mi
i
2
: Modal frequency

The modal damping value can be input differently according to the frequency and in this case, the modal
damping value is added to the modal damping matrix ΦT CΦ , which is composed of other general damping
values such as mass-proportional damping, stiffness-proportional damping etc. Hence, modal separation of
the modal equilibrium equation is possible when the modal damping matrix ΦT CΦ is a diagonal matrix, and
this is applicable when the proportional damping coefficient and structural damping is constant for each
element and when damping elements (spring, damper) do not exist. If not, the coupling between equilibrium
equations of each mode, due to the un-diagonalized modal damping term, needs to be considered.

Enforced motion
When enforced motion is given in the mode superposition method, it cannot be applied directly to the modal
equilibrium equation. FEA NX uses the following processes to apply enforced motion.
Firstly, the equilibrium equation (5.3.4) is separated into the DOFs with and without enforced motion.

 M11 M12  u1  C11 C12  u1   K11 K12   u1   f1 

M M  u   C C  u   K K  u   f  (5.3.10)
 21 22   2   21 22   2   21 22   2   2
u1 : Displacement of unconfined DOF
u2 : Displacement of DOF confined by enforced motion
f1 : Load acting on unconfined DOF
f2 : Confining force of DOF confined by enforced motion

Separating the unconfined DOF displacement u1 into the following quasi-static displacement u1qs and dynamic
relative displacement y is as follows:

u1  u1qs  y
1
(5.3.11)
u1qs  K 11 K 12u 2

314 | Section 3. Effective Mass and Mode Superposition

ANALYSIS REFERENCE Chapter 5. Algorithm

Rearranging for the dynamic relative displacement y ,

M11y  C11y  K11y  f1   M11K11

1
K12  M12  u2 (5.3.12)

The damping related terms on the right hand side of the equation were ignored. Applying the mode
superposition method and expressing the equation using the modal relative displacement y t   Φ11x t  is
as follows:

Φ11
T
M11Φ11  x  Φ11
T
C11Φ11  x  Φ11
T
K11Φ11  x  Φ11
T
f1   M11K11
1
K 12  M12  u 2 
 (5.3.13)
1
u1  u1qs  y  K11 K12u 2  Φ11x

If the K11 has singularity because of an existing rigid-body mode in the structure, the singularity can be
removed by appropriate shifting using the stiffness matrix K11 and mass matrix M11 .

Residual vector
As explained above, errors can occur due to high order modes that are not included in the eigenmode shape Φ
when using mode superposition. To reduce such errors, FEA NX uses the residual vector R , which is composed
perpendicular to the existing eigenmode, for the mass matrix M and stiffness matrix K .

R  K 1 (I  MΦΦT )F (5.3.14)

Here, F is generally composed of the load vector and the damping force is included when a damping element
exists.
3
FEA NX uses the method suggested by Dickens etc. to find the augmented mode shapes perpendicular to the
residual vector R . This is added to the existing eigenmode shape Φ for applying mode superposition.

3
J.M. Dickens, J.M. Nakagawa, and M.J. Wittbrodt, “A Critique of Mode Acceleration and Modal Truncation
Augmentation Methods for Modal Response Analysis” Computers & Structures, Vol 62, No. 6, 1997, pp. 985-998

Section 3. Effective Mass and Mode Superposition | 315

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section44
Section Dynamic Response
FEA NX uses direct time integration and mode superposition to obtain the transient response of the linear
4.1 equation of motion shown in equation (5.3.4). For the direct time integration of linear problems, the implicit
Time Integration method is used.

Implicit direct integration

4
FEA NX uses the  method ( HHT   ) suggested by Hilber, Hughes, Taylor for implicit direct integration.
5
The HHT   method is a general form of the Newmark method and has a controllable numerical damping
effect. Using this, the high frequency noise can be controlled and it has a 2 order accuracy for time steps, just
like the Newmark method. The HHT   method uses the following modified dynamic equilibrium equation:

Man1  (1   H ) Cvn1  f int,n1  f ext ,n1    H Cvn  f int,n  f ext ,n   0 (5.4.1)

Here, a n1 and v n1 each represent the acceleration and velocity vector of the n 1 th time step and
 H [1 3,0] is the coefficient that determines the numerical damping effect. When considering the effects
of non-mechanical strain, such as thermal expansion of the material, and the internal forces due to in-situ
stress and pore pressure, the internal forces of linear analysis can be expressed as the following equation
including the product of stiffness matrix and DOF.

f int,n1  Ku n1  f nonmech,n1  f int,0 (5.4.2)

Introducing the time step equation from the Newmark method, the velocity, displacement and acceleration at
time steps n, n  1 can be expressed using the following relationship:

v n1  v n  t  an1  (1   )an 

1
(5.4.3)
u n1  u n  tv n  t 2  2 an1  (1  2 )an 
2

Recomposing the equilibrium equation (5.4.1) using equations (5.4.2) and (5.4.3), the following simultaneous
equation with the displacement at time n  1 as a variable can be obtained as follows:

4
H.M Hilber, T.J.R. Hughes, and R.L. Taylor, “Improved Numerical Dissipation for Time Integration Algorithms in
Structural Dynamics,” Earthquake Engineering and Structural Dynamics, Vol 5, No. 3, 1977, pp. 283-292
5
M. Newmark, “A Method of Computation for Structural Dynamics,” ASCE Journal of the Engineering Mechanics
Division, Vol. 5, No. EM3, 1959, pp. 67-94

316 | Section 4. Dynamic Response

ANALYSIS REFERENCE Chapter 5. Algorithm

K eff u n1  f eff

1 (1   H )
K eff  M C  (1   H )K ,
t 2 t
f eff  f int,0  (1   H ) f ext ,n1  f nonmech ,n1    H f ext ,n  f nonmech ,n   (5.4.4)
 1 n 1 n  1  
M u  v   1 a n  
 t t  2  
2

 (1   H ) n  (1   H )  n    
C u   1 v  t (1   H )   1 a n    H Ku n
 t     2  

eff
The right hand side f from equation (5.4.4) is determined by the internal force and calculated displacement,
velocity, acceleration at time step n . When the right hand side is determined, the displacement vector un1 at
n 1 can be calculated using the simultaneous equation solver explained in the section above. The velocity and

acceleration at n 1 can be obtained by substituting this calculated displacement into the Newmark time step
equation (5.4.3). The transient response of the structure can be calculated by the time integration that repeats
the processes outlined above.
eff
The effective stiffness matrix ( K ) in the left hand side of equation (5.4.4) reuses the once decomposed
matrix when the time step is kept constant, allowing effective analysis by only repeating the front-back
substitution process.
HHT   time integration has unconditional stability when   (1  2 H ) / 2 ,   (1   H ) / 4 and when
2

 H  0 . It is specialized in to the Newmark method that uses the average acceleration. FEA NX uses a
default value of  H  0.05 .

Damping effect
FEA NX considers two types of damping: mass-proportional damping and stiffness-proportional damping.
There is also mode damping, which is only applied for mode superposition as mentioned in section 5.3.2. The
damping effects in linear time history analysis are applied to the damping matrix C in the following form:

C   ej Mej   jeK ej  B (5.4.5)

 ej : Mass proportional damping coefficient for j th element
 e
j : Stiffness proportional damping coefficient for j th element
th
M ej : Mass matrix of j element
th
K ej : Stiffness matrix of j element

B : Damping matrix due to damping element (damper)

Section 4. Dynamic Response | 317

 ANALYSIS REFERENCE Chapter 5. Algorithm

Application of mode superposition

To use time integration using mode superposition, the mass in the mode equilibrium equation (5.3.6) is set to
'1' and rewritten as follows:

pi
i (t )  C ij j (t )  i2i (t )  pi (t )  pi (t  t )  t
t (5.4.6)
C ij  [C]ij  [ΦT CΦ]ij

Time integration using mode superposition can be classified into two types, depending on the ductile state of
the mode damping matrix C ij :

► Uncoupled system
If the mode damping matrix C ij is diagonalized and the ductility is removed, the response is analyzed
independently for each mode and the displacement and velocity of each time step is determined from the
displacement and velocity of the previous time step using the following equation. The modal integral
coefficients a
i
, b
i
at the i th mode can be obtained by finding the particular solution and homogeneous
solution of (5.4.6) and applying it to the initial condition (displacement and velocity of the previous time step).

in1   a11i a12i  in   b11i b12i   pin 

 n1    i i  n
 i i   n 1 
(5.4.7)
i   a21 a22  i  b21 b22   pi 

► Coupled system
If ductility is not removed from the mode damping matrix, the ductility between modes needs to be
considered and modal analysis cannot be performed independently. In this case, FEA NX separates the mode
damping matrix into the following diagonal component ( Cdiag ) and off-diagonal component ( Coff ) and treats
the damping force of the off-diagonal component as an external force for analysis.

C  Cdiag  Coff (5.4.8)

In this case, all displacements are independent and the mode velocity is softened to compose the following
simultaneous equation. If the time step is fixed, it can be solved without extra matrix decomposition, just like
direct time integration.

I B12 Coff  ξ n1   A11

T
A12  ξ n   B11 B12  p n  CToff ξ n 
      
0 I  B CToff  ξ n1   A 21 A 22  ξ n  B 21 B 22   p n1  (5.4.9)
 22 
A  diag (a ), B  diag (b )
i i

318 | Section 4. Dynamic Response

ANALYSIS REFERENCE Chapter 5. Algorithm

Initial condition of mode superposition

When initial displacement and initial velocity are given, the initial displacement i0 and initial velocity i0 in the
modal coordinate system is defined as follows. Using all modes gives an equation, and using partial modes
gives an approximate relationship.

1 T
i0  i Mu0
mi (5.4.10)
1 T
i0  i Mv 0
mi

i : i th eigenmode shape
u0 : Initial displacement
v0 : Initial velocity

Frequency response analysis calculates the structural response under a load vibrating at a uniform frequency.
4.2 All loads in frequency response analysis are defined in the frequency domain and expressed as a function of
Frequency Response excitation frequency. In other words, the load in frequency response analysis can be expressed using the
following complex harmonic function when the angular excitation frequency is  .

f  t   f   eit (5.4.11)

The response can also be expressed in the same form.

u  t   u   eit (5.4.12)

Using this, the equation of motion is expressed in the following form:

 2M  iC  K  u    f   (5.4.13)

Here, both the load and displacement are expressed as complex numbers. When expressing the complex value
using magnitude/phase angle, the magnitude represents the maximum load or displacement within the
vibration period and the phase angle is the position (angle) at which this maximum value occurs. On the other
hand, when expressing the complex value using real component/imaginary component, the real component is
the load or displacement magnitude at the starting point of the vibration period and the imaginary component
is the load or displacement after 1/4 period (  / 2 ). Hence, the imaginary component changes with the
vibration period. The relationship between magnitude/phase angle and real component/imaginary component
is as follows:

Section 4. Dynamic Response | 319

 ANALYSIS REFERENCE Chapter 5. Algorithm

u  ur2  ui2 : Magnitude

  tan (ui / ur )
1
: Phase angle
ur  u cos : Real component
ui  u sin  : Imaginary component

Direct frequency response analysis

When using the direct solver for direct frequency response analysis, solving the simultaneous equation (5.4.13)
gives the frequency response u   . If there is no damping, equation (5.4.13) is a real number simultaneous
equation. But if there is damping, it is a complex number simultaneous equation. The solution can be found
accurately using the direct method, but calculation is very inefficient for large problems or when many
frequencies exist because the simultaneous equation needs to be recomposed and solved for each frequency.

Response spectrum analysis is a method of evaluating the structural response due to base motion (uniform
4.3 shaking of nodes confined by the boundary condition), especially earthquakes, and is the most generalized
Response Spectrum method for seismic design. This method assumes a linear system response and only evaluates the maximum
response. Hence, analysis using time integration outlined in sections 5.4.1 and 5.7 is appropriate for problems
that have dominant nonlinearity or when results considering the simultaneity of a particular time step are
important.
The maximum response is evaluated as a mode combination, which reflects the mode participation rate on the
modal response corresponding to the predefined spectrum function. Here, because simultaneity of the modal
maximum response is not considered and the response itself is calculated as a combination, the response
spectrum analysis results can be seen as an approximate solution for time integration. Hence, if the spectrum
function is defined for a particular acceleration or particular seismic wave, the response spectrum analysis
result obtains an approximate maximum value of the linear transient response analysis result for the input
acceleration. However, the analysis results for seismic design are more generally obtained using the design
response spectrum, made from the statistical historical seismic waves in a particular region or country.

Modal spectrum response

The static equilibrium equation for response spectrum analysis is shown in equation (5.3.6), and the maximum
modal response can be expressed using the spectrum data as follows:

imax  max i (t )  i S D i ,  i 

imax  max i (t )   i SV i ,  i  (5.4.14)
imax  max i (t )   i S A i ,  i 

SD i ,  i  : Displacement spectrum data

SV i ,  i  : Velocity spectrum data

320 | Section 4. Dynamic Response

 ANALYSIS REFERENCE Chapter 5. Algorithm

S A i ,  i  : Acceleration spectrum data

th
i : Participation factor of the i mode

Substituting equation (5.4.14) into equation (5.3.5) can express the contribution of the maximum modal
displacement, velocity, acceleration as an equation of spectral data.

uimax  i i S D i ,  i   i i S A i ,  i  / i2

vimax  i i SV i ,  i   i i S A i ,  i  / i (5.4.15)
amax
i  i i S A i ,  i 

A point in the spectrum data is defined as the absolute maximum modal response value of the natural period
(natural frequency), and the effects of the modal damping ratio is included. Because the maximum response of
each period is very diverse for the response spectrum of a particular acceleration history, it is expressed as a
very complex graph form. However for the design response spectrum, a simple line combination in log scale as
shown in figure 5.4.1 is generally used:

500
Figure 5.4.1 Example of
Pseudo Acceleration, S A

acceleration response spectrum 300

200

  0.02
100
  0.05
70   0.10
50

30
102 101 100 101
Period(sec)

Modal combination method

The maximum physical quantities for each mode (maximum value of each displacement, stress, member force,
reaction force etc. component) is called Rimax If the actual maximum physical quantity is assumed as the sum
.
of the maximum values of each mode, simply adding the maximum values of each mode is sufficient. However,
because the maximum values of each mode cannot be guaranteed to occur at the same time step, the
maximum value cannot be found using only simple linear superposition.

Section 4. Dynamic Response | 321

 ANALYSIS REFERENCE Chapter 5. Algorithm

N
Rmax   Rimax (5.4.16)
i 1

Hence, a modal combination method needs to be introduced to evaluate the maximum value approximately.
Various modal combination methods that consider the superposition characteristics or damping effects have
been introduced, but because there is no definite method that gives an appropriate value for all cases, the
characteristics of each modal combination method need to be understood.

► Summation of the absolute value (ABS)

N
Rmax   Rimax (5.4.17)
i 1

This method assumes that all modal responses have the same phase and judges all absolute maximum modal
values to occur at the same time. Hence, it provides the largest value.

► Square root of the summation of the squares (SRSS)

R 
N
max 2
Rmax  i (5.4.18)
i 1

This method provides appropriate results when each mode is sufficiently separated:

► Naval research laboratory method (NRL)

 R 
N
max 2
Rmax  Rmmax  i (5.4.19)
i 1,i  m

This method removes one mode ( m ) that has the maximum absolute value from the SRSS method, and like
the SRSS method, this method provides appropriate results when each mode is sufficiently separated.

Because these methods above are effective only when the modes are sufficiently separated and not adjacent,
the US Nuclear Regulatory Commission (NRC) regulatory guide 1.92(1976) suggests appropriate evaluation
methods for maximum values when multiple modes are adjacent.

► Ten percent method (TENP)

322 | Section 4. Dynamic Response

ANALYSIS REFERENCE Chapter 5. Algorithm

N  i 1 
Rmax   R i
2
 2 Ri R j  (5.4.20)
i 1  j 1 

This method includes effects of all adjacent frequency modes within 10% of the SRSS. Here, the frequencies of
two modes i, j ( j  i) are judged to be adjacent within 10% frequency if the following condition is satisfied:

i   j
 0.1 (5.4.21)
i

► Complete quadratic combination method (CQC)

N i
Rmax   R  R
i 1 j 1
i ij j (5.4.22)

Here,  ij is the cross-correlation coefficient, which is defined as follows:

8  i j ( i  rij m )rij 3 2
ij  (5.4.23)
(1  rij )  4 i j rij (1  rij 2 )  4( i2   j2 )rij 2
2 2

rij : Frequency ratio (  j / i ),  j  i

If i j in equation (5.4.23), ij  1 regardless of the damping ratio. If the damping ratio is '0', ij  1 for all
nodes and the results are the same as the SRSS results. When the damping ratios of two modes are identical, it
can be simplified to equation (5.2.24)

8 2 (1  rij )rij 3 2
ij  ( i   j   ) (5.4.24)
(1  rij 2 ) 2  4 2 rij (1  rij )2

Sign of modal combination result

Because modal combination methods are displayed as absolute values of the mode results, all response
spectrum results always have a positive (+) value. However, for directional results such as reaction force or
deformed shape, appropriate signs need to be applied. The most general method for determining the sign of
the combined results is following the sign of the major mode. The major mode is defined as the mode, out of
the modes that have the largest mass participation rate for each directional component, that has the most
closest direction to the defined spectrum direction (load direction).

Section 4. Dynamic Response | 323

 ANALYSIS REFERENCE Chapter 5. Algorithm

Spectrum data correction

Spectrum data is a function form for the natural frequency and modal damping ratio, as shown in equation
(5.4.14). However, because the user cannot know the frequency before analysis, the spectrum data is defined
as a table with a constant interval. Hence, the interpolation is used when reading the spectrum value of the
applicable frequency or period of the structure and linear interpolation on a logarithmic scale, which expresses
the spectrum response for natural period change, is most generally used. When entering the spectrum data for
multiple damping ratios, linear interpolation on a logarithmic scale is performed in the same way for the
structural modal damping ratio.
However when spectrum data is available for only one damping ratio, there is no data for interpolation and a
special interpolation method is needed for that single damping ratio. The Japan specifications for highway
bridges (2002) suggest the following correction factor for the damping ratio:

1.5
CD     0.5 (5.4.25)
40  1

Figure 5.4.2 Correction factor for 2.0

damping ratio
Correction factor, CD

1.5 1.5
CD     0.5
40  1
A  0.05, 1.0 
1.0

0.5

0.0
0.0 0.1 0.2 0.3 0.4 0.5
Damping ratio, 

When the damping ratio is '0.05', CD  1 (Point A ) and equation (5.4.25) connotes the correction factor when
the damping ratio of the spectrum data is '0.05'. Hence, when the damping ratio (  spectrum ) of the spectrum
data is not '0.05', the ratio of correction factors corresponding to each damping is applied as the final damping
correction factor, as shown in equation (5.4.26).

CD  i 
Rimax  Rimax (5.4.26)
CD  spectrum 

324 | Section 4. Dynamic Response

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 5 Nonlinear Finite Element Solution

Nonlinear Finite Element Solution is a method of converging the accumulated incremental solution from
iterative calculations to the correct solution, and it is processed as shown in figure 5.5.1.

Figure 5.5.1 Accumulated

incremental solution and f
nonlinear finite element
convergence t t
fext
K i 1
gi
f int,i

t
f ext u i  ui 1

ui1

t t  t
u u u

In the figure, t fext and t t fext each represent the external forces at time t and time t   t , and the solution
and incremental solution between time t and time t  t can be expressed as the following relationship:

t t
u t u  u (5.5.1)
u : Incremental solution occurring at time increment  t

If iterative calculation is performed for nonlinear analysis in the time increment  t , the accumulated
incremental solution is as follows:

n
u    ui or ui1  ui   ui 1 (5.5.2)
i 1

ui : Accumulated incremental solution up to i th iterative calculation

 ui1 : Incremental solution occurring at i  1 th iterative calculation

 ui1 is calculated from the following linear simultaneous equation using the tangential stiffness matrix K i1 .

Section 5. Nonlinear Finite Element Solution |325

 ANALYSIS REFERENCE Chapter 5. Algorithm

 ui1  K i11 g i (5.5.3)

gi : Residual force, unbalanced force

t t
The unbalanced force gi is expressed as the following difference between external force fext and internal

force fint ,i .

gi  t t fext  fint ,i (5.5.4)

Equations (5.5.2)-(5.5.4) are iterated until it satisfies the user specified convergence criteria, and the
convergence criteria judges using the change in member force, displacement or energy etc.

Line search
FEA NX provides the line search function to improve the performance of the basic iterative solutions explained
above. The fundamental concept of line search is the introduction of a scalar value  during the process of
adding the calculated incremental solution  ui1 to the accumulated incremental solution for improved
accuracy. In this case, the accumulated incremental solution is calculated as follows:

ui1  ui   ui1 (5.5.5)

Assuming that the calculated ui1 above satisfies the equilibrium state and uses the principal of stationary
total potential energy, the line search problem results in finding the  at which the derivative of the total
potential energy for  is '0'.

s( )   ui 1T g( )  0 (5.5.6)

326 | Section 5. Nonlinear Finite Element Solution

 ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.5.2 Conceptual diagram

of line search algorithm

Potential energy
Acceptable range


tan (s( ))
1

Exact solution
 uTi1g( )  0

Assuming linear change for the energy derivative s( ) about  , the  that satisfies equation (5.5.6) is
calculated as follows:

s(  0)
 (5.5.7)
s(  1)  s(  0)

Here, the slopes at which  is '0' or '1' can be expressed as follows:

s (  0)   u i 1T g i
(5.5.8)
s (  1)   u i 1T g i 1

Because the assumptions made for the line search algorithm are not accurately satisfied for the real case, the
s( ) calculated from equation (5.5.7) is generally not '0'. In FEA NX the processes outlined above are
repeated until the s( j ) / s(  0) value is below the user-specified constant value.

Initial stiffness, Newton Raphson, Modified Newton Raphson

The iterative methods in nonlinear analysis can be classified into the Initial stiffness, Newton Raphson,
Modified Newton Raphson methods depending on the calculation point of the tangential stiffness. The Initial
stiffness method continuously maintains the tangential stiffness calculated at the start point of analysis. The
Newton Raphson method recalculates the tangential stiffness for each iterative calculation. The Modified
Newton Raphson calculates the tangential stiffness where a change in external force occurs. Because
tangential stiffness, matrix calculation, and matrix decomposition requires a long calculation time, using the

Section 5. Nonlinear Finite Element Solution |327

 ANALYSIS REFERENCE Chapter 5. Algorithm

Initial stiffness and Modified Newton Raphson methods is faster than the Newton Raphson method when
problems do not occur in the convergence process. FEA NX does not classify the Initial stiffness and Modified
Newton Raphson methods explicitly. Defining the tangential stiffness recalculation point can give the effects
of all iterative methods.

Automatic stiffness matrix recalculation

It is important to select an appropriate calculation point for tangential stiffness, depending on the
characteristics of the target analysis model such as nonlinearity, evenness of the converging solution, etc. FEA
NX provides the automatic tangential stiffness update, which judges an appropriate recalculation point by
considering the overall characteristics of the nonlinear problem such as convergence characteristics or
determination of divergence etc., as a nonlinear finite element solution. Tangential stiffness update is
performed when the following conditions are satisfied:

► When the expected number of iterative calculations is larger than the user defined maximum number
► When the solution is determined to diverge

Convergence condition
The convergence of the iterative solution is judged using the force norm, displacement norm and energy norm.

giT gi
Force norm ratio  (5.5.9)
fint,iT fint,i

 uiT  ui
Displacement norm ratio  (5.5.10)
uiT ui

 u iT g i
Energy magnitude ratio  (5.5.11)
u i T fint,i

For general nonlinear systems, all convergence norms decrease simultaneously as the system converges.
Particularly, the force norm represents the size of the unbalanced force and has the closest relationship with
the degree of satisfaction of the nonlinear equation. On the other hand, the displacement norm represents the
size of the incremental solution and is not appropriate as a single convergence norm for problems with a very
large local stiffness, such as for systems using the penalty method.

FEA NX compares a single or multiple norms out of these three norms to the user-provided tolerance to
determine convergence.

Divergence determination and load bisection

Determining the divergence for a solution is an important criterion used in automatic tangential stiffness
update, and the divergence rate Ei is determined fundamentally.

328 | Section 5. Nonlinear Finite Element Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

 ui T g i
Ei  (5.5.12)
 ui T gi1

When the absolute value of the divergence rate is larger than '1' ( Ei  1 ), the nonlinear analysis solution is
judged to have possible divergence and necessary measures are taken on the algorithm, such as recalculation
of the stiffness matrix or load bisection.

Load bisection is applied when the increment of the current load step is too large to obtain a converging
solution, such as when the solution diverges or when the number of required iterative calculations is larger
than the user defined maximum number etc. By restarting the iterative calculation through bisecting the
current load increment, an inappropriate load increment size can be dealt with flexibly. FEA NX performs load
bisection automatically until the user defined maximum bisection level is reached.

Automatic time increment adjustment

To increase the efficiency of nonlinear analysis, FEA NX includes a function that automatically adjusts the time
increment size as the base of the nonlinear analysis astringency. The fundamental time increment size and
maximum increment size is determined by the user input. When using the automatic time increment
adjustment function in nonlinear analysis, the time increment size of a particular time step increases or
decreases based on the number of iterative calculations needed for convergence in the previous increment
step.

t i1  ns t i (1  ns  ns ,max ) (5.5.13)

Here, the increment adjustment factor ( ns ) is limited to natural numbers to obtain the maximum number of
nonlinear solutions in the user intended point or load size. The increment adjustment factor has a range from
the minimum value '1' representing the initial increment, and the maximum value ( ns ,max ) provided by the
user.

Quasi-Newton method
The quasi-Newton method is a type of nonlinear solution is the generalized form of the secant method. It
maintains the advantages of the Modified Newton Raphson method, which recomposes the stiffness matrix
only when a load increment is present, and improves the problem of low astringency. In other words, costs do
not occur for recomposition during iterative calculation of the stiffness matrix and effective calculation using
the decomposed stiffness is possible. At the same time, this method can be used to improve the astringency
and general performance.

Section 5. Nonlinear Finite Element Solution |329

 ANALYSIS REFERENCE Chapter 5. Algorithm

6
FEA NX uses the BFGS (Broyden-Fletcher-Goldfarb-Shanno) method , a type of quasi-Newton method. The
inverse matrix of the stiffness matrix from iterative calculations in nonlinear finite element analysis is adjusted
by the following BFGS update process.

K j 1  Γ jT K j 11Γ j  z j δ j δ jT (5.5.14)

Here, j represents the BFGS update index, and matrix Γ j and scalar z j can de expressed as follows.

Γ j  I  z j γ jδ jT
1 (5.5.15)
zj 
δ jT γ j

Also, the quasi-Newton vectors δ j and γ j are expressed using the increment solution  ui , which applies
the line search factor  from the i th iterative calculation, and the difference between unbalanced forces
during iterative calculation as follows.

δ j  u i  u i 1   u i
(5.5.16)
γ j  g i  g i 1

The i th incremental solution during iterative calculations is calculated using the j th BFGS updated stiffness
matrix and unbalanced forces, as shown below:

 ui  K j 1gi1  Γ jT K j11Γ j gi1  z j δ j δ jT gi1 (5.5.17)

The inverse matrix of the stiffness matrix is not actually modified by the BFGS update process; the incremental
solution is calculated during iterative calculations using a recursive method. In other words, it maintains the
decomposed form of the initial stiffness matrix with no BFGS updates. The incremental solution can be found
using simple recursive vector operations. The quasi-Newton vector is saved for these operations. The saved
vector is erased when incremental analysis converges and the stiffness matrix is recomposed.

Arc-length method
Figure 5.5.3 displays the various displacement load paths, including the unstable equilibrium path. When
performing static nonlinear analysis for these phenomena, analysis of the unstable static equilibrium state
after the limit point cannot be performed when the general load controlled nonlinear solution is used. In other
words, general nonlinear solutions cannot find the converging solution after the limit point. When using the

6
Matthies, H. and Strang, G., “The solution of nonlinear finite element equations,” International Journal for Numerical
Methods in Engineering, Vol. 14, Issue 11, pp. 1613-1626, 1979

330 | Section 5. Nonlinear Finite Element Solution

 ANALYSIS REFERENCE Chapter 5. Algorithm

displacement controlled method, the analyzable region increases locally, but this is not a general solution and
tracing is impossible for the snap-back phenomenon. In this case, the arc-length method can be used. The arc-
length method can successfully trace the equilibrium path even when the static equilibrium state includes an
unstable region.

Figure 5.5.3 Various unstable

equilibrium paths
Snap-back
Snap-through

Load
Load

Displacement Displacement

Brittle failure
Ductile failure
Load

Load

Displacement Displacement

The external forces in the arc-length method are assumed to be proportional to the load parameter  , which is
an independent scalar variable. Hence, the arc-length method can be seen to increase the DOF of the
fundamental finite element problem by '1'. However, because the algorithm is composed such that the
parameter  and accumulated incremental solution satisfy the arc-length constraint, the final number of DOF
is maintained. The unbalanced forces including the load parameter can be expressed as follows:

gi (ui , i )  ifext  fint ,i (5.5.18)

Here, linearizing the condition that the unbalanced forces that occur at the i  1 th iterative calculations due to
the incremental solution  ui1 and incremental load parameter i1 is 0. The relationship between the
incremental solution and incremental load parameter can be obtained as follows:

 ui1  K i11 (gi  i1 fext ) (5.5.19)

Using this, the accumulated incremental solution at the i  1 th iterative calculation is as follows:

Section 5. Nonlinear Finite Element Solution |331

 ANALYSIS REFERENCE Chapter 5. Algorithm

ui1  ui   u  i1uT (5.5.20)

 u  K i11gi : Incremental solution for unbalanced force
1
uT  K i 1 fext : Displacement generated for total external force

7
FEA NX uses the Crisfield, Riks, or Modified Riks method arc-length constraints. The Crisfield method is used
as the default arc-length constraint:

ui1T ui1  l 2 (5.5.21)

l : Arc-length

The incremental load parameter i1 can be calculated from the equation above, and substituting this can
calculate the i  1 th iterative calculation solution. Like the general nonlinear solution, this process is repeated
until the user specified convergence criteria is satisfied and the convergence criteria are the same as that of
the general nonlinear solution. In other words, convergence is judged using the change in member force,
displacement, or energy.

When using the arc-length method, accurate load state calculation can be difficult because the load increment
is determined by the arc-length constraint condition and cannot be controlled by the user. Hence, the
applicable range of the arc-length method is limited to problems that need tracing of the unstable equilibrium
state. No additional advantages exist for general nonlinear problems.

Over-relaxation method

The over-relaxation method is one of the methods to improve the convergence rate by multiplying the
estimated unbalance force by the coefficient at the iterative calculation.

Although it is a basic approach rather than a line search method, it is a method that is very similar to the initial
stiffness method because the formula is very simple and unlike the line search method, the additional analysis
time is required in the iterative calculation.

Ke ui  Ke ui1  r  ui 1  (5.5.22)

Where ω is the excess relaxation coefficient.

7
Crisfield, M.A., “An arc-length method including line searches and accelerations,” International Journal for Numerical
Methods in Engineering, Vol. 19, Issue 9, pp 1269-1289, 1983

332 | Section 5. Nonlinear Finite Element Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

fn1

Figure 5.5.4 Over-relaxation fn

method
  1.0

The initial relaxation coefficient is directly input by the user, 1.2 is defined as the default value, and should not
exceed 2.0 at maximum.

Enhanced predictor

The initial displacement estimation method is a method of predicting the initial displacement at the present
stage using the load factor ratio taken from the present stage divided by previous stage and multiplied by the
displacement result of the previous stage as shown in the following equation.

n
u npredictor  u n1
n1 (5.5.23)

The estimated displacements do not exactly coincide with those of the current step, but they are useful for
iterative calculations because they predict closer results than the displacements estimated by elastic stiffness.

Particularly, it is more effective when the material model has a large plasticity. However, because it is less
accurate than the estimated tangent stiffness by Newton-Raphson, it is recommended to use it with the initial
stiffness method.

In general, the initial stiffness method is a stable method for solving the problem.

Section 5. Nonlinear Finite Element Solution |333

 ANALYSIS REFERENCE Chapter 5. Algorithm

f
Figure 5.5.5 Estimation of Initial Initial stiffness
Deformation (enhanced fn1
predictor)

fn Enhanced Predictor stiffness

f n1

334 | Section 5. Nonlinear Finite Element Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

Section 6 Strain/Stress Measurement

Considering Large Deformation
For geometric linear analysis, the strain and stress are defined without considering the shape difference before
and after deformation. The strain in geometric linear analysis is generally defined as follows:

1 u u T
ε [ ( ) ] (5.6.1)
2 X X

u : Displacement

X : Coordinates before or after deformation

For geometric nonlinear analysis that considers large deformation, the strain can be defined using various
methods and a corresponding stress exists for each strain to define virtual work.

Definition of strain
Strains that consider large deformations include Green strain, Green-Lagrange strain and rate of deformation
or strain rate. The Green strain tensor E is defined as follows:

ds 2  dS 2  2dX  E  dX (5.6.2)
X : Coordinates of a particular position on the structure before deformation

The Green strain can be seen as the difference between the squared value of the differential length before
deformation dS and after deformation ds . The Green strain tensor can be defined using the deformation
gradient as follows:

1 T
E (F  F  I ) (5.6.3)
2

If only rigid motion occurs, Green strain does not occur and it is appropriate as a measurement of deformation.

The rate of deformation D is defined by the velocity gradient as follows:

1 v
D (L  LT )  sym[ ] (5.6.4)
2 x

Section 6. Strain/Stress Measurement Considering Large Deformation | 335

 ANALYSIS REFERENCE Chapter 5. Algorithm

v : Velocity vector

X : Coordinates of a particular position on the structure before deformation

In other words, it corresponds to the symmetric part of the velocity gradient tensor. The rate of deformation
can be seen as a value for the squared differential length.

ds 2 (5.6.5)
 2dx  D  dx
t

The rate of deformation also does not occur when only rigid motion exists, and it has the following relationship
with the Green strain:

D  FT  E  F1 (5.6.6)

Because the rate of deformation is a rate of change with time, it is generally time integrated and used as a
strain. If analysis that considers geometric nonlinearity of a material is performed in FEA NX the strain is
computed by time integrating the rate of deformation. FEA NX uses the rate of deformation or strain rate.

Definition of stress
When geometric deformation is large, the stress can also be defined using various methods. FEA NX uses the
Cauchy stress.

Because Cauchy stress is a value that satisfies the equilibrium equation of the current shape, it is also known as
true stress σ and is defined as follows:

n  σ  d  df  td (5.6.7)

The Cauchy stress of the shape before deformation can be converted to the 2nd PK stress (Piola-Kirchhoff: S )
as follows:

S  JF 1  σ  F T (5.6.8)

336 | Section 6. Strain/Stress Measurement considering large deformation

 ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.6.1 Shape before/after n

deformation and force direction n0
for defining stress df
df

F  1d f
d
d0

Reference configuration Current configuration

The physical meaning of 2nd PK stress is not clear, but it is useful in describing the equation of motion when
coupled with the Green strain and so it is often used to define the behavior or materials with an energy
potential such as rubber. FEA NX uses the following stress and strain integration method for all materials:

Stress rate and strain rate integration

Elasto-plastic materials, viscoelastic materials etc. do not have energy potential, but use a constitutive
relationship consisting of the strain rate and objective stress rate. The Jaumann stress rate used is defined as
follows:

σ J  σ  w  σ  σ  wT (5.6.9)

Strain rate and objective stress rate have the following relationship from the constitutive equation of the
material:

σJ  C : D (5.6.10)

Reflecting the central difference and considering the structural rotation in equation (5.6.10), the equation that
calculates the stress at step n 1 using the stress and strain increments calculated at step n:

σ n1  R  σ n  RT  C : ε (5.6.11)

Section 6. Strain/Stress Measurement Considering Large Deformation | 337

 ANALYSIS REFERENCE Chapter 5. Algorithm

8
The rotation amount increment R is calculated as follows to satisfy the incrementally objective stress
condition.

1 1
R  ( I  W ) 1 (I  W ) (5.6.12)
2 2

Particularly, the strain increment and incremental spin calculation is performed on the n  1 / 2 shape.

1 u u T 1 u u T

ε  ( [ ] ), w  ( [ ] ) (5.6.13)
2 x n1/2 x n1/2 2 x n1/2 x n1/2

The strain rate integral is also calculated using equation (5.6.11) and uses the structural rotation amount
increment R .

8
Hughes, T.J.R. and Winget, J., “Finite rotation effects in numerical integration of rate constitutive equations arising in
large deformation analysis,” International Journal for Numerical Methods in Engineering, Vol. 15, 1980

338 | Section 6. Strain/Stress Measurement considering large deformation

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 7 Nonlinear Dynamic Response

FEA NX supports nonlinear time history analysis that includes geometric, material nonlinearity and it is based
on implicit time integration.

The dynamic equilibrium equation in nonlinear time history analysis uses the HHT   method as implicit time
7.1 integration, just like for linear time history analysis, and uses the following modified equilibrium equation.
Implicit Time

Integration t
 Mv n1   (1   H ) Cv n1  f int,n1  f ext ,n1    H Cv n  f int,n  f ext ,n   0 (5.7.1)

In nonlinear time history analysis, the effects of the mass matrix rotation due to geometric nonlinearity are
considered. The rotational inertia part of the mass matrix is modified for each iterative calculation, according
to the finite rotation of the nodes, and the inertial force generated from the rate of change of the mass matrix
is considered in analysis.

Nonlinear time history analysis calculates the convergence solution for each time step using the nonlinear
finite element solution in section 5.5. The unbalanced forces are expressed from equation (5.7.1) as follows:


g n1 
t
 M n1v n1   (1   H ) Cn1v n1  f int,n1  f ext ,n1    H Cn1v n  f int,n  f ext ,n  (5.7.2)

The tangential stiffness matrix can be found by applying the time step equation (5.4.3), from the Newmark
method for velocity and acceleration, onto the unbalanced force and differentiating for the displacement DOF,
as follows:

1  (1   H ) n1 f int,n1 f ext ,n1

A M n1  M n1  C  (1   H )  (1   H ) (5.7.3)
t 2
t t u u

Here,  H  0.05 is used as the default value, just like linear time history analysis. Also, to secure
unconditional stability, the following values are used :   (1  2H ) / 2 ,   (1   H )2 / 4

Angular velocity and angular acceleration

When considering geometric nonlinearity in nonlinear time history analysis, the angular velocity and angular
acceleration need to be updated by reflecting the effects of body axis system rotation. Defining the body axis

Section 7. Nonlinear Dynamic Response | 339

 ANALYSIS REFERENCE Chapter 5. Algorithm

system where finite rotation occurs as  , the Newmark time step equation in the coordinate system of the
body axis system is as follows:

ωn1  ωn  t  αn1  (1   )αn  (5.7.4)

ω , α : Angular velocity and angular acceleration about the body axis system 

Using the base vector perpendicular to the body axis system e , the equation above can be expressed for the
GCS.

ωn1  t αn1   en1en  ωn  t (1   )αn  (5.7.5)

The product of the perpendicular base vectors in equation (5.7.5) is the same as the incremental rotation
matrix found below:

C  exp(ˆ) (5.7.6)

ˆ : Skew symmetric matrix for the rotation amount increment

The rotation amount increment can be expressed by the Newmark method using the incremental rotation
matrix, as shown below.

 1  
θ  t 2  α n1  C tωn  t 2     α n  (5.7.7)
 2  

Rearranging equation (5.7.7) for angular velocity and angular acceleration respectively and substituting into
equation (5.7.5) gives the following updating equations for both angular velocity and angular acceleration.

       n
ω n1  θ  C 1   ω n  t 1  α 
t      2  
(5.7.8)
1  1 n  1  n
α n1  θ  C  ω  1  α 
t 2   t   2  

340 | Section 7. Nonlinear Dynamic Response

ANALYSIS REFERENCE Chapter 5. Algorithm

Damping effect
Mass proportional damping and stiffness proportional damping are also considered in nonlinear time history
analysis, just like linear time history analysis. In this case, the damping matrix is composed similarly to
equation (5.4.5). The mass matrix, which is used to calculate the damping matrix in nonlinear time history
analysis, considers the rotational effects due to finite rotations and the stiffness matrix only uses the stiffness
matrix due to material nonlinearity.

C   ej Mej   jeKemat , j  B (5.7.9)

Kmat : Stiffness matrix due to material nonlinearity

Section 7. Nonlinear Dynamic Response | 341

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 8 Contact Condition

Contact analysis fundamentally assumes that two objects in a space can be in contact, but cannot penetrate
each other (non-penetration condition), and is nonlinear in behavior or in condition from a physical point of
view. The type of contacts are general contact (considers the impact and impact friction between two objects
in analysis) and rough contact (does not consider sliding) shown in figure 5.8.1, and welded contact (two
objects are welded from the start of analysis) shown in figure 5.8.2. Here, the welded contact is assigned
depending on the position of two objects at the start of analysis and can be seen as linear.

Figure 5.8.1 Concept of general

contact and rough contact

Rough contact
General contact

Figure 5.8.2 Concept of welded

contact

Welded contact

342 | Section 8. Contact Condition

ANALYSIS REFERENCE Chapter 5. Algorithm

Relationship between contact condition and analysis type

The contact condition can be used for initially adjacent bodies in structural analysis, consolidation analysis and
seepage analysis.

The contact can be classified with the node-to-surface contact, or surface-to-surface contact. Node-to-surface
contact takes less time, but the solution accuracy is relatively low because the nodes of the main object tend
to penetrate through the sub object. On the other hand, surface-to-surface contact takes longer but the non-
penetrating conditions are satisfied relatively accurately, allowing the accurate simulation of the structural
behavior. FEA NX supports the surface-to-surface contact.

The general contact can be used in nonlinear structural analysis (static, dynamic) and consolidation analysis.
The general contact corresponds to the nonlinear condition, and its behavior is different according to the
geometric nonlinearity consideration in terms of analysis techniques. In case of considering geometric
nonlinearity, the possibility of contact is considered for all master segments under the assumption of large
displacement. On the other hand, only the contact closed to within initial user-defined distance between
master segment and slave node is considered.

Contact plane search

Contact search uses the slave node/master segment algorithm. This algorithm determines contact by the
adjacency between the slave node and master segment, or how much the slave node penetrates the master
segment. Generally, the order of the slave node defined object and master segment defined object does not
matter. But from a numerical point of view, the master segment needs to be defined on the object with a
relatively larger stiffness, or relatively element-sparse object, to obtain more accurate analysis results.

To determine the actual contact of the slave node and master segment, the global search process is conducted.
Global search is the process that determines the preliminary slave nodes where objects or segments can
collide in space. Contact search is performed locally for the slave node and master segment sets determined
by global search.

To determine whether nodes and planes are actually in contact, the slave nodes need to be projected
orthogonally on the master segment as shown in figure 5.8.3. Defining vector r from the origin to the
projected point (A), and vector xs from the origin to the slave node gives the following equation:

r c ,c 
  x s  r c ,c   0

(5.8.1)
r c ,c 
  x s  r c ,c   0


Section 8. Contact Condition | 343

 ANALYSIS REFERENCE Chapter 5. Algorithm

Here, c ,c  is the position of contact point (A) on the master segment, expressed in the natural coordinate
system. The c ,c  that satisfies the equation above can be calculated numerically using the Newton-
Raphson method. The coordinate increment  c , c  for applying the Newton-Raphson method is as
follows:

 r   r 
    r r       
 
 r   
  
     r 
r c ,c   xs  (5.8.2)
     
   

Using the initial condition as c ,c    0,0 , the equation above converges easily when the position of
contact point (A) or slave node (B) is not far from the master segment. If the next slave node is checked and
determined to have penetrated the contact plane, the force (contact force) proportional to the penetration
depth is added to the slave node and contact plane.

Figure 5.8.3 Normal relationship

between slave node and master
segment


nA

Calculation of contact force

The displacement relationship between the slave node and master segment, which are determined to be in
contact, is confined using the penalty method. In FEA NX, the gap and contact force is defined using the
following equations (5.8.3), (5.8.4):

g N   xB  x A   n A (5.8.3)
f C  k N g N if g N  0 (5.8.4)
kn : Penalty coefficient

x A , xB : Position vectors of point (A) on master segment and slave node (B)

344 | Section 8. Contact Condition

ANALYSIS REFERENCE Chapter 5. Algorithm

nA : Normal vector of point (A) on master segment

The penalty coefficient k N has the effect of applying an elastic stiffness between the master segment and
slave node, and the non-penetrated condition is satisfied approximately depending on the size. For the welded
condition or the sliding contact condition, the contact force is assigned even when g N is positive and the
initially adjacent segment and slave nodes are not separated and its effects reflected. FEA NX automatically
calculates the penalty coefficient using the following equation (5.8.5):

fs M i
plane/shell elements : ki  Asi
h (5.8.5)
fK
solid elements : ki  s i Asi
h
fs : Proportionality coefficient

Ki : Bulk modulus

Mi : Coefficient of expansion

Ai : Area

Vi : Volume
1 n Vmi 1 n
h : Length of master segment solid : 
n i Ami
plane/shell :  Ami1 2
n i

The proportionality coefficient f s above is determined differently depending on the analysis type or contact
condition type, and the user can modify it to effectively satisfy the non-penetrated condition or welded
condition. For the welded condition, a resistant force against sliding in the lateral direction of the contact
plane is assigned and its size is as follows:

 u B  u A   t xA 
 
gT    (5.8.6)
 u  u   t y 
B A A

f T  kT gT (5.8.7)
kT : Penalty coefficient

t xA , t Ay : Lateral vector of point (A) on master segment

For the welded condition, the penalty coefficient kT uses the same value as k N .

Section 8. Contact Condition | 345

 ANALYSIS REFERENCE Chapter 5. Algorithm

As equation (5.8.4.) and (5.8.7), the contact force of linear elasticity or resistance force for sliding is not
suitable for the general contact used in nonlinear analysis. In order to obtain the convergence, much oscillation
can be occurred because the stiffness of general contact rapidly changes from 0 to kN according to the sign
of gap. To compensate for this, FEA NX uses the modified contact force as follows:

fC 0 if g N  d1
f0 g g
fC  [( N  1)(exp( N  1)  1)] if  d 2  g N  d1
(exp(1)  1) d1 d1
f0 d d
fC  [( 2  1)(exp( 2  1)  1)]  k N ( g N  d 2 ) if g N   d 2
(exp(1)  1) d1 d1
(5.8.8)

Figure 5.8.4 Relation between gap

and modified contact force fC

fC 0 f C  exponential f C  linear

 d1 d2 gN

There exists more complex nonlinearity in the resistance force against horizontal sliding because this force
occurs only if the vertical force is applied. As a result, the equation (5.8.7) is modified and applied as follows:

kT f C
fT  gT (5.8.9)
kN d2

In the above equation, the discontinuity of force which may occur at the moment of contact suddenly
disappeared can be minimized by proportion to the vertical contact force and horizontal force. In case of
general contact, the friction can be considered additionally.

f  fT   f C
0 (5.8.10)

346 | Section 8. Contact Condition

ANALYSIS REFERENCE Chapter 5. Algorithm

 : Friction coefficient

The behavior which satisfied with the above equation becomes elastic movement, and if it does not satisfy the
above equation due to the large horizontal movement, sliding occurs. When you compare (5.8.9) and (5.8.10),
the horizontal relative displacement within certain distance d 2 shows elastic movement, and the more
relative displacement can be seen that the slip is assumed. As the general contact considering friction creates
an asymmetric stiffness matrix, numerical calculated efficiency is significantly inhibited. Also, the large friction
coefficient (over 0.3~0.4) or horizontal elastic modulus kT can be the factors causing convergence problem.

Section 8. Contact Condition | 347

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 9 Slope Stability Solution

Slope stability for an embankment or excavation is one of the most frequently dealt problems in geotechnical
engineering. The slope always has a self-weight potential energy due to gravity and if external forces such as
pore pressure, applied load, earthquake, wave force etc. act on the slope, its stability is greatly affected. Here,
slope failure can occur when the active force of the slope is greater than the resistant force of the soil. Slope
stability analysis evaluates the stability against failure using the relationship between the active force and
resistant force of the slope.

The widely used limit equilibrium method can only evaluate the stability for the basic given conditions.
However, the actual collapse of the ground generates a large local deformation and fails at the limit. Hence,
establishing an analysis method that can trace the deformed shape continuously from the initial deformation
to collapse is important for stability analysis of the ground. Recently, active research is ongoing to apply finite
element method as a method of evaluating stability such as slope failure, using the strengths of finite
elements such as the ease in checking the deformed shape, even at various load and boundary conditions.

Figure 5.9.1 Slope failure process STEP #1 : Stable Ground Status STEP #2 : Input the external Force

pore water pressure, applied

Load, earthquake load, wave load

Potential energy (W) Potential energy (W)

(Self Weight) (Self Weight)

STEP #3 : Shear Stress within Slope STEP #4 : Compute Factor of Safety

W Resisting force
Factor of Safety =
T Driving Force

348 | Section 9. Slope Stability Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

The following slope stability analysis methods have been suggested.

► Mass procedure and slice method according to the limit equilibrium theory
► Limit theory according to the rigid-plastic theory
► Finite element method according to the elasto-plastic theory

In FEA NX, the usable slope stability analysis methods that use the finite element method are the strength
reduction method and the stress analysis method based on the limit equilibrium theory.

Slope stability analysis using the finite element method are detailed approximate solutions that satisfy all the
9.1 equilibrium force conditions, compatibility conditions, constitutive equations and boundary conditions of each
Strength Reduction point on the slope. This numerical analysis method can simulate nearly actual failure shapes, reflect the field
conditions better and can analyze the minimum safety factor and failure behavior of the slope in detail.
Method Particularly, the failure process is automatically simulated without any assumptions made to the failure plane
9
of the slope.

The strength reduction method (a slope stability analysis method based on the finite element method)
gradually decreases the shear strength and performs analysis until the point where the calculation does not
converge. This point is considered to be the failure point of the slope, and the maximum strength reduction
ratio at this point is thought of as the minimum safety factor of the slope. This method is costly because it
requires multiple nonlinear analyses, but it can provide more accurate results in a reasonable time for
improved data processing speeds. Also, the strength reduction method can verify the deformation process
from the initial slope to failure without any required assumptions for the failure plane.

Strength reduction theory

To simulate slope failure using the strength reduction method, the safety factor is computed at an arbitrary
point where the Mohr circle is in contact with the failure envelope, as shown in the figure below. The stress
state at this point can be determined as the failure state and when this failure point increases, overall slope
collapse occurs. The finite element analysis at this limit state diverges, and the safety factor at this point is
defined as the minimum safety factor.

9
Griffiths, D.V. and Lane, P.A. (1999). Slope Stability Analysis by Finite Elements, Geotechnique, 49(3), 387-403

Section 9. Slope Stability Solution | 349

 ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.9.2 Strength reduction

method

Critical Status Line Critical Status Line /
W A Factor of Safety

Mohr’s Circle for A point

Calculation of minimum safety factor

The material models used in the strength reduction method are Mohr Coulomb, Drucker Prager and Modified
Mohr Coulomb. For the input variables used here, all variables are assumed to have a constant value except
cohesion, friction angle and dilatency angle, which determine shear failure. The cohesion, friction angle and
dilatency angle corresponding to ground elements (plane strain, axisymetric, solid) are gradually decreased
and the safety factor Fs at slope failure is computed.


Fs  (5.9.1)
f

Here,  is the shear strength of the slope material, and can be expressed using the Mohr-Coulomb criteria as
follows:

  c   n tan  (5.9.2)

Also,  f is the shear stress of the active plane and can be calculated as follows:

 f  c f   n tan  f (5.9.3)
c
cf  : Shear strength factor (Cohesion)
SRF
 tan  
 f  tan 1   : Shear strength factor (Friction angle)
 SRF 

SRF : Strength reduction factor

350 | Section 9. Slope Stability Solution

 ANALYSIS REFERENCE Chapter 5. Algorithm

For the strength reduction method, the SRF value just before the non-convergence is evaluated as the safety
factor. Hence, the safety factor can be slightly different depending on the user input number of convergence
and convergence criteria.
Other material models can be included for slope safety analysis, but the strength reduction is not applied to
these models.

Strength reduction using the arc-length method

The main difference between the existing strength reduction method and the method using the arc-length
method is the method of increasing/decreasing the safety factor, which is the standard for strength reduction.
The existing method computes the safety factor of the next step by controlling the safety factor of the current
step by the user defined increment. Hence, ineffective calculation is performed for very stable models or
unstable models without the engineer's judgment because the uniform safety factor is incremented. However,
using the arc-length method, the arc length is computed by the convergence speed of the previous step and
thus a more appropriate safety factor increment can be obtained.
Defining the projected stress on the new failure plane caused by the strength decrease of the in-situ stress
state σ0 as σ1 , and introducing an additional arc-length parameter  to apply the arc-length method, the
stress within the element can be assumed as follows:

σ  (1   )σ0   σ1  σ (5.9.4)

Here, σ is the stress component corresponding to  that is needed to maintain the equilibrium with the
external forces on the failure plane.

Figure 5.9.3 Arc-length load  fos    0 

vector and stress flow diagram
fos  0    1

σ0
σ σ new fos    1
f AL

σ1

Using equation (5.9.4), the unbalanced forces of the i th iterative calculation for nonlinear analysis can be
expressed as a function of the DOF vector and arc-length parameter:

gi (ui , i )  fext    BT (1   )σ 0   σ1  σ  d  (5.9.5)

e
e

Section 9. Slope Stability Solution | 351

 ANALYSIS REFERENCE Chapter 5. Algorithm

The external force including self weight ( fext ) does not change with  . This method is different from the
general arc-length method for unstable equilibrium state analysis in that the internal forces change with the
stress assumption in equation (5.9.4). To apply the Newton-Raphson based nonlinear finite element analysis,
equation (5.9.4) can be expanded for the incremental solution and increment parameters as shown below:

gi1  gi  Ki1 ui1  i1f AL  0 (5.9.6)

g
K i1   : Tangential stiffness matrix
u
g
f AL     BT σ 0 d     BT σ1d  : Arc-length force vector
 e
e
e
e

Using this, the accumulated incremental solution that occurs at the i  1 th iterative calculation is as follows:

ui1  ui   u  i1uT (5.9.7)

 u  K i11gi : Incremental solution for unbalanced force
1
uT  K i 1 f AL : Displacement generated for arc-length load vector

Here, i1 is calculated using the arc-length constraint condition. The related information is already
explained in the nonlinear solution section and is hence omitted. This process is repeated until the user defined
convergence criteria is satisfied. The safety factor of each step can be calculated through this process. This
subsequent process is also repeated until the rate of change for the safety factor is within a certain condition.

Figure 5.9.4 Safety factor

increment for arc-length fosn
force

fosn 1
fext

Arc-length

displacement

352 | Section 9. Slope Stability Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

The limit equilibrium method is one of the most often used slope stability methods for actual design. However,
9.2 this method cannot find the stress history of the actual slope or the change in ground behavior. On the other
Stress Analysis hand, slope stability analysis using the finite element method can consider the slope formation process and
other ground characteristics, but it requires a longer analysis time because it performs multiple nonlinear
Method analyses.
Recently, much research has been done in this area as to the strengths of the using limit equilibrium method
and finite element based slope stability analysis simultaneously. FEA NX provides the slope stability analysis
method that uses the finite element stress analysis results. This method is based on the virtual sliding surface
of the limit equilibrium method and the stress results of stress analysis.
This method computes the safety factor for multiple assumed virtual sliding surfaces using the stress results of
finite element analysis, and the minimum safety factor and corresponding critical section is computed. The
provided ground material models are Mohr Coulomb, Drucker Prager and Modified Mohr Coulomb, just like
the strength reduction method.

Calculation of minimum safety factor

The safety factor used in the finite element method is defined as follows:

Fs 

S
f d
(5.9.8)

S m d

Here,  m is the generated shear stress,  f the shear strength, and for Mohr Coulomb materials, it can be
expressed as follows:

 f  c   n tan 
(5.9.9)
m 
1
2
 y   x  sin 2   xy cos 2

Here, the directional stress normal to the sliding surface  n is as follows:

 n   x sin 2    y cos2    xy sin 2 (5.9.10)

c : Cohesion

 : Internal friction angle of the material

 : Angle between horizontal plane and sliding surface

x ,  y : Normal stress in the x direction and y direction respectively

 xy : Shear stress

Section 9. Slope Stability Solution | 353

 ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.9.5 Stress components of Grid

the slope
Radius
n sliding surface

m
 xy
 x


 xy
y

Stress integration along virtual sliding plane

To calculate the safety factor, the line integral of stresses along the virtual sliding plane need to be performed.
For this, the stress value at an arbitrary position is required. The stress is calculated using the inner product of
the nodal stress and the shape function at that position.

node
σ   Ni σ inode (5.9.11)
i1

Ni : Shape function at node i

 inode : Nodal stress at node i

 : Stress at arbitrary point within the element

Here, the nodal stress is calculated through the stress recovery technique, which uses the nodal average
method. In other words, FEA NX calculates the nodal stress through extrapolation of the integral stresses of
each node sharing element, and the final nodal stress for shared nodes is applied as the average value of the
calculated nodal stresses.

The stress integral along the virtual sliding plane in the 2D GCS is transformed into an integral form in 1D local
coordinate system and calculated using the following equation:

L 1 L nint
    TWi i 
n2
n1  x , y d  
2 1
T   d  
2 i1
(5.9.12)

 : Coordinate variable in local coordinate system

Wi : Constant of integration at integral point i

T : Transformation matrix that transforms stress in the element coordinate system to GCS

354 | Section 9. Slope Stability Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

L : Element length

 : Shear stress  m or shear strength  f on virtual sliding plane

The finalized total safety factor for the given virtual sliding plane is as follows:

nel


n2
 f d
n1
Fs  i 1
nel
(5.9.13)

n2
 md 
n1
i 1

nel : Number of elements passing the virtual sliding plane

n1 : Start point of the virtual sliding plane within the element

n2 : End point of the virtual sliding plane within the element

The stress analysis method based on the limit equilibrium method uses the stress field computed from the
finite element method and the virtual sliding surface of the limit equilibrium method. Hence, it has a stress
distribution and deformed shape from finite element analysis, and is optimized for obtaining the critical
section from the limit equilibrium method. Compared to the strength reduction method, it requires a much
shorter analysis time and can accurately compute various ground or reinforcement members without any
particular assumptions.

Section 9. Slope Stability Solution | 355

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 10 Equivalent Linear Solution

Equivalent linear analysis simulates nonlinearly behaving ground according to strain using linear analysis. This
method simplifies the complex nonlinear ground properties into linear equivalent properties for linear analysis.
FEA NX supports free field analysis for analyzing the in-situ ground behavior before construction. FEA NX also
supports 2D equivalent linear analysis for ground-structure coupled analysis.

Free field analysis finds the ground response for an input load on the in-situ ground state before any
10.1 construction of structures. Free field analysis is used for ground surface vibration prediction to determine the
Free Field Analysis design response spectrum, computation of dynamic stress and strain to evaluate liquefaction, and
determination of earthquake load that causes ground or structural instability.
Method
Free field analysis finds the ground response due to the vertically transmitted shear waves that pass the linear
viscoelastic region. The analysis ground consists of multiple strata that are infinite in the horizontal direction
and a semi-infinite bottom layer as shown in figure 5.10.1. Each stratum is homogeneous and assumed to have
isotropic material properties. The vibrations in the analysis model are caused by shear waves, which penetrate
and reflect the model ground vertically, and displacement only occurs in the horizontal direction. Hence, the
wave equation (5.10.1) needs to be satisfied for all strata.

 2u  2u  3u
  G 2  2 G 2 (5.10.1)
t 2 x x t

u : Horizontal displacement
 : Mass density

G : Shear modulus

 : Hysteretic damping ratio

356 | Section 10. Equivalent Linear Solution

 ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.10.1 Free field analysis Layer Coordinate Propagation

Properties Thickness
model No. System Direction
u1

1 x1 G1 D1 1 h1
u2

x2
... um

m xm Gm D m  m hm
um+1

m+1 xm+1 Gm 1 Dm 1 m 1 hm+1

um+2

xm+2
...

uN

N GN D N  N
xN
hN=∞
Particle motion
Incident wave
Reflected wave

Representing the displacement function as a harmonic function like equation (5.10.2) and transforming
equation (5.10.1) into the frequency domain gives the governing equation (5.10.3), and the stress-
displacement relationship is equation (5.10.4).

u( x, t )  u( x, )eit (5.10.2)
 2
G* u ( x,  )   2u ( x,  )  0 (5.10.3)
x 2

 ( x,  )  G * u ( x,  ) (5.10.4)
x

 ( x,  ) : Shear stress in frequency domain

*
G : Complex shear modulus

10
The complex shear modulus above uses the following equation suggested by Udaka.

G*  G(1  2 2  i2 1   2 ) (5.10.5)

10
Udaka, Takekazu. (1975). Analysis of Response of Large Embankments to Traveling Base Motions, Department of Civil
and Environmental Engineering. Berkeley: University of California, p346.

Section 10. Equivalent Linear Solution | 357

 ANALYSIS REFERENCE Chapter 5. Algorithm

The free field ground model is generally expressed as figure 5.10.1 to find the solution of the 1D wave
transmission equation. The layer boundary number is assigned from the ground surface and expressing the
response of the m th layer as um , the response can be represented using the following function of the depth
xm from the top part of the m th layer:

um ( xm , )  Am ()eikmxm  Bm ()eikmxm (5.10.6)

* *

 m ( xm , )  ik G ( Am ()e  Bm ()e ikm

(5.10.7)
* *
* * ikm xm xm
m m )
*
k m : Wave number of the m th layer
Am : Layer response coefficient of elastic waves transferred upwards in the m th layer

Bm : Layer response coefficient of elastic waves transferred downwards in the m th layer

The following compatibility condition and force equilibrium condition need to be satisfied at the adjacent layer
boundaries:

um ( xm  hm )  um1 ( xm1  0) 

  m  1, 2,...,( N  1) (5.10.8)
 m ( xm  hm )   m1 ( xm1  0) 

Substituting equations (5.10.6) and (5.10.7) into equation (5.10.8) can derive the relationship between the
coefficients as follows:

 Am1  Bm1  Am eikmhm  Bme  ikmhm 

* *

 
 * *
kmGm  ikm* hm 
, m  1, 2,..., ( N  1) (5.10.9)
 m1
A  Bm1  * *
( Am e ikm* hm
 Bm e ) 
 km1Gm1 

Rearranging this to derive the relationship between the response coefficients of adjacent ground layers gives
the following recurrence relationship shown in equation (5.10.10):

1 1
Am1  Am (1   m* )eikmhm  Bm (1   m* )e  ikmhm
* *

2 2 (5.10.10)
1 1
Bm1  Am (1   m* )eikmhm  Bm (1   m* )e  ikmhm
* *

2 2
hm : Thickness of m th layer
 *
m : Dynamic stiffness ratio between adjacent layers

Here, the dynamic stiffness ratio between adjacent layers can be expressed as follows:

358 | Section 10. Equivalent Linear Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

* *
kG
 m*  *
m m
* (5.10.11)
k G
m1 m1

Because shear stress is always '0' at the ground surface, A1  B1 can be known from equation (5.10.7). Hence,
the response coefficient of the m th layer can be found by applying equation (5.10.10) in order from the first
layer.

Am ( )  am ( ) A1 ( )
(5.10.12)
Bm ( )  bm ( ) B1 ( )

Here, a1  b1  1 can be known. The transfer function Hij ( ) between layer boundary i and layer boundary
j is as follows:

ui ( ) ai ( )  bi ( )
H ij ( )   (5.10.13)
u j ( ) a j ( )  b j ( )

If the transfer function H ij ( ) is determined and the response u j ( ) at layer boundary j is given, the
response ui ( ) at layer boundary i can be found using the following equation:

ui ( )  H ij ( )u j ( ) (5.10.14)

Also, frequency domain response can be inverted to the time domain using the FFT (fast Fourier transform)
method.

The main difference between ground-structure interaction problems and general structural dynamics
10.2 problems is the radiation damping phenomena due to the infinite ground domain. Whilst general damping
2D Equivalent Linear properties damp the structural motion through material friction, radiation damping releases wave energy into
the infinite domain of the ground, which contributes to the damping phenomena of structural kinetic energy.
Analysis
Radiation damping is included in the damping term of the equation of motion, and its size is determined by
the wave form of the externally transmitted wave. The wave form can be easily modeled in the frequency
domain and it is efficient to consider this using frequency domain solvers. General ground materials are
fundamentally heterogeneous and the nonlinearity of its mechanical behavior is very severe.

To correctly analyze the ground-structure interaction problem, the radiation damping phenomena above and
important nonlinearity characteristics need to be considered simultaneously. Hence, the frequency domain

Section 10. Equivalent Linear Solution | 359

 ANALYSIS REFERENCE Chapter 5. Algorithm

solver is used for easy modeling of radiation damping and the equivalent linear method is used to analyze
material nonlinearity.

The analysis process uses the frequency domain analysis process that uses FFT (fast Fourier transform) and the
equation of motion composes a combined ground-structural system that has the free field response as the
input motion. Hence, the ground response and structural response are obtained together for one computation.
Interpolation is used to reduce the number of frequencies where the solution to the equation of motion needs
to be found, and this method uses the interpolation of the transfer function solutions of two continuous
frequencies. Hence, the setting of basic frequencies is important and a sufficiently high limit frequency needs
to be defined to secure validity of the ground motion analysis.

As mentioned above, the input motion in the time domain can be converted to the frequency domain using
FFT to calculate the structural response under a vibration load with a constant frequency. All loads in
frequency response analysis are defined in the frequency domain and are expressed as functions of the
assigned frequency. In other words, when the angular excitation frequency is  , the load in frequency
response analysis can be expressed as the following complex harmonic function:

f  t   f   eit (5.10.15)

The corresponding response can also be expressed in the same form:

u  t   u   eit (5.10.16)

Using this, the equation of motion can be expressed in the following form:

 2M  iC  K  u    f   (5.10.17)

Using equation (5.10.5) suggested by Udaka, the equation above can be modified and rearranges into the
following equation:

 M  K  u    f   (5.10.18)
2 *

K*  K(1  2 2  i2 1   2 ) : Complex stiffness coefficient

Energy transmitting boundary

It is difficult to model the nearly infinite ground accurately using the 2D model used for ground-structure
analysis. Hence, the model boundary needs to be set at an engineering appropriate position and the set
boundaries need to be processed to simulate actual site conditions.

360 | Section 10. Equivalent Linear Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

FEA NX uses the energy transmitting boundary suggested by Lysmerand Wass to express the load f ( ) in
11

equation (5.10.18) as follows:

f ()  K*f u f ()  (R  L)(u()  u f ()) (5.10.19)

u f ( ) : Free field analysis displacement of the transmit boundary
*
K f : Complex stiffness matrix of the transmit boundary

R : Stiffness matrix of the right transmit boundary

L : Stiffness matrix of the left transmit boundary

Direct frequency response analysis

When using the direct solver for direct frequency response analysis, solving the simultaneous equation (5.10.18)
gives the frequency response u   . The solution can be found accurately using the direct method, but
calculation is very inefficient for large problems or when many frequencies exist because the simultaneous
equation needs to be recomposed and solved for each frequency. To supplement this, efficient analysis can be
performed by interpolation using the transfer function.

Enforced motion
The input motion of equivalent linear analysis is generally earthquake loading. FEA NX performs analysis that
use enforced motion. First, the equilibrium equation (5.10.17) is separated into DOFs with and without
enforced motion.

 M11 M12  u1  C11 C12  u1   K11 K12   u1   f1 

M M  u   C C  u   K K  u   f  (5.10.20)
 21 22   2   21 22   2   21 22   2   2

u1 : Displacement of unconfined DOF

u2 : Displacement of DOF confined by enforced motion
f1 : Load acting on unconfined DOF
f2 : Confining force of DOF confined by enforced motion

Separating the unconfined DOF displacement u1 into the following quasi-static displacement u1qs and
dynamic relative displacement y is as follows:

11
Lysmer, J. and Wass, G. (1972). Shear waves in plane infinite structures. Proc. ASCE, Vol. 98, EM1, pp. 85-105.

Section 10. Equivalent Linear Solution | 361

 ANALYSIS REFERENCE Chapter 5. Algorithm

u1  u1qs  y
u1  u1qs  y (5.10.21)
u1  u1qs  y

The quasi-static displacement, velocity, and acceleration can be calculated using the following equations:

u1qs  K111K12u 2
u1qs  K111K12u 2 (5.10.22)
1
u  K K12u 2
qs
1 11

Dynamic relative displacement, relative velocity, and relative acceleration are expressed as follows:

y  u1  u1qs
y  u1  u1qs (5.10.23)
y  u1  u1qs

Also, frequency domain analysis results can be found by inverting to the time domain using the FFT (fast
Fourier transform) method.

362 | Section 10. Equivalent Linear Solution

ANALYSIS REFERENCE Chapter 5. Algorithm

Section 11 Linear Buckling Analysis

Linear buckling analysis is used to determine critical load factors of a structure and the corresponding buckling
mode shapes. For linear buckling analysis, the equilibrium equation for a structure considering its geometric
stiffness due to stresses is as follows:

Ku  K Gu  p (7.11.1)

K : Elastic stiffness matrix

KG : Geometric stiffness matrix due to stresses
u : Total displacement of the structure
p : Loads acting on the structure

In linear analysis, stresses in a structure are proportional to loads, and the geometric stiffness matrix is
proportional to the stresses. Therefore, if we assume that the load, p is proportional to a reference load, p ,
the geometric stiffness, K G can be expressed as,

KG  KG
(7.11.2)
p  p

p : Reference load
KG : Geometric stiffness matrix corresponding to the reference load
 : Load factor

Substituting equation (7.11.2) for equation (7.11.1), we obtain,

Ku  KGu  p (7.11.3)

The equilibrium state as in equation (7.11.3) can be either stable or unstable depending on the magnitude of
the load factor,  . In order to assess its stability, the perturbation,  u is added to u at an equilibrium state.

Section 11. Linear Buckling Analysis | 363

 ANALYSIS REFERENCE Chapter 5. Algorithm

K(u   u)  KG (u   u)  p (7.11.4)

The non-perturbation terms in equation (7.11.4) are eliminated using the equilibrium equation (7.11.3), and
then we obtain an eigenvalue problem as follows:

(K   KG ) u  0 (7.11.5)

At this point, the stability of the equilibrium state can be assessed by the following matrix equation

K   KG  0 : Stable status
K   KG  0 : Unstable status

So the eigenvalue,  , which satisfies equation (7.11.5), can be referred to as the critical load factor at which
instability of the equilibrium begins. And the corresponding eigen mode,
 (   u ) represents the buckling

shape of the structure. The critical load causing buckling in the structure can be expressed as
p considering
the critical load factor and the reference load. Fig. 7.1.1 shows a column subjected to a compression load in its
equilibrium state and the buckling shape due to the critical load.

364 | Section 11. Linear Buckling Analysis

ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.11.1
Buckling of Column under
Compression

In linear buckling analysis, the geometric stiffness matrix, needs to be set up and the eigenvalue problem of
equation (8.6) must be solved. MIDAS calculates the geometric stiffness making use of the linear stresses or
internal forces in elements and calculates the eigenvalue problem by the Lanczos iteration method. G K
Expression (7.11.5), which is an eigenvalue problem in linear buckling analysis, can be simply expressed in the
following form.

(K  mKG )m  0 (7.11.6)

m : Critical load factor

m : Buckling mode shape

FEA NX uses the Lanczos iteration method to solve eigenvalue problems of equation (7.11.6). Equation (7.11.6)
is similar to the eigenvalue problem, but unlike mass matrix, geometric stiffness matrix K G is not positive
definite. Therefore, the linear buckling analysis replaces it with m   m / (1  m ) and uses the Lanczos
iteration method to which Shift-invert technique is applied. Calculated critical load factors are output
sequentially from low absolute values.

In the linear buckling analysis, the critical load factor  has a "+" or "-" symbol depending on the direction of
the reference load p . When finite element models and loads are complex, the symbols of the critical load

Section 11. Linear Buckling Analysis | 365

 ANALYSIS REFERENCE Chapter 5. Algorithm

factor may change with the buckling mode. Thus, FEA NX provides a function to calculate only critical loads
with "+" symbols if necessary.

Permanent loads, such as self-weight, generally remain constant unlike variable loads, but cause geometric
stiffness because they produce stresses. FEA NX provides a function that can be used as a variable load or a
constant load by selecting the type of load. Loads selected for constant loads have a constant value
independent of the load factor  , resulting in only geometric stiffness. Linear buckling analysis involving
constant loads becomes a problem of eigenvalues in the same form as equation (7.11.7).

(K  K *G  mK G )m  0 (7.11.7)

K *G : Geometric stiffness for stresses of constant loads

366 | Section 11. Linear Buckling Analysis

 ANALYSIS REFERENCE Chapter 5. Algorithm

Section 12 Fatigue Analysis

Fatigue failure is a phenomenon in which a member fails when loading smaller than the member’s yield
strength is repeatedly applied. Fatigue analysis can be performed on the basis of stress (stress-life method)
and strain (strain-life method). The stress-life method is mainly used for fatigue analysis under the condition
of a low stress level relative to the yield stress. The strain-life method on the other hand is used in fatigue
analysis under the condition of a high stress level at which plastic deformations take place. Generally, the life
obtained through the strain-life method is relatively more accurate than that obtained through the stress-life
method. Especially when the stress level is relatively high as in a stress concentration zone, fatigue analysis
using the strain-life method is appropriate. This manual summarizes the stress-life method, which uses an S-N
curve and the strain-life method, which uses an E-N curve.

The stress-life method predicts the extent of fatigue under a given loading history using the relationship
12.1 between the number of loading cycles (N) and stress amplitude (S) at the time when failure occurs from a
Stress-Life Method cyclically applied constant load. The S-N curve shows the relationship between the stress amplitude (S)
occurring from a cyclic loading of constant amplitude applied to the structure and the number of cycles to
failure (N) when the stress of the corresponding amplitude is repeated. For a fatigue analysis using the stress-
life method, a linear elastic analysis on the structure is performed first, and then equivalent stresses such as
von Mises stresses or principal stresses are obtained. These stresses are then applied to an S-N curve to
predict the number of loading cycles required to reach the fatigue failure.

Each material retains different S-N curves based on various conditions. But since fatigue tests cannot be
performed for all the predictable cases, correction factors are applied to the S-N curve, which has been
obtained from a standard uniaxial fatigue test. Also, the structures in general are subjected to cyclic loadings
of variable amplitudes. Individual stress amplitudes are extracted from the cyclic stresses of variable
amplitudes using the rainflow-counting technique and are applied to the S-N curve.

The advantages and disadvantages of the stress-life method are as follows:

►Fatigue analysis can be performed by a relatively simple algorithm.
►Calculations are simple, and the analysis time is short.
►Since elastic deformations are dealt with, it is not valid for the cases of large plastic deformations.

Section 12. Fatigue Analysis | 367

 ANALYSIS REFERENCE Chapter 5. Algorithm

Cyclic Loading
a
When a stress of constant amplitude is regularly applied as shown in Figure 5.12.1, the stress amplitude, ,
m
and the average stress, , can be calculated as follows:

 max   min
a 
2 (5.12.1)
 max   min
m 
2 (5.12.2)
Figure 5.12.1
Cyclic loading amplitude and
average stress

When a stress of constant amplitude with the average stress being zero is regularly repeated, the fatigue life is
calculated using an S-N curve in the stress-life method. An S-N curve typically retains the shape as shown in
Figure 5.12.2, connecting two points, P and Q by a straight line. Point P is obtained by the number of cycles
( N 0 ) causing fatigue failure when the stress amplitude in the magnitude of 90% of the maximum allowable
stress ( Su ) is applied. Point Q is calculated using the fatigue limit stress amplitude ( Se ), which induces no
6
fatigue failure, or the stress amplitude, which induces failure upon applying cyclic loading repeated 10 times.
Point R is located on the extended line from the line connecting P and Q, and represents failure stress ( S f )
when a load is applied once. Generally for steel, N0  1000 and Se  0.5Su can be assumed.

368 | Section 12. Fatigue Analysis

ANALYSIS REFERENCE Chapter 5. Algorithm

Figure 5.12.2
Example of typical S-N curve

By defining b for the slope of the straight line on the S-N curve shown in Figure 5.12.2, the number of cycles
to fatigue failure, N , pertaining to a cyclically applied specific stress amplitude, S , can be calculated as
Equation (5.12.4).

b
 log S  log Se 
log N e  log N
(5.12.3)
1
S  b
N  Ne  
 Se  (5.12.4)

Miner’s cumulative damage theory

When a number of stress amplitudes exist, the degree of damage of a material is calculated by accumulating
individual damages from each of the stress amplitudes. Using the number of cycles corresponding to a specific
stress amplitude, n i , and fatigue life, N i , the cumulative damage is calculated as,

ni
Damage   (5.12.5)
i Ni

The fatigue life of a structure is the reciprocal of the damage.

1
Life  (5.12.6)
Damage

Section 12. Fatigue Analysis | 369

 ANALYSIS REFERENCE Chapter 5. Algorithm

Average stress effects

If average stresses,  m , are different, the fatigue lives are also different as shown in Figure 5.12.1, even if the
stress amplitudes,  a , are identical. In order to account for the effects of average stresses, Goodman and
Gerber each proposed the following equations.

►Goodman (England, 1899)

a m
 1 (5.12.7)
Se Su

►Gerber (Germany, 1874)

2
a  
  m  1 (5.12.8)
S e  Su 

►Soderberg (USA, 1930)

a m
 1 (5.12.9)
Se Sy

►Morrow (USA, 1960s)

a m
 1 (5.12.10)
Se Sf

►SWT (Smith-Watson-Topper, 1970)

Se   a max (5.12.11)

It must be cautioned here that Se is a modified value for considering the average stress. Using the modified
Se , the S-N curve is reconstructed to reflect the effects of the average stress.

The number of cycles causing failure due to an increase in stress resulting in plastic deformations is drastically
12.2 reduced in comparison with a case in which no plastic deformations take place. In order to consider such a
Strain-Life Method phenomenon, the degree of fatigue failure is examined from the strain-life relationship, which is referred to as
the strain-life method. In order to calculate plastic strains at the state of cyclic loading, the history of strains
obtained through a nonlinear analysis can be directly used. Alternatively, the Neuber’s Law can be applied to
12
the stress results of a linear elastic analysis to obtain the plastic strains at the state of cyclic loading. Linear
elastic analysis in conjunction with the calculation of plastic strains by the Neuber’s Law was widely used in the
past due to relatively less demand for calculations. For more accurate strain calculations however, the

370 | Section 12. Fatigue Analysis

ANALYSIS REFERENCE Chapter 5. Algorithm

nonlinear history analysis method has been the preference lately. Similar to the stress-life method, the strain-
life method is also suitable for cases with plastic deformations taking place in certain stress concentration
zones rather than the cases of entire structures developing plastic strains.

Linear elastic analysis and strain calculation using the Neuber’s Law
The magnitude of a cyclic strain under the action of cyclic constant loading is obtained by summing the elastic
deformation and the plastic deformation as expressed in Equation (5.12.12). The magnitude of the elastic
deformation (  e ) is calculated from the change in stress (  ) and the modulus of elasticity ( E ). The
magnitude of the plastic deformation (  p ) is calculated using the change in stress, cyclic deformation
strength coefficient (  ' ) and the coefficient of cyclic strain hardening ( n ' ). Since the change in stress used
here is also a stress, which reflected plastic deformation as shown in Figure 5.12.3, it must be calculated by
using the Neuber’s Law expressed in Equation (5.12.13). The stress obtained from the results of linear elastic
analysis corresponds to K j S in Equation (5.12.13).
The change in strain in Equation (5.12.12) is obtained by an iterative method because the magnitude of change
in stress includes the magnitude of change in strain. The magnitude of change in stress ( S ) due to cyclic
loading of variable amplitudes is calculated using the rainflow-counting technique in the same manner as the
stress-life method.

1
  e  p     n '
    (5.12.12)
2 2 2 2 E  2 ' 
K f  S e     
1/2 1/2
(5.12.13)
Kf : Stress concentration factor

S : Magnitude of change in nominal stress (magnitude of change in stress

at a location where stress concentration effects do not exist)

Section 12. Fatigue Analysis | 371

 ANALYSIS REFERENCE Chapter 5. Algorithm

 e  p
Figure 5.12.3
Typical stress-strain curve
 B
when plastic deformation
occurs.
A

2

O 

  
C 2
2 2

The Neuber’s Law differently defines the relationship between the elastic stress and strain as shown in
Equations (5.12.14) and (5.12.15) depending on the types of equivalent stress to be applied.

►von Mises equivalent stress

 K f S   
2
S S
e   , (5.12.14)
3G  3E  3G
 
 2 1    

►Maximum shear stress or maximum principal stress

S  K f S 
2

e  ,   (5.12.15)
E E

Strain calculation using nonlinear history analysis

In case a strain history used in the strain-life method is directly calculated through nonlinear analysis, the cyclic
deformation strength coefficient (  ' ) and the coefficient of cyclic strain hardening ( n ' ) are not used. Instead
the nonlinear material properties are specified for nonlinear analysis, and the nonlinear analysis results are
used. The strain is the sum of elastic strain and plastic strain. Each strain can be calculated using the principal
 , ,
strains ( 1 2 3 ) as shown in Equations (5.12.16) – (5.12.18), depending on the type of equivalent strains.
Equivalent stress or nonlinear analysis results are used to calculate the effects of the average stress. In order to

372 | Section 12. Fatigue Analysis

ANALYSIS REFERENCE Chapter 5. Algorithm

consider the changes in stress and strain simultaneously, the multi components rainflow-counting technique,
which uses a number of components simultaneously, is applied.

►von Mises equivalent strain

1/2
 e 1   e  e 2   e  e  2   e  e  2 
  1  2    2  3    3  1  
2 2 1  v   2 2   2 2   2 2  

1/2
(5.12.16)
2  1p  2p    2p  3p    3p 1p  
2 2 2
 p
          
2 3  2 2   2 2   2 2  
 

►Maximum shear strain

 e  1e   2e  2e   3e  3e  1e 
 max  or or 
2  2 2 2 
 
(5.12.17)
 p  1   2
p p
 2   3
p p
 3  1p
p

 max  or or 
2  2 2 2 
 

►Maximum principal strain

 e 1e

2 2 (5.12.18)
 p
 p
 1
2 2

E-N curve
Figure 5.12.4 represents an E-N curve, which shows the relationship between the strain and the number of
cycles causing fatigue failure as expressed in Equation (5.12.19). By applying a number of individual strain
amplitudes to the E-N curve, which have been extracted through rainflow-counting, each number of cycles
and the corresponding individual damage level are obtained. As was done in the S-N technique, the Miner’s
cumulative damage theory is used to linearly accumulate individual damage levels and predict the final
damage and life.

  f
'

  2 N    f'  2 N 
b c
(5.12.19)
2 E

The advantages and disadvantages of the strain-life method are as follows:

►Using the results of linear elastic analysis, the method can be applied to fatigue analysis for a case in which
partial plastic deformations occur.

Section 12. Fatigue Analysis | 373

 ANALYSIS REFERENCE Chapter 5. Algorithm

►It is possible to calculate using the results of nonlinear history analysis.

►Calculations are complex and the analysis time is long.

 f
Figure 5.12.4 100
Relationship between Total strain
strain (  ) and the   f
c
 2 N    f  2 N 


1 b c

2
10
number of cycles to fatigue 2 E
failure ( N ) Strain amplitude
102
b
103
Elastic strain
 f
104 E
Plastic strain

105
100 101 102 103 104 105 106 107
Reversals to failure, 2N

Average stress effects

If the average stress,  m , is different as shown in Figure 5.12.1, the fatigue life is also different even if the
magnitude of change in strain in a structure is identical. In order to consider the effects of average stress in the
strain-life method, the following proposed equations are used:

►SWT (Smith-Watson-Topper, 1970)

  f 
' 2

 2 N    f'  f'  2 N 
b c
 max 
2b

2 E (5.12.20)

 max  m
2

►Morrow (1968)
  f   m 
'

  2 N    f'  2 N 
b c
(5.12.21)
2 E

►Manson-Halford (1981)

  f   m 
c
'
  f'   m  b
  2 N    f'   2 N 
b c
 (5.12.22)
2 E   '
f 

374 | Section 12. Fatigue Analysis