Problème de barre 1D en MEF (méthode par éléments finis) - méthode explicite

[invite67c38467]

Bonjour,

Avec la méthode des éléments finis, je dois écrire un programme Scilab permettant de résoudre par méthode explicite un problème de barre à 1 dimension. J'ai déjà écris quelque chose mais suis surpris pas les résultats que j'obtiens. Je poste donc ici pour vous demander votre avis, et éventuellement corriger ce qui ne va pas. Je joins ici également le sujet : Task_2.pdf (pour l'instant je cherche juste à écrire le programme, pas répondre aux questions ^^ ).

Explication du programme

Je considère que chaque élément de ma barre est un ressort de raideur où est le module de Young, la section de la barre, le nombre de nœuds de ma barre (et donc le nombre d'éléments de ma barre), et la longueur totale de ma barre. La matrice de raideur est donc :

Je pose l'hypothèse que ma barre est fixe à une extrémité, ma matrice devient donc :

Et donc pour trouver le déplacement de chaque nœud à l'instant je cherche à résoudre l'équation suivante :

où est la matrice des masse (matrice diagonale, ou chaque terme vaut la masse divisé par le nombre de nœuds), la matrices des force extérieures, pour mon problème, pour l'instant, il s'agit juste d'une force constante à l'extrémité.

Pour résoudre cette équation j'utilise la méthode explicite, comme vous pouvez le voir ci-dessous dans l'extrait de code Scilab :

D = zeros(n, nt); //d : déplacement
Dp = zeros(n, nt); // d point : vitesse
Dpp = zeros(n, nt); // d point point : accélération
//valeurs initiales du déplacement, vitesse, accé
Dinit = zeros(n,1);
Dpinit = zeros(n,1);
Dppinit = zeros(n,1);

Rint= zeros(n, nt);
Rint(:,1) = K * Dinit;
D(:, 1) = Dinit + deltat * Dpinit  + 1/2*deltat^2*inv(M)*(P(:,1)-Rint(:,1));
Dp(:,1) = Dpinit + 1/2*deltat*Dppinit;

for i=1:(nt-1)
    
    Dpp(:,i) = inv(M)*(P(:,i)-Rint(:,i));
    D(:,i+1) = D(:,i) + deltat * Dp(:,i) + deltat^2 * inv(M)*(P(:,i)-Rint(:,i));
    Rint(:,i+1) = K * D(:, i+1);

    Dp(:,i+1) = Dp(:,i) + deltat*Dpp(:,i);
end

Ci-dessous le code complet commenté :

 Cliquez pour afficher

clear all
clc

l = 1000; // bar length (mm)
e = 70000; // aluminium Young's modulus (MPa)
density = 2.7; // aluminium
a = 10*15; // area (mm^2)
n = 10;
n=n-1
tmax = 5; // (s)
deltat = 0.01; // (s)

X=0:(l/(n)):(l-l/n); // vecteur de coordonnées des nœuds
X=X';
T=0:deltat:(tmax-deltat); // vecteur des time steps
T=T';

// Fabriquer la matrice K :
k = a*e*(n-1)/l; // ATTENTION
Diag1 = ones(n,1);
Diag2 = ones(n-1,1);
K = 2*diag(Diag1) - diag(Diag2, 1) - diag(Diag2, -1);
K(1,1)=1;
//K(n,n)=1;
K = k*K;

M = density*10^(-3)*a*l/n * diag(ones(n,1)); // Fabriquer la matrice es masses M

nt = tmax / deltat; // nombre de times steps

// Création d'une force constante
P = zeros(n,nt);//nt
P(1, :)= 1000;//sin(T'*2*%pi()/(l/3*sqrt(density/e))); // (N)

D = zeros(n, nt); //d : déplacement
Dp = zeros(n, nt); // d point : vitesse
Dpp = zeros(n, nt); // d point point : accélération
//valeurs initiales du déplacement, vitesse, accé
Dinit = zeros(n,1);
Dpinit = zeros(n,1);
Dppinit = zeros(n,1);

Rint= zeros(n, nt);
Rint(:,1) = K * Dinit;
D(:, 1) = Dinit + deltat * Dpinit  + 1/2*deltat^2*inv(M)*(P(:,1)-Rint(:,1));
Dp(:,1) = Dpinit + 1/2*deltat*Dppinit;

for i=1:(nt-1)
    
    Dpp(:,i) = inv(M)*(P(:,i)-Rint(:,i));
    D(:,i+1) = D(:,i) + deltat * Dp(:,i) + deltat^2 * inv(M)*(P(:,i)-Rint(:,i));
    Rint(:,i+1) = K * D(:, i+1);

    Dp(:,i+1) = Dp(:,i) + deltat*Dpp(:,i);
end


D= [D; zeros(1,nt)];
X = [X;l]
plot3d(X,1000*T,1000*D);

Résultat obtenu
J'affiche donc le résultat sous forme 3D : axe X ce sont les coordonnées de la barre, axe Y le temps, et axe Z le déplacement des nœuds :
deplcament.png
Comme vous pouvez le voir, j'obtiens un déplacement des nœuds périodique alors que la force que j'applique est constante. Je trouve cela très bizarre étant donné que si j'applique une force constante c'est "comme si" je faisais un essai de traction... Or lors d'un essai de traction la barre est constamment en extension... Ensuite je me dis que j'ai modélisé tout ça par des ressorts, est-ce une erreur ? Je ne vois pas comment faire autrement et je ne sais pas si je dois croire ce le résultat que j'obtiens... Quel est donc votre avis ? Merci.
-----

[invite178a8847]

Bonjour,

C'est un peu tard pour te répondre, mais ça pourra servir aux autres. L'équation que tu résous par la méthode des éléments finis correspond à un système masse plus ressort sans frottements, il n'est donc pas étonnant que tu observes une oscillation. Il faut que tu ajoutes un terme dissipatif si tu veux converger vers l'équilibre statique.