24小时热门版块排行榜    

查看: 1178  |  回复: 5

springer_

木虫 (著名写手)


[交流] 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
回复此楼

» 猜你喜欢

» 抢金币啦!回帖就可以得到:

查看全部散金贴

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
简单回复
dsctg2楼
2017-02-24 11:48   回复  
springer_(金币+1): 谢谢参与
2017-02-24 14:30   回复  
springer_(金币+1): 谢谢参与
发自小木虫IOS客户端
2017-11-10 23:54   回复  
springer_(金币+1): 谢谢参与
发自小木虫Android客户端
1401022405楼
2018-01-11 19:19   回复  
springer_(金币+1): 谢谢参与
发自小木虫Android客户端
hxdtj20146楼
2018-01-27 19:42   回复  
springer_(金币+1): 谢谢参与
发自小木虫IOS客户端
相关版块跳转 我要订阅楼主 springer_ 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 今天系统多次维护,明天很可能放榜! +5 zju2000 2026-08-16 5/250 2026-08-16 23:23 by Haru815
[基金申请] 时间戳又变了8-15 +13 archvillain 2026-08-15 25/1250 2026-08-16 20:13 by zhaosm1982
[基金申请] 2027广东省杰青 +3 奶牛小黑 2026-08-15 6/300 2026-08-16 20:07 by 奶牛小黑
[基金申请] 咱们一起用铁证分析2026国家社科基金中标与否 +7 启萌科技 2026-08-12 26/1300 2026-08-16 12:35 by 启萌科技
[基金申请] 我的国基提前知道中了,可是同事的操作让我实在接受不了,怎么会有这样的人 +11 家与远方 2026-08-10 16/800 2026-08-16 10:28 by ray43
[基金申请] 有时候,自然基金真的不能太认真 (我的申报经验) +10 majunge000 2026-08-11 12/600 2026-08-16 08:18 by xli1984
[基金申请] filecode=后面第一个是大写字母 +6 wangze12014 2026-08-14 7/350 2026-08-15 20:12 by gltch
[基金申请] 小木虫上这么多卖论文的,真有人买论文么?感觉没必要啊 +12 Tide man 2026-08-10 13/650 2026-08-15 16:34 by 氺木
[基金申请] 各位道友,我要去昆明玩几天,回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
[基金申请] 是这周出结果还是下周出结果? +4 yuleib84 2026-08-11 4/200 2026-08-14 23:05 by lfy8008
[基金申请] 奇怪,两个人的filecode固定段从头到尾一模一样 +8 布布和一二 2026-08-10 11/550 2026-08-14 14:58 by Equinoxhua
[硕博家园] 读博的好处 +4 lnee 2026-08-11 4/200 2026-08-14 10:20 by ahsoarli
[基金申请] filecode +15 documentary 2026-08-10 17/850 2026-08-14 10:08 by kissu88
[基金申请] FileCode能看出啥? +10 要乐观耀哥 2026-08-10 32/1600 2026-08-14 09:37 by 要乐观耀哥
[基金申请] 重要来源:本周末出结果 +10 瞬息宇宙 2026-08-12 10/500 2026-08-13 15:46 by likettle
[基金申请] 不应该看fileCode +7 且听虎啸 2026-08-12 9/450 2026-08-13 14:27 by flydreamws
[基金申请] 结合人工智能,周易传统文化,filecode打分制来了,3分以上希望很大。 +3 Tide man 2026-08-12 4/200 2026-08-13 08:35 by ZJTJZ
[基金申请] 为什么网上很多人说本周 12号出结果 +6 瞬息宇宙 2026-08-10 7/350 2026-08-11 19:25 by Tide man
[基金申请] 国自然结果 +4 Vierhys 2026-08-10 8/400 2026-08-10 15:06 by Vierhys
[基金申请] 据悉今年马上要出结果了 +7 瞬息宇宙 2026-08-10 8/400 2026-08-10 12:42 by Vivilian
信息提示
请填处理意见