Nonlinear finite elements/Solution procedure

The finite element system of equations is

${\displaystyle \mathbf {M} ~{\ddot {\mathbf {u} }}=\mathbf {f} _{\text{ext}}-\mathbf {f} _{\text{int}}=\mathbf {f} ~.}$

Let's assume that the mass matrix is diagonal. Let us also assume that we would like to solve this problem using an explicit method that uses central differencing.

Let the time step size be ${\displaystyle \Delta t}$. Let us also assume that the time step is constant. Then, at time step ${\displaystyle n}$, the time is ${\displaystyle t_{n}=n\Delta t}$.

Let

${\displaystyle \mathbf {u} _{n}=\mathbf {u} (t_{n})~.}$

Let us take a half step to compute the velocity at the middle of the timestep. Then,

${\displaystyle \mathbf {v} ^{n+1/2}={\cfrac {\mathbf {u} ^{n+1}-\mathbf {u} ^{n}}{\Delta t}}~.}$

We will use this half-step velocity to compute the acceleration

${\displaystyle \mathbf {a} ^{n}={\ddot {\mathbf {u} }}^{n}={\cfrac {\mathbf {v} ^{n+1/2}-\mathbf {v} ^{n-1/2}}{\Delta t}}~.}$

Now,

${\displaystyle {\ddot {\mathbf {u} }}^{n}=\mathbf {M} ^{-1}~\mathbf {f} ^{n}~.}$

Therefore,

${\displaystyle \mathbf {M} ^{-1}~\mathbf {f} ^{n}={\cfrac {\mathbf {v} ^{n+1/2}-\mathbf {v} ^{n-1/2}}{\Delta t}}\qquad \implies \qquad \mathbf {v} ^{n+1/2}=\mathbf {v} ^{n-1/2}+\Delta t~\mathbf {M} ^{-1}~\mathbf {f} ^{n}~.}$

1. Initialize:
   1. Set ${\displaystyle t=0}$ and ${\displaystyle n=0}$.
   2. Set the initial velocities (${\displaystyle \mathbf {v} ^{0}}$) at the nodes.
   3. Set the initial accelerations (${\displaystyle \mathbf {a} ^{0}}$) at the nodes.
   4. Set the initial stresses (${\displaystyle \sigma ^{0}}$) at the element Gauss points.
   5. Compute the lumped mass ${\displaystyle \mathbf {M} }$ at the nodes.
2. Set the displacements and velocities at the first time step:
   1. Displacement: ${\displaystyle \mathbf {u} ^{0}=0~.}$
   2. Velocity: ${\displaystyle \mathbf {v} ^{1/2}=\mathbf {v} ^{0}+{\cfrac {\Delta t}{2}}~\mathbf {a} ^{0}~.}$
3. Apply essential BCs: ${\displaystyle [\mathbf {v} ^{1/2}]_{\Gamma _{v}}={\bar {\mathbf {v} }}(\mathbf {x} ,\Delta t/2)~.}$
4. Update the nodal displacements: ${\displaystyle \mathbf {u} ^{1}=\mathbf {u} ^{0}+\Delta t~\mathbf {v} ^{1/2}~.}$
5. Set ${\displaystyle n=1}$ and ${\displaystyle t=\Delta t}$.
6. Loop through the following steps until ${\displaystyle t=t_{\text{max}}}$.
7. For each element (at the Gauss points):
   1. Compute the strain measure: ${\displaystyle \varepsilon ^{n}=u_{,X}~~{\text{or}}~~D^{n}=v_{,x}~.}$
   2. Compute the stress: ${\displaystyle P=E^{PF}~\varepsilon ^{n}~~{\text{or}}~~{\dot {\sigma }}=E^{\sigma D}~D^{n}~.}$
8. For each node:
   1. Compute the internal force ${\displaystyle (f_{i}^{n})_{\text{int}}}$.
   2. Compute the external force ${\displaystyle (f_{i}^{n})_{\text{ext}}}$.
   3. Compute the total force ${\displaystyle f_{i}^{n}}$.
   4. Compute the acceleration ${\displaystyle a_{i}^{n}}$: ${\displaystyle a_{i}^{n}=M_{ii}^{-1}~f_{i}^{n}~.}$
   5. Update the velocity ${\displaystyle v_{i}^{n+1/2}}$. ${\displaystyle v_{i}^{n+1/2}=v_{i}^{n-1/2}+\Delta t~a_{i}^{n}~.}$
   6. Apply the essential boundary conditions.
   7. Update the displacement ${\displaystyle u_{i}^{n+1}}$. ${\displaystyle u_{i}^{n+1}=u_{i}^{n}+\Delta t~v_{i}^{n+1/2}~.}$
9. Update the counters: ${\displaystyle n=n+1}$, ${\displaystyle t=t+\Delta t}$.

If the time step ${\displaystyle \Delta t}$ is too large, the algorithm may not be stable and may give unreasonable results. To provide a check on the time step size, we use the CFL (Courant-Friedrichs-Lewy) condition to determine ${\displaystyle \Delta t}$. This condition (in 1-D) states that

${\displaystyle {\Delta t_{\text{crit}}={\cfrac {l_{0}}{c_{0}}}}}$

where ${\displaystyle l_{0}}$ is the initial length of the element, and ${\displaystyle c_{0}}$ is speed of sound in the material (wavespeed) given by

${\displaystyle {c_{0}={\cfrac {E^{PF}}{\rho _{0}}}~.}}$