一種平面問(wèn)題的加權(quán)Nitsche間斷伽遼金有限元法

鄧小蔚 張健飛 王明威

摘 要：以經(jīng)典間斷伽遼金有限元法求解彈性力學(xué)界面問(wèn)題，存在著由于穩(wěn)定系數(shù)取值不當(dāng)引起的數(shù)值不穩(wěn)定問(wèn)題，而加權(quán)Nitsche間斷伽遼金有限元法可以緩解這種問(wèn)題，但僅應(yīng)用于常量單元離散的情況。為解決上述問(wèn)題，基于加權(quán)Nitsche間斷伽遼金有限元法，針對(duì)平面彈性力學(xué)問(wèn)題，推導(dǎo)了四節(jié)點(diǎn)四邊形單元離散情況下的加權(quán)系數(shù)和穩(wěn)定參數(shù)的計(jì)算公式，建立了權(quán)重與穩(wěn)定參數(shù)間的定性依賴(lài)關(guān)系。通過(guò)建立和求解廣義特征值問(wèn)題，實(shí)現(xiàn)了加權(quán)系數(shù)和穩(wěn)定參數(shù)的自動(dòng)計(jì)算，使得高階單元的使用成為可能。通過(guò)數(shù)值試驗(yàn)檢驗(yàn)了方法的收斂性和穩(wěn)定性。結(jié)果表明：在求解均勻或材料分區(qū)不均勻介質(zhì)問(wèn)題時(shí)，加權(quán)Nitsche間斷伽遼金有限元法均表現(xiàn)出良好的穩(wěn)定性，且計(jì)算結(jié)果具有較高的精度。所提出的方法在一定程度上無(wú)須人工干預(yù)，具有高效率、高精度和良好的穩(wěn)定性，可以應(yīng)用于復(fù)雜界面問(wèn)題。

關(guān)鍵詞：彈性力學(xué);間斷Galerkin有限元法;加權(quán)Nitsche法;穩(wěn)定系數(shù);界面問(wèn)題;高階單元

中圖分類(lèi)號(hào)：O343.1?文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1008-1542（2018）06-0567-10

間斷伽遼金（discontinuous Galerkin， DG）有限元方法是有限元法（finite element method， FEM）[1]和有限體積法（finite volume method， FVM）[2]的推廣，它采用完全間斷的分離多項(xiàng)式空間對(duì)試函數(shù)和近似解進(jìn)行離散，具有更好的靈活性。對(duì)于不連續(xù)問(wèn)題[3-4]求解，DG有限元法弱遵循邊界條件，允許在單元界面處出現(xiàn)間斷的基函數(shù)，為處理不連續(xù)場(chǎng)提供了一種天然的計(jì)算框架。

DG法由于其具有局部保守性、穩(wěn)定性、高精度等優(yōu)勢(shì)，且易于處理復(fù)雜邊界、具有懸掛節(jié)點(diǎn)的不規(guī)則網(wǎng)格，已被用于處理多類(lèi)問(wèn)題，如反應(yīng)擴(kuò)散方程[5]、對(duì)流擴(kuò)散方程[6]、裂紋擴(kuò)展問(wèn)題[7-9]、復(fù)合材料[10]等。常見(jiàn)的DG法可按對(duì)稱(chēng)性分為兩類(lèi)，其中不對(duì)稱(chēng)DG法包括：LDG（local discontinuous Galerkin）[6]，常用于解決對(duì)流擴(kuò)散問(wèn)題，該類(lèi)方法中除主變量外，通量也是未知的;OBB（oden babuska baumann）[11]在單元上表現(xiàn)出局部質(zhì)量守恒的性質(zhì)，利于解決對(duì)流和擴(kuò)散問(wèn)題;NIPG（non-symmetric interior penalty Galerkin）[12]，與OBB法相比，可通過(guò)調(diào)整懲罰參數(shù)來(lái)提高準(zhǔn)確率;IIPG（interior penalty discontinuous Galerkin）[13]，在流體問(wèn)題中具有較好的穩(wěn)定性和準(zhǔn)確性;SIPG（symmetric interior penalty Galerkin）[14]是一種對(duì)稱(chēng)DG法，常用于解決不具有嚴(yán)格性穩(wěn)定要求的問(wèn)題。

經(jīng)典DG法在求解平面彈性界面問(wèn)題時(shí)，一般按算術(shù)平均計(jì)算界面上的通量，依靠數(shù)值經(jīng)驗(yàn)選擇穩(wěn)定參數(shù)，帶有較大的主觀性，對(duì)穩(wěn)定性有較大影響[15]，而且還影響計(jì)算精度。而精確的界面和內(nèi)部應(yīng)力計(jì)算結(jié)果是進(jìn)行開(kāi)裂等界面非線性行為分析的基礎(chǔ)。采用加權(quán)Nitsche法，通過(guò)加權(quán)平均計(jì)算界面上的通量，可以提高DG法的穩(wěn)定性。目前，關(guān)于加權(quán)DG有限元法，常應(yīng)變?nèi)切螁卧x散情況的研究較多，推導(dǎo)了詳細(xì)的加權(quán)系數(shù)和穩(wěn)定參數(shù)的計(jì)算公式，獲得了較為穩(wěn)定的結(jié)果[16-19]。筆者基于加權(quán)Nitsche間斷伽遼金有限元法，針對(duì)平面彈性力學(xué)中的界面問(wèn)題和四節(jié)點(diǎn)四邊形單元離散情況，推導(dǎo)了加權(quán)系數(shù)和穩(wěn)定參數(shù)的計(jì)算公式，編制了相應(yīng)的程序，并通過(guò)數(shù)值試驗(yàn)檢驗(yàn)了該方法的收斂性和穩(wěn)定性。

1?彈性力學(xué)界面問(wèn)題控制方程

3?空間離散

將Ω離散化為四節(jié)點(diǎn)四邊形實(shí)體單元和界面單元，如圖2所示?？蓪⑴c界面單元相鄰的實(shí)體單元分別稱(chēng)為左單元和右單元。局部坐標(biāo)系中，界面單元的節(jié)點(diǎn)編號(hào)、面號(hào)以及積分點(diǎn)編號(hào)如圖3所示。

4?數(shù)值分析

4.1?簡(jiǎn)支梁?jiǎn)栴}

4.1.1?收斂性測(cè)試

工況1

在界面Γ*上，各個(gè)應(yīng)力分量的相對(duì)誤差范數(shù)隨單元尺寸的變化如圖5所示。對(duì)于均勻介質(zhì)情況，以彈性力學(xué)解析解為參照標(biāo)準(zhǔn)。由圖5可知，應(yīng)力相對(duì)誤差范數(shù)隨單元尺寸呈基本線性的變化趨勢(shì)，當(dāng)單元尺寸足夠小時(shí)，計(jì)算結(jié)果趨向于解析解。

工況2

界面上各應(yīng)力分量的相對(duì)誤差范數(shù)隨單元尺寸的變化如圖6所示。對(duì)于雙材料分區(qū)不均勻介質(zhì)情況，以有限元軟件ANSYS數(shù)值解為參照標(biāo)準(zhǔn)。從圖6可以看出，隨單元尺寸減小，應(yīng)力相對(duì)誤差范數(shù)也相應(yīng)減小，且逐漸趨向于零?？梢?jiàn)，當(dāng)單元尺寸足夠小時(shí)，本文方法計(jì)算結(jié)果與ANSYS數(shù)值解趨于一致。

4.1.2?穩(wěn)定性測(cè)試

本文方法與經(jīng)典DG有限元法相比，主要不同在加權(quán)系數(shù)γ與穩(wěn)定系數(shù)α的選取上，其中本文方法根據(jù)式（26）、式（27）選取，而經(jīng)典DG有限元法采用算術(shù)平均計(jì)算界面通量，即γ1=γ2=0.5，依靠數(shù)值經(jīng)驗(yàn)選擇穩(wěn)定參數(shù)α，一般取α=δpG/|s|，其中G為剪切模量，|s|表示界面的長(zhǎng)度，δp為經(jīng)驗(yàn)參數(shù)。對(duì)于SIPG法，δp取值范圍是[20，+∞）[15]。

工況1

改變?chǔ)膒的值，觀察經(jīng)典DG有限元法求解結(jié)果的變化情況，并將其結(jié)果與本文方法結(jié)果比較，兩者的相對(duì)誤差范數(shù)如表1所示。結(jié)合表1可知，對(duì)于均勻介質(zhì)情況，大致在δp∈[20，1×108）內(nèi)，本文方法的精度與經(jīng)典SIPG法很接近;約在δp∈[103，108）時(shí)，經(jīng)典SIPG法的求解結(jié)果精度較高;約δp>1×108時(shí)，應(yīng)力結(jié)果出現(xiàn)不穩(wěn)定現(xiàn)象。

以彈性解析解為對(duì)照，對(duì)比δp=1×106，δp=1×1012時(shí)界面的應(yīng)力值，以及本文方法結(jié)果如圖7所示。可以看出，δp=1×1012時(shí)的界面應(yīng)力出現(xiàn)明顯震蕩現(xiàn)象。

工況2

與均質(zhì)問(wèn)題類(lèi)似，選取幾組不同的δp值，求出經(jīng)典DG有限元法與本文方法結(jié)果的誤差范數(shù)，如表2所示。結(jié)合表2得出：雙材料分區(qū)不均勻介質(zhì)情況，δp在[500，6×104）內(nèi)時(shí)，經(jīng)典SIPG法的求解結(jié)果精度較高;δp>6×104時(shí)，應(yīng)力結(jié)果出現(xiàn)不穩(wěn)定現(xiàn)象。

δp=1×104，δp=1×108時(shí)界面的應(yīng)力值，本文方法結(jié)果以及ANSYS數(shù)值解答如圖8所示?？梢钥闯?，本文方法結(jié)果與δp=1×104所得結(jié)果精度相當(dāng)，δp=1×108時(shí)，應(yīng)力結(jié)果出現(xiàn)明顯震蕩。

理論上隨著δp增大，穩(wěn)定系數(shù)α相應(yīng)增大，經(jīng)典SIPG法的精度應(yīng)該隨之提高。但從上述結(jié)果可知，當(dāng)δp取值過(guò)小時(shí)，計(jì)算精度不夠;當(dāng)δp取值過(guò)大時(shí)，計(jì)算結(jié)果會(huì)出現(xiàn)震蕩現(xiàn)象;而本文方法無(wú)需經(jīng)驗(yàn)確定穩(wěn)定系數(shù)，避免了數(shù)值不穩(wěn)定性，同時(shí)也具有較高的精度。

4.2?多材料分區(qū)不均勻界面問(wèn)題

如圖9所示正方形區(qū)域，邊長(zhǎng)為5 m，不計(jì)自重。多個(gè)界面將整個(gè)區(qū)域Ω劃分為不同的子區(qū)域Ωm，m=1，2，…，14。荷載分為2種情況，即頂面均勻受壓和頂面均勻受剪，集度q=10 N/m2。

子區(qū)域均為不同材料，Ωm中彈性模量如表3所示，泊松比ν均為0.3。

將本文方法所得應(yīng)力結(jié)果與ANSYS數(shù)值解答對(duì)比，如圖10和圖11所示。

通過(guò)以上的應(yīng)力云圖比較可知：本文方法在求解多材料分區(qū)不均勻界面問(wèn)題時(shí)，不僅具有與連續(xù)伽遼金有限元法相同的精度，而且各個(gè)界面上的穩(wěn)定參數(shù)和加權(quán)系數(shù)都是自動(dòng)生成，無(wú)須人工干預(yù)，提高了算法的效率和穩(wěn)定性，從而可以應(yīng)用于復(fù)雜界面問(wèn)題。

5?結(jié)?論

基于加權(quán)Nitsche間斷伽遼金有限元法，針對(duì)平面彈性力學(xué)界面問(wèn)題，推導(dǎo)了四節(jié)點(diǎn)四邊形單元離散下的加權(quán)系數(shù)和穩(wěn)定參數(shù)計(jì)算公式，這些公式同樣也可以推廣到其他高階單元的情況，使得高階單元的使用成為可能。數(shù)值試驗(yàn)結(jié)果表明，該方法在求解均勻或材料分區(qū)不均勻介質(zhì)問(wèn)題時(shí)具有良好的收斂性和穩(wěn)定性，計(jì)算精度與連續(xù)間斷伽遼金有限元法相同。與經(jīng)典間斷伽遼金有限元法相比，該法不需要人為設(shè)定穩(wěn)定參數(shù)，避免了穩(wěn)定參數(shù)取值過(guò)小引起的精度問(wèn)題和取值過(guò)大引起的數(shù)值不穩(wěn)定問(wèn)題，可獲得更為精確的界面和內(nèi)部應(yīng)力計(jì)算結(jié)果，從而為進(jìn)一步分析開(kāi)裂和滑動(dòng)等界面非線性行為提供可靠的基礎(chǔ)。

目前筆者只研究了平面彈性力學(xué)界面問(wèn)題的穩(wěn)定加權(quán)Nitsche間斷伽遼金有限元法，未來(lái)需進(jìn)一步將其推廣應(yīng)用于解決空間問(wèn)題和復(fù)雜非線性界面問(wèn)題，并利用ANSYS，ABAQUS等商用軟件進(jìn)行二次開(kāi)發(fā)，拓展相關(guān)功能，逐步實(shí)現(xiàn)大規(guī)模工程應(yīng)用。

參考文獻(xiàn)/References：

[1]?王勖成. 有限單元法基本原理和數(shù)值方法[M]. 北京： 清華大學(xué)出版社， 1997.

[2]?李人憲. 有限體積法基礎(chǔ)[M]. 北京：國(guó)防工業(yè)出版社， 2008.

[3]?王永婧， 張冬雯， 于健騏. 基于T-S模型的一類(lèi)時(shí)滯非線性系統(tǒng)模型預(yù)測(cè)控制[J]. 河北工業(yè)科技， 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]?歐陽(yáng)志， 鞠興華. 基于非線性破壞準(zhǔn)則的地基臨界荷載計(jì)算[J]. 河北工業(yè)科技， 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]?張榮培. 一類(lèi)非線性擴(kuò)散方程的間斷有限元方法研究[D]. 北京： 中國(guó)工程物理研究院， 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]?琚宏昌， 陳國(guó)榮. 混凝土梁裂紋擴(kuò)展的內(nèi)聚區(qū)模型數(shù)值模擬[J]. 鄭州大學(xué)學(xué)報(bào)（工學(xué)版）， 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]王明威， 張健飛， 鄧小蔚. 彈性力學(xué)問(wèn)題的加權(quán)Nitsche間斷伽遼金有限元法[J]. 河南科技大學(xué)學(xué)報(bào)（自然科學(xué)版）， 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.