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

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

Chapitre 3

19 format(A,x,I1,2x,A,x,$)

do j=1,3

write(20,21)c_noeud(j,i),"|"

21 format(F8.4,x,A,$)

enddo

write(20,22);22 format(2/:)

enddo

!!!!!!!!!!!! écriture la matrice Ke

write(20,23)"Ke =",E/vol,"*"

23 format(A,x,F18.4,x,A)

do i=1,12

write(20,11)"|"

11 format(x,A,x,$)

do j=1,12

write(20,12)ke(i,j)

12 format(F8.4,x,"|",$)

enddo

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

enddo

close(20)

call system('start c:\matrice_K_B\metrice_K.txt')

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end program

!

!procédure pour la table de localisation pour les matrices Bi

!cette sous routine permet d'insérer la matrice identité à la

!place de Bi

!

subroutine table_Bi(col,tab)

integer i,j

real::col

real, dimension(1:3,1:12)::tab

tab=0;j=3*col-2

do i=1,3

tab(i,j)=1;j=j+1

enddo

return

end

!procédure pour calculer les termes ji la matrice jacobienne inverse

!

! procédure pour calculer le déterminant de la matrice

! de la matrice jacobienne de transformation

!warning: cette procédure est pour ce cas spécificique

!il ne peut être utilisé pour calculer le déterminant d'1e autre matrice

function 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)

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,det_jab

det_jab=((x_2-x_1)*(y_3-y_1)*(z_4-z_1)+(y_2-y_1)*(z_3-z_1)*(x_4-x_1)+&

&(z_2-z_1)*(x_3-x_1)*(y_4-y_1)-(x_2-x_1)*(z_3-z_1)*(y_4-y_1)-(y_2-&

&y_1)*(x_3-x_1)*(z_4-z_1)-(z_2-z_1)*(y_3-y_1)*(x_4-x_1))

return

end function det_jab

!

subroutine produit(mat1,mat2,n,c,n1,c1,res)

integer::i,j,k,n,c,n1,c1

real, dimension(1:n,1:c1)::res

real, dimension(1:n,1:c)::mat1

real, dimension(1:n1,1:c1)::mat2

do i=1,n

do j=1,c1

res(i,j)=0.

do k=1,c

res(i,j)=res(i,j)+mat1(i,k)*mat2(k,j)

enddo

enddo