# Numerical integration in time for finite elements

I am trying to solve $$M\ddot{u}=-Ku+F_\text{ext}$$ for a 2D linear elastic model with $$M$$ be the mass matrix,$$K$$ the stiffness matrix and $$F_\text{ext}$$ the external load vector coming from a uniformly distributed load acting on one edge of the model.(Note: $$F_\text{ext}$$ is not time-dependent). An explicit time-scheme is used and more specific Forward-Euler scheme. The steps that I follow are:

1. Initial conditions $$\dot{u}_0=0$$ $$u_0=0$$
2. Solve $$M\ddot{u}_n=-Ku_n+F_\text{ext}$$ using an iterative solver
3. Update $$u_\text{n+1}=u_n+dt\dot{u}_n$$
4. Update $$\dot{u}_\text{n+1}=\dot{u}_n+dt\ddot{u}_n$$
5. Go back to 2 for next time step

Based on this implementation I noticed that the output values (valocity,displacment,acceleration) go to infinity.What is the main issue that can cause this problematic behaviour?I want to note that the used time-step is small $$10^{-6}$$ so I don't think is a stability issue. Here is the main routine:

for(int i=0;i<2*NN;i++){
RHS[i]=0;;
}

for(int i=0;i<2*NN;i++){
double sum=0;
for(int j=0;j<2*NN;j++){
sum+=K_global[i][j]*displ[j];
}
RHS[i]=Fext[i]-sum;

}

BoundaryCondForRHS(NN,NEy,dbc,RHS);//rows connected with BC are set to zero

/*update*/
for(int i=0;i<2*NN;i++){
displ[i]=dt*veloc[i]+displ[i]; //displ at t->n+1
veloc[i]=dt*accel[i]+veloc[i]; //veloc at t->n+1
}

• Do you mean Backward Euler? Forward Euler is not implicit. Oct 12 '20 at 22:02
• How are you coupling $\ddot u$ and $u$? It seems like you are actually solving $M \ddot u_{n+1}=K u_n+ F_{ext}$ Oct 12 '20 at 22:11
• @StevenRoberts I meant explicit toy are right Oct 12 '20 at 22:12
• @CharlieS I am using $\ddot{u}_n$ to update velocities.Is my implementation wrong? Oct 12 '20 at 22:14
• Why not use backward Euler? You'll have to do a linear solve either way. Might as well get better stability. Oct 12 '20 at 22:17

Your step two is to solve the original ODE, which doesn't make sense. I'll write out the steps for applying Forward Euler to your second order ODE. Forward Euler solves the first order ODE $$M \dot{y} = f(y)$$ with the steps \begin{align} M k_1 &= f(y_n) \\ y_{n+1} &= y_n + dt \, k_1 \end{align}

Let $$v = \dot{u}$$. You finite element problem can be rewritten as $$\begin{bmatrix} I & 0 \\ 0 & M \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}' = \begin{bmatrix} v \\ K u + F_{ext} \end{bmatrix}.$$

When we apply backward Euler to this we get \begin{align} \begin{bmatrix} k_1 \\ M \ell_1 \end{bmatrix} &= \begin{bmatrix} v_n \\ K u_n + F_{ext} \end{bmatrix}, \\ \begin{bmatrix} u_{n+1} \\ v_{n+1} \end{bmatrix} &= \begin{bmatrix} u_n + dt \, k_1 \\ v_n + dt \, \ell_1 \end{bmatrix} \end{align}

So the steps are

1. Initial conditions $$u_0$$ and $$\dot{u}_0$$
2. Solve the linear system $$M \ell_1 = K u_n + F_{ext}$$ for $$\ell_1$$
3. $$u_{n+1} = u_n + dt \, \dot{u}_n$$
4. $$\dot{u}_{n+1} = \dot{u}_n + dt \, \ell_1$$
5. Go back to two

Edit: Now that the code is posted, it looks like you are following these steps. I suspect if you look at the eigenvalues of the Jacobain $$\begin{bmatrix} 0 & I \\ M^{-1} K & 0 \end{bmatrix},$$ there will be some (if not all) along the imaginary axis, which is outside the stability region. I would recommend trying Backward Euler or the Newmark method.

• I think you wrote exactly what I did ,with $l1 = \ddot{u}_n$ Oct 12 '20 at 22:34
• So to be crystal clear ,the problem lies to the time-integration method and not to implementation? Oct 12 '20 at 22:49
• Yes, it looks like the steps are the same. I got a bit confused with the $\ddot{u}_n$ notation. I don't see any issues with the code. For a short timespan, I would expect forward Euler to converge, even with eigenvalues outside the stability region. If possible, I would try to construct the matrix above and check the eigenvalues. Oct 12 '20 at 23:08
• Is it possible, that the time step is not small enough?For example may a time step of $10^{-9}$ or something like that is needed Oct 12 '20 at 23:33
• Your linear elastic problem will oscillate and thus have imaginary eigenvalues. As stated above, for purely imaginary eigenvalues, explicit Euler is unconditionally unstable: that is, no matter how small your time step, your solution will diverge. If you make the time step small, it will grow slower but it will still grow. You really need to use a different integration method like Newmark, trapezoidal rule or something which is stable for imaginary eigenvalues. Oct 13 '20 at 7:05