一种平面问题的加权Nitsche间断伽辽金有限元法

邓小蔚 张健飞 王明威

摘 要：以经典间断伽辽金有限元法求解弹性力学界面问题，存在着由于稳定系数取值不当引起的数值不稳定问题，而加权Nitsche间断伽辽金有限元法可以缓解这种问题，但仅应用于常量单元离散的情况。为解决上述问题，基于加权Nitsche间断伽辽金有限元法，针对平面弹性力学问题，推导了四节点四边形单元离散情况下的加权系数和稳定参数的计算公式，建立了权重与稳定参数间的定性依赖关系。通过建立和求解广义特征值问题，实现了加权系数和稳定参数的自动计算，使得高阶单元的使用成为可能。通过数值试验检验了方法的收敛性和稳定性。结果表明：在求解均匀或材料分区不均匀介质问题时，加权Nitsche间断伽辽金有限元法均表现出良好的稳定性，且计算结果具有较高的精度。所提出的方法在一定程度上无须人工干预，具有高效率、高精度和良好的稳定性，可以应用于复杂界面问题。

关键词：弹性力学;间断Galerkin有限元法;加权Nitsche法;稳定系数;界面问题;高阶单元

中图分类号：O343.1 文献标志码：A

文章编号：1008-1542（2018）06-0567-10

间断伽辽金（discontinuous Galerkin， DG）有限元方法是有限元法（finite element method， FEM）[1]和有限体积法（finite volume method， FVM）[2]的推广，它采用完全间断的分离多项式空间对试函数和近似解进行离散，具有更好的灵活性。对于不连续问题[3-4]求解，DG有限元法弱遵循边界条件，允许在单元界面处出现间断的基函数，为处理不连续场提供了一种天然的计算框架。

DG法由于其具有局部保守性、稳定性、高精度等优势，且易于处理复杂边界、具有悬挂节点的不规则网格，已被用于处理多类问题，如反应扩散方程[5]、对流扩散方程[6]、裂纹扩展问题[7-9]、复合材料[10]等。常见的DG法可按对称性分为两类，其中不对称DG法包括：LDG（local discontinuous Galerkin）[6]，常用于解决对流扩散问题，该类方法中除主变量外，通量也是未知的;OBB（oden babuska baumann）[11]在单元上表现出局部质量守恒的性质，利于解决对流和扩散问题;NIPG（non-symmetric interior penalty Galerkin）[12]，与OBB法相比，可通过调整惩罚参数来提高准确率;IIPG（interior penalty discontinuous Galerkin）[13]，在流体问题中具有较好的稳定性和准确性;SIPG（symmetric interior penalty Galerkin）[14]是一種对称DG法，常用于解决不具有严格性稳定要求的问题。

经典DG法在求解平面弹性界面问题时，一般按算术平均计算界面上的通量，依靠数值经验选择稳定参数，带有较大的主观性，对稳定性有较大影响[15]，而且还影响计算精度。而精确的界面和内部应力计算结果是进行开裂等界面非线性行为分析的基础。采用加权Nitsche法，通过加权平均计算界面上的通量，可以提高DG法的稳定性。目前，关于加权DG有限元法，常应变三角形单元离散情况的研究较多，推导了详细的加权系数和稳定参数的计算公式，获得了较为稳定的结果[16-19]。笔者基于加权Nitsche间断伽辽金有限元法，针对平面弹性力学中的界面问题和四节点四边形单元离散情况，推导了加权系数和稳定参数的计算公式，编制了相应的程序，并通过数值试验检验了该方法的收敛性和稳定性。

1 弹性力学界面问题控制方程

3 空间离散

将Ω离散化为四节点四边形实体单元和界面单元，如图2所示。可将与界面单元相邻的实体单元分别称为左单元和右单元。局部坐标系中，界面单元的节点编号、面号以及积分点编号如图3所示。

4 数值分析

4.1 简支梁问题

4.1.1 收敛性测试

工况1

在界面Γ*上，各个应力分量的相对误差范数随单元尺寸的变化如图5所示。对于均匀介质情况，以弹性力学解析解为参照标准。由图5可知，应力相对误差范数随单元尺寸呈基本线性的变化趋势，当单元尺寸足够小时，计算结果趋向于解析解。

工况2

界面上各应力分量的相对误差范数随单元尺寸的变化如图6所示。对于双材料分区不均匀介质情况，以有限元软件ANSYS数值解为参照标准。从图6可以看出，随单元尺寸减小，应力相对误差范数也相应减小，且逐渐趋向于零。可见，当单元尺寸足够小时，本文方法计算结果与ANSYS数值解趋于一致。

4.1.2 稳定性测试

本文方法与经典DG有限元法相比，主要不同在加权系数γ与稳定系数α的选取上，其中本文方法根据式（26）、式（27）选取，而经典DG有限元法采用算术平均计算界面通量，即γ1=γ2=0.5，依靠数值经验选择稳定参数α，一般取α=δpG/|s|，其中G为剪切模量，|s|表示界面的长度，δp为经验参数。对于SIPG法，δp取值范围是[20，+∞）[15]。

工况1

改变δp的值，观察经典DG有限元法求解结果的变化情况，并将其结果与本文方法结果比较，两者的相对误差范数如表1所示。结合表1可知，对于均匀介质情况，大致在δp∈[20，1×108）内，本文方法的精度与经典SIPG法很接近;约在δp∈[103，108）时，经典SIPG法的求解结果精度较高;约δp>1×108时，应力结果出现不稳定现象。

以弹性解析解为对照，对比δp=1×106，δp=1×1012时界面的应力值，以及本文方法结果如图7所示。可以看出，δp=1×1012时的界面应力出现明显震荡现象。

工况2

与均质问题类似，选取几组不同的δp值，求出经典DG有限元法与本文方法结果的误差范数，如表2所示。结合表2得出：双材料分区不均匀介质情况，δp在[500，6×104）内时，经典SIPG法的求解结果精度较高;δp>6×104时，应力结果出现不稳定现象。

δp=1×104，δp=1×108时界面的应力值，本文方法结果以及ANSYS数值解答如图8所示。可以看出，本文方法结果与δp=1×104所得结果精度相当，δp=1×108时，应力结果出现明显震荡。

理论上随着δp增大，稳定系数α相应增大，经典SIPG法的精度应该随之提高。但从上述结果可知，当δp取值过小时，计算精度不够;当δp取值过大时，计算结果会出现震荡现象;而本文方法无需经验确定稳定系数，避免了数值不稳定性，同时也具有较高的精度。

4.2 多材料分区不均匀界面问题

如图9所示正方形区域，边长为5 m，不计自重。多个界面将整个区域Ω划分为不同的子区域Ωm，m=1，2，…，14。荷载分为2种情况，即顶面均匀受压和顶面均匀受剪，集度q=10 N/m2。

子区域均为不同材料，Ωm中弹性模量如表3所示，泊松比ν均为0.3。

将本文方法所得应力结果与ANSYS数值解答对比，如图10和图11所示。

通过以上的应力云图比较可知：本文方法在求解多材料分区不均匀界面问题时，不仅具有与连续伽辽金有限元法相同的精度，而且各个界面上的稳定参数和加权系数都是自动生成，无须人工干预，提高了算法的效率和稳定性，从而可以应用于复杂界面问题。

5 结 论

基于加权Nitsche间断伽辽金有限元法，针对平面弹性力学界面问题，推导了四节点四边形单元离散下的加权系数和稳定参数计算公式，这些公式同样也可以推广到其他高阶单元的情况，使得高阶单元的使用成为可能。数值试验结果表明，该方法在求解均匀或材料分区不均匀介质问题时具有良好的收敛性和稳定性，计算精度与连续间断伽辽金有限元法相同。与经典间断伽辽金有限元法相比，该法不需要人为设定稳定参数，避免了稳定参数取值过小引起的精度问题和取值过大引起的数值不稳定问题，可获得更为精确的界面和内部应力计算结果，从而为进一步分析开裂和滑动等界面非线性行为提供可靠的基础。

目前笔者只研究了平面弹性力学界面问题的稳定加权Nitsche间断伽辽金有限元法，未来需进一步将其推广应用于解决空间问题和复杂非线性界面问题，并利用ANSYS，ABAQUS等商用软件进行二次开发，拓展相关功能，逐步实现大规模工程应用。

参考文献/References：

[1] 王勖成. 有限单元法基本原理和数值方法[M]. 北京： 清华大学出版社， 1997.

[2] 李人宪. 有限体积法基础[M]. 北京：国防工业出版社， 2008.

[3] 王永婧， 张冬雯， 于健骐. 基于T-S模型的一类时滞非线性系统模型预测控制[J]. 河北工业科技， 2018， 35（1）：37-42.

WANG Yongjing， ZHANG Dongwen， YU Jianqi. Model predictive control for a class of nonlinear systems with time-delay based on T-S model[J]. Hebei Journal of Industrial Science and Technology， 2018，35（1）：37-42.

[4] 欧阳志， 鞠兴华. 基于非线性破坏准则的地基临界荷载计算[J]. 河北工业科技， 2017， 34（4）：259-264.

OUYANG Zhi， JU Xinghua. Critical load calculation of foundation based on nonlinear failure criterion[J]. Hebei Journal of Industrial Science and Technology， 2017， 34（4）：259-264.

[5] 张荣培. 一类非线性扩散方程的间断有限元方法研究[D]. 北京： 中国工程物理研究院， 2012.

[6] COCKBURN B， DAWSON C. Some extensions of the local discontinuous galerkin method for convection-diffusion equations in multidimensions[J]. Mathematics of Finite Elements and Applications， 1999：225-238.

[7] MERGHEIM J， KUHL E， STEINMANN P. A hybrid discontinuous Galerkin/interface method for the computational modelling of failure[J]. International Journal for Numerical Methods in Biomedical Engineering， 2004， 20（7）：511-519.

[8] 琚宏昌， 陳国荣. 混凝土梁裂纹扩展的内聚区模型数值模拟[J]. 郑州大学学报（工学版）， 2008， 29（3）：10-14.

QU Hongchang， CHEN Guorong. Numerical simulation of crack propagating using cohesive zone model[J]. Journal of Zhengzhou University（Engineering Science）， 2008， 29（3）：10-14.

[9] NGUYEN V P. Discontinuous Galerkin/extrinsic cohesive zone modeling： Implementation caveats and applications in computational fracture mechanics[J]. Engineering Fracture Mechanics， 2014， 128：37-68.

[10]van der MEER F P， SLUYS L J， HALLETT S R， et al. Computational modeling of complex failure mechanisms in laminates[J]. Journal of Composite Materials， 2011， 46（5）：603-623.

[11]ODEN J T， BABUSKA I， BAUMANN C E. A discontinuous hp finite element method for diffusion problems[J]. Journal of Computational Physics， 1998， 146（2）：491-519.

[12]RIVIERE B， SHAW S， WHEELER M F， et al. Discontinuous Galerkin finite element methods for linear elasticity and quasistatic linear viscoelasticity[J]. Journal Numerische Mathematik， 2003， 95（2）：347-376.

[13]DAWSON C， SUN S， WHEELER M F. Compatible algorithms for coupled flow and transport[J]. Computer Methods in Applied Mechanics and Engineering， 2004， 193：2565-2580.

[14]ARNOLD D N. An interior penalty finite element method with discontinuous elements[J]. Siam Journal on Numerical Analysis， 1982， 19（4）：742-760.

[15]LIU R， WHEELER M F， DAWSON C N. A three-dimensional nodal-based implementation of a family of discontinuous Galerkin methods for elasticity problems[J]. Computers and Structures， 2009， 87（3/4）：141-150.

[16]HARARI  I， DOLBOW  J. Analysis of an efficient finite element method for embedded interface problems[J]. Computational Mechanics， 2010， 46（1）： 205-211.

[17]ANNAVARAPU C， HAUTEFEUILLE M， DOLBOW J E. Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods[J]. International Journal for Numerical Methods in Engineering， 2012， 92（2）：206-228.

[18]ANNAVARAPU C， HAUTEFEUILLE M， DOLBOW J E. A robust Nitsches formulation for interface problems[J]. Computer Methods in Applied Mechanics and Engineering， 2012， 225/228：44-54.

[19]王明威， 張健飞， 邓小蔚. 弹性力学问题的加权Nitsche间断伽辽金有限元法[J]. 河南科技大学学报（自然科学版）， 2018，39（3）：83-88.

WANG Mingwei， ZHANG Jianfei， DENG Xiaowei.  Weighted nitsche discontinuous Galerkin finite element method for elasticity problems[J]. Journal of  Henan University of Science and Technology（Natural Sciences）， 2018，39（3）：83-88.

[20]JIANG W， ANNAVARAPU C， DOLBOW J E， et al. A robust Nitsches formulation for interface problems with spline-based finite elements[J]. International Journal for Numerical Methods in Engineering， 2015， 104（7）：676-696.

[21]EMBAR A， DOLBOW J， HARARI I. Imposing Dirichlet boundary conditions with Nitsches method and spline-based finite elements[J]. International Journal for Numerical Methods in Engineering， 2010， 83（7）：877-898.