Intégration numérique dans le temps pour les éléments finis

Question

J'essaye de résoudre $M\ddot{u}=-Ku+F_\text{ext}$ pour un modèle élastique linéaire 2D avec $M$ être la matrice de masse,$K$ la matrice de rigidité et $F_\text{ext}$ le vecteur de charge externe provenant d'une charge uniformément répartie agissant sur un bord du modèle (Remarque: $F_\text{ext}$ne dépend pas du temps). Un schéma temporel explicite est utilisé et un schéma Forward-Euler plus spécifique. Les étapes que je suis sont:

1. Conditions initiales $\dot{u}_0=0$ $u_0=0$
2. Résoudre $M\ddot{u}_n=-Ku_n+F_\text{ext}$ en utilisant un solveur itératif
3. Mettre à jour $u_\text{n+1}=u_n+dt\dot{u}_n$
4. Mettre à jour $\dot{u}_\text{n+1}=\dot{u}_n+dt\ddot{u}_n$
5. Revenir à 2 pour la prochaine étape de temps

Sur la base de cette implémentation, j'ai remarqué que les valeurs de sortie (valeur, déplacement, accélération) vont à l'infini Quel est le principal problème qui peut causer ce comportement problématique? Je veux noter que le pas de temps utilisé est petit $10^{-6}$donc je ne pense pas que ce soit un problème de stabilité. Voici la routine principale:

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

Votre deuxième étape consiste à résoudre l'ODE d'origine, ce qui n'a aucun sens. Je vais écrire les étapes pour appliquer Forward Euler à votre deuxième ordre ODE. Forward Euler résout le premier ordre ODE$$ M \dot{y} = f(y) $$ avec les étapes $$ \begin{align} M k_1 &= f(y_n) \\ y_{n+1} &= y_n + dt \, k_1 \end{align} $$

Laisser $v = \dot{u}$. Votre problème d'éléments finis peut être réécrit comme$$ \begin{bmatrix} I & 0 \\ 0 & M \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}' = \begin{bmatrix} v \\ K u + F_{ext} \end{bmatrix}. $$

Lorsque nous appliquons Euler en arrière à cela, nous $$ \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} $$

Donc les étapes sont

1. Conditions initiales $u_0$ et $\dot{u}_0$
2. Résolvez le système linéaire $M \ell_1 = K u_n + F_{ext}$ pour $\ell_1$
3. $u_{n+1} = u_n + dt \, \dot{u}_n$
4. $\dot{u}_{n+1} = \dot{u}_n + dt \, \ell_1$
5. Revenir à deux

Modifier: maintenant que le code est publié, il semble que vous suivez ces étapes. Je soupçonne si vous regardez les valeurs propres du Jacobain$$ \begin{bmatrix} 0 & I \\ M^{-1} K & 0 \end{bmatrix}, $$il y en aura quelques-uns (sinon tous) le long de l'axe imaginaire, qui est en dehors de la région de stabilité. Je recommanderais d'essayer Backward Euler ou la méthode Newmark.