冶金过程数值模拟

1、上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*本科生课程本科生课程冶金过程数值模拟冶金过程数值模拟上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 授课内容授课内容上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类NoYes上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金

2、数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类连续性方程100组元i质量衡算ci(wi)DiRi质量浓度(质量分数)运动方程ueffFb- peff=+t热量方程cpTeffqeff=+t湍流动能eff/PrG-Pr=1.0湍流动能耗散速度eff/PrC1G/-C2/Pr=1.3; C1=1.44; C2=1.92divdivgradS uSxyzxxyyzz uuu上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类先考虑忽略非线性的通用二阶偏微分方程（PDE）：22

3、2220abcdefgxx yyxy 式中a、b、c、d、e、f、g皆是坐标的函数，而不是本身的函数。可根据下式对上式进行性质判断：24bac 如果0，PDE为双曲型方程（hyperbolic partial differential equation），如波动方程：22222uuax上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类椭圆型偏微分方程描述一类稳态问题，其变量与时间无关，要求变量在一个闭区域上的解。这类问题是边值问题。稳态导热问题、有回流的流动及对流换热等问题，其控制方程都属于椭

4、圆型问题。如考虑厚板内的一维稳态导热问题，其控制方程和边界条件分别是：000 LT xTTxxT xLT其解为： 00LTTT xTxL于是得到椭圆型PDE所具有的重要性质：(1). 区域内任意点的温度将受到两边界上温度的影响；(2). 忽略源项时，T(x)只取决于边界温度，且一定处于两边界温度之间。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类椭圆型PDE描述一类因变量与时间有关的问题，又称初值问题。其求解区域是一个开放区，计算时由初值出发，沿时间行进的方向推进，依次获得相应于给定边界条

5、件的解，其求解方法又称步进算法。考虑厚板内的一维非稳态导热，其控制方程及初始条件、边界条件分别为： 220,0 0,;,iLT xT xTTaxTT T LT其解析解为：22021,sinexpnnn xanT xTBLL式中： 002sin, 1,2,3,.Lnin xBT xTdxnLL于是，说明如下：(1). 与椭圆型PDE类似，边界温度T0、TL影响区域内任意点的温度T(x)；(2). 只要求初始条件(即=0时的条件)便能求解，不要求终了条件(=)；(3). 初始条件只影响其之后的温度，不影响其之前的温度，即影响区域内任意点的所有其后时间的温度，影响的程度随时间而减弱，且对于不同空间点

6、，影响程度也不同；(4). 稳态是=时到达的状态，这时的解变得与Ti(x)无关，于是又回到了其椭圆型PDE的空间性质一维稳态导热问题；(5). 在忽略源项时，温度取决于其初始条件和边界条件；(6). 这类问题说明，待求变量受时间影响是“一侧”单项影响，受空间影响则是“两侧”双向影响，因缘于时间影响，所以是抛物型。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类考虑如图所示的渠道中的一维流体流动。流速u为一常数，且u0。对于0，入口上有流体温度稳定为T0，控制方程和边界条件、初始条件分别为：0

7、0 ,0;0,ppic Tc uTT xT TTx其解无论是关于还是x都是双曲型。该问题的解为：0 ,0 ; or , iTx uT xTxuT xTx u说明如下：(1). 入口边界条件以一个有限的速度u向前传播；(2). 入口边界条件在x点前没有变化，直到x/u为止；(3). 区域上有的边界条件(x=0)将影响区域内的解，而区域下游的边界条件不会影响区域内的解。x=x2/u=x1/ux2x1=0T0Ti上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 偏微分方程的数学分类偏微分方程的数学分类通用标量传输方程是一个偏微分方程，由这

8、个通用的传输方程可以导出各种形式的偏微分方程。如果体系为稳态，且没有流动，控制方程主要决定于扩散项，则得到椭圆型方程，如稳态导热、稳态扩散问题；如果体系为非稳态，方程将表现为抛物型方程，如非稳态导热、非稳态扩散等；如果存在对流项，则标量传输的对流项会表现为双曲线型；对大多数工程问题，方程表现为混合型行为。Sxyzxxyyzz uuu上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法上海大学上海大

9、学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法针对求解区域，划分，原则上网格形式随意，但一般按照平行于坐标轴的方向划分；网格的交点即，我们就是用这些节点位置的通量函数值ij代替连续的通量函数(x,y)；在节点之间的中间位置寻找分界线(分界面)，这些分界线(或分界面)统称为；由界面围出的面积(体积)就称之为。区域离散化的四要素：沿坐标轴方向联结相邻节点而形成的曲线簇；：要求解未知通量的几何位置；：应用控制方程或者守恒定律的最小几何单位，控制容积的通量由节点通量代表；：规定了与各节点相对应的控制容积的分界面位置。(i,j)xi

10、yj上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法(i,j)先节点后界面；界面位于相邻两节点的正中间；代表控制容积的节点位置不一定在控制容积的正中心。先界面后节点；节点位置在控制容积的正中心；界面不一定位于相邻两节点的正中间位置。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法边界节点所代表的控制体不同。外节点法位于非角顶上的边界节点代表了半个控制体；而内节点法则可看成厚度为零的控制体的代表。x外节点法内节点法在非均匀离散中，

11、外节点法的界面是处于节点之间的中间位置，因此所计算的通量值精度较高，但节点位置不处于控制容积的正中心，所以用节点位置的值代替控制容积的值，会有一定误差；内节点法则正好相反。求解区域内如果有物性突变，则内节点法处理这种突变界面比较容易，因为内节点法容易保证控制单元内的物性均一，比如模拟相变过程。离散过程以尽量均匀为宗旨，应注意相邻控制体的厚度变化不宜过大，同一控制体各个方向的尺寸一般不宜相差太远，这些都容易降低计算准确度。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法为了书写和编程方便，必须对节点加上编号，原则

12、上编号顺序是任意的，但习惯采用与坐标轴正方向相一致的顺序编号。如用i表示x方向的节点的编号，j表示y轴方向，k表示z轴方向，因此，节点坐标为(xi,yj)的点可以直接写成(i,j)。由于i,j和k的值都随坐标轴的前进而增加，所以有下列关系111111; ; ; iiiiiijjjjjjxxxxxxyyyyyy在节点(xi,yj)上的某传输量的值，如温度，也可以写成T(i,j)或Tij。在非稳态传递过程中，还必须对时间域进行离散化处理，把连续时间离散成k。如时间步长为，则也有k+ k= k+1，式中k表示k时刻的时间步长。当然，时间步长也可以是不均匀的，把时间和空间统一起来表示，可以写成诸如Ti

13、jk或更普遍形式的ijk。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法考虑一个无非稳态项和对流项、且具有恒定扩散系数的一维标量传输方程：xxi-1ii+1x将i-1和i+1分别在节点i附件做泰勒级数展开，可得到2222331122 and 22iiiiiiiixxddddxOxxOxdxdxdxdx如果只取至1阶项，其余高阶项均舍弃，则有11 and iiiiiiddOxOxd

14、xxdxx上式中的O(x)就是用差商代替微商过程中的舍入误差，并且有0lim0 xx 220dSdx如果原式两两相减，则得到所谓一阶中心差商（之前的差分则是一阶向后差分和一阶向前差分）2112iiidOxdxx上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法考虑一个无非稳态项和对流项、且具有恒定扩散系数的一维标量传输方程：xxi-1ii+1x将i-1和i+1分别在节点i附件做泰勒级数展开，可得到2222331122 and 22iiiiiiiixxddddxOxxOxdxdxdxdx220dSdx如果两两相加，

15、则有二阶中心差商2211222iiiidOxdxx将上式代入控制方程的定义式，同时略去高阶极小项21122iiiixS这就是上述一维标量控制方程定义式在点i处的离散方程，显然，最主要的变化就是微分运算转变成为代数方程四则运算。如果空间离散x足够小，差商代替微商的误差就足够小。如果对空间离散的网格所有节点都采用类似的表达，则可以得到一组关于值的离散代数方程。：将这个无非稳态项和对流项、且具有恒定扩散系数的一维标量传输方程的离散过程扩展到二维体系，推导其二维传输方程的离散表达形式。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化

16、方法离散化方法上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法仍然考虑一个无非稳态项和对流项的一维标量传输方程：0ddSdxdxxWPEwexwxe考虑如图所示的一维网格，存储值的节点用W、P和E表示。控制容积界面用w和e表示，假设界面面积单位为1。则可以在整个控制容积P上对上式进行积分0VVdddVSdVdxdx根据高数中的奥氏公式，有wVVeddddddVndAAAdxdxdxdxdx 代入，则有10ewVddSdVdxdxA 上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金

17、数值冶金数值 数值方法数值方法 离散化方法离散化方法xWPEwexwxe对上式中的一阶导数用差分形式代替。同时，对于一维问题，假设单元的y、z方向的尺寸可取为单位1，单元的体积为x11=x，则上式可以写成10ewVddSdVdxdxA 0eEPwPWewS xxx 对上式展开，按不同节点处的值合并同类项，可得到一个通用的表达形式;PPWWEEWWwEEePWEaaabaxaxaaabS x 对于求解域内的所有控制容积都可以导出上上式，于是得到一系列代数方程。除边界节点处的离散化方程外，该方程组是一个非常特殊的三对角线性方程组，可以用各种十分有效的直接法或迭代法求解。所以为了保持方程组的线性性质

18、，往往需要对边界节点离散方程中的一些非线性性质项进行线性化处理。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法离散化过程是从对控制方程在整个控制容积上进行积分开始的，这与有限差分直接从微分方程推导不同。并且在控制容积上始终满足物理上的通量守恒，不论网格尺寸的大小。这保证了方程中的各项具有明确的物理意义。这也是有限体积法优于有限差分法和有限元法的地方。物理上的守恒不一定保证计算的精确。的解可能是不精确的，但一定满足守恒性。(d/dx)e为e界面上的扩散通量。控制容积的通量平衡是按照控制容积的界面来书写的，因此的

19、梯度d/dx必须在单元的界面上进行计算。和S在节点之间分布的假设不要求一定要设成相同。对于传输过程的通用方程中各项的离散形式取决于待求量的近似处理方法，即取决于所取的值在各节点间的分布函数关系(型线)，即所采用的插值方法。习惯上按下述方法选取。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法需要选定值随x变化的插值方法。一般采用阶梯式，即同一控制容积内各处的值相同，等于其节点的值P，因此：sPPwdxx需要选定值随的变化规律。一般采用阶梯显式，即在整个时间步长内， 值取为当前时刻的值，只有当达到下一时刻(+ )

20、时才跃迁为，于是：ewewuuduu的一阶导数值随时间的变化一般选为显示阶跃式的变化，则有ewewdxxxx上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法ewewuuduuewewdxxxx对于上两式中的部分项，我们可以表达为;22PWPWPEEPewewewuuuuuuxxxx一般假设源项S对时间和空间x均呈阶梯式变化，则有：ewSdxdSx 上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法xyPSWNEiji+1i-1j+

21、1j-1swneQnQsQwQe如图，考虑二维稳态导热问题，图中所示为求解区域内的任意微元体。假设该微元体自身无热源，从热平衡的原理来看，稳态的温度场的含意是：从邻近节点W(西)、E(东)、N(北)、S(南)(即(i-1,j)、(i+1,j)、(i,j+1)、(i,j-1)4个节点通过微元体界面w、e、n、s传给P节点(即(i,j)节点)所在微元体的热流量总和应该为零，即0wensQQQQ首先，考虑从W节点通过界面w传热Qw到节点P，根据傅立叶导热定律wwwxTQAx 对上式中的一阶微分进行差分处理，且Aw=yz=y1=y，则,1,1,1.1,i jijiji jwwijijTTTTQAyxx

22、 其中，为导热系数。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法xyPSWNEiji+1i-1j+1j-1swneQnQsQwQe同理，导出Qe、Qn、Qs，并全部代入总热流为零公式，得到1,1,1,122220iji jiji ji ji jTTTTTTxy合并同类项后得到,1,1,1,1,1,1,1,1i ji jijijijiji ji ji ji ja TaTaTaTaT其中1,1,1,1,1,1,1,1,1,1,1,1,1;1;1;1; where ;ijijiji ji ji ji ji ji

23、jijiji ji jm nmnm nm nm nm nax xax xay yay yaaaaaxxxyyy 1,1,1,1,1,1,0iji jiji ji ji ji ji jiji ji ji jTTTTTTTTx xx xy yy y 特殊地，对均匀离散网络，上式可简化为上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方法离散化方法xyPSWNEiji+1i-1j+1j-1swneQnQsQwQeWPEwexWPEWPEwexWPE分段线性分布阶梯分布控制容积法的变量（包括参量）的分布包括分段线性分布和阶梯分布（如

24、图），其中，阶梯分布因无法体现交界面上的通量梯度，一般只用于源项和物性参数在时间和空间域上的分布，而分段线性分布则主要用来计算变量的梯度，也可用于时域上的分布。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方程求解离散化方程求解控制方程的离散使得微分形式控制方程转变为离散化的代数方程组，必须对这组代数方程组进行求解才能获得的离散解，并用这套离散解代替隐含的解析解。代数方程组线性非线性方程中的系数待求函数的

25、函数方程中的系数待求函数的函数直接法迭代法上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方程求解离散化方程求解用前面的方法得到的代数方程组可以写成矩阵的形式：A = B式中，A是系数矩阵，=1,2,T为由的离散值组成的矢量，B是由源项产生的矢量。对上述矩阵公式，如果满足克拉默法则，则可以直接求逆1 = A B如果A-1存在，或者|A|0，则必有解。对于所要求解的离散化方程组，系数矩阵A是稀疏矩阵，且对于结构网格是带状矩阵。对于给定类型的方程组，例如纯扩散，矩阵是对称的，为三对角型线性方程组。可用求三对角型线性方程组的三对角

26、矩阵算法（tri-diagonal matrix algorithm, TDMA）追赶法求解。Numerical Recipes in CNumerical Recipes in Fortran77Numerical Recipes in Fortran90上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 离散化方程求解离散化方程求解在各种计算传热和计算流体力学中，迭代法是最广泛适用的求解方法。这种方法的逻辑是“预估校正”（guess and correct）原理，反复调用离散方程组来不断改进预设解的方式，并不断逼近问题的真解。比如

27、一个简单的迭代方法高斯赛德尔法（Gauss-Seidel method），其过程可简化为：YesNo预设求解区域内所有网格点处的的离散值；依次访问每一个网格点，用下式对各格点值进行校正：EEWWPPaaba遍历求解区域直到覆盖所有网格点，至此完成一次迭代；判断是否满足一个合适的收敛标准，比如网格点中值的最大变化值小于10-6，满足则停止，否则转入步骤。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 目录目录上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 差分方程的特性

28、判断差分方程的特性判断解析解是精确解，离散解是近似解，所以存在精度问题；精度的判断首先是截断误差O(x)n)，并由此来定义离散方法的阶数n；从截断误差的形式判断，阶数越高越精确，网格划分越密越精确，这有一定道理；但实际情况是“阶数越高越精确，网格划分越密越精确”这句话经常出错，主要因为：(1). 虽然整个数值解的精度取决于求解区域上的各节点离散方程的截断误差，但邻近边界的内节点，往往难于得到高截断误差的表达式；(2). 对于对流换热问题，除了上述数学上的一些误差外，更要顾及差分格式在物理特性上的表现，而这不是都能与截断误差联系起来的；(3). 过分细密的网格，会是计算机的运算次数大大增加，且误

29、差的传播会使舍入误差被过分放大。所以现在的计算采用二阶精度截断误差的差分格式就可以了。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 差分方程的特性判断差分方程的特性判断差分方程相容的数值方法是一种随着网格变得越来越细密，截断误差接近于零的方法。所谓差分格式的相容性（consistency），是指差分方程应该是微分方程组的某种近似。对于非稳态问题，无论空间和时间，都必须考虑截断误差。如果截断误差是空间步长x（或时间步长t）的某次幂，那么必定相容。有时我们会遇到方法的截断误差是O(x/t)，这时，除非x降低得比t快，否则就不能保证数

30、值方法的相容性。相容性是一个非常重要的性质，没有它就不能保证改进网格能改进数值解。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科生课程 吴永全吴永全*冶金数值冶金数值 数值方法数值方法 差分方程的特性判断差分方程的特性判断差分方程的稳定性是指差分格式在实际运算时，所得的近似解能否任意逼近差分方程的准确解。“截断误差”、“舍入误差”、“试验误差”都将影响差分方程的近似解对准确解的偏离，而差分格式的稳定性（stability）要求所得的近似解能够任意逼近准确解，即要求把误差的传播控制在可以接受的范围内。需要指出的是，稳定性是一个差分格式的固有属性，凡是稳定的差分格式，任何一个信息扰动在计算过程中被放大的程度都是有限的；而不稳定的差分格式，无论什么误差都会在计算中被不断放大，当计算时间足够长时，误差的传播会让得到的解毫无意义。求解非稳态问题的稳定性更加重要，通常稳定性分析允许我们确定时间行进时在其解中的误差残留是否有限，这有时需要依赖CFD专业人士的经验和直觉。上海大学上海大学冶金工程冶金工程专业本科生课程专业本科