程，且兼顾考虑邻点速度修正值对所计算点速度修正值的影响及源项与速度场之间的同步性。通过对方腔顶盖驱动流的数值模拟验证了算法及程序的可靠性和有效性，特别当网格加密到128×128 以后，SIMPLEXT 算法仍然收敛的突出优越性是SIMPLE等算法无法比拟的。
在SIMPLE 的众多改进算法中，应用较为广泛的是SIMPLER和SIMPLEC。SIMPLEC 算法考虑了求解速度修正方程的“协调性”问题，而且压力修正值不作亚松弛。SIMPLER 则是Patankar 提出的，它对压力的初值及更新采用独立求解压力方程的方法来得到，压力不作亚松弛，由压力修正方程求出的压力修正值仅用于改进速度。目前，SIMPLE 系列算法不仅在不可压缩流体Navier-Stokes 方程的数值求解中应用得非常广泛，并且也被成功地应用于可压缩流体流场的数值计算中。
3 前期的研究结果与理论分析工作
3.1电磁场的有限元计算
3.1.1 ANSYS电磁场计算的基本理论
电磁场的基本理论由麦克斯韦方程组描述。微分形式的麦克斯韦方程组如下：
?D (3-1) ??H?J??t?B ??E?? (3-2)
?t ??D?? (3-3) ??B?0 (3-4) 式(3-1)~(3-4)分别为安培全电流定律、法拉第电磁感应定律、电场的高斯定律和磁场的高斯定律，式中ρ是电荷密度，J为电流密度向量，t为时间。此外控制各个电磁参数的本构方程为：
B??H (3-5) D??E (3-6) J??E (3-7) 式(3-5)~(3-7)分别为介质的磁特性本构关系、介质的电特性本构关系和介质的的欧姆定律：式中， H 为磁场强度矢量， J 为总电流密度矢量，D 为电通密度矢量， E 为电场强度矢量， B 为磁感应强度矢量，ρ为电荷体密度。
麦克斯韦方程组中，电磁变量是相互交织在一起的，数学求解困难较大，因此实际求解电磁场问题时是引入标量电势、标量磁势和矢量磁势，将电场变量和磁场变量分离开，形成磁场偏微分方程和电场偏微分方程，可分别写为：
?2A ?A???2???J (磁场偏微分方程) (3-8)
?t2 21
?2???????2?? (电场偏微分方程) (3-9)
?t?2式中的μ和ε分别为介质的磁导率和介电常数，A 和Φ分别是定义的矢量
磁势和标量电势。采用有限元法进行数值求解，解得磁势和电势的场分布值，然后再经过转化可得到电磁场的各种物理量。
有限元法利用在每一个单元内假设的近似函数来分片地表示全求解域上待求的未知场函数。单元内的近似函数由未知场函数在单元的各个节点的数值和其插值函数来表达。这样一来，一个问题的有限元分析中，未知场函数在各个结点上的数值就称为新的未知量，从而使一个连续的无限自由度问题变成离散的有限自由度问题。一经求解出这些未知量，就可以通过插值函数计算出各个单元内场函数的近似值，从而得到整个求解域上的近似解。
3.1.2本文先期建立的二维与三维ANSYS-FEM电磁场计算模型
作为将凝固传输统一模型与三维电磁场耦合计算而进行的第一阶段工作，为确定电磁计算结果和数据转换的有效性，同时对冷坩埚电磁定向凝固系统、电磁熔炼和搅拌系统和连铸系统进行了建模和电磁场计算。其模型分别见图3-1，模型比例为1:1。 (a) (b)
(c) (b)
(d)
图3-1 电磁场计算模型 (a) 2D电磁连铸模型; (b) 3D电磁连铸模型;
(c) 3D电磁定向凝固系统; (d) 3D电磁加热和搅拌系统
22
3.2二维/三维电磁场有限元计算结果的转换
3.2.1 二维/三维电磁结果FEM?FDM/FVM转换的线性计算 3.2.1.1 ANSYS有限元计算结果的存储形式
ANSYS节点法求解电磁场得到的是节点解，静磁场计算结果的插值转换需要输出结点坐标信息(*c.lis)、单元-材质-节点信息(*mn.lis)、节点形式的磁感应强度结果(*b.lis)、有单元信息的磁感应强度结果(*e.lis)四个文件。对谐波和行波磁场，则还需焦耳热(*j.lis)和电流密度(*d.lis)两个结果文件。实际模型中由于形状复杂，无法全部用六面体单元剖分网格，在某些区域会以退化的单元形式出现。通常情况下，8节点空间等参数单元合并其中1个或几个节点，便可以退化为4~7节点的单元。在ANSYS模型中，六面体形单元只以楔形、金字塔形和四面体形3种退化单元形状来协调网格划分，如图3-2所示，不出现7节点6面体单元(即图2a中六面体合并任意相邻两点)。故一般情况下，电磁场的计算值存储于4种不同形状和不同有效节点数目的单元或单元节点上。二维问题与三维情况下类似，只是二维情况下单元只有四边形和由四边形退化成的三角形两种形状。
(a) 5 4 1 8 7 (b) 6 5 7, 8 (c) 5、6, 7, 8 (d) 3 5, 6, 7, 8 4 3, 4 3, 4 1 2 3 6 1 2 2 2 1 图3-2 ANSYS三维电磁有限元计算所用的4种单元形状 (a) 六面体单元；(b) 楔形单元；
(c) 金字塔形单元；(d) 四面体形单元
3.2.1.2 有限元计算结果的有限差转换过程
2维的FEM剖分单元为大小不均且形状不规则的三或四边形，见图3-3a)；而FDM的剖分单元为规则的矩形，节点/中心代表点位于矩形的中心，见图3-3b)；转换程序要利用FEM的结果计算出FDM中心代表点上的电磁场值，包括磁感应强度BY/BZ，焦耳热qh和感应电流密度JG。于是要判明每一个FDM中心代表点要位于哪一个FEM单元上，参见图3-3c)。判明后利用该单元的节点值及节点坐标构造插值函数，计算出该FDM中心代表点上的相应电磁场量值。
23
(a) (b) (c) 图3-3 (a) 铸件区域的有限元剖分；(b)同一区域的有限差剖分 (c) 有限差节点与有限元剖分的单元/节点之间的几何关系
空间位置关系判断及插值计算FDM中心代表点在二维空间上位于哪一个FEM单元中的判断过程如图3-4(b)、(c)所示。
A P2 P1 P2 B B C D A D P1 P3 C (a) (b) (c) 图3-4 判定FDM点与FEM单元空间关系的示意图
平面三角形、四边形的面积S3，S4计算公式如下：
S3?y2 z2y3 z3y4 z41?y1 z11???y2 z2 1 , S4??y4 z4y1 z12?y2 z2y3 z32y3 z3 1y1 z1 1?? (3-10) ?运用上面的公式时，要求三角形与四边形的顶点按逆时针方向排列，这样计算出的面积值才不为负。
对于如图3-2(a)所示的三维一次8节点六面体单元，设A为其内部任意一点，8个节点按右手定则顺序依次对其编号为i、j、k、l、m、n、p、q。根据有限元线性插值理论，电磁场量节点值在单元内沿三个坐标轴(x，y，z)方向上将呈线性变化，x方向分量u?u(x,y,z)可写成
u?a1?a2x?a3y?a4z?a5xy?a6xz?a7yz?a8xyz (3-11)
式中，a1、a2?a8为待定系数。六面体单元内任意一点A在x方向上的电磁场分量u可写成：
24
??iuA??N?????????式中，
1xiyiyjykylymynypyq1xj1xk??1xl1xm1xn1xp1xq?j??k??l??m??n??p??q????? (3-12) ???zizjzkzlzmznzpzqxiyixjyjxkykxlylxmymxnynxpypxqyqxizixjzjxkzkxlzlxmzmxnznxpzpxqzqyiziyjzjykzkylzlymzmynznypzpyqzqxiyizixjyjzjxkykzkxlylzlxmymzmxnynznxpypzpxqyqzq (3-13)
?N?为形函数，是坐标x、y、z的函数，而???是六面体单元8个节点在x方向
的电磁场分量构成的列向量，由ANSYS计算并导出。?i??q分别表示将?内的坐标值xi、yi、zi?xq、yq、zq替换为xA、yA、zA。同理，y、z方向的电磁场分量v、w的插值计算公式可类似地写出。
对于在如图3-2(d)所示的四面体一次单元，同样设A为四面体内部任意一点，4个节点按右手定则顺序依次对其编号为i、j、k、l。则A点处某一电磁场量值为：
uA??N?[uivA??N?[viwA??N?[wiujvjwjukvkwkul]T (3-14a) vl]T (3-14b) wl]T (3-14c)
式中，
?N????Ni?1?1?1V?6?1??1xixjxkxlNjyiyjykylNk?Vi?Nl????VVjVVkVVl?? (3-15) V?zi?zj?? (3-16) zk??zl?V是由i、j、k、l四点构成的四面体体积，Vi、Vj、Vk、Vl为四面体内的
A点分别取代i、j、k、l构成的四面体体积。
25