一種平面問題的加權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.