| 查看: 1198 | 回复: 5 | |||
[交流]
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 |
回复此楼» 猜你喜欢
40+中年人失业后，还能找到好工作吗？
已经有13人回复
-
用NBS溴化时溶剂除了CCl4外还有什么
已经有6人回复
-
申请科研项目纯靠实力能成功吗？
已经有11人回复
-
极限拉扯，如何抉择
已经有15人回复
-
二审4个多月了还是Awaiting Reviewer Scores
已经有4人回复
-
OpenAI又攻克数学难题了，离AI做有机合成还有多远？
已经有13人回复
-
课题组会帮你通过考核么？
已经有7人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有9人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有8人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有10人回复
简单回复
dsctg2楼
2017-02-24 11:48
回复
springer_(金币+1): 谢谢参与
2017-02-24 14:30
回复
springer_(金币+1): 谢谢参与
是 发自小木虫IOS客户端
2017-11-10 23:54
回复
2018-01-11 19:19
回复
springer_(金币+1): 谢谢参与
是 发自小木虫Android客户端
2018-01-27 19:42
回复
springer_(金币+1): 谢谢参与
顶 发自小木虫IOS客户端
200