冶金过程数值模拟 Numerical Modelling of Metallurgical Processing

1. 本科生课程——冶金过程数值模拟 冶金过程数值模拟Numerical Modelling of Metallurgical Processing 主讲：吴永全 上海大学，材料科学与工程学院，材料工程系，冶金工程教研室

2. 冶金数值—— 授课内容 授课内容 绪论 数学描述 数值求解方法 冶金过程数值模拟 绪论 数学描述 数值求解方法 冶金过程数值模拟 流场计算简介 导热问题数值模拟 导热问题数值模拟 对流与扩散数值模拟 对流与扩散数值模拟

3. 一维稳态对流与扩散 一维稳态对流与扩散 1 1 对流项的其他离散格式 2 多维对流和扩散 3 虚假扩散 4 冶金数值—— 对流与换热—— 目录

4. 冶金数值—— 对流与换热—— 一维稳态对流与扩散 积累项 对流项 扩散项 源项 对流产生于流体流动，本次课的任务就是在已知流场（即速度和密度）的情况下，求得φ的解，即求由于对流和扩散同时存在时某自变量φ（如温度或浓度）的解。计算流场部分最后表达。对流项和扩散项之间具有不可分割的联系，因此需要把这两项处理成一个单位。 从数学角度，对流项不过是一阶导数项，其离散毫无问题。但从物理角度来看，这是最难离散的导数项。这主要与对流作用带有强烈的方向性有关。对流项离散方程的构建是否合适强烈地影响到数值解的准确性、稳定性和经济性。

5. (Δx)w (Δx)e w e uw ue W(i-1) P(i) E(i+1) (Δx) 冶金数值—— 对流与换热—— 一维稳态对流与扩散 考虑无源一维稳态问题，控制方程为： 其中，u代表x方向的速度。同时连续方程为： 该式对可压缩流体和不可压缩流体都成立。换言之，单位时间单位面积上流过的质量流量为常数。当流体为不可压缩流体时，ρ=const.，则有各断面上速度处处相等。 区域离散如图。 针对如图所示的控制体定解区域对控制方程积分：

6. 冶金数值—— 对流与换热—— 一维稳态对流与扩散 对于物理量φ的一阶导数用中心差商代替，于是有 式中，系数1/2为内插因子，表示假设界面处于中心。对于不同的界面位置则要采用其它的内插因子。把上式代入顶式得到： 为了把方程写的更简洁，我们定义两个新的符号： 式中，F称为对流（或流动）强度，其物理意义为单位时间单位面积上流过的流体质量，F越大，流过的质量越多，对流强度越大。F值可正可负，仅取决于流动的方向。D称为扩散传导系数（diffusion conductance），D值具有与F相同的量纲，表示扩散状况，但D值永远为正。当Γ=0时，D值也等于0。 2014年9月6日5时24分

7. 冶金数值—— 对流与换热—— 一维稳态对流与扩散 这两个变量可以组成贝克来数Pe，它是以网格间距定义的无量纲特征数： 当Re=0，意味着没有对流只有扩散；Pe=∞，表示扩散作用可以忽略，仅有对流；Pe=1，表示对流和扩散的作用相当；当Pe>>1，表示对流作用远大于扩散作用。 引入上述两个变量后，控制方程可以离散表达为： 一个显然的结果是，对流的引入并未改变最终的离散化方程形式，只是使系数计算和导热问题有些差别，所以前面涉及的离散化方程的求解方法依然适用。 2014年9月6日5时24分

8. 冶金数值—— 对流与换热—— 一维稳态对流与扩散 • 问题的引出 • 连续方程ρu=const表示Fe=Fw，这时aP=aE+aW，说明满足连续方程就确保了系数之和规则成立。但在迭代过程中，中间过程的迭代值可能不保证满足连续方程，所以上式中不能直接省掉Fe和Fw。根据差分格式稳定性的4个基本准则之一：各系数为正，发现必须满足D-(1/2)F≥0，也就是D≥|F/2|，或者|Pe|≤2。 • 也就是说，要用中心差分格式，就要保证上面的贝克来数关系，只能选取较小的网格间距，同时还违背了斯卡巴勒准则。 • 另外，中心差分还不能处理Γ=0的情况，因为满足连续方程后，aP=0，导致方程无法进行计算。 寻求其它方法或者离散格式 2014年9月6日5时24分

9. 一维稳态对流与扩散 一维稳态对流与扩散 1 1 对流项的其他离散格式 对流项的其他离散格式 2 2 多维对流和扩散 3 虚假扩散 4 冶金数值—— 对流与换热—— 目录

10. 冶金数值—— 对流与换热—— 对流项的其他离散格式 对流扩散问题的严格解（解析解） 对于对流强度和扩散系数都已知的对流扩散问题，我们可以得到严格解。如求解区域为0≤x≤L，边界条件为 其严格解为 其中，这里的Pe=(ρuL)/Γ。 由图中的精确解可以知道，当Pe=0，即纯扩散时，φ与x的关系式是线性的；对于小Pe而言，变化关系偏离线性不大。因此，在这种情况下采用中心差分格式是可行的。但随着|Pe|的增加，节点间的φ值越来越为上游节点所影响，当|Pe|>>1时，即强制对流时，节点间的大部分区域的φ值几乎就是上游节点的φ值。显然，这个时候如果再用线性分布来近似就会使所得离散方程违反正系数规则，从而导致物理上的不真实解。 2014年9月6日5时24分

11. (Δx)w (Δx)e w e uw ue W(i-1) P(i) E(i+1) (Δx) 冶金数值—— 对流与换热—— 对流项的其他离散格式 上风格式 针对中心差分的缺陷，考虑到大Pe时节点间的大部分区域的φ值基本上等于上游节点的φ值，于是在扩散项仍采用中心差分格式的同时，对流项中的φe与φw值均取上游节点的φ值，如 于是就有 这种格式的系数可以看到不会有负系数问题，从而保证了任何Pe都能得到物理上真实的解。但一个明显缺陷是：无论小Pe还是大Pe，计算误差都较大。 2014年9月6日5时24分

12. 冶金数值—— 对流与换热—— 对流项的其他离散格式 指数格式 根据严格解的构造可以定义 式中，J表示由对流流量ρuφ和扩散流量 -Γdφ/dx组成的总流量。于是，控制方程可以改写成dJ/dx=0。 将一维无源稳定严格解代入J的定义式，得到 对于控制容积面w和e，分别用φW和φP以及φP和φE替代φ0和φL，并用hW替代L，于是 代入J的控制方程离散形式得到 这个方法可以满足正系数规则，也可以保证一维无源稳态问题的严格解。但这个严格解只对一维无源稳态问题成立，对一般的对流和扩散问题（多维、有源等），将失去严格解的含义。 2014年9月6日5时24分

13. 冶金数值—— 对流与换热—— 对流项的其他离散格式 混合格式 做指数格式中系数aE/De与Pe的函数关系如右图。 发现在Pe分别趋向于-∞、 ∞、0的时候， aE/De分别趋向于-Pe、0、1-(Pe)/2。 混合格式实际上构成了对严格解的包络线，它代表了一种合理近似。混合格式就是由这三条线组成，即： 代入控制方程并整理可以得到 与指数格式一样，混合格式没有要求控制容积面一定要在两节点的正中间，可以在两节点间的任意位置。 2014年9月6日5时24分

14. u φ=1 φ=0 x=0 x=L Δx Δx n 1 2 3 i-1 i i+1 n-1 x=0 x=L Δx 冶金数值—— 对流与换热—— 对流项的其他离散格式 举例四 设场变量φ经过对流扩散过程从一维区域的x=0点传输到x=L点，如图。流体密度为ρ=1.0kg/m3，L=1.0m，扩散系数Γ=0.1kg/(m·s)。求流速分别为0.1m/s和2.5m/s时，将区域离散成10个节点网格时，φ的分布。 Step 1: 解析解 此问题的解析解为： 可采用中心差商离散。 Step 2: 求解区域离散 采用外节点法将区域离散成10等分，每个控制容 积长度为0.1m，相应的F=ρu，D=Γ/Δx，Fe=Fw=F，De=Dw=D，对所有控制容积成立。 2014年9月6日5时24分

15. u φ=1 φ=0 x=0 x=L 冶金数值—— 对流与换热—— 对流项的其他离散格式 举例四 设场变量φ经过对流扩散过程从一维区域的x=0点传输到x=L点，如图。流体密度为ρ=1.0kg/m3，L=1.0m，扩散系数Γ=0.1kg/(m·s)。求流速分别为0.1m/s和2.5m/s时，将区域离散成10个节点网格时，φ的分布。 Step 3: 控制方程离散 用中心差商离散，其离散方程为： 针对节点1和10，采用控制容积积分得到（中间过程略）： 追赶法 2014年9月6日5时24分

16. u φ=1 φ=0 x=0 x=L 冶金数值—— 对流与换热—— 对流项的其他离散格式 举例四 设场变量φ经过对流扩散过程从一维区域的x=0点传输到x=L点，如图。流体密度为ρ=1.0kg/m3，L=1.0m，扩散系数Γ=0.1kg/(m·s)。求流速分别为0.1m/s和2.5m/s时，将区域离散成10个节点网格时，φ的分布。 5点离散 Step 3: 控制方程离散 用中心差商离散，其离散方程为： 针对节点1和10，采用控制容积积分得到（中间过程略）： 20点离散 2014年9月6日5时24分

17. u φ=1 φ=0 x=0 x=L 冶金数值—— 对流与换热—— 对流项的其他离散格式 举例四 设场变量φ经过对流扩散过程从一维区域的x=0点传输到x=L点，如图。流体密度为ρ=1.0kg/m3，L=1.0m，扩散系数Γ=0.1kg/(m·s)。求流速分别为0.1m/s和2.5m/s时，将区域离散成10个节点网格时，φ的分布。 上风格式 Step 3: 控制方程离散 用中心差商离散，其离散方程为： 针对节点1和10，采用控制容积积分得到（中间过程略）： 混合格式 2014年9月6日5时24分

18. 冶金数值—— 对流与换热—— 对流项的其他离散格式 小结：几种格式的比较 • 中心差分格式的解，在小|Pe|时和严格解很符合，但当|Pe|=2的附近，误差明显增大。当|Pe|>2后，结果均已超出边界值，违反了物理上的真实性。 • 上风格式对任何Pe都能得到物理上真实解，但在整个Pe范围内都有显著误差。这是由于在小Pe时，对流项用了上游节点值计算；而在大Pe时，扩散项仍用中心差分格式离散所致。 • 混合格式所得结果在|Pe|≤2时和中心差分格式的相同，|Pe|>2时比上风格式有很大改进，因为此时将扩散项取成了零，因此和严格解符合很好。 2014年9月6日5时24分

19. 一维稳态对流与扩散 一维稳态对流与扩散 1 1 对流项的其他离散格式 对流项的其他离散格式 2 2 多维对流和扩散 多维对流和扩散 3 3 虚假扩散 4 冶金数值—— 对流与换热—— 目录

20. N n (Δx)n- P W w e E Δy (Δx)s+ s S Jn (Δx)e- (Δx)w+ Δx Jw Je Js 冶金数值—— 对流与换热—— 多维对流和扩散 二维对流和扩散问题 对于二维对流和扩散问题，通用方程的二维形式为： 写成总流量形式 将上式在如图所示的控制容积上积分，得到

21. N n (Δx)n- P W w e E Δy (Δx)s+ s S Jn (Δx)e- (Δx)w+ Δx Jw Je Js 冶金数值—— 对流与换热—— 多维对流和扩散 二维对流和扩散问题 对于二维对流和扩散问题，通用方程的二维形式为： 写成总流量形式 将上式在如图所示的控制容积上积分，得到 式中,源项已经线性化，S=SC+SPφP；并假定控制容积内的ρ和φ均可用P点的ρP和φP表示；J是控制容积面上的总流量，如Je=∫Jdy。 2014年9月6日5时24分

22. N n (Δx)n- P W w e E Δy (Δx)s+ s S Jn (Δx)e- (Δx)w+ Δx Jw Je Js 冶金数值—— 对流与换热—— 多维对流和扩散 二维对流和扩散问题 针对连续方程： 也在控制容积上进行积分，得到： 式中，F是通过控制容积面上的流量，且假定e、w、n、s的流量分别代表其各自的F（对流强度），即： 式中“0”表示时间τ的上一时间层的已知值。各对流项F的定义已知，相应的扩散传导系数D则由下式定义： 将离散的连续方程与控制方程合并，得到： 采用混合格式代入，得到： 2014年9月6日5时24分

23. 冶金数值—— 对流与换热—— 多维对流和扩散 二维对流和扩散问题 上述方程是采用混合格式进行推导得到的，相应的系数表达是混合格式。如果写成更一般的格式，可以得到： 式中的A(|Pe|)称为贝克来函数，它代表着对流扩散差分格式的形式。下表给出不同格式下的计算形式。 2014年9月6日5时24分

24. 冶金数值—— 对流与换热—— 多维对流和扩散 三维对流和扩散问题 二维问题的推导同样适用于三维情况，这里直接给出三维对流扩散问题的离散化方程： 式中 流率（对流强度）F及扩散传导系数D分别定义为 贝克来数仍然定义为：Pe=F/D。

25. 一维稳态对流与扩散 一维稳态对流与扩散 1 1 对流项的其他离散格式 对流项的其他离散格式 2 2 多维对流和扩散 多维对流和扩散 3 3 虚假扩散 虚假扩散 4 4 冶金数值—— 对流与换热—— 目录

26. 冶金数值—— 对流与换热—— 虚假扩散 • 在处理有流动存在的问题时，经常会遇到所谓的虚假扩散（或人工扩散）问题，比如因为上风格式中的系数aW和aE等系数比中心差分格式相应的系数大|F|/2，这相当于上风格式在真实的广义扩散系数中增加了一个大小为|F|/2的虚假的扩散系数。 • 虚假扩散一般产生于下列三种情况： • 离散化时，非稳态项或对流项采用一阶截断误差的格式，为此可以采用二阶上风格式（又称为QUICK格式（quadratic upstream interpolation for convertive kinematics））等高阶离散化格式加以避免； • 流动方向与网格线呈倾斜交叉状的多维对流扩散问题； • 建立差分格式时没有考虑非常数的源项的影响。

27. WW W w P e E EE i-2 i-1 i i+1 i+2 uw ue 冶金数值—— 对流与换热—— 虚假扩散 QUICK格式 QUICK格式是一种对流项的二阶上风格式，是英国Leonard于1979年提出的用于计算控制容积界面值的二阶插值格式。它利用控制容积界面两侧的紧挨着的邻近点和处于上风侧的一个远邻近节点，共三个节点的值来进行插值计算。 如图所示。 认为，界面w处的对流项中的φw值除邻近节点P、W影响外，还受其上游的远邻近节点WW影响，而对于反向流动的情况，则受节点E的影响。于是有 而扩散项中的可采用上述三点构造的拟合曲线在界面处的斜率计算，也可采用中心差商格式计算。

28. WW W w P e E EE i-2 i-1 i i+1 i+2 uw ue 冶金数值—— 对流与换热—— 虚假扩散 QUICK格式 以一维稳态对流问题为例。针对一维对流扩散问题控制容积积分方程 其中，当uw≥0,ue≥0时，式中对流项采用QUICK格式计算，有 扩散项采用二阶中心差分格式，整理得到最后的通式 2014年9月6日5时24分

29. 冶金数值—— 对流与换热—— 虚假扩散 举例五 利用QUICK格式计算举例四中当u=0.2m/s时的一维对流扩散问题。 当u=0.2m/s时，F=ρu=Fe=Fw=1.0×0.2=0.2kg/(m2·s)，D=Γ/Δx=De=Dw=0.1/0.2=0.5kg/(m2·s)，(Pe)e=(Pe)w=F/D=0.4。例四中内节点3、4的离散可以直接应用上述QUICK格式，但边界点1、2、5需要特别处理，以适应节点3、4的QUICK格式。 节点1： 控制容积西侧界面φw=φA。但计算控制容积节点1时需要用到界面西侧节点，为此我们在距离外边界Δx/2处补充一个外部镜像点O，且满足 于是，在计算边界节点1右侧的节点值可以计算成 由节点P、E、O构造的曲线在边界处的斜率为(9 φP-8φA- φE)/(3 Δx)，从而节点1的离散方程 2014年9月6日5时24分

30. 冶金数值—— 对流与换热—— 虚假扩散 举例五 利用QUICK格式计算举例四中当u=0.2m/s时的一维对流扩散问题。 当u=0.2m/s时，F=ρu=Fe=Fw=1.0×0.2=0.2kg/(m2·s)，D=Γ/Δx=De=Dw=0.1/0.2=0.5kg/(m2·s)，(Pe)e=(Pe)w=F/D=0.4。例四中内节点3、4的离散可以直接应用上述QUICK格式，但边界点1、2、5需要特别处理，以适应节点3、4的QUICK格式。 节点5： 节点2： 最后得到的矩阵方程为 2014年9月6日5时24分

31. 冶金过程数值模拟 祝大家好运！