一个自己编的采用有限单元法计算温度场的例子

ilxy

发一个自己编的采用有限单元法计算温度场的例子
温度描述：一个杆件，半径为2in，长度为9in，左边的温度为200F,周围环境的温度为0F,杆件的导热率KXX=3 Btu/(h-in-F),杆件同周围环境的换热系数为h=1.0 Btu/(h-in2-F).
所有的几何条件可以作为变量输入，有限元的划分份数也可以自己控制。结果文件在output文件中。

program heatone
use imsl
implicit none
real tempend,tempsur,kxx,htran,radis
real len,elelen,danlen,area,xishuqian,kyoubu(2,2),forceyou(2),hpl
real eledanyuan(2,2),forcedanyuan(2)
integer num,i,j
real,allocatable::tempzong(,output(:,
real,allocatable::eleduo1(:,,eleduo2(:,,eleduo(:,,forceduo(
real,allocatable::elezong(:,,forcezong(
integer,allocatable::dingwei(:,

!在这个例子中，长度的单位是in
len=9.0
radis=2.0
area=(radis)**2.0*3.14

kxx=3.0
htran=1.0
tempsur=0.0
tempend=200.0

num=9
allocate(tempzong(num+1),eleduo1(2*num,2),eleduo2(2*num,2),forceduo(2*num),dingwei(num,2))
allocate(elezong(num+1,num+1),forcezong(num+1),output(num+1,2),eleduo(2*num,2))

elelen=len/(num*1.0)

danlen=len/num
xishuqian=area*kxx/danlen
hpl=htran*2*3.14*(radis)*danlen/6.0

do i=1,num
eleduo1(2*i-1,1)=xishuqian
eleduo1(2*i-1,2)=xishuqian*-1
eleduo1(2*i,1)=xishuqian*-1
eleduo1(2*i,2)=xishuqian

eleduo2(2*i-1,1)=hpl*2.0
eleduo2(2*i-1,2)=hpl
eleduo2(2*i,1)=hpl
eleduo2(2*i,2)=hpl*2.0

eleduo(2*i-1,1)=eleduo1(2*i-1,1)+eleduo2(2*i-1,1)
eleduo(2*i-1,2)=eleduo1(2*i-1,2)+eleduo2(2*i-1,2)
eleduo(2*i,1)=eleduo1(2*i,1)+eleduo2(2*i,1)
eleduo(2*i,2)=eleduo1(2*i,2)+eleduo2(2*i,2)

enddo

do i=1,2
do j=1,2
kyoubu(i,j)=0.0
enddo
enddo


!下面的部分是将kyoubi集成到总刚里面
kyoubu(2,2)=htran*area

eleduo(2*num-1,1)=eleduo(2*num-1,1)+kyoubu(1,1)
eleduo(2*num-1,2)=eleduo(2*num-1,2)+kyoubu(1,2)
eleduo(2*num,1)=eleduo(2*num,1)+kyoubu(2,1)
eleduo(2*num,2)=eleduo(2*num,2)+kyoubu(2,2)


do i=1,num
forceduo(2*i-1)=0.0
forceduo(2*i)=0.0
enddo


forceyou(1)=0.0
forceyou(2)=htran*area*tempsur

!下面的部分是将forceyou集成到总单元荷载里面
forceduo(2*num-1)=forceduo(2*num-1)+forceyou(1)
forceduo(2*num)=forceduo(2*num)+forceyou(2)


do i=1,num,1
dingwei(i,1)=i
dingwei(i,2)=i+1
enddo

do i=1,num+1,1
do j=1,num+1,1
  elezong(i,j)=0.0
  forcezong(i)=0.0
  tempzong(i)=0.0
  enddo
  enddo


do i=1,num,1

  eledanyuan(1,1)=eleduo(2*i-1,1)
  eledanyuan(1,2)=eleduo(2*i-1,2)
  eledanyuan(2,1)=eleduo(2*i,1)
  eledanyuan(2,2)=eleduo(2*i,2)

  forcedanyuan(1)=forceduo(2*i-1)
  forcedanyuan(2)=forceduo(2*i)

  elezong(dingwei(i,,dingwei(i,)= elezong(dingwei(i,,dingwei(i,)+eledanyuan
  forcezong(dingwei(i,)= forcezong(dingwei(i,)+forcedanyuan
enddo

!下面开始求解方程
write(*,100) elezong
100 format(4f10.6)
pause 1

write(*,200) forcezong
200 format(1f10.5)
pause 1

elezong(1,1)=10**5*elezong(1,1)
!注意不能重复乘第二次了
forcezong(1)=elezong(1,1)*tempend


tempzong=(elezong).ix. (forcezong)

do i=1,1+num,1
output(i,1)=(i-1.0)*danlen
output(i,2)=tempzong(i)
enddo

open(4,file='c:/nongdu.txt')
write(4,50) transpose(output)
close(4)
50 format(2f10.4)

write(*,300) tempzong
300 format(1f10.4)
pause 1

end program

ilxy

程序见附件。

敦诚

上传程序吧，我给你加分！这样无法复制

shiwencailx

顶一个！------------------------3Q

TBE_Legend

支持一下，ilxy在哈尔滨？

敦诚

我也在啊！！哈哈哈

sp600

F语言？？？

liuph54

。我真想学习一下，谢谢楼主

liuph54

希望有机会交流，我的qq：281282002，期待中。

southme007

都是哈尔滨的弟兄~~呵呵 支持一下

xtfbnu

楼主好厉害，顶一个

jinghemin

楼主好厉害

xuyong0209

学习中，顶

qingzi039

谢谢楼主  收藏了！

200802172

这个我看不懂啊，我的ＱＱ：３０４８２０８３２，有会的教教我啊

vab123

什么语言编的。。。

nlpyt

支持一下，不容易啊！！

fzj3893811

谢谢啦，我也在搞有限元，向你学习！！

dbhj2

无私的楼主，学习了

火影忍者一号

强烈支持，顶