技术摘要：
本发明提供的一种基于稳定广义有限元的界面问题模拟方法，包括以下步骤：根据界面问题GFEM/XFEM构建具有集合富集策略的GFEM；对GFEM富集策略中的富集函数进行修改，构造新型的SGFEM基本函数；引入单边距离函数，对SGFEM基本函数进行调整；基于调整后的SGFEM基本函数对  全部
背景技术：
在许多工程应用和物理现象中都会出现界面问题，譬如流固耦合、双材料、多相 流、多孔介质、三维纤维复合材料和自由边界问题等等[1]Brenner  H .Interfacial transport  processes  and  rheology .Elsevier，2013 .[2]Li  Z，Ito  K .The  immersed  interface method：numerical  solutions  of  PDEs  involving  interfaces  and  irregular  domains.SIAM，  2006.这些问题的解在界面处通常是不连续的。对于复杂界面 问题，数值模拟技术是不可或缺的。以二维椭圆和抛物界面问题为例进行本文的说明，三维 问题可以用类似的思想解决。 令 是一个有界区域， 是其边界。界面Γ将Ω分割成两个区域Ω0和Ω1，当 处理光滑界面问题，即Γ是光滑的，不含有角点，如图1所示。考虑椭圆界面问题： 该问题的解满足界面条件： [u]Γ＝0，在Ω上       (2) 在 上 其中 是边界 的单位外法向量， 是界面处的单位法向量，[u]Γ＝u1-u0表示两 个区域物理量在界面处的跳跃。a(P)为分片常数函数， 其中a0，a1为界面常数，满足0＜ζ0≤ai≤ζ1＜∞，i＝0，1和 显然a(P)在界 面Γ上不连续。条件(2)表明真解在界面处连续，但是 在界面处是不连续的，因此真解的 光滑性很低，给数值模拟方法带来挑战。 传统有限元方法(finite  element  method，FEM)的逼近空间 的如下： 其中Ih为网格节点的指标集，φi为有限元中形函数，如果网格是三角形， 的 函数为分片线性函数，如果是四边形， 的函数为分片双线性函数。传统FEM由于多项式 函数在描述不连续特性上的不足，面对诸如此类界面问题时，必须沿界面进行网格匹配或 者细化网格才能有效地模拟问题[3]Lin  T，Lin  Y，Zhang  X .Partially penalized  immersed  finite  element  methods  for  elliptic  interface  problems .SIAM  J .  6 CN 111611738 A 说　明　书 2/14 页 Numer.Anal.，2015，53：1121-1144.，特别在界面随时间演化的复杂工程问题中，网格划分 负担很重，难以适应工程计算的需要。 近十几年来，基于非拟合网格的有限元法获得了重要进展，这种方法的网格划分 简单且网格是固定的，其网格所在位置与界面无关。广义或扩展有限元方法  [4]Babuska  I，Banerjee  U，Osborn  J，Survey  of  meshless  and  generalized  finite  element methods：a  unified  approach，Acta  Numer，2003，12∶1-125.[5]FriesT-P，Belytschko  T.  The  extended/generalized  finite  element  method：an  overview  of  the  method  and  its applications.Internat.J.Numer.Methods  Engrg .，2010，84：253-304.[6]Babuska  I ，  M e l e n k  J .T h e  p a r t i t i o n  o f  u n i t y  f i n i t e  e l e m e n t  m e t h o d ， Internat.J .Numer.Methods  Engrg .  1997，40：727-758.(generalized/extended  FEM， GFEM/XFEM)属于非拟合有限元法的一种，它在固定规则网格基础上，耦合能够表达真解局 部信息的特殊函数，从而提高对复杂非光滑问题的逼近精度或满足对特定问题的特殊逼近 要求，这些特殊函数称为富集函数。GFEM/XFEM在过去20多年来引起了广泛的兴趣，且成功 应用于诸多工程问题，譬如材料模型、裂纹扩展和固液相互作用等等[4，5]。下面的叙述中 为了简单起见，使用GFEM来表示GFEM/XFEM。基于GFEM对单元形函数构造理论的深入研究， 具有任意内部特征(如裂纹、空洞等)及外部特征(如凹角、棱边等)的复杂问题，都可以在简 单且与区域无关的有限元网格上进行求解[4，5]。GFEM  很早就应用于求解固定界面和移动 界面问题[7]Kergrene  K，Babuska  I，Banerjee  U.  Stable  Generalized  Finite  Element  Method  and  associated  iterative  schemes：application to  interface  problems.Comput.Methods  Appl.Mech.Engrg.，2016，305：1-36.[8]Chessa J，Belytschko  T.An  extended  finite  element  method  for  two-phase  fluids.J.Appl.Mech.，  2003， 70：10-17.[9]   N，Cloirec  M，Cartraud  P，Remacle  J .A  computational  approach  to  handle  complex  microstructure  geometries.Comput.Methods  Appl.Mech.Engrg，  2003，192 ：3163-3177.[10]Sauerland  H，Fries  T-P.The  stable  XFEM  for two- phase  flows.Comput.&Fluids，2013，87：41-49.[11]Cheng  K，Fries  T-P，  Higher-order  XFEM  for  curved  strong  and  weak  discontinuities，Internat.J .Numer .  Methods  Engrg .，2010，82：564-590 .[12]Zhang  Q，Babuska  I .A  stable  generalized finite  element  method  (SGFEM)  of  degree  two  for  interface  problems .Comput.Methods Appl.Mech.Engrg.，2020，363：112889。 当前常用的两种常规GFEM：拓扑广义有限元法(topological  GFEM)和修正广义有 限元法(corrected  GFEM)。拓扑GFEM的逼近空间为 其中富集函数D为关于界面Γ的距离函数，定义为D(P)＝dist(P，Γ)， 是 富集节点的指标集，定义为 es代表网格单元。[11]  中分 析拓扑GFEM中由于混合单元(blending  elements)会导致较差的逼近结构，他们在[11]中 提出了修正GFEM，其逼近空间为。 7 CN 111611738 A 说　明　书 3/14 页 其中 定义 ξ为斜坡函数 ，定义为 当前这些GFEM在求解界面问题时有如下几个主要难点：第一，GFEM可能缺少稳定 性，其刚度矩阵的条件数比FEM的条件数大得多[13]Babuska  I  and  U .  Banerjee  U.Stable  generalized  finite  element  method，Comput.Methods  Appl.Mech.  Engrg.， 2011(201-204)：91-111.；第二，GFEM可能缺少鲁棒性，当界面接近单元边界时，条件数会变 得病态[14]Zhang  Q，Babuska  I，Banerjee  U .Robustness  in  stable  generalized  finite  element  methods(SGFEM)applied  to  Poisson  problems  with  crack  singularities，Comput.Methods  Appl.Mech.Engrg.，2016，311：476-502。稳定性和鲁棒性 的缺失使得GFEM的刚度矩阵的条件数很大，从而导致消元法中灾难性的舍入误差，或者迭 代求解线性方程组时收敛速度会变得非常缓慢。此外，由于引进了有限元形函数之外的富 集函数，导致很多有限元的标准操作不能适用于GFEM，例如质量集中[15]Thomee  V.Galerkin  finite  element  methods  for  parabolic  problems.Springer，1984.[16] Menouillard  T，RéthoréJ，Combescure  A，et  al.Efficient  explicit  time  stepping  for  the  eXtended  Finite  Element  Method(X-FEM) .Internat.J.Numer.Methods  Engrg，  2006，68(9)：911-939.[17]Elguedj  T，Gravouil  A，Maigre  H .An  explicit  dynamics  extended  finite  element  method .Part  1：Mass  lumping  for  arbitrary  enrichment  functions.  Comput.Methods  Appl.Mech.Engrg，2009，198(30-32)：2297-2317.
技术实现要素：
本发明为克服现有的GFEM在求解界面问题时存在的稳定性和鲁棒性缺失从而导 致消元法中灾难性的舍入误差，或者迭代求解线性方程组时收敛速度会变得非常缓慢的技 术缺陷，提供一种基于稳定广义有限元的界面问题模拟方法。 为解决上述技术问题，本发明的技术方案如下： 基于稳定广义有限元的界面问题模拟方法，包括以下步骤： S1：根据界面问题GFEM/XFEM构建具有集合富集策略的GFEM； S2：对GFEM富集策略中的富集函数进行修改，构造新型的SGFEM基本函数； S3：引入单边距离函数，对SGFEM基本函数进行调整； S4：基于调整后的SGFEM基本函数对界面问题进行求解，完成界面问题的模拟。 上述方案中，本发明针对界面问题提出了一种SGFEM，既提高了传统GFEM的精度， 还大幅降低刚度矩阵的条件数，使得界面问题的数值模拟更加精确和鲁棒。本方法适用于 任意相对位置的网格和界面。 上述方案中，本发明开发的新型SGFEM可以改进传统GFEM/XFEM的精度和条件数， 使得界面问题的数值模拟更加精确和鲁棒，由于界面问题广泛存在于工程应用和物理现象 中，所以本发明所述的方法结果对解决工程应用和物理现象中的问题具有很高的应用价 值。 其中，所述步骤S1具体为： 令 是一个有界区域， 是其边界，界面Γ将Ω分割成两个区域Ω0和Ω1，处 8 CN 111611738 A 说　明　书 4/14 页 理光滑界面问题时，即Γ是光滑的，不含有角点；考虑椭圆界面问题： 在 上 该问题的解满足界面条件： [u]Γ＝0，在Ω上       (2) 其中 是边界 的单位外法向量， 是界面处的单位法向量，[u]Γ＝u1-u0表示两 个区域物理量在界面处的跳跃；a(P)为分片常数函数， 其中a0，a1为界面常数，满足0＜ζ0≤ai≤ζ1＜∞，i＝0，1和 显然a(P)在界 面Γ上不连续；条件(2)表明真解在界面处连续，但是 在界面处是不连续的，因此真解的 光滑性很低，给数值模拟方法带来挑战；一般有限元方法FEM的逼近空间 表示如下： 其中Ih为网格节点的指标集，φi为有限元中形函数，如果网格是三角形， 的 函数为分片线性函数，如果是四边形， 的函数为分片双线性函数；而广义或扩展有限元 方法GFEM/XFEM属于非拟合有限元法的一种，它在固定规则网格基础上，耦合能够表达真解 局部信息的特殊函数，从而提高对复杂非光滑问题的逼近精度或满足对特定问题的特殊逼 近要求，这些特殊函数称为富集函数；因此，构建具有集合富集策略的拓扑广义有限元法 GFEM。 其中，所述拓扑GFEM的逼近空间为 其中富集函数D为关于界面Γ的距离函数，定义为D(P)＝dist(P，Γ)， 是富集节点的指标集，定义为 es代表网格单元；对拓扑 GFEM进行修正，修正后的GFEM逼近空间为： 其中 定义 ξ为斜坡函数 ，定义为 其中，所述步骤S2具体为通过增加一个简单的减去插值项的局部处理步骤对  GFEM的富集函数D进行替换，从而构造新型的SGFEM基本函数。 其中，在所述步骤S2中，增加一个简单的减去插值项的局部处理步骤具体为： 对于光滑界面问题，用 代替D，其中 表示连续函数f基于区域Ω的有限 元插值函数；基于这种改进后的GFEM称为SGFEM[13，14]；SGFEM的逼近空间定义如下： 9 CN 111611738 A 说　明　书 5/14 页 其中 与拓扑GFEM(4)中的相同；对于富集函数减掉它们的有限元插值，这是减 轻其与有限元函数线性相关性的举措之一[12-14]。 其中，所述步骤S3具体为： 不同于传统GFEM中的距离函数D(P)：＝dist(P，Γ)(Γ表示界面)，引入新的单边 距离函数 具体表示为： 其在Γ的一侧为D(P)，另一侧为0；这种单边距离函数比距离函数D的表达式更简 单，因此数值积分和算法都得到了简化。而且，该单边距离函数使得SGFEM最佳收敛性的证 明更加方便[12]。采用单边距离函数的SGFEM的逼近空间定义如下： 调整后的SGFEM基本函数可以保持最佳的收敛性，还可以使得刚度矩阵的尺度条 件数与限元法同阶。 其中，在所述步骤S3中，为了便于展示结果，对GFEM的刚度矩阵尺度条件数  SCN进 行定义： 令A为FEM，GFEM，或SGFEM的刚度矩阵，D为一个对角阵，定义如为 则A的尺度条件数定义如下： K：＝κ2(DAD)        (9) 其中κ2为对称矩阵的条件数，为该对称矩阵的最大特征值与最小非零特征值之 比；FEM的刚度矩阵条件数K -2FEM＝O(h )；式子(8)的刚度矩阵尺度条件数KSGFEM跟有限元同 阶，且具备鲁棒性，这解决了常规GFEM的稳定性了鲁棒性困难。 其中，所述方法还包括以下步骤： S5：基于质量集中策略对调整后的SGFEM基本函数进行处理，达到减少质量矩阵数 和简化算法的目的。 其中，在所述步骤S5中，所述质量集中策略具体为： 采取典型抛物界面问题来说明质量集中策略： 满足界面条件： 参考界面问题(1)的离散化过程对抛物界面问题(9)进行处理，其中时间离散采用 10 CN 111611738 A 说　明　书 6/14 页 向后差分格式：即用节点O＝t0＜t1＜...＜tN＝T代替时间区间[0，T]，称 为 时间步长，令{tn：tn＝nτ，0≤n≤N}， 因此抛物界面问题(9)基于有限维空 间 的全离散格式如下： 向后Euler-Galerkin格式： 找到 使得 Crank-Nicolson  Galerkin格式： 找到 使得 其中 将(17)和(18)改写成线性 方程组，分别为 和 其中 为刚度矩阵， 为质量矩 阵， 为 的基底； 对形函数的阶进行重排列，使得SGFEM的质量矩阵M有如下形式： 其中M11与标准有限元FE部分有关，M22与SGFEM的富集部分有关，M12与标准FE和富 集部分的内积有关；M 211的矩阵规模为n ，M22的矩阵规模为Cn，n2Cn，均为常数；所谓质量集中 策略就是将M用一个对角矩阵 去近似代替，其中： 换言之， 是由将M中每行元素集中起来置于该行的对角线位置上；本方法采取只 集中M中有限元相关部分的质量矩阵，即只对M中的M11做集中，其它部分保持不变，集中后的 质量矩阵为： 11 CN 111611738 A 说　明　书 7/14 页 其中， 是一个对角矩阵，对角线元素为M11相应行的行和。 其中，在实际应用中，因为M11的矩阵规模为n2，M22的矩阵规模为Cn，C为常数，在离 散化规模很大时，有n2＞＞Cn，即富集距离函数的部分，也就是界面附近的富集节点个数与 有限元节点个数相比是少一维的。 与现有技术相比，本发明技术方案的有益效果是： 本发明提供的一种基于稳定广义有限元的界面问题模拟方法，针对界面问题提出 了一种SGFEM，既提高了传统GFEM的精度，还大幅降低刚度矩阵的条件数，使得界面问题的 数值模拟更加精确和鲁棒。本方法适用于任意相对位置的网格和界面。 附图说明 图1为界面问题的区域示意图； 图2为本发明所述方法流程示意图； 图3为不同种GFEM的富集节点格式示意图； 图4与界面相交的单元示例和每个三角形上的高斯积分规则示意图； 图5为不同界面系数比的能量相对误差EE，左：(ρ＝10)，右：(ρ＝100)； 图6为不同界面系数比的尺度条件数SCN，左：(ρ＝10)，右：(ρ＝100)； 图7为当δ＝1.88×10-2时不同种GFEM的富集方案示意图； 图8为刚度矩阵尺度条件数； 图9为质量集中策略：左：集中前，12849个非零点；右：集中后，4529个非零点。