《有限差分法-导热模拟》由会员分享,可在线阅读,更多相关《有限差分法-导热模拟(5页珍藏版)》请在金锄头文库上搜索。
1、有限差分法-导热模拟有限差分法(Finite Differential Method)是基于差分原理的一种数值计算法。其基 本思想:将场域离散为许多小网格,应用差分原理,将求解连续函数P的泊松方程的问题转 换为求解网格节点上P的差分方程组的问题。、利用有限差分法离散三维傅立叶热传导微分方程:QTdtV 2T jQx 2=aV 2T(IT)z f6(i,j,k)/#2(i+1,j,k)3(i,j-1町冬Y/0(i,j,k)4(i,j+1,k)1(i-1,j,k)密x/5(i,j,k-1)图1三维节点图解:将三维温度场域划分为足够小的正方体网格,网格之间距离为h,图一显示为节0(i,j,k)及其周
2、围的节点 l(i-l,j,k)、2(i+l,j,k)、3(i,j-l,k)、4(i,j+l,k)、5(i,j,k-l)、 6(i,j,k+l)。节点上的电位分别用T,T,T,T,T,T,T表示0123456由有限差分法得:(12)a2T、T - 2T+ TT(i -1,j,k) - 2T(i,j,k)+ T(i +1,j,k)()沁 1o亠=dx 2 x=x0h 2h 2同理:(13)月2T、T -2T + TT(i,j-1,k)-2T(i,j, k) + T(i,j +1,k)() U 304 =ay 2 y=yoh 2h 2(巴)U T5 一 2T0 + T6 = T(i,j,k 一 D
3、一 2T(i,j,k)+ T(i,j,k + D (14) dz 2 z=z0h 2h 2将时间t划分为足够小的时间段,时间节点之间的距离为g,则采用有限差分法的后向 差分法得:dTdt(15)将式(1-2)、(1-3)、(1-4)、(1-5)代入式(1-1)得:Tn (i, j, k)二TnM,j k)旦Tn (i _ 1, j, k) + (i + 1, j, k) + gh2Tn (i, j -1, k) + Tn (i, j +1, k) + Tn (i, j, k -1) + T” (i, j, k +1) - 6Tn (i, j, k) 上式整理可推出三维热传导差分公式n(L +
4、6r)Tn (i, j, k) - rT” (i -1, j, k) - rT” (i +1, j, k) - rT” (i, j -1, k)-rTn (i, j +1, k) - rT”(i, j, k -1) - rT”(i, j, k +1) = T”-1 (i, j, k)其中:r環求解完毕。二、用C或C+编写一个计算二维传热的有限差分程序。1、条件假设:初始热量密度大小均匀分布,且为1000的平面置于边界温度为0的环境 中,仿真其降温过程。平面大小为长x宽=100x50,仿真时长为10(可变)。2、网格划分:在平面上,以长为i向,以宽为j向建立坐标系,设网格间距h = 1将平面划分
5、为100 x 50个网格。在时间上,假设时间增量g = 1。3、二维热传导差分公式:参照上面三维热传导差分公式可得二维热传导差分公式:G + 4r)Tn(i, j)- rTn(i-1, j)- rTn(i +1, j)- rTn(i, j-1)- rTn(i, j +1) = Tn-1(i, j)(2-1)ag 其中:r =瓷在公式中,Tn-l(i, j)为前一时间节点下当前点的热量密度。为后续迭代,可将公式整理为:. 1 T ” (i, j)=1 + 4rag其中:r = /h 2(2-2)n-1(i, j) + rTn (i -1, j) + rTn (i +1, j) + rTn (i,
6、 j -1) + rTn (i, j + 1Jn-1(i, j) + rT n (i -1, j) + rTn (i +1, j) + rT “ (i, j -1) + rT n (i, j + 1Jk+1kk+1k4、差分方程组的求解方法:在此选择高斯一一赛德尔迭代法。对于每一个温度节点, 迭代公式:1T n (i, j)=-k+11 + 4r(2-3)ag 其中:r = Lh 2式中k为迭代次数。迭代停止条件为:maxr (i, j) - T (i, j) I i = 1,2,.100、j = 1,2,.5O、k = 1,2,3.k+1k5、算法实现:agD初始化:输入r =厉,定义实型数
7、TIN为初始热量值,实型数PERd为迭代精度。定义数组 Tdate(TIME , LENGTHI , LENGTHJ) , TIME 二 10 为离散时间长度, LENGTHI = 102 (加上两个边界条件)为I向离散量,LENGTHLJ = 52 (加上两个边 界条件)为J向离散量。三维数组每个组元存放时间TIME时刻坐标(LENGTHI, LENGTHJ) 处的热量密度。边界条件:Tdate(n, 0,j)=0ln=0,1,2,3,.j=0,1,2,3,.51Tdate(n, 101,j)=0|n=0,1,2,3,.j=0,1,2,3,.51Tdate(n,i,0)=0ln=0,1,2,
8、3,i=0,1,2,3,.101Tdate(n,i,0)=0|n=0,1,2,3,.j=0,1,2,3,.101初始条件:Tdate(0,i,j)=TIN|i=1,2,3,100, j=1,2,3,.50 定义Maxv存放迭代时的停止条件变量:max*T (I, j) - T (I, j)| I = 1,2,.100、j = 1,2,.50、k = 1,2,3.K+1k2)迭代:For t=1: TIME:以上一时间点的值作为下一个时间点迭代的初始值。Tdate(n, i, j)= Tdate(n-1, i, j)li=1,2,3,100, j=1,2,3,. .50定义中间变量TFor j=
9、1:50:For i=1:100:T = -l 1 + 4rWhile MaxvPERCI*PERCIN-1(1,j)+ rTn (I - 1,j) + rTn (I + 1,j) + rTn (I,j - 1) + rTn (I,j + 1)k +1kk+1kMaxv = max T 一 Tn (I, j)龙 11 = 1,2,3.100 j = 1,2,3,.50Tn (i, k)= TEnd forEnd forDo whileEnd for3)图形输出将C+语言编写程序迭代计算出的Tdate1010252数组输出到matlab/Tdate.txt中, 用matab导入Tdate.txt
10、中和数据,经处理后以灰度位图显示。数值计算代码见文件夹thermal中的thermal.cpp文件,图形显示代码见文件夹matlab 中的 picture.m 中。文件结构detbugthermal.neb有卩嗟分法-导热 ffi.docmatlab令picture mdat 乩 txtther mal令thermal cppthermal.sinthermal,suoPicture.m, data.txt, thermal.cpp 作为独立的三个文档发表为 picture (改为.m) .txt, tdate(产 生数据).txt, thermal(改为.CPP).txt。下载改名后创建相应文件夹,然后运行,thermal.cpp在 visual stdio 2005环境编译,picture.m在matlab7.0以上版本运行!你可以看到动态效果,good luck!展开效果nu Dnu o