6.4 Assemblage et calcul des modes propres

Le calcul des matrices de masse et de raideur utilise la technique d'assemblage décrite dans les chapitres précédents.

On calcule tout d'abord une matrice de masse et de raideur élémentaire sur l'élément de référence, ce qui permet d'obtenir les matrices élémentaires de masse $\mathbf{M}^{k}$ et de raideur $\mathbf{K}^{k}$ sur un élément $e_{k}$ de longueur $h$:

$$\mathbf{K}^{k}=\frac{c_{0}^{2}}{h}\left[\begin{array}{cc} 1 ... ...3} & \frac{1}{6}\\ \frac{1}{6} & \frac{1}{3}\end{array}\right]$$

On effectue ensuite un assemblage élément par élément pour calculer les matrices globales de raideur $\mathbf{K}$ et de masse $\mathbf{M}$, puis on applique les conditions aux limites.

Le programme Matlab correspondant à cet algorithme d'assemblage est décrit ci-dessous (programme 7.4):

programme Matlab 7.4: Assemblage de la matrice de masse et de rigidité

function [K,M]=assemblage(Ne,c0,L)
% assemblage de la matrice de masse et de rigidite
% pour une corde (C) M. BUFFAT
% maillage
h=1/Ne; N=Ne+1;
% matrice elementaire
Kh=[1 -1; -1 1];
Mh=[1/3 1/6; 1/6 1/3];
% assemblage
M=zeros(N,N);K=zeros(N,N);
for k=1:Ne
    M(k:k+1,k:k+1)=M(k:k+1,k:k+1)+h*Mh;
    K(k:k+1,k:k+1)=K(k:k+1,k:k+1)+c0^2/h*Kh;
end;
% C.L. on elimine les nds 1 et NN
M=M(2:N-1,2:N-1);
K=K(2:N-1,2:N-1);
% fin
return

On note que la prise en compte des conditions aux limites a été effectue en éliminant les lignes et les colonnes des 2 noeuds frontières (noeud 1 et noeud N).

Pour calculer les valeurs propres et vecteurs propres de la matrice $\mathbf{A}=\mathbf{M}^{-1}\mathbf{K}$, on résout le problème équivalent de valeurs propres généralisés:

$$\mathbf{M}x=\lambda\mathbf{K}x$$

en utilisant la fonction de calcul de valeurs propres eig de Matlab en spécifiant que les matrices sont symétriques définie positives.

Le programme Matlab correspondant est décrit ci-dessous (programme 7.4):

programme Matlab 7.4: Modes propres de vibration d une corde

% mode de vibration d'une corde (C) M. BUFFAT
L=1; c0=1;
% maillage
Ne=10; h=1/Ne; X=[0:h:L];
N=Ne+1;
% assemblage
[K,M]=assemblage(Ne,c0,L);
% calcul des valeurs propres D et vecteurs propres V 
% du probleme  M.X=lambda*K.X
[V,D]=eig(K,M,'chol');
% trace des premiers modes propres
X=[0:h:L]'; 
Xe=[0:h/10:L]';
figure(1);k=1; Vex=sin(k*pi*X);
plot(X,[0;V(:,k);0],'r-*',Xe,sin(k*pi*Xe)*max(V(:,k))/max(Vex));
title('comparaison calcul (Ne=10) theorie: mode k=1');
xlabel('x'); ylabel('Amplitude');
figure(2);k=6; Vex=sin(k*pi*X);
plot(X,[0;V(:,k);0],'r-*',Xe,sin(k*pi*Xe)*max(V(:,k))/max(Vex));
title('comparaison calcul (Ne=10) theorie: mode k=6');
figure(3);k=8; Vex=sin(k*pi*X);
plot(X,[0;V(:,k);0],'r-*',Xe,sin(k*pi*Xe)*max(V(:,k))/max(Vex));
title('comparaison calcul (Ne=10) theorie: mode k=8');
figure(4); I=[1:N-2];
plot(I,sqrt(diag(D))/2/pi,'r-*',I,I*pi/2/pi);
title('comparaison calcul (Ne=10) theorie: frequence propre');
xlabel('mode'); ylabel('frequence');