有限差分法-导热模拟有限差分法(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)及其周围的节点 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 一 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 + 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将平面划分为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, 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 =厉，定义实型数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,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=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中和数据，经处理后以灰度位图显示。数值计算代码见文件夹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