有限元剛度矩陣求解程序（fortran）
該程序用于求解《有限元基礎教程》P43的例題。
module Maxtri_K
contains
!---------------------------------------------------------
!-----------purpose:求解剛度矩陣
subroutine solve_KM(n_num,coods_all,E,nu,t,KM)
!-----------定義子程序------------------------------------
!-----------purpose:求解所有單元的剛度矩陣----------------
implicit none
real::E,nu,t
real::KM(12,12)
real::K(6,6)
real::coods(6)
real::coods_all(12)
integer::n_num(4,3)
integer::i,j,m
integer::temp_i,temp_j
integer::temp_n_i
integer::temp_n_j
integer::temp_row
integer::temp_col
real::temp_data
KM=0
coods(:)=coods_all(1:6)
call solve_single_K(coods,E,nu,t,K)
do m=1,4,1
do i=1,6,1
if(i/2.0<=1) then
temp_i=1
else if(i/2.0>1.and.i/2.0<=2) then
temp_i=2
else
temp_i=3
end if
temp_n_i=n_num(m,temp_i) !得到行對應的節點號
!write(*,*) temp_n_i
if(mod(i,2)==0) then
temp_row=2*temp_n_i
else
temp_row=2*(temp_n_i-1)+1
endif !得到在整體剛度矩陣中的行temp_row
write(*,"('i= ',I3,' row=',I3)") i,temp_row
!-------------------------------------------------------------------
do j=1,6,1
if(j/2.0<=1) then
temp_j=1
else if(j/2.0>1.and.j/2.0<=2) then
temp_j=2
else
temp_j=3
end if
temp_n_j=n_num(m,temp_j) !得到列對應的節點號
if(mod(j,2)==0) then
temp_col=2*temp_n_j
else
temp_col=2*(temp_n_j-1)+1
endif !得到在整體剛度矩陣中的行temp_col
temp_data=K(i,j)
KM(temp_row,temp_col)=KM(temp_row,temp_col)+temp_data
write(*,"('row= ',I3,' col= ',I3)") temp_row,temp_col
end do
end do
end do
open(11,file='fout.txt')
do i=1,12,1
write(11,"(12F5.0)") KM(i,:)*4.0/E/t
end do
write(*,*) KM(1,1)*4.0/E/t
end subroutine
subroutine solve_single_K(coods,E,nu,t,K)
implicit none
real::coods(6) !單元三個節點的坐標,i,j,m
real::E,nu,t !彈性模量和泊松比，厚度
real::K(6,6) !單元剛度矩陣
real::A !單元面積
real::B(3,6) !應變轉換矩陣
real::D(3,3) !彈性矩陣
real::S(3,6) !應力轉換矩陣
integer::i,j
call A_area(coods,A)
call B_maxtri(coods,A,B)
call D_maxtri(E,nu,D)
call S_maxtri(D,B,S)
call K_maxtri(B,S,K)
do i=1,6,1
do j=1,6,1
K(i,j)=t*A*K(i,j)
end do
end do
end subroutine
subroutine A_area(coods,A)
implicit none
real::coods(6)
real::A
A=coods(3)*coods(6)+coods(1)*coods(4)+coods(5)*coods(2)-coods(3)*coods(2)-coods(5)*coods(4)-coods(1)*coods(6)
A=A/2.0
end subroutine
subroutine B_maxtri(coods,A,B)
implicit none
real::coods(6)
real::B(3,6)
real::B_ijm(3,2) !B矩陣的分量
real::temp_i(2),temp_j(2),temp_m(2)
real::ai,bi,ci
real::A
integer::i,j
do i=1,3,1
do j=1,2,1
B_ijm(i,j)=0
end do
end do
do i=1,3,1
ai=coods(3)*coods(6)-coods(5)*coods(3)
bi=coods(4)-coods(6)
ci=coods(5)-coods(3)
temp_i=coods(1:2)
temp_j=coods(3:4)
temp_m=coods(5:6)
coods(1:2)=temp_j
coods(3:4)=temp_m
coods(5:6)=temp_i
B_ijm(1,1)=bi/2.0/A
B_ijm(2,2)=ci/2.0/A
B_ijm(3,1)=ci/2.0/A
B_ijm(3,2)=bi/2.0/A
B(1:3,(2*i-1):(2*i))=B_ijm
end do
end subroutine
subroutine D_maxtri(E,nu,D)
implicit none
real::E
real::nu
real::D(3,3)
integer::i,j
do i=1,3,1
do j=1,3,1
if(j==i) then
D(i,j)=1.0
else
D(i,j)=0
endif
end do
end do
D(1,2)=nu
D(2,1)=nu
D(3,3)=(1-nu)/2.0
do i=1,3,1
do j=1,3,1
D(i,j)=E*D(i,j)/(1-nu**2)
end do
end do
end subroutine
subroutine S_maxtri(D,B,S)
implicit none
real::D(3,3)
real::B(3,6)
real::S(3,6)
integer::i,j,m
do i=1,3,1
do j=1,6,1
S(i,j)=0
do m=1,3,1
S(i,j)=S(i,j)+D(i,m)*B(m,j)
end do
end do
end do
end subroutine
subroutine K_maxtri(B,S,K)
implicit none
real::B(3,6)
real::S(3,6)
real::K(6,6)
real::BT(6,3)
real::i,j,m
do i=1,3,1
do j=1,6,1
BT(j,i)=B(i,j)
end do
end do
do i=1,6,1
do j=1,6,1
K(i,j)=0
do m=1,3,1
K(i,j)=K(i,j)+BT(i,m)*S(m,j)
end do
end do
end do
end subroutine
end module
program main
use Maxtri_K
implicit none
real::coods_all(12)
real::E,nu,t
real::KM(12,12)
integer::n_num(4,3)
E=2.0e5
nu=0
t=10
n_num=reshape((/1,2,5,3,2,4,3,5,3,5,2,6/),(/4,3/))
coods_all=(/0,2000,0,1000,1000,1000,0,0,1000,0,2000,0/)
call solve_KM(n_num,coods_all,E,nu,t,KM)
end
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
