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

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

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

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

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

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

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

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

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

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

1 三場變分法

1.1 數(shù)學描述

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

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

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

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

式中：usα為實體面Γα上的位移。Γ mα為實體面和媒介面之間經(jīng)過投影之后所產(chǎn)生的重疊區(qū)域，具體投影方法將在2.2 節(jié)詳細介紹。

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

1.2 弱形式

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

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

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

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

2 三場對偶mortar 元

2.1 空間離散

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

式中：Nα和Nm分別為子區(qū)域和媒介面的單元形函數(shù)矩陣；dα和dm為相應(yīng)的離散位移。

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

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

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

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

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

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

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

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

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

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

2.2 界面投影及數(shù)值積分

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

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

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

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

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

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

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

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

2.3 線性方程組

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

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

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

2.4 乘子凝聚

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

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

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

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

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

3 數(shù)值算例

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

3.1 接觸分片試驗

接觸分片試驗常用于檢驗接觸數(shù)值算法處理非協(xié)調(diào)網(wǎng)格的能力，一般僅考慮兩個子區(qū)域。為檢驗處理約束交叉問題的能力，此處將2 m×2 m× 2 m 的計算區(qū)域考慮為5 個子區(qū)域。實體網(wǎng)格如圖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。媒介面網(wǎng)格如圖4(b)所示，單元為各向等長，長度為1/3 m。

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

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

計算結(jié)果如圖5 所示。在圖5(a)中，黑色線條表示位移等值線。可見，計算所得最大豎向位移為-2×10-5m，與理論解吻合，且圖中未見明顯位移不連續(xù)現(xiàn)象。在圖5(b)中，黑色線條表示網(wǎng)格線?？梢姡嬎闼秘Q向應(yīng)力為均勻分布，圖示的數(shù)值在-1±10-6Pa 之間。為反映計算準確性，定義位移、應(yīng)力及能量的相對誤差如下：

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

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

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

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

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

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

3.2 帶圓孔平板

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

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

根據(jù)圖6(b)所示意的10 個子區(qū)域，剖分了如圖7 所示的疏密相間三維計算網(wǎng)格，以檢驗數(shù)值算法處理復雜子區(qū)域劃分的能力。計算區(qū)域平面尺寸為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 帶圓孔平板計算網(wǎng)格 Fig.7 Computational mesh of the plate with circular hole

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

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

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

圖9 帶圓孔平板的協(xié)調(diào)有限元網(wǎng)格 Fig.9 Conforming mesh of the plate with circular hole

圖10 數(shù)值計算結(jié)果與解析解的比較 Fig.10 Comparison of the numerical and analytical results

4 關(guān)于媒介面網(wǎng)格的討論

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

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

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

5 結(jié)論

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

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

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

(3) 相比Nistche 方法，無需增加額外的穩(wěn)定項，也無需推導本構(gòu)關(guān)系的變分。

數(shù)值算例驗證了三場對偶mortar 元的有效性，表明該方法可實現(xiàn)約束交叉問題的高精度求解，可用于解決非協(xié)調(diào)網(wǎng)格的計算需求。