24小时热门版块排行榜

返回列表

springer_

木虫 (著名写手)

应助: 2 (幼儿园)
金币: 3760.7
帖子: 1163
在线: 164.9小时
虫号: 2748954

[交流] Fortran平面桁架有限元程序

program  truss_2D
use prep
use solve

implicit none

integer :: i,j,m,n,nel,nne,nn,nodof,edof,gdof
integer :: e2s(4)
real :: L,FN
real :: kl(4,4),kg(4,4),T(4,4)

integer,allocatable  :: elemNodes (:,:） ,nf(:,:）

real ,allocatable  :: Coords (:,:）,prop(:,:）, &
         KK(:,:）,loads(:,:）

real ,allocatable  :: nodedisp(:,:）,F(:） ,edg(:）, &
            fg(:）,fl(:）,delta(:）


open ( 10,file = 'data.txt' )
open ( 11,file = 'out.txt'  )

read (10,*)  nel    ! Number of elements
read (10,*)  nne    ! Number of nodes per element

allocate ( elemNodes(nel,nne),prop(nel,nne) )

read (10,*)  nn       ! Number of nodes
read (10,*)  nodof    ! Number of degrees of freedom per node
edof = nodof * nne

allocate ( Coords(nn,nodof),nf(nn,nodof),             &
            loads(nn,nodof),nodedisp(nn,nodof),edg(edof), &
            fg(edof),fl(edof) )

read (10,*) ( (elemNodes(i,j), j=1,nne), i=1,nel )
read (10,*) ( (prop(i,j),    j=1,nne), i=1,nel )
read (10,*) ( (Coords(i,j), j=1,nodof), i=1,nn )
read (10,*) ( (nf(i,j),       j=1,nodof), i=1,nn )
read (10,*) ( (loads(i,j),    j=1,nodof), i=1,nn )

gdof = 0
do i = 1,nn
      do j = 1,nodof
         if ( nf(i,j)/=0 ) then
            gdof = gdof + 1
            nf(i,j) = gdof
         end if
      end do
end do

allocate ( KK(gdof, gdof), F(gdof),delta(gdof) )

F = 0.
call truss_F( m,n,nn,nodof,nf,gdof,loads,F)

e2s = 0
KK  = 0.
do i = 1,nel
      call truss_T(i,nel,nne,nn,nodof,elemNodes,Coords,L,T)
      call truss_kl (i,nel,nne,L,prop,kl)
      call truss_kg (T,kl,kg)

      call truss_e2s(i,j,nn,nel,nne,nodof,nf,elemNodes,e2s)

      call form_KK (m,n,edof,gdof,kg,e2s,KK)

end do

call fem_Solver(KK,F,gdof,delta)  ! 求解

nodedisp = 0.
forall ( i = 1:nn,j = 1:nodof,nf(i,j)/=0 )
      nodedisp(i,j) = delta( nf(i,j) )
end forall

do i = 1,nel
      call truss_T(i,nel,nne,nn,nodof,elemNodes,Coords,L,T)

      call truss_kl (i,nel,nne,L,prop,kl)
      call truss_kg (T,kl,kg)

      call truss_e2s(i,j,nn,nel,nne,nodof,nf,elemNodes,e2s)

      edg = 0.
      do j = 1,edof
         if ( e2s(j) /= 0 ) then
            edg(j) = delta( e2s(j) )
         end if
      end do
      fg = matmul( kg,edg )
      fl = matmul( T,fg )
      FN = fl(3)
      write (11,100) i
      write (11,200) FN
   end do
100  format (/,T10,'单元',I2 )
200  format (T10,'轴力=',F18.4)

end program truss_2D

回复此楼