上海大学冶金工程专业本科生课程

1、上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*本科生课程本科生课程冶金过程数值模拟冶金过程数值模拟上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 授课内容授课内容上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热先考虑忽略对流换热的导热微分方程：pTcTS式中，为介质密度，kg/m3；

2、cp为质量横压热容，J/(kg)；S为单位体积内产生的热量（内热源，产生为正，消耗为负），W/m3；为导热系数，J/(m2ms)。不同坐标系下的具体表达形式：pTTTTcSxxyyzz211pTTTTcrSrrrrzz22222111sinsinsinpTTTTcrSrrrrr其中，对于柱坐标系：x=rcos，y=rsin，z=z； 对于球坐标系：x=rsincos，y=rsincos，z=rcos。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热第一类边界条件：直接给定边界上的温度值或随时

3、间变化的函数。 000 or xxTTTT第二类边界条件：直接给定边界上的温度的导数值或导数值随时间变化的函数。 000 or xxTTDDxx第三类边界条件：直接给定边界上的温度的导数值与温度之间的函数。Lfx LTh TTx440sLwx LTTTx 44Lf0sLwx LTh TTTTx 式中，TL, Tf, Tw分别为界面、流体和环境温度，K；h为表面换热系数，W/(m2K)；0为斯芯藩玻尔兹曼常数，5.6710-8W/(m2K4)；s为表面黑度（发射率），其值取01。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题

4、数值模拟 一维稳态导热一维稳态导热x=0123i-1 i i+1nn+1x=Lxx/2x/2WEPi-1ii+1we(x)w(x)ex将求解区域n等分，得到(n+1)个节点，相应地共得到(n+1)个控制容积，从左到右依次编号。其中内节点共(n-1)个，从左到右编号分别是2到n，相应的控制容积宽度为x。边界节点共2个，均为半控制容积，左边界节点编号为1，右边界节点编号为n+1，相应控制容积宽度为x/2。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热x=0123i-1 i i+1nn+1x=L

5、xx/2x/2WEPi-1ii+1we(x)w(x)ex首先针对稳态导热问题。显然，稳态导热中温度对时间的偏微分为零，于是220d TSdx于是用差分格式代替微分格式211222iiiTTTd Tdxx代入并整理211122iiiiSTTTx这就是内节点i的离散化方程，也称内节点i的差分方程。采用同样的方法对边界节点进行分析，可以得到：恒温对流换热绝热辐射换热1wnTT212222h xSh xnnfTTTx212SnnTTx0024411f22ssxxSnnnTTTxT 上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值

6、模拟 一维稳态导热一维稳态导热x=0123i-1 i i+1nn+1x=Lxx/2x/2WEPi-1ii+1we(x)w(x)ex根据区域离散化，得到(n+1)个节点，相应地共有(n+1)个待求的未知温度，显然需要建立(n+1)个方程才能得到确定的解。对于(n-1)个内节点，我们构建的(n-1)个控制微分方程对应。再加上左右两个边界节点的差分方程就共有(n+1)个方程了。联立这(n+1)个方程组成封闭的方程组，有唯一解。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热图中一等截面直肋，处于温

7、度为80的流体中。肋表面与流体之间的表面传热系数为45W/(m2K)，肋基处温度300 。肋片端部绝热。肋片由铝合金制成，热导率为110W/(m)。肋片的厚度0.01m，高度为H=0.1m。h, Tfh, TfxHL=1Twx由于肋片的横向毕渥数Bi=h/=0.0041=eps x=G*x0+f; %迭代公式 n=n+1; tol=norm(x-x0); x0=x; if(n=M) disp(Warning:迭代次数太多，可能不收敛！); return; end fprintf(%3in,n); fprintf(%3.4fn,x); fprintf(%3.5fnnn,tol);end%在前面具

8、有gauseidel函数的基础上，在Matlab中求解的指令集：A=1 -0.332 0 0 0 0 0 0 0 0;-0.498 1 -0.498 0 0 0 0 0 0 0;0 -0.498 1 -0.498 0 0 0 0 0 0;0 0 -0.498 1 -0.498 0 0 0 0 0;0 0 0 -0.498 1 -0.498 0 0 0 0;0 0 0 0 -0.498 1 -0.498 0 0 0;0 0 0 0 0 -0.498 1 -0.498 0 0;0 0 0 0 0 0 -0.498 1 -0.498 0;0 0 0 0 0 0 0 -0.498 1 -0.498;0

9、 0 0 0 0 0 0 0 -0.992 1;b=199.674;0.326;0.326;0.326;0.326;0.326;0.326;0.326;0.326;0.649;x0=300;290;280;270;260;250;240;230;220;210;x,n=gauseidel(A,b,x0,1.0e-4)上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热对于三对角线元素组成的线性方程组：1111222221111100NNNNNNNNNbcTdabcTdabcTdabTd 当i=1

10、，由第一个方程得到：11121 2111cdTTPTQbb 当i=2，由第二个方程得到：2221232 322 122 12cda QTTPTQa Pba Pb 对任意i，有：11111iiiiiiiiiiiiiiicdaQTTPTQa Pba Pb 特殊地，对于TN有：11NNNNNNNNda QTQa Pb追追赶赶上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热雅克比法是一种迭代法，所谓迭代就是按次序逐点地由老值（前次迭代的值，带(k)上标的温度值）计算新值（后一次迭代的值，用(k+1)

11、上标表示的温度值），当全部节点的新值都计算出来后，一起用新值替代老值，作为下一轮计算的老值，这是一种逐批更新的迭代方法。 1 where ,kkkkPWWEEPPWECPPTa Ta TbaaaabSxS Tx 迭代过程可以起步于整个温度场的初始试探值，在计入边界条件后，经过反复迭代直到所有节点都满足 1kkPPTT就得到问题的收敛解。初始试探值对最终解的正确性没有影响。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维稳态导热一维稳态导热该方法也是按次序逐点由老值计算新值，但与雅克比方法相比，它在任何节点只要已

12、经计算出来新值（方程右边带(k+1)上标的温度值）便立即取代老值（方程右边带(k)上标的温度值）。所以逐点更新是高斯-赛德尔法的特点。比如，如果计算次序是由右边界向左边界推进，则： 11 where ,kkkPWWEEPPWEPCTa Ta TbaaaaSx bSx高斯-赛德尔法的优点是边界信息深入内部的速度比较快，因而收敛比较快。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一

13、维非稳态导热一维非稳态导热常物性、无内热源的一维非稳态导热控制方程：2222pTTTacxx其中，a=/(cp)为热扩散率。xii+1n-1n n+1(i+1,n-1)(i+1,n)(i+1,n+1)(i,n)对于一维非稳态导热问题，时间和空间坐标可用作平面坐标系（如图）表示。任一节点(i,n)上的温度值Tni表示空间第n个节点在i时刻的温度值。时间一阶微商采用向前差商：1,iinnn iTTT空间二阶微商采用二阶中心差商：211222iiinnnTTTTxxO2Ox代入并整理：111221 2iiiinnnnaaTTTTxx 2Ox 上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课

14、程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热差分格式图解差分方程显示差分截断误差：O()+(x)2稳定性条件：Fo=a/ (x)21/2隐式差分截断误差：O()+(x)2稳定性条件：Fo=a/ (x)2，绝对稳定六点差分（克兰克-尼克尔森格式）截断误差：O()2+(x)2稳定性条件：Fo=a/ (x)2，绝对稳定但一般也要求满足Fo0T=0oCt00 12345 6因为采用显示格式，所以从稳定性考虑确定时间步长的极限值为2260.00410 108s22 10pxc 同时考虑计算精度，我们选取=2s。上海大学上海大学冶金工程冶金工程专业本科

15、生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。板厚L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：pTTcxx初始条件：o0s,200 CT边界条件：o0s,0,00s,0,0 CxTxxT 对于内部节点2、3、4的显式差分格式为L=2cmx=0. 4cm绝热t0T=0oCt00 12345 6对于节点1

16、，采用元体平衡法时，因为西侧板面绝热，所以只有来自东面的传热1111 2iiiinnnnTFo TTFo T12111iiiipTTTTcxx12111iiiiFo TTTT11121iiiiTTFo TT对于节点5，同样采用元体平衡法时，东西侧均有传热14565552iiiiiipTTTTTTcxxx154561 32iiiiTFoTFo TFoT上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面

17、保持绝热。板厚L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：pTTcxx初始条件：o0s,200 CT边界条件：o0s,0,00s,0,0 CxTxxT 将各种数据代入得到L=2cmx=0. 4cm绝热t0T=0oCt00 12345 6同时代入T6=0，于是2272102110 0.0048pFoaxcx 111212312134231453415458786868685iiiiiiiiiiiiiiiiiiTTTTTTTTTTTTTTTTTT11112213314415

18、57 81 80001 86 81 80001 86 81 80001 86 81 80001 85 8iiiiiiiiiiTTTTTTTTTT显然，这种显示差分得到的结果在计算过程中非常方便，无需迭代，每次只需要知道前一次的计算结果，就能够代入直接计算下一时刻的结果了。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热步数时刻s(0)x=0cm1x=0.2cm2x=0.6cm3x=1.0cm4x=1.4cm5x=1.8cm(6)x=2cm0020020020020020020020012

19、200200200200200150024200200200200193.75118.75036200200200199.22185.1698.44048200200199.9197.56176.0784.670510199.99199.99199.62195.17167.3374.930612199.94199.94199.11192.24159.2667.750714199.84199.84198.36188.98151.9562.250816199.65199.65197.37185.52145.3657.90918199.37199.37196.17181.98139.4554.360

20、1020198.97198.97194.8178.44134.1351.40-0.6 -0.4 -0.2 0.00.20.40.60.81.01.21.41.61.82.0050100150200T, oCx, cm 0s 2s 4s 6s 8s 10s 50s 100s 200s 300s 400s 500s 600s 700s 800s上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。

21、板厚L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：pTTcxx初始条件：o0s,200 CT边界条件：o0s,0,00s,0,0 CxTxxT 对于内部节点2、3、4的隐式差分格式为L=2cmx=0. 4cm绝热t0T=0oCt00 12345 6对于节点1，采用元体平衡法时，因为西侧板面绝热，所以只有来自东面的传热1112111iiiipTTTTcxx1112111iiiiFo TTTT111121iiiTTFoTFo对于节点5，同样采用元体平衡法时，东西侧均有传热11

22、1114565552iiiiiipTTTTTTcxxx1115456213iiiiTFoTTFoTFo1111112iiiinnnnTFo TTTFo上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。板厚L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：pTTcxx初始条

23、件：o0s,200 CT边界条件：o0s,0,00s,0,0 CxTxxT 将各种数据代入得到L=2cmx=0. 4cm绝热t0T=0oCt00 12345 6同时代入T6=0，于是2272102110 0.0048pFoaxcx 111211112312111342311145341115465981081081081128iiiiiiiiiiiiiiiiiiiTTTTTTTTTTTTTTTTTTT1111221331441559 81 80001 810 81 80001 810 81 80001 810 81 80001 811 8iiiiiiiiiiTTTTTTTTTT 对于隐式差分，

24、上一时间层的结果直接作为下一时间层的高斯-赛德尔迭代的常数向量b，进行计算。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 一维非稳态导热一维非稳态导热步数时刻s(0)x=0cm1x=0.2cm2x=0.6cm3x=1.0cm4x=1.4cm5x=1.8cm(6)x=2cm0020020020020020020020012199.9958199.9958199.9621199.6255196.2925163.2993024199.9786199.9786199.8412198.7361190.5159136.082

25、8036199.9364199.9364199.5986197.3199183.7117115.6703048199.8554199.8554199.2073195.4296176.5295100.1720510199.7217199.7217198.6526193.1452169.363188.2490612199.5226199.5226197.9296190.5532162.440678.94840714199.247199.247197.0419187.7353155.884971.58840816198.886198.886195.9985184.7633149.75265.6781

26、0918198.4334198.4334194.8119181.6975144.057660.86201020197.8848197.8848193.4967178.587138.792856.88080-0.6 -0.4 -0.2 0.00.20.40.60.81.01.21.41.61.82.0050100150200T, oCx, cm 0s 2s 4s 6s 8s 10s 50s 100s 200s 300s 400s 500s 600s 700s 800s-0.6 -0.4 -0.2 0.00.20.40.60.81.01.21.41.61.82.0050100150200 2s 4

27、s 6s 8s 10s 300s 2s 4s 6s 8s 10s 300sT, oCx, cmblack: explicitred: implicit上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 二维导热二维导热PNSEWnswe(x)w+(x)e-x(x)n-(x)s+y二维问题的网格划分如图所示。x方向的步长为x，y方向的步长为y。控制容积的形状可以是正方形（均匀网格）、矩形

28、、三角形或者六边形等。网格的步长可以是均分，也可以不均分，是具体需求而定。如果是非稳态问题，除了空间坐标的离散，还包括时间坐标的离散，于是在原有的二维中需要增加一维，从而变为三维的离散空间。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 二维导热二维导热对于二维常物性控制方程（无热源）：2222pTTTcxy对时间的一阶微商用一阶向前差商代替，对坐标的二阶微商用二阶中心差商代替：1, ,kki ji ji j kTTT21,1,22, ,2kkkiji jiji j kTTTTxx2,1,122, ,2kkki j

29、i ji ji j kTTTTyy代入控制方程得到内部节点的显式差分方程：1,1,1,1,1,2222 where kkkkkkkki jiji jiji ji ji ji jpaaTTTTTTTTacxy1,1,1,1,1,2214 where kkkkkki jijiji ji ji jpaTFo TTTTFo TFocxx如果是均匀离散x=y，则有：上两式就是显式差分的一般性表达，很明显，每个时间层k+1的格点(i,j)温度都表达成了上一时间层k的以该格点(i,j)为中心的5个格点温度的显函数。该格式的稳定性条件是：Fo1/4。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程

30、 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 二维导热二维导热节点差 分 方 程稳定性条件内节点显式Fo1/4隐式无条件平直绝热边界节点显式Fo1/4隐式无条件平直对流边界节点显式Fo1/(4+2Bi)隐式无条件平直辐射边界节点显式Fo1/2隐式无条件二维非稳态导热问题节点差分方程（无内热源，x=y）1,1,1,1,1,1 4kkkkkki jijiji ji ji jTFo TTTTFo T11111,1,1,1,1,14kkkkkki jijiji ji ji jTFo TTTTTFo1,1,11,21 4kkkkki ji ji jiji jTFo TTTFo T1

31、111,1,11,21 4kkkkki ji ji jiji jTFo TTTTFo1,1,11,221 22kkkkki ji ji jijfi jTFo TTTFoBiTFoBiFo T1111,1,11,22124kkkkki ji ji jijfi jTFo TTTFoBiTTBiFo04421,1,11,21 2spkkkkkki ji ji jiji ji jfcxTFo TTTFo TTT 00442211111,1,11,1 22ssppkkkkkki ji ji ji jiji jfcxcxFo TTFo TTTTT 上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课

32、程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 二维导热二维导热如图所示一长400mm，宽200mm的矩形板，导热率是常数，无内热源。设板的周围各边界的温度给定，且不随时间变化，厚度方向的温度变化忽略。试求该矩形板内的温度分布。L=400mmH=200mmxy00100860100400100 xx Lyy WTTTx LT将求解区域（矩形板）的x方向8等分（m=8），y方向4等分（n=4），x=y=50mm，因各边界温度均已知或可求，故只需要计算内节点的温度即可。内节点的控制方程离散后可直接写成：,1,1,1,14i jijiji ji jTTTTT上海大学上海大学冶金

33、工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 导热问题数值模拟导热问题数值模拟 二维导热二维导热如图所示一长400mm，宽200mm的矩形板，导热率是常数，无内热源。设板的周围各边界的温度给定，且不随时间变化，厚度方向的温度变化忽略。试求该矩形板内的温度分布。L=400mmH=200mmxy00100860100400100 xx Lyy WTTTx LT1,23,21,11,31,12,11,22,21,23,22,12,11,13,12,2,11,11,1,26,15,17,16,27,16,17,20.25250.2562.50.250.25500.25100

34、5040.251000.25327.5iiiiTTTTTTTTTTTTTTTTTTTiTTTTTTT1,32,31,22,32,31,33,32,2,1,1,1,1,31,31,3,26,36,25,27,26,16,37,26,27,17,30.25500.25250.250.25250.250.250.25215i jijiji ji jiiiiTTTTTTTTTTTTTTTTTTTTTTTTTTTT5,37,36,27,36,37,2250.25240TTTTT10.25000000.2500000000000000.2510.2500.25000.2510.2500.25000.2510.2500.25000.2510.2500.25000.2510.250.250000000.2510.2500.25010.25000000.25000.250.2510.2500.25000.2500.2510.2500.25000.2500.2510.2500.25000.2500.2510.2500.