Nonlinear finite elements/Solution procedure

Solution procedure

The finite element system of equations is

$𝐌\ddot{𝐮}={𝐟}_{\text{ext}}-{𝐟}_{\text{int}}=𝐟.$

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 $\mathrm{\Delta }t$. Let us also assume that the time step is constant. Then, at time step $n$, the time is $t_n=n\mathrm{\Delta }t$.

Let

${𝐮}_n=𝐮(t_n).$

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

${𝐯}^{n+1/2}=\frac{{𝐮}^{n+1}-{𝐮}^n}{\mathrm{\Delta }t}.$

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

${𝐚}^n={\ddot{𝐮}}^n=\frac{{𝐯}^{n+1/2}-{𝐯}^{n-1/2}}{\mathrm{\Delta }t}.$

Now,

${\ddot{𝐮}}^n={𝐌}^{-1}{𝐟}^n.$

Therefore,

${𝐌}^{-1}{𝐟}^n=\frac{{𝐯}^{n+1/2}-{𝐯}^{n-1/2}}{\mathrm{\Delta }t}⟹{𝐯}^{n+1/2}={𝐯}^{n-1/2}+\mathrm{\Delta }t{𝐌}^{-1}{𝐟}^n.$

Algorithm

1. Initialize:
   1. Set $t=0$ and $n=0$.
   2. Set the initial velocities (${𝐯}^0$) at the nodes.
   3. Set the initial accelerations (${𝐚}^0$) at the nodes.
   4. Set the initial stresses (${\sigma }^0$) at the element Gauss points.
   5. Compute the lumped mass $𝐌$ at the nodes.
2. Set the displacements and velocities at the first time step:
   1. Displacement: ${𝐮}^0=0.$
   2. Velocity: ${𝐯}^{1/2}={𝐯}^0+\frac{\mathrm{\Delta }t}{2}{𝐚}^0.$
3. Apply essential BCs: $[{𝐯}^{1/2}{]}_{{\mathrm{\Gamma }}_v}=\overline{𝐯}(𝐱,\mathrm{\Delta }t/2).$
4. Update the nodal displacements: ${𝐮}^1={𝐮}^0+\mathrm{\Delta }t{𝐯}^{1/2}.$
5. Set $n=1$ and $t=\mathrm{\Delta }t$.
6. Loop through the following steps until $t=t_{\text{max}}$.
7. For each element (at the Gauss points):
   1. Compute the strain measure: ${\epsilon }^n=u_{,X}\text{or}{D}^n=v_{,x}.$
   2. Compute the stress: $P=E^{PF}{\epsilon }^n\text{or}\dot{\sigma }=E^{\sigma D}{D}^n.$
8. For each node:
   1. Compute the internal force $(f_i^n{)}_{\text{int}}$.
   2. Compute the external force $(f_i^n{)}_{\text{ext}}$.
   3. Compute the total force $f_i^n$.
   4. Compute the acceleration $a_i^n$: $a_i^n=M_{ii}^{-1}{f}_i^n.$
   5. Update the velocity $v_i^{n+1/2}$. $v_i^{n+1/2}=v_i^{n-1/2}+\mathrm{\Delta }t{a}_i^n.$
   6. Apply the essential boundary conditions.
   7. Update the displacement $u_i^{n+1}$. $u_i^{n+1}=u_i^n+\mathrm{\Delta }t{v}_i^{n+1/2}.$
9. Update the counters: $n=n+1$, $t=t+\mathrm{\Delta }t$.

Stability of the explicit algorithm

If the time step $\mathrm{\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 $\mathrm{\Delta }t$. This condition (in 1-D) states that

$\mathrm{\Delta }{t}_{\text{crit}}=\frac{l_0}{c_0}$

where $l_0$ is the initial length of the element, and $c_0$ is speed of sound in the material (wavespeed) given by

$c_0=\frac{E^{PF}}{{\rho }_0}.$