一种基于动态因子驱动的递进式多场耦合算法

一种基于动态因子驱动的递进式多场耦合算法∗何再兴，张树有，伊国栋浙江大学机械与能源工程学院，杭州(310027)E-mail ：2251hezaix@摘 要：以热——结构耦合为模型，提出了一种基于动态因子驱动的递进式多场耦合算法。

该算法首先根据实际物理量变化的大小将模型初步分裂成显式积分区域和隐式积分区域，对不同区域分别进行物理场分割计算，然后检验显式积分区域边界场量的变化大小是否超过临界值，对超过临界值的显式积分区域，根据其移动增量控制动态因子对该区域其进行二次分裂，并对划分后的隐式积分子区域进行进一步的积分运算，从而得到符合精度要求的运算结果。

本文的算法在低压断路器过载保护元件——双金属片在通电流情况下的分析中得到实际应用验证。

关键词：多场耦合，分裂，分割，动态因子，二次划分 中图分类号：TH7021. 引言由于工程实际中存在大量的多物理场的耦合作用，多场耦合(multiphysics problem)已经成为CAE 的研究热点之一[1,2]。

耦合场仿真使用的数值方法大多数以牛顿法和松弛法为基础，多场耦合算法在构造上主要有两类方法：区域分裂法与物理场分割算法。

区域分裂法主要是根据物理模型不同区域的场变化大小，分别采用隐式和显式积分，以达到在保证计算精度的前提下减小计算量的目的。

直接积分方法是多场耦合常用的数值方法，该方法主要有显式积分法和隐式积分法两种。

隐式方法的优点在于数值稳定性和能量耗散性能较好，精度高；显式积分法最大的优点是计算效率高，但是它是一种条件稳定的算法，精度低。

为了把隐式算法的数值稳定性和高精度特性与显式算法的计算效率结合起来，许多学者在这方面做了大量的研究。

Hughes [3,4]基于单元分离法提出了一种显隐式混合算法， Chiang [5]通过设计方程组并行求解器来实现显隐式混合并行计算的功能。

物理场分割算法( partitioned solution)主要是将不同物理场分割开来，在一个计算步中只含有一个物理场量。

用预测值替代某一场变量的真实值,时间积分在仅考虑一个场变量的子系统中进行,而相互作用项作为外力来考虑,从而实现不同物理场的解耦。

分割算法在不同物理场的求解分开进行,降低了计算复杂度，提高了计算效率。

主要的算法有Huang [6],[7]提出的无条件稳定的交叉迭代法，以及多次迭代方法[8]等分割算法。

分裂算法和分割算法分别通过不同的方式来减小计算量，分裂算法通过显隐式积分的有效混合来提高运算效率，但是求解耦合场方程组的运算量仍然很大；分割算法若是全部采用隐式示积分，则运算精度虽然可以保证，但是计算效率会比较低。

随着工程实际中越来越复杂的耦合场问题的出现，迫切需要一种高效与高精度的多场耦合算法。

鉴于此，本文提出了一种综合区域分裂与物理场分割算法优点的多场耦合算法，并通过实例进行了验证。

2. 热—结构耦合模型在一般三维问题中，瞬态温度的场变量()t z y x ,,,φ在直角坐标中应满足的微分方程是（边界条件应该满足：123Γ+Γ+Γ=Γ，Γ是Ω域的全部边界）：∗本课题由国家863项目（2006AA04Z114）和高等学校博士学科点专项科研基金（20040335060）资助。

0x x x ck k k Q t t x y y z z φφφφρρ⎛⎞∂∂∂∂∂∂∂⎛⎞⎛⎞−−−−=⎜⎟⎜⎟⎜⎟∂∂∂∂∂∂∂⎝⎠⎝⎠⎝⎠ (1) 边界条件：φφ= （在1Γ边界上） (2) x x y y z z k n k n k n q x y z φφφ∂∂∂++=∂∂∂ （在2Γ边界上） (3) ()a b xx y y z z k n k n k n h x y zφφφφφ∂∂∂++=−∂∂∂ （在3Γ边界上） (4) x k ,y k ,z k 分别是材料沿x ，y ，z 方向的热传导系数()()K m W ⋅；()t z y x Q Q ,,,=为物体内部热源密度()W kg ；x n ，y n ，z n 是外界法线的方向余弦。

其有限元方程：Q KT T M =+•+=P b P （5） 式中，•T 是节点温度向量，向量b P 由节点上的给定温度产生，M 是热容矩阵，K 是热传导矩阵，P 是热源、热流、对流组成的等效右端项。

结构的瞬态动力响应分析，可表述为：()t f u K u C u M m mm=++••• （6）初始条件为初始时刻给定的()0u 和()0•u 。

式中mM 为结构质量阵，mC 为结构阻尼阵。

温度-结构耦合场控制方程的离散形式可以表示为：00000000t t Mu u Ku F C K Q T T T ••••••⎧⎫⎧⎫⎧⎫⎡⎤⎡⎤⎡⎤⎧⎫⎪⎪⎪⎪⎪⎪++=⎨⎬⎨⎬⎨⎬⎨⎬⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎩⎭⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭⎩⎭（7） 其中，u 为节点位移向量；T 为节点温度向量；,,ttC K H 分别为热容矩阵和热传导矩阵；而结构场的节点荷载向量{}{}{}prnd th F F F F =++，其中{}thF 为温度应力产生的节点荷载。

3. 多场耦合算法在结构有限元分析时, 结构中总有一部分区域由于变形或应力梯度较大,而有些区域的变形和应力梯度则相对较小。

在这种情况下,显隐式混合积分法在前一部分区域作隐式积分计算, 而在后一部分区域作显式积分计算。

条件稳定的交叉迭代法是分割算法中的一种比较成熟的算法,该算法时间积分分别在仅考虑一个场变量的子系统中进行, 将每个子系统中的相互作用项作为外力来考虑。

基于区域分裂的显隐式混合积分法和物理场分割法都有比较成熟的算法，如何很好的划分显式积分和隐式积分区域是综合两种算法优点的一个难点。

本文采用先划分后修正的方法：在模型初始划分好之后，对不同区域分别进行显、隐式积分运算，对不符合精度要求的显式积分区域，进行二次划分。

3. 1 隐式积分区域二次划分在初始划分好显式和隐式积分区域，并分别进行运算求解之后，需要通过检查显式积分△是否分别大于临界值U m和T m。

如果有一项成立，则必须将区域的应变U△和温度变化T原显式积分区域进行重新划分，并且对划分后的隐式积分子区域进行重新积分运算。

△等于T m的两条临界线，理想的二次划分应该是先找出应变U△等于U m和温度变化T取能使划分后隐式积分子区域较大的一条，然后将临界线向原显式积分区域的内部推进一段距离（方向为每个节点对应的曲线内法向），边界上每一点向内推进的距离，通过公式（8）来计算（其中Pt为原显式积分区域边界BL上任意一点，Pt//为新的显隐式积分子区域的分界线BL〞上与Pt对应的点，Pt′为BL′上Pt的对应点）：L＝|| Pt′- Pt//|| =a*|| Pt – Pt′|| （8）如图1中的（a）、（b）和（c）所示，BL为原显式积分区域的边界线，其内部为显式积△等于U m和温分区域，外部为与其相邻的隐式积分区域，BL′即为Um或者Tm线（应变U △等于T m的两条临界线中的包含面积较小者），BL〞为通过公式（8）中的点运算度变化T后得到的新的分界线，参数a可取0.1到0.3之间的某一个小数。

以新的分界线BL〞边界，将原显式积分区域划分为BL和BL〞之间的隐式积分子区域和BL〞内部的显式积分子区域。

在下一步计算中对隐式积分子区域重新进行隐式积分计算，这样就解决了由于初始划分粗糙不准确造成的本应进行隐式积分的区域被划分到显式积分区域引起造成的计算精度不够的问题。

图1 三种不同显式积分区域的二次划分示意图3. 2动态因子及显式积分区域二次划分算法改进显式积分区域二次划分的理想划分可以有效解决初始划分不准确的问题，但是它的运算量是非常大的：1、需要找出BL′；2、对BL〞上的每一个节点都进行公式（8）的运算；3、对BL和BL〞之间的隐式积分子区域进行新的隐式积分。

由于第二组运算的存在，使得运算量总和非常大，因为它需要对边界线上每一点都进行运算，因而使得计算效率很低，在极端的情况下甚至会低于全隐式积分算法，所以必须进行简化，找到新的运算方式，从而提高计算效率。

显式积分区域二次划分的理想划分算法改进：对于原显式积分区域的边界上的节点，不再对每一点都进行公式（8）的运算，而是通过提取一个根据边界上的U △（或者T △）值情况而变化的动态因子，通过它来确定将BL ′向内（使隐式积分子区域增大的方向）移动的增量值的大小，通过移动BL ′得到新的边界线，将这条新的边界线作为显、隐式积分子区域的分界线。

将原来的每个点不同的移动量L(公式（8）)修改为一个统一的平移量（公式9）。

max min max min max min '**2L b c U U L L U U ⎧⎫−+⎧⎫⎪⎪⎪⎪=+⎨⎬⎨⎬+⎪⎪⎪⎪⎩⎭⎩⎭V V V V V V或者 max min max min max min '**2L b c T T L L T T ⎧⎫−+⎧⎫⎪⎪⎪⎪=+⎨⎬⎨⎬+⎪⎪⎪⎪⎩⎭⎩⎭V V V V V V (9)令max minmax minU U U U θ−=+V V V V 或者max minmax minT T T T θ−=+V V V V （10）则{}max min '**2L b c L L θ+⎧⎫⎪⎪=+⎨⎬⎪⎪⎩⎭V V （11）公式（10）中定义的θ被称为移动量控制动态因子，将边界BL ′上的每个节点沿法向移动'L 距离，就得到了显隐式积分子区域的分界线BL 〞。

图1中的（d ）、（e ）、（f ）为采用了改进算法后的二次划分边界BL 〞确定示意图。

从公式（10）可以看出，动态因子θ的大小主要受到原显式区域边界上应变的最大值maxUV 和最小值minUV（或者温度变化的最大值maxTV和最小值minTV）的影响，当BL 是一条等应变（等温度变化）曲线时，这时候θ为零，'L 最小（max min '*2L b L L +⎧⎫⎪⎪=⎨⎬⎪⎪⎩⎭V V ），可见初始划分的区域边界的形状对于显式积分区域的二次划分有着十分显著的影响，在其它条件不变的情况下，如果初始边界的形状与等应力线（或者等温度变化线）相同，可以使得在二次划分的时候BL ′到BL 〞的移动量达到最小；相反的，形状相差得越大，则移动量越大，因此，初始分割虽然是一次粗糙的划分，但是其划分的质量好坏，直接影响到计算机的处理效率。

利用动态因子，可对二次划分的算法做进一步的改进：令()1*c m U U θ=− (或()1*c m T T θ=−) （12） 不需要寻找U=U △m （或者T=T △m ）的等应变线（或者等温度变化线）BL ′，而是直接找出U=U △c （或者T=T △c ）的等应变线（或者等温度变化线），作为显式区域二次划分的分割线BL 〞。