Numerical integration in time for finite elements

Question

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

ConjugateGradient(2*NN,M_global,RHS,accel);//find acceleration at t->n


/*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
}

Answer 1

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.