冶金工程专业本科生课程课件

1、*冶金过程数值模拟Numerical Modelling of Metallurgical Processing主讲：吴永全上海大学，材料科学与工程学院，材料工程系，冶金工程教研室本科生课程冶金过程数值模拟*冶金过程数值模拟Numerical Modelling *冶金数值 授课内容授课内容绪论数值求解方法冶金过程数值模拟数学描述绪论数值求解方法冶金过程数值模拟数学描述导热问题数值模拟对流与扩散数值模拟流场计算简介导热问题数值模拟*冶金数值 授课内容授课内容绪论数值求解方法冶金过程数值*冶金数值 导热问题数值模拟 目录一维稳态导热1一维非稳态导热2二维导热3一维稳态导热1三维导热4*冶金数值

2、导热问题数值模拟 目录一维稳态导热1一维*冶金数值 导热问题数值模拟 一维稳态导热先考虑忽略对流换热的导热微分方程：式中，为介质密度，kg/m3；cp为质量横压热容，J/(kg)；S为单位体积内产生的热量（内热源，产生为正，消耗为负），W/m3；为导热系数，J/(m2ms)。不同坐标系下的具体表达形式：直角坐标柱坐标球坐标其中，对于柱坐标系：x=rcos，y=rsin，z=z； 对于球坐标系：x=rsincos，y=rsincos，z=rcos。控制方程*冶金数值 导热问题数值模拟 一维稳态导热先考虑忽略*冶金数值 导热问题数值模拟 一维稳态导热第一类边界条件：直接给定边界上的温度值或随时间变

3、化的函数。边界条件第二类边界条件：直接给定边界上的温度的导数值或导数值随时间变化的函数。第三类边界条件：直接给定边界上的温度的导数值与温度之间的函数。对流换热辐射换热对流辐射换热式中，TL, Tf, Tw分别为界面、流体和环境温度，K；h为表面换热系数，W/(m2K)；0为斯芯藩玻尔兹曼常数，5.6710-8W/(m2K4)；s为表面黑度（发射率），其值取01。*冶金数值 导热问题数值模拟 一维稳态导热第一类边界*冶金数值 导热问题数值模拟 一维稳态导热一维导热问题的数值求解第一步：求解区域离散化x=0123i-1ii+1nn+1x=Lxx/2x/2WEPi-1ii+1we(x)w(x)ex将

4、求解区域n等分，得到(n+1)个节点，相应地共得到(n+1)个控制容积，从左到右依次编号。其中内节点共(n-1)个，从左到右编号分别是2到n，相应的控制容积宽度为x。边界节点共2个，均为半控制容积，左边界节点编号为1，右边界节点编号为n+1，相应控制容积宽度为x/2。*冶金数值 导热问题数值模拟 一维稳态导热一维导热问*冶金数值 导热问题数值模拟 一维稳态导热一维导热问题的数值求解第二步：控制方程离散化x=0123i-1ii+1nn+1x=Lxx/2x/2WEPi-1ii+1we(x)w(x)ex首先针对稳态导热问题。显然，稳态导热中温度对时间的偏微分为零，于是于是用差分格式代替微分格式代入并

5、整理这就是内节点i的离散化方程，也称内节点i的差分方程。采用同样的方法对边界节点进行分析，可以得到：恒温对流换热绝热辐射换热非线性*冶金数值 导热问题数值模拟 一维稳态导热一维导热问*冶金数值 导热问题数值模拟 一维稳态导热一维导热问题的数值求解第三步：离散化方程的求解x=0123i-1ii+1nn+1x=Lxx/2x/2WEPi-1ii+1we(x)w(x)ex根据区域离散化，得到(n+1)个节点，相应地共有(n+1)个待求的未知温度，显然需要建立(n+1)个方程才能得到确定的解。对于(n-1)个内节点，我们构建的(n-1)个控制微分方程对应。再加上左右两个边界节点的差分方程就共有(n+1)

6、个方程了。联立这(n+1)个方程组成封闭的方程组，有唯一解。*冶金数值 导热问题数值模拟 一维稳态导热一维导热问*冶金数值 导热问题数值模拟 一维稳态导热一维导热问题的数值求解举例一：肋片稳态导热图中一等截面直肋，处于温度为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);

7、 x0=x; if(n=M) disp(Warning:迭代次数太多，可能不收敛！); return; end fprintf(%3in,n); fprintf(%3.4fn,x); fprintf(%3.5fnnn,tol);end%在前面具有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

8、;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.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)*冶金数值 导热问题数值模拟 一维稳态

9、导热一维导热问*冶金数值 导热问题数值模拟 一维稳态导热线性方程组的求解TDMA法(TriDiagonal Matrix Algorithm)对于三对角线元素组成的线性方程组：当i=1，由第一个方程得到：当i=2，由第二个方程得到：对任意i，有：特殊地，对于TN有：用上一步的Pi和Qi代入到下一步求出下一步的Pi+1和Qi+1以上一步系数Pi和Qi为前提，因为TN=QN，直接往上代入即求出所有的Ti追赶*冶金数值 导热问题数值模拟 一维稳态导热线性方程组*冶金数值 导热问题数值模拟 一维稳态导热线性方程组的求解雅克比(Jacobi)法雅克比法是一种迭代法，所谓迭代就是按次序逐点地由老值（前次迭

10、代的值，带(k)上标的温度值）计算新值（后一次迭代的值，用(k+1)上标表示的温度值），当全部节点的新值都计算出来后，一起用新值替代老值，作为下一轮计算的老值，这是一种逐批更新的迭代方法。迭代过程可以起步于整个温度场的初始试探值，在计入边界条件后，经过反复迭代直到所有节点都满足就得到问题的收敛解。初始试探值对最终解的正确性没有影响。*冶金数值 导热问题数值模拟 一维稳态导热线性方程组*冶金数值 导热问题数值模拟 一维稳态导热线性方程组的求解高斯-赛德尔(Gauss-Seidel)迭代法该方法也是按次序逐点由老值计算新值，但与雅克比方法相比，它在任何节点只要已经计算出来新值（方程右边带(k+1)

11、上标的温度值）便立即取代老值（方程右边带(k)上标的温度值）。所以逐点更新是高斯-赛德尔法的特点。比如，如果计算次序是由右边界向左边界推进，则：高斯-赛德尔法的优点是边界信息深入内部的速度比较快，因而收敛比较快。*冶金数值 导热问题数值模拟 一维稳态导热线性方程组*冶金数值 导热问题数值模拟 目录一维稳态导热1一维非稳态导热2二维导热3一维稳态导热1一维非稳态导热2三维导热4*冶金数值 导热问题数值模拟 目录一维稳态导热1一维*冶金数值 导热问题数值模拟 一维非稳态导热常物性、无内热源的一维非稳态导热控制方程：控制方程其中，a=/(cp)为热扩散率。区域离散化xii+1n-1nn+1(i+1,

12、n-1)(i+1,n)(i+1,n+1)(i,n)对于一维非稳态导热问题，时间和空间坐标可用作平面坐标系（如图）表示。任一节点(i,n)上的温度值Tni表示空间第n个节点在i时刻的温度值。控制方程离散化时间一阶微商采用向前差商：空间二阶微商采用二阶中心差商：代入并整理：以x为特征长度的傅里叶数，称为网格傅里叶数，用Fo表示。*冶金数值 导热问题数值模拟 一维非稳态导热常物性、*冶金数值 导热问题数值模拟 一维非稳态导热控制方程差分格式差分格式图解差分方程显示差分截断误差：O()+(x)2稳定性条件：Fo=a/ (x)21/2隐式差分截断误差：O()+(x)2稳定性条件：Fo=a/ (x)2，绝

13、对稳定六点差分（克兰克-尼克尔森格式）截断误差：O()2+(x)2稳定性条件：Fo=a/ (x)2，绝对稳定但一般也要求满足Fo1x(i+1,n-1)(i+1,n)(i+1,n+1)(i,n)x(i,n-1)(i,n)(i,n+1)(i+1,n)x(i,n-1)(i,n)(i,n+1)(i+1,n-1)(i+1,n)(i+1,n+1)显式差分：【优点】是Tni+1表达式的右端只涉及i时间层的温度，Tn-1i、Tni、Tn+1i，因此只要知道上层节点温度，就可以逐一求解差分方程；【致命缺点】是计算式中的Fo必需不大于1/2，否则计算会出现计算格式的稳定性问题，即不收敛或计算数据出现振荡。隐式差分

14、：【优点】是绝对（无条件）稳定；【缺点】是每一步都需要方程组联立求解，计算量大。六点差分：数学上绝对（无条件）稳定，物理意义上还需要满足Fo0T=0oCt00123456因为采用显示格式，所以从稳定性考虑确定时间步长的极限值为同时考虑计算精度，我们选取=2s。*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 一维非稳态导热数值求解举例二：无限大薄板非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。板厚L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一

15、合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：初始条件：边界条件：Step 2: 控制方程离散化有限差分法对于内部节点2、3、4的显式差分格式为L=2cmx=0. 4cm绝热t0T=0oCt00123456对于节点1，采用元体平衡法时，因为西侧板面绝热，所以只有来自东面的传热对于节点5，同样采用元体平衡法时，东西侧均有传热*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 一维非稳态导热数值求解举例二：无限大薄板非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。板厚L=2cm，热导

16、率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：初始条件：边界条件：Step 3: 线性方程组的获得及其求解将各种数据代入得到L=2cmx=0. 4cm绝热t0T=0oCt00123456同时代入T6=0，于是显然，这种显示差分得到的结果在计算过程中非常方便，无需迭代，每次只需要知道前一次的计算结果，就能够代入直接计算下一时刻的结果了。*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 一维非稳态导热数值求解举例二：无限大薄板非稳态导热步数时刻s(0)x=0c

17、m1x=0.2cm2x=0.6cm3x=1.0cm4x=1.4cm5x=1.8cm(6)x=2cm0020020020020020020020012200200200200200150024200200200200193.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.250816

18、199.65199.65197.37185.52145.3657.90918199.37199.37196.17181.98139.4554.3601020198.97198.97194.8178.44134.1351.40*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 一维非稳态导热数值求解举例二：无限大薄板非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。板厚L=2cm，热导率为10W/(mK)，cp=10106J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s

19、等时刻板内的温度分布。控制方程：初始条件：边界条件：Step 2: 控制方程离散化隐式差分法对于内部节点2、3、4的隐式差分格式为L=2cmx=0. 4cm绝热t0T=0oCt00123456对于节点1，采用元体平衡法时，因为西侧板面绝热，所以只有来自东面的传热对于节点5，同样采用元体平衡法时，东西侧均有传热*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 一维非稳态导热数值求解举例二：无限大薄板非稳态导热一无限大薄板，初始状态为均匀温度200oC在某一时刻=0，板东侧面温度突然降到0oC，西侧面保持绝热。板厚L=2cm，热导率为10W/(mK)，cp=1010

20、6J/(m3K)。采用有限差分法的显示格式，选用一合理的时间步长，计算=10s、20s等时刻板内的温度分布。控制方程：初始条件：边界条件：Step 3: 线性方程组的获得及其求解将各种数据代入得到L=2cmx=0. 4cm绝热t0T=0oCt00123456同时代入T6=0，于是对于隐式差分，上一时间层的结果直接作为下一时间层的高斯-赛德尔迭代的常数向量b，进行计算。*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 一维非稳态导热数值求解举例二：无限大薄板非稳态导热步数时刻s(0)x=0cm1x=0.2cm2x=0.6cm3x=1.0cm4x=1.4cm5x=1

21、.8cm(6)x=2cm0020020020020020020020012199.9958199.9958199.9621199.6255196.2925163.2993024199.9786199.9786199.8412198.7361190.5159136.0828036199.9364199.9364199.5986197.3199183.7117115.6703048199.8554199.8554199.2073195.4296176.5295100.1720510199.7217199.7217198.6526193.1452169.363188.2490612199.52261

22、99.5226197.9296190.5532162.440678.94840714199.247199.247197.0419187.7353155.884971.58840816198.886198.886195.9985184.7633149.75265.67810918198.4334198.4334194.8119181.6975144.057660.86201020197.8848197.8848193.4967178.587138.792856.88080*冶金数值 导热问题数值模拟 一维非稳态导热数值求解*冶金数值 导热问题数值模拟 目录一维稳态导热1一维非稳态导热2二维导热3

23、一维稳态导热1一维非稳态导热2三维导热4二维导热3*冶金数值 导热问题数值模拟 目录一维稳态导热1一维*冶金数值 导热问题数值模拟 二维导热1). 区域离散化PNSEWnswe(x)w+(x)e-x(x)n-(x)s+y二维问题的网格划分如图所示。x方向的步长为x，y方向的步长为y。控制容积的形状可以是正方形（均匀网格）、矩形、三角形或者六边形等。网格的步长可以是均分，也可以不均分，是具体需求而定。如果是非稳态问题，除了空间坐标的离散，还包括时间坐标的离散，于是在原有的二维中需要增加一维，从而变为三维的离散空间。*冶金数值 导热问题数值模拟 二维导热1). 区域离*冶金数值 导热问题数值模拟

24、二维导热2). 控制方程的离散有限差分法对于二维常物性控制方程（无热源）：对时间的一阶微商用一阶向前差商代替，对坐标的二阶微商用二阶中心差商代替：代入控制方程得到内部节点的显式差分方程：如果是均匀离散x=y，则有：上两式就是显式差分的一般性表达，很明显，每个时间层k+1的格点(i,j)温度都表达成了上一时间层k的以该格点(i,j)为中心的5个格点温度的显函数。该格式的稳定性条件是：Fo1/4。改成全隐式格式该怎么办？*冶金数值 导热问题数值模拟 二维导热2). 控制方*冶金数值 导热问题数值模拟 二维导热3). 边界条件的离散节点差 分 方 程稳定性条件内节点显式Fo1/4隐式无条件平直绝热边

25、界节点显式Fo1/4隐式无条件平直对流边界节点显式Fo1/(4+2Bi)隐式无条件平直辐射边界节点显式Fo1/2隐式无条件二维非稳态导热问题节点差分方程（无内热源，x=y）*冶金数值 导热问题数值模拟 二维导热3). 边界条*冶金数值 导热问题数值模拟 二维导热4). 举例三如图所示一长400mm，宽200mm的矩形板，导热率是常数，无内热源。设板的周围各边界的温度给定，且不随时间变化，厚度方向的温度变化忽略。试求该矩形板内的温度分布。L=400mmH=200mmxyStep 1: 生成计算网格将求解区域（矩形板）的x方向8等分（m=8），y方向4等分（n=4），x=y=50mm，因各边界温度均已知或可求，故只需要计算内节点的温度即可。Step 2