Modélisation et simulation par éléments finis : cas d'un tablier de pont.

Chapitre 3

3.3.5. Programme (matrice_K) de Calcul de la matrice de rigidité élémentaire pour un élément tétraédrique à 4 noeuds

3.3.5.1. Structure du programme

? Données du problème (Input)

- Module de Young E ;

- Coefficient de Poisson y ;

- Coordonnées des noeuds de l'élément réel (xi, yi, zi).

? Résultat (Output):

- Matrice de rigidité élémentaire Ke

3.3.5.2. Code source en FORTRAN

Program matrice_k

implicit none

integer, parameter :: n=3 ! dimension du tableau

real, dimension(1:3,1:12) ::B1,B2,B3,B4

real, dimension(1:3,1:3):: jacb_inv

integer :: i,j,lin,col,lin1,col1

real, dimension(1:12,1:12)::Ke

real, dimension(1:6,1:12)::B,res1,res2,res3,res4

real, dimension(1:12,1:6)::Bt,res5

real, dimension(1:6,1:3)::B_1,B_2,B_3,B_4

real, dimension(1:6,1:6)::H !!!!!!! matrice des constantes d'élasticté

real::x_1,y_1,z_1,x_2,y_2,z_2,x_3,y_3,z_3,x_4,y_4,z_4,vol,det_jab,nu,E

real,dimension(1:3,1:4)::c_noeud

lin=6;col=12;;lin1=6;col1=6

!!! Entrer des donnees

26 write(*,*)"Entrer le module de young E et le coefficient de POISSON v"

read(*,*,err=26)E,nu

!!!!!! Verification de la valeur de E et nu

if ((E.le.0).or.((nu.lt.0).or.(nu.gt.0.4999))) then

write(*,28)

write(*,*)" Veuillez entrer des valeurs convenables pour E et v (E>0)

(0<v<0.5) "

write(*,28)

goto 26

endif

do j=1,4

27 write(*,24)"Entrer les coordonnees du noeud ",j 24 format(x,A,2x,I1)

read(*,*,err=27)c_noeud(1,j),c_noeud(2,j),c_noeud(3,j) enddo

write(*,28)

!!! Affectation des coordonnées

x_1=c_noeud(1,1);y_1=c_noeud(2,1);z_1=c_noeud(3,1) x_2=c_noeud(1,2);y_2=c_noeud(2,2);z_2=c_noeud(3,2) x_3=c_noeud(1,3);y_3=c_noeud(2,3);z_3=c_noeud(3,3) x_4=c_noeud(1,4);y_4=c_noeud(2,4);z_4=c_noeud(3,4) !!!!!!! vérification de la singularité de J if ((det_jab(x_1,y_1,z_1,x_2,y_2,z_2,x_3,y_3,z_3,x_4,y_4,z_4)).eq.0) then write(*,*)" la matrice jacobienne de transformation est singuliere " write(*,*)" cela peut etre du a une grande distorsion de l'element " write(*,*)" Veuillez changer svp les coordonnees des noeuds"

write(*,28);28 format(/:)

86 /176

Chapitre 3

goto 27 endif

11111111111111111111111111

jacb_inv(1,1)=(y_3-y_1)*(z_4-z_1)-(z_3-z_1)*(y_4-y_1) jacb_inv(1,2)=(z_2-z_1)*(y_4-y_1)-(z_4-z_1)*(y_2-y_1) jacb_inv(1,3)=(y_2-y_1)*(z_3-z_1)-(y_3-y_1)*(z_2-z_1) jacb_inv(2,1)=(z_3-z_1)*(x_4-x_1)-(z_4-z_1)*(x_3-x_1) jacb_inv(2,2)=(x_2-x_1)*(z_4-z_1)-(x_4-x_1)*(z_2-z_1) jacb_inv(2,3)=(z_2-z_1)*(x_3-x_1)-(x_2-x_1)*(z_3-z_1) jacb_inv(3,1)=(x_3-x_1)*(y_4-y_1)-(y_3-y_1)*(x_4-x_1) jacb_inv(3,2)=(y_2-y_1)*(x_4-x_1)-(x_2-x_1)*(y_4-y_1) jacb_inv(3,3)=(x_2-x_1)*(y_3-y_1)-(y_2-y_1)*(x_3-x_1)

1111111111111111111111111111111

B_1=reshape((/-jacb_inv(1,1)-jacb_inv(1,2)-jacb_inv(1,3),0.,0.,0.,& &-jacb_inv(2,1)-jacb_inv(2,2)-jacb_inv(2,3),0.,0.,0.,& &-jacb_inv(3,1)-jacb_inv(3,2)-jacb_inv(3,3)& &,0.,-jacb_inv(3,1)-jacb_inv(3,2)-jacb_inv(3,3),-jacb_inv(2,1)-jacb_inv(2,2)& &-jacb_inv(2,3),-jacb_inv(3,1)-jacb_inv(3,2)-jacb_inv(3,3),0.,-jacb_inv(1,1)-& &jacb_inv(1,2)-jacb_inv(1,3),-jacb_inv(2,1)-jacb_inv(2,2)-jacb_inv(2,3),& &-jacb_inv(1,1)-jacb_inv(1,2)-jacb_inv(1,3),0./),(/6,3/),order=(/2,1/)) B_2=reshape((/jacb_inv(1,1),0.,0.,0.,jacb_inv(2,1),0.,0.,0.,jacb_inv(3,1),0 .,jacb_inv(3,1)& &,jacb_inv(2,1),jacb_inv(3,1),0.,jacb_inv(1,1),jacb_inv(2,1),jacb_inv(1,1), 0./),(/6,3/),order=(/2,1/)) B_3=reshape((/jacb_inv(1,2),0.,0.,0.,jacb_inv(2,2),0.,0.,0.,jacb_inv(3,2),0 .,jacb_inv(3,2)& &,jacb_inv(2,2),jacb_inv(3,2),0.,jacb_inv(1,2),jacb_inv(2,2),jacb_inv(1,2), 0./),(/6,3/),order=(/2,1/)) B_4=reshape((/jacb_inv(1,3),0.,0.,0.,jacb_inv(2,3),0.,0.,0.,jacb_inv(3,3),0 .,jacb_inv(3,3)& &,jacb_inv(2,3),jacb_inv(3,3),0.,jacb_inv(1,3),jacb_inv(2,3),jacb_inv(1,3), 0./),(/6,3/),order=(/2,1/)) H=reshape((/1-nu,nu,nu,0.,0.,0.,nu,1-nu,nu,0.,0.,0.,nu,nu,1-nu,0.,0.,0.,0.,0.,0.,(1-2*nu)/2,0.,0.,0.,0.,& &0.,0.,(1-2*nu)/2,0.,0.,0.,0.,0.,0.,(1-2*nu)/2/),(/6,6/),order=(/2,1/)) 111 function pour calculer le coefficient multiplicateur de Ke vol=36*((1/6.)*det_jab(x_1,y_1,z_1,x_2,y_2,z_2,x_3,y_3,z_3,x_4,y_4,z_4))*(1 -2*nu)*(1+nu)

11111111111111 appel des procedures

call table_Bi(1.,B1);call table_Bi(2.,B2);call table_Bi(3.,B3);call
table_Bi(4.,B4)

call produit(B_1,B1,lin,3,3,col,res1);call produit(B_2,B2,lin,3,3,col,res2) call produit(B_3,B3,lin,3,3,col,res3);call produit(B_4,B4,lin,3,3,col,res4) B=res1+res2+res3+res4;Bt=transpose(B)

call produit(Bt,H,12,6,6,6,res5);call produit(res5,B,col,lin1,lin,col,ke) call system('mkdir c:\matrice_K_B')

open(unit=20,file=' c:\matrice_K_B\metrice_K.txt',status='unknown')

write(20,15)"+

+"

write(20,15)"| Matrice de rigidité élémentaire d'un élément tétraédrique à 4 noeuds|"

write(20,15)"+

+"

15 format(A)

write(20,16);16 format(/:)

11111111111 écriture des noeuds

write(20,18)"+ +"
write(20,18)"| Coordonnées des noeuds (x, y, z) |"

write(20,25)"+ +"
18 format(x,A);25 format(x,A,/:)

do i=1,4

write(20,19)"Noeud ",i,"|"

87 /176