Equation Solving Methods for Nonlinear Calculations

Introduction

The finite element method is used whenever mechanical problems cannot be solved analytically. Often, nonlinear effects such as failure under pressure (geometric nonlinearity), plastic yielding (material nonlinearity), and contact or kinematic degrees of freedom are also taken into account. These effects, especially for nonlinear material models, can be considered by using an iterative calculation method.

FEM Setting

Basic steps for FEM setting (further information can be found in [1]):

1. Weak form of equilibrium
   $${\int }_B\underset{\mathrm{nonlinear} \mathrm{component}}{\underbrace{grad{\delta }_u\cdots \sigma dV}}={\int }_B{\delta }_u·\rho bdV+{\int }_{\partial B}{\delta }_u·t_{o}dA$$

   δu	Virtual (test) displacement
   t0	Initial load factor
   σdV	Internal forces
   ρbdV	Solid forces
   B	Integrated area

   Weak Form
2. Conversion to Voigt notation with Tensor 4 Step
   $${\int }_{\mathrm{B}}\mathrm{\delta \epsilon }\cdots \mathrm{ℂ}\cdots \mathrm{\epsilon dV} = {\int }_{\mathrm{B}}\underset{\mathrm{solid} \mathrm{forces}}{\underbrace{{\mathrm{\delta u}}^{\mathrm{T}}\mathrm{bdV}}} + {\int }_{\partial \mathrm{B}}\underset{\mathrm{surface} \mathrm{area} \mathrm{forces}}{\underbrace{{\mathrm{\delta u}}^{\mathrm{T}} · {\mathrm{t}}_0\mathrm{dA}}}$$

   C	Stiffness matrix
   δε	Variation of the strain state
   B	Integrated area
   ε	Strain
   U	Deformation

   Weak Form in Voigt Notation
3. This notation is used in the following text to solve the approximate solution of the FE approach for nonlinear materials.
4. For this, the displacement field is multiplied by the approach functions.
   $$\mathrm{u}(\mathrm{x},\mathrm{t}) = \mathrm{H}\left(\mathrm{x}\right)\hat{\mathrm{u}}\left(\mathrm{t}\right)$$

   u(x,t)	Displacement over time (load increment)
   H	Form function
   û	Nodal displacement

   Displacement Field via Trial Function

Inserting the derivative of displacement into the weak form. Numerical integration is used to calculate the nodal displacement, stresses, and strains in post-processing by means of the material rule.

Newton-Raphson Iteration Sequence

Due to the nonlinear material behavior, the material matrix C in Equation 2 above changes with each expansion step. The standard calculation method for solving this problem is the Newton-Raphson iteration. It is used to linearize the function in a starting point. The stiffness matrix C of the preliminary step is always used in the iteration. In the linearized iteration step, a tangent is placed at the zero of the function.

The equations belonging to the flowchart in the image above are the following:

1. Dividing the load into load steps.
   $${}^{\mathrm{t} + ∆\mathrm{t}}\mathrm{f}_{\mathrm{ext}} = {}^{\mathrm{t}}\mathrm{f}_{\mathrm{ext}} + ∆\mathrm{f}$$

   fext	External load
   t	Time/load step
   Δf	New load increment

   Load in Load Steps
2. Predictor Step
   $${}_0{}^{\mathrm{t}}\mathrm{K}{∆}_0\mathrm{\phi } = {}^{\mathrm{t}+∆\mathrm{t}}\mathrm{f}_{\mathrm{ext}} - {}_0{}^{\mathrm{t}}\mathrm{f}_{\mathrm{int}}$$

   K	Stiffness matrix of the previous time step
   t+Δtfext	External forces increased by a further load step
   0tfint	Internal forces of the previous time step
   ϕ	Strain

   Predictor Step
3. In Point 3, iteration of the flowchart, the total distortion reduced by the plastic distortion is calculated (corrector step). The objective of the iterative calculation is always that the sum of the loads is zero. However, this is impossible numerically. Therefore, a break-off limit ε is defined at which the calculation is interrupted as sufficiently accurate.
   $$\mathrm{R} = \left|{}^{\mathrm{t} + ∆\mathrm{t}}\mathrm{f}_{\mathrm{ext}} - {}_{\mathrm{n}}{}^{\mathrm{t} + ∆\mathrm{t}}\mathrm{f}_{\mathrm{int}}\right| < \mathrm{\epsilon }$$

   R	Break-off limit
   fext	External forces
   fint	Internal forces
   ε	Epsilon break-off limit
   t	Time step

   Break-Off Limit

In the program, the break-off limit can be set among the calculation parameters.

The following image shows the flow of a Newton-Raphson iteration. In the first iteration

the break-off R or ε is not reached. The tolerance limit is not reached in the second iteration (red) either. Only in the third iteration is the distance of the tangent stiffness so small that convergence is achieved.

As already mentioned, the deformation is continuously summed up during the iteration.

Conclusion

The Newton-Raphson iteration has the consistency or convergence order 2. The number of "correct" locations in the iteration doubles with each step. Thus, a Newton-Raphson iteration converges quadratically, and the accuracy increases with each iteration when the method converges. However, if the method does not converge, the error goes to infinity and the calculation is stopped.

Error causes are, for example, a slope of the load-deformation curve that is too steep, and a curve gradient in the plastic area that is too flat. If the load-deformation curve in the image above showed a break in the second iteration step that was too strong, the material tangent and thus the stiffness matrix would not correctly display the slope of the elastic area. In this case, the slope of the root would be wrong for the plastic area. This is one of the reasons that an increase in the load steps is accompanied by an improved convergence.