基于三場變分原理的對偶mortar 有限元法

周墨臻，張丙印，張頂立，方黃城

(1. 北京交通大學城市地下工程教育部重點實驗室，北京 100044；2. 清華大學水沙科學與水利水電工程國家重點實驗室，北京 100084)

對計算域剖分整體協調網格是傳統有限元方法對前處理建模的基本要求。然而，對復雜工程問題的大規模計算，這一要求可能成為巨大挑戰[1]。在工程應用中，使用非協調網格的需求十分普遍，例如，疏密網格過渡或自適應網格加密[2―4]、多尺度模擬[5―6]、非線性接觸[7―9]、基于區域分解的并行計算[10―11]、邊界元或有限體積法等與有限元的耦合[12―14]、多物理場耦合[15―16]等。目前，面向有限元的非協調網格處理技術主要有點-面法[7,17]、面-面法[14,16,18―19]、Nistche 法[20]、特殊界面連接單 元[3,6,10,21]、三場變分法[2,22―23]、四場變分法[13,15]、多面體單元[24―25]、比例邊界有限元[26]、局部無網格法[1,5]等。

在這些方法中，mortar 元[18―19]作為最具代表性的一種面-面方法，由于求解精度高且具有較好靈活性而得到廣泛應用，并成功嵌入 ABAQUS 和ANSYS 等商業軟件[27―29]。然而，隨著工程應用對有限元計算在復雜化和精細化等方面的要求不斷提高，mortar 元表現出了以下不足：

1) 約束交叉問題。如圖1 所示，當兩個以上子區域的界面在同一處發生重疊，二維問題會出現交叉點，三維問題則出現交叉線，從而導致界面約束條件之間發生交叉，需進行特殊處理。

圖1 交叉點或交叉線問題 Fig.1 Cross-point or cross-line problems

2) 主從偏見問題。Mortar 元將非協調網格兩側的界面分為主面和從面，其主要目的是將從面用于離散Lagrange 乘子(以下均簡稱為乘子)。這種主從關系是由數值算法所引入的非物理概念，需要人為指定，且主要依賴于計算經驗。尤其對圖1 所示的約束交叉情形，主從關系的選擇更為困難。

3) 求解效率問題。Mortar 元一般采用罰函數法、乘子法或增廣Lagrange 法施加界面約束，可引起矩陣病態、矩陣非正定或額外迭代層等問題。這些問題在中小規模計算時并不突出，但對大規模計算則會嚴重影響整體求解效率。

上述問題雖然均可在計算力學領域內找到各自的解決途徑，但若要同時解決，且還保留mortar元的固有優勢，則仍需深入研究。例如，Nistche法[20]可解決問題1)和2)，但通常需要增加穩定項以避免矩陣病態，而對非線性材料，還需推導本構關系的變分，處理困難；再如，使用分片常量函數插值乘子[14,30]可解決問題1)和3)，但代價是乘子不連續，精度較低；此外，采用對偶基函數插值乘子，也就是對偶mortar 元[16,31]，可解決問題3)。由于該特點，對偶mortar 元得到了較快發展，也是目前等幾何分析中的熱點問題之一[32―38]。但對偶mortar元對問題1)也需要特殊處理，一般是在約束交叉處縮減乘子[32―39]，但這一處理需檢測約束交叉面片，復雜工程應用存在不便。

本文基于mortar 元的理論框架，采用三場變分原理，提出了引入媒介面的對偶mortar 有限元法。該方法可同時解決約束交叉、主從偏見和求解效率問題，無需縮減乘子，對材料非線性問題也不引入額外困難，為非協調網格處理提供了一種新途徑。

1 三場變分法

1.1 數學描述

圖2 以三個子區域Ωα為例(上標α用以代表不同的子區域)，示意了本文所涉及的各類邊界條件。Dirichlet 和Neumann 邊界分別記為和。假定在每個Ωα內各自剖分協調網格，而不同子區域之間的網格相互獨立，由此導致的網格非協調表面記為Γα，并簡稱為實體面。

圖2 子區域界面及媒介面示意 Fig.2 Schematic of subdomain interfaces and the medium face frame

區別于mortar 元，本文引入一個虛擬的媒介面Γm，其網格剖分完全獨立于各個子區域的實體網格。由于各子區域之間以及子區域與媒介面之間均不滿足有限元計算的連續性要求，因此，需要構造連續性條件。為敘述簡便，以下均以基于位移求解的應力變形有限元計算為例。

將實體區域和媒介面內任意物質點的位形分別記為Xa和Xm(如圖2 示)，位移分別記為αu和um。界面間的位移連續性條件可描述為：

式中：usα為實體面Γα上的位移。Γ mα為實體面和媒介面之間經過投影之后所產生的重疊區域，具體投影方法將在2.2 節詳細介紹。

在mortar 元中，界面連續性條件定義為從面與主面之間的相對位移為零，顯然從面和主面均對應于本文的實體面，因而導致主從偏見問題。本文的連續性條件定義為實體面和媒介面之間的相對位移為零，不同的實體面之間無需定義連續性條件，從而可直接解決約束交叉和主從偏見問題。

1.2 弱形式

采用Lagrange 乘子法施加式(1)所示的耦合約束條件，則約束變分原理給出：

式中：δΠ為泛函變分；為各個子區域的虛功，由常規有限元理論確定，此處略去其具體表達；和分別為界面虛功和界面連續性條件；αλ表示乘子，物理含義是界面力。可見，上述變分原理共涉及三個待求解的未知場變量：αu、um和αλ。因此，式(2)也稱為三場變分法[22―23]。

若將位移進一步分解為剛體位移和變形，則三場變分法將演變為四場變分法[15]。四場變分法的主要特點是利用剛體平衡等條件，布置媒介面的網格結點。本文采用的三場變分法只需滿足式(2c)所示的連續性條件，而無需其他額外方程，且媒介面的網格剖分完全獨立于實體區域。

mortar 元由于不涉及媒介面的位移場mu，因此，對應二場變分原理。本文以下將mortar 元的理論框架由二場變分原理擴展到用于三場變分原理，并通過采用對偶mortar 元解決求解效率方面的 問題。

2 三場對偶mortar 元

2.1 空間離散

各子區域的位移αu采用線性實體等參單元，媒介面的位移um采用線性表面等參單元：

式中：Nα和Nm分別為子區域和媒介面的單元形函數矩陣；dα和dm為相應的離散位移。

乘子λα離散在實體面Γα上：

式中：zα表示離散后的乘子；ψα為乘子的插值函數矩陣。mortar 元一般直接取ψα=Nα，而在對偶mortar 元中，ψ α則采用滿足如下雙正交特性的對偶基函數[28]：

式中，δjk為Kronecker 函數。為確定ψα，可將積分域Γ m,α分解到每個實體面的表面單元上：

式中，e mα表示媒介面經投影后與實體表面單元之間的重疊區域，其計算方法將在2.2 節具體介紹。

假定ψα為Nα的線性組合：

式中：jka為待定系數；Ae為系數矩陣。將式(7)代入式(6)可得：

由式(8)可求解得到待定系數jka，從而確定對偶基函數αψ的具體形式。值得指出的是，上式中矩陣eM的大小僅取決于單個表面單元的結點個數，因而其求逆計算的代價很小。

將式(3)和式(4)分別代入式(2b)和式(2c)，可分別得到界面虛功和連續性條件的離散形式：

式中，dsα表示usα的空間離散形式。矩陣Dα和Mα分別為：

式中，3I表示3×3 的單位矩陣，因為三維問題中位移共有3 個自由度。由式(10a)可見，矩陣αD為嚴格對角矩陣，這是由式(5)所確定的。這一特性對提高求解效率具有關鍵作用，2.4 節將具體介紹。

2.2 界面投影及數值積分

在早期mortar 元中，界面積分域e mα直接取為整個從面單元。研究表明，該種做法可能導致較大數值誤差。本文將mortar 元的投影與積分改進研究[18]推廣用于三場對偶mortar 元。主要步驟如下：

1) 如圖3(a)示，將實體表面單元eα和媒介面單元em均按eα中心點處的法向進行投影，從而將三維面單元降維成二維平面單元。

2) 如圖3(b)示，對兩個平面單元使用多邊形裁剪算法，得到二者之間的多邊形重疊區域e mα。

3) 如圖3(c)示，對e mα進一步作三角剖分。圖3(c)示意了剖分得到的7 個三角形子片段et。

4) 如圖3(d)示，將相關的積分項在三角形子片段et上作積分，本文采用7 點Hammer 積分。

圖3 投影方法和積分方案 Fig.3 Projection method and numerical integration

由此，式(8b)～式(8c)和式(10a)～式(10b)可分別計算如下：

式中：Hi為第i個積分點的積分權系數；表示Jacobian 行列式。需注意的是，上式中的插值函數、、和定義在eα或me上，而積分點則定義在三角形子片段上。因此，需要根據積分點的笛卡爾坐標，在eα和me上作等參逆變換，在得到等參坐標之后，再計算插值函數。

2.3 線性方程組

由式(2)和式(9)，可寫出總體線性方程組如下：

式中：上標s為實體面Γα上的自由度子集；上標R則表示實體區域中除實體面Γα之外的自由度子集。K、Fext、ds、dR、D、M和z表示由各子區域Ωα的相應子塊所構成的矩陣或列向量：

式中：n為子區域的數目；Kα和分別為第α個子區域的總體剛度矩陣和外力，由確定，屬于有限元的常規計算內容，此處略去具體表達。D α和Mα由式(11c)～式(11d)計算后集成得到。

2.4 乘子凝聚

在式(12)中，矩陣第3、4 行的對角元素均為零，使得矩陣失去正定性，導致鞍點問題。其次，乘子自由度z為待求解的未知量，增加了矩陣規模。為此，基于通用矩陣變換[40]，對乘子進行凝聚：

式中，ds表示實體面與媒介面的相對位移，可通過矩陣C表示如下：

將式(15)代入式(12)的第4 行可得：

式(16)將耦合約束變換成為完全解耦形式。也即，對式(14)而言，式(16)等價為Dirichlet 邊界條件，可采用直接消去法施加。這表明：雖然媒介面引入了額外自由度，但實體面的自由度可全部消去，矩陣規模可基本不受影響。同時，經過變換之后，式(14)不再包含乘子自由度z，且矩陣恢復了正定性。至此，只需求解滿足式(16)的式(14)，而無需再求解式(12)，從而解決了求解效率方面的問題。

需要指出的是，若使用常規mortar 元，則矩陣D不是對角陣，而式(15)需求解1-D，計算代價很高。三場對偶mortar 元由于使用式(5)所示的雙正交條件，使得D為嚴格對角矩陣，其求逆代價可忽略不計，式(15)可很方便地計算得到。

3 數值算例

基于所提出的三場對偶mortar 有限元，采用C++語言自主編制了相應的三維計算程序。以下采用兩個數值算例，對該計算程序進行驗證。兩個算例均采用楊氏模量E= 1×105Pa、泊松比ν= 0 的線彈性材料。此處，選取泊松比ν= 0 的主要目的是簡化理論解，雖與實際工程計算存在差別，但屬于檢驗非協調網格處理效果的一般做法[26,41]。

3.1 接觸分片試驗

接觸分片試驗常用于檢驗接觸數值算法處理非協調網格的能力，一般僅考慮兩個子區域。為檢驗處理約束交叉問題的能力，此處將2 m×2 m× 2 m 的計算區域考慮為5 個子區域。實體網格如圖4(a)所示，其中，1Ω～4Ω剖分各向等長的六面體單元，單元長度分別為1/4 m、1/5 m、1/6 m 和1/7 m，Ω5的單元尺寸為2/7 m×2/7 m×1/3 m。媒介面網格如圖4(b)所示，單元為各向等長，長度為1/3 m。

圖4 接觸分片試驗計算網格 Fig.4 Computational mesh of the contact patch test

對圖4(b)所示的約束交叉線，mortar 元需檢測相應面片并縮減乘子，本文方法則無需任何特殊處 理，也無需檢測。同時，計算時完全無需對Ω1～Ω5的實體面區分主從關系，很好地避免了mortar 元的主從偏見問題。計算的邊界條件為：Ω1～Ω4的上表面施加大小為1 Pa 的法向均布壓力，Ω5的下表面所有結點均取為固支約束。顯然，該問題的理論解是單向應力狀態，且各處豎向應力均為-1 Pa。

計算結果如圖5 所示。在圖5(a)中，黑色線條表示位移等值線。可見，計算所得最大豎向位移為-2×10-5m，與理論解吻合，且圖中未見明顯位移不連續現象。在圖5(b)中，黑色線條表示網格線。可見，計算所得豎向應力為均勻分布，圖示的數值在-1±10-6Pa 之間。為反映計算準確性，定義位移、應力及能量的相對誤差如下：

式中，ue、σe和εe分別為位移、應力及應變的理論解。

圖5 接觸分片試驗的數值計算結果 Fig.5 Numerical results of the contact patch test

統計得到的位移、應力及能量誤差分別為1.83× 10-8、4.17×10-8和4.63×10-8。考慮到機器精度，可認為計算結果與理論值吻合。若直接使用二場mortar 元，而不在交叉點處縮減乘子，則上述誤差分別為1.06×10-4、1.34×10-3和1.44×10-3。可見，本文所提出的三場對偶mortar 元無需處理交叉點，也可對約束交叉問題具有很高的求解精度。

為檢驗求解效率，采用預處理共軛梯度法PCG(preconditioned conjugate gradients)求解線性方程組，分別對三場對偶mortar 元和罰函數法進行測試。表1 統計了PCG 方法在使用不同預處理時的迭代次數。可見，在三種預處理情況下，三場對偶mortar 元的迭代次數均遠少于罰函數法，驗證了該方法在求解效率方面的優勢。

表1 線性方程組求解迭代次數 Table 1 Iterations of solving the linear system

需指出的是，對Lagrange 乘子法所導致的鞍點問題，PCG 方法無法求解，需使用穩定雙共軛梯度法 BICG-STAB(Bi-conjugate gradients stabilized method)等特殊方法，且求解效率較低。而增廣Lagrange 法則由于引入額外迭代層，即便對線性問題也可能需要多次求解線性方程組，整體效率受較大影響。對這兩種情況，本文不再一一驗算。

3.2 帶圓孔平板

為進一步驗證對約束交叉問題的處理能力，考慮如圖6(a)所示的彈性力學經典平面應力問題：帶圓孔無限平板的單向均勻拉伸。該問題的位移及應力解析解可參見文獻[20]。

鑒于對稱性，可取圖6(a)的1/4 進行計算。而對無限域，則截取圓孔周圍一定范圍內的平板進行分析，如圖6(b)所示。為使得應力變形狀態與圖6(a)一致，對圖6(b)所示模型，施加x向和y向對稱約束，并根據應力的解析解，計算左側及上部邊界的面力，將其作為面力邊界施加在相應表面[20]。

根據圖6(b)所示意的10 個子區域，剖分了如圖7 所示的疏密相間三維計算網格，以檢驗數值算法處理復雜子區域劃分的能力。計算區域平面尺寸為4 m×4 m，厚度為1 m，無限遠處的拉力p= 1 Pa，圓孔半徑取為1 m。在圖6(b)和圖7(a)中，分別對笛卡爾坐標系和柱坐標系進行了示意。

圖6 帶圓孔無限大平板示意 Fig.6 Schematic of the infinite plate with circular hole

圖7 帶圓孔平板計算網格 Fig.7 Computational mesh of the plate with circular hole

計算結果如圖8 所示。在圖8(a)～圖8(b)中，黑色線條表示位移等值線。可見，計算所得x向和y向位移分布規律與理論解一致，且未出現明顯位移不連續現象。在圖8(c)中，黑色線條表示網格線。可見，計算所得應力σxx在圓孔頂部出現明顯應力集中現象，最大應力值為2.987 Pa，與解析解所給出的3 Pa 十分接近。為進一步對比分析，本文另補充了針對圖9 所示協調網格的計算。

圖10 給出了計算所得應力σxx沿圓孔的分布及其與解析解的對比情況，柱坐標θ的定義如圖7(a)所示。由圖可見，非協調網格得到了與協調網格幾乎一致的計算結果，且均與解析解十分接近，驗證了三場對偶mortar 元處理非協調網格的有效性。

圖8 帶圓孔平板的數值計算結果 Fig.8 Numerical results of the plate with circular hole

圖9 帶圓孔平板的協調有限元網格 Fig.9 Conforming mesh of the plate with circular hole

圖10 數值計算結果與解析解的比較 Fig.10 Comparison of the numerical and analytical results

4 關于媒介面網格的討論

在上述的兩個算例中，媒介面的網格均具有兩個特點：1) 自身協調；2) 密度低于實體網格。

在實際工程應用中，滿足第一個特點并不會引入很大困難。這是因為媒介面相比于實體區域而言，只需剖分降維后的網格。相當于以剖分一個協調的面網格的代價，實現了不同子區域的網格各自獨立剖分，且子區域的劃分形式不限。

第二個特點，更準確地說，是要求媒介面的網格不能比實體網格中較密的那一側更密。若不滿足這一要求，則會導致式(14)中第三行存在線性相關性，也即矩陣病態。實際應用通常可方便地滿足這一要求，更多討論細節可參見文獻[22]。考慮到工程問題的復雜性，允許媒介面網格任意剖分而無需格外顧及網格密度的限制，同樣具有研究意義。關于該問題的解決，限于篇幅將另文討論。

5 結論

本文將mortar元從二場變分原理推廣為引入獨立媒介面的三場變分原理，并采用對偶基函數實現了乘子的凝聚，提出了三場對偶mortar 有限元。相比于已有的數值算法，該方法具有如下優勢：

(1) 由于媒介面的引入，自動解決了mortar 元的約束交叉問題和主從偏見問題，無需縮減乘子等其他特殊處理。

(2) 采用對偶基函數插值乘子，求解效率高于罰函數法、Lagrange 乘子法和增廣Lagrange 法。

(3) 相比Nistche 方法，無需增加額外的穩定項，也無需推導本構關系的變分。

數值算例驗證了三場對偶mortar 元的有效性，表明該方法可實現約束交叉問題的高精度求解，可用于解決非協調網格的計算需求。