1、热传导方程的数值解法及应用 主讲人: 陈鹏主要内容1.热传导方程的建立2.用有限差分法建立热差分模型3.双层玻璃中的一维热传导4.利用PDE工具箱设计面包烤盘5.利用差分模型研究浴缸水温的变化规律问题引入 思考这样一个问题： 对于体积相同形状不同的同种物体，什么形状 最能减少热量损失的？热传导方程的建立推导物体的热传导方程时，需要利用能量守恒定律和关于热传导的Fourier定律：热传导的Fourier定律定律(用自己的语言组织)： 时间内，沿某面积元 的外法线方向流过的热量 与该面积元两侧的温度变化率 成正比，比例系数为 .自然条件下温度趋于减少，所以等式右边有个负号.即： 又叫导热系数(单位
2、: )， 为该面元的外法向单位向量.dt ds dqun d d d d dnuq k e s t k u s tn k ne k2W/m对于一个封闭的体积元 ,在 时间内其内部的热量的变化为为 .通过对体积元的闭合面积分，得到：进一步对时间积分，我们可以得到从 到 时刻流入体积元内部的热量 ，即：又由高斯公式： dt dQd d d = d dQ q s k u s t 1t 2t2 21 11 ( d )d d dt tt tQ Q t k u s t 1Q 2 2 21 1 1 21 d d ( )d d d d d dt t tt t tQ k u s t k u v t k u x
3、y z t 根据初中所学的热力学公式：物体吸收的热量等于该物体的比热容、质量与温度增量的乘积.我们得到:其中： 为物体的密度.对 进一步变形，可以得到： 2 2 1= ( , , , ) ( , , , )d d dQ cm u c u x y zt u x y zt x y z 2Q212 d d d dtt uQ c x y z tt 根据热量守恒，即：我们得到：即：简写为： 1 2Q Q2 21 1d d d d = d d d dt tt t uk u x y z t c x y z tt = uk u c t 2 2 22 2 2=u k u u ut c x y z 2u a u
4、若物体内部有热源， 时间内，在 处的体积元内所产生的热量为 ，同样容易得到含热源的热传导方程：其中： .( , , , )F x y zt 2 +u a u f = Ff c dt ( , , )x y z如果时间足够长，温度不再变化，此时 ,得到稳定场方程. 无热源条件下得到Laplace 方程：有热源条件下得到Poisson 方程：但是，在绝热的边界条件下，而内部有不灭热源是无法达到热平衡的。这样，我们的热传导方程便全部建立起来了.=0u2 =0u2 =u f通过热传导方程的推导过程，我们还可以得到一个有意思 的结论： 对于体积相同形状不同的同种物体，球形是最能减少热量损失的。 因为对于体
5、积相同的同种物体，用内部热源将其从室温升高到相同的某一温度，其所需的热量是一致的，但是根据傅里叶积分公式，如果物体表面积的增加，那么它的面积分量也会增加，从而导致其散热更快，在相同时间内，其热量损失越多。这可以解释自然界的许多现象。有限差分法使得到热传导数值解成为可能有限差分法就是将带求解的区域划分为无数多个微小的网格(或称为元胞)，网格上承载着位置和温度信息，用网格上的温度近似代替物体所在位置真实的温度.网格是一个为了便于分析和理解的数据结构，在求解的过程中并不存在.有限差分法不仅可以用在求解热传导方程中，在所有的(偏)微分方程都可以用它来求解，包括复杂的电磁场矢量方程.为了简单起见，我们建
6、立二维矩形网格，将温度在二维平面对空间和时间进行差分.将 对时间进行向后差分，得： 将 对空间进行中心差分，得：u 1 10=lim n n n nij ij ij ijt u u u uut t t u 2 1, , 1,22 22 , 1 , , 1 22 22 ( ) (*)2 ( )i j i j i jx i xy j y i j i j i jx i xy j y u u uu O xx xu u uu O yy y 上式误差之所以为 的高阶无穷小可以通过泰勒公式来证明。 泰勒公式展开为佩亚诺余项形式： 同理：两项相加：便得稍微整理一下便得（*）式 2x 2 2 21, , 21=
7、 + ( )2!i j i j x i x x i xy j y y j yu uu u x x O xx x 2 2 21, , 21= ( ) ( ) ( )2!i j i j x i x x i xy j y y j yu uu u x x O xx x 2 2 21, 1, , 2+ =2 ( )i j i j i j uu u u x O xx 最后结合热传导方程，就可以得到差分形式，即下一时刻位于(i,j)号网格中，所有位置的温度 与上一时刻周围网格温度的关系为：式中的未知参数是为了书写方便和便于编程所设计的，其中：显然，对于三维空间的的网格，我们也可以类似的得到：1, , 1,
8、1, , 1 , 1(1 2 2 ) ( ) ( )n n n n n ni j i j x y x i j i j y i j i ju u r r r u u r u u 1,ni ju 2 2 2 2/ , /x yr a t x r a t y 1, , , , 1, , 1, , , 1, , 1, , , 1 , , 1(1 2 2 2 ) ( ) ( ) ( )n n n n n n n ni jk i jk x y z x i jk i jk y i j k i j k z i jk i jku u r r r r u u r u u r u u 这样的处理还没有完，由于边界的
9、情况未知，所以我们需要对边界进行特殊处理。边界条件一般分为三类：边界温度已知、边界温度的法向梯度已知、两者的线性组合已知。 第一种最简单，只要设定一个初始温度 ，之后的每一次迭代过程中，都保持边界温度只随时间有关(已知)，与网格无关，但是边界内部的网格温度需要根据边界上的温度而求得。0,i ju 对于第二种边界条件，需要将边界的每一个网格与其周围网格的平均值近似的代替该点由于网格内部热传导所形成的变化温度，如果考虑边界绝热，那么只需要考虑这一项就可以了；如果边界非绝热，我们还需要加上一项从外界流入某边界网格的热量除以物体比热和密度的乘积。 对于第三种边界条件，则处理方法是：在一次迭代过程中，依
10、次执行第一种边界条件的处理方法和第二种边界条件的处理方法。双层玻璃的功效这是一个一维热传导方程的问题。问题很简单，就是有两个其他条件完全相同的房屋，它们的差异在于窗户。其中一间是用双层玻璃，一间用双层玻璃合起来那么厚的单层玻璃，研究一下双层玻璃的隔热性能。d d墙l室内 T1 室外 T2 2d墙室内 T1 室外 T2利用前面的傅里叶定律，对于单层玻璃，单位时间流出室内的热量为：k1玻璃的导热系数k2空气的导热系数 2d墙室内 T1 室外 T2SdTTkQ 2 1211 Q1对于双层玻璃，根据热量守恒(连续性方程)得到.k1指玻璃的导热系数k2指空气的导热系数 d d墙l室内 T1 室外 T2S
11、dTTkSlTTkSdTTkQ bbaa 212112 Ta指内层玻璃的外侧临界温度Tb指外层玻璃的内侧临界温度 Q2化简得：我们对双层玻璃和单层玻璃的热量损失做一个比值： 1 22 1 1 22T TQ k Sd kl k d 1 21 1 221 21 1 212 222T Tk Sd kl k dQ T TQ kl k dk Sd 从结果可以看出，当双层玻璃的间距 时， ，即利用双层玻璃的热量损失确实比单层少，所以说双层玻璃确实能够起到比单层玻璃更好的隔热效果.为了进一步量化我们的结论，我们可以查找相关资料，代入具体数值.21 1 222QQ kl k d 0l 2 1Q Q通过查找资料
12、得知，在室温条件下k1=410-3810-3,k2=2.510-4, k1/k2=1632.保守起见，取k1/k2 =16.取 h=l/d=4, 则 Q1/Q2=0.03即双层玻璃窗与同样多材料的单层玻璃窗相比，可减少97%的热量损失。 l/dQ1/Q2 4200.060.030.02 6设计最优化烤盘在烤焙面包的过程中，其实也是有讲究的，如何精确设计烤焙面包的烤盘的形状，使得在一个烤箱中一次能烤更多的面包,但是烤焙面包不会出现局部烤焦.很显然矩形烤盘是最节省空间的，但是它的四个角最容易烤焦，而圆形烤盘是受热最均衡的，但是它无法形成镶嵌图形，从而造成空间的浪费。我们需要将其中的规律提炼出来，并
13、且用一个易于控制的参数来调配空间利用率与受热均衡度之间的关系，从而可以烤焙出满足不同人口味的面包.很显然，当烤箱温度一到达设定温度时，其温度达到恒定，所以该问题需要考虑二维平面上的热稳态泊松方程，其中的热源是根据温度来设定的。 ( ) extendk u kf h T T采用诺伊曼条件（Neumann boundary condition），即第二类边界条件.流入面包各面的法向热流密度是一个常数，即： 这样我们得到了二维稳态的热传递方程模型：0 Tk qn2 22 20( ) ( ) extendT Tk h T Tx y Tk qn对于这个例子我们采用matlab提供的工具箱来求解，在下一个
14、例子我会完全通过自己编程来求解.利用matlab提供的pdetool工具箱，我很可以得到稳态时刻，烤盘上的面包温度分布情况.其结果如图所示.从图中可以看出多边形烤盘的热量分布从边缘向中间逐渐递减，且边角温度变化异常明显。随着边数地增加，边界的平均温度逐渐趋于平稳，当边界为圆形时，热方差达到最小值。regularpolygon 弧度 边长 半径 热均值 热方差 热标准差 标幺值tetragon 1.5708 1.0000 0.7071 140.0062 0.0117 0.1080 1.0000pentagon 1.2566 0.7624 0.6485 138.1211 0.0084 0.0918
15、 0.4414hexagon 1.0472 0.6204 0.6204 137.2175 0.0075 0.0865 0.2586heptagon 0.8976 0.5246 0.6045 136.7134 0.0070 0.0839 0.1690octagon 0.7854 0.4551 0.5946 136.3997 0.0068 0.0823 0.1138nonagon 0.6981 0.4022 0.5880 136.1955 0.0069 0.0830 0.1379decagon 0.6283 0.3605 0.5833 136.0491 0.0066 0.0815 0.0862he
16、ndecagon 0.5712 0.3268 0.5799 135.9419 0.0064 0.0800 0.0345dodecagon 0.5236 0.2989 0.5774 135.8649 0.0067 0.0819 0.1000tridecagon 0.4833 0.2754 0.5754 135.8050 0.0068 0.0822 0.1103tetradecagon 0.4488 0.2554 0.5738 135.7590 0.0069 0.0831 0.1414pentadecagon 0.4189 0.2381 0.5725 135.7179 0.0073 0.0854
17、0.2207hexadecagon 0.3927 0.2230 0.5715 135.6789 0.0069 0.0832 0.1448heptadecagon 0.3696 0.2097 0.5707 135.6552 0.0068 0.0822 0.1103.circle 0 0 0.5642 135.4438 0.0062 0.0790 0上面是通过正多边形做出来的结果.但其空间占有率到底如何呢.其中用到了运筹学的优化算法，由于篇幅有限，这里不做介绍.最终我们对我们的烤盘形状做了一个改进，即把多边形的顶角用圆角代替.但是如何控制圆角的半径呢？我们知道圆角半径越大，空间利用率就越小，我们用
18、面包面积相对剩余率 来衡量面包剩余的程度与圆角半径的关系，对于不同多边形，为了统一比较，我们对 也进行标幺化处理. nn2 2 2 22 2 22 21tan 2 2 tan1 cosmin 2 1max min 1 cos 2 tannn n n n n nnn nn nn n nnr n r nr rR nA AR A n n 最后绘制不同多边形的标幺化后的面包面积相对剩余率曲线.0 0.1 0.2 0.3 0.4 0.5 0.600.10.20.30.40.50.60.70.80.91Radius Of CornerCustomer SatisfactionSatisfaction Curve 645812cricle从图中可以看出，随着圆角半径的增加， 逐渐减少，圆角半径在00.1变化不是很明显，随着圆角半径的继续增大， 急速下降。因此需要选择一个合理的圆角半径，使得 损失不大，因此我们选择0.1的圆角半径进行进一步的优化。 *n *n*n
