基于三場(chǎng)變分原理的對(duì)偶mortar 有限元法

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

(1. 北京交通大學(xué)城市地下工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 100044；2. 清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

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

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

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

圖1 交叉點(diǎn)或交叉線問(wèn)題 Fig.1 Cross-point or cross-line problems

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

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

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

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

1 三場(chǎng)變分法

1.1 數(shù)學(xué)描述

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

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

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

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

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

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

1.2 弱形式

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

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

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

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

2 三場(chǎng)對(duì)偶mortar 元

2.1 空間離散

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

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

乘子λα離散在實(shí)體面Γα上：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 線性方程組

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

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

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

2.4 乘子凝聚

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

式中，ds表示實(shí)體面與媒介面的相對(duì)位移，可通過(guò)矩陣C表示如下：

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

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

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

3 數(shù)值算例

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

3.1 接觸分片試驗(yàn)

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

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

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

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

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

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

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

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

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

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

3.2 帶圓孔平板

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

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

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

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

圖7 帶圓孔平板計(jì)算網(wǎng)格 Fig.7 Computational mesh of the plate with circular hole

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

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

圖8 帶圓孔平板的數(shù)值計(jì)算結(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ù)值計(jì)算結(jié)果與解析解的比較 Fig.10 Comparison of the numerical and analytical results

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

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

在實(shí)際工程應(yīng)用中，滿足第一個(gè)特點(diǎn)并不會(huì)引入很大困難。這是因?yàn)槊浇槊嫦啾扔趯?shí)體區(qū)域而言，只需剖分降維后的網(wǎng)格。相當(dāng)于以剖分一個(gè)協(xié)調(diào)的面網(wǎng)格的代價(jià)，實(shí)現(xiàn)了不同子區(qū)域的網(wǎng)格各自獨(dú)立剖分，且子區(qū)域的劃分形式不限。

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

5 結(jié)論

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

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

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

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

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