劉傳奇 許廣濤 ,2 魏宇杰 ,2,3,*

1 中國科學(xué)院力學(xué)研究所非線性力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,北京 100190

2 中國科學(xué)院大學(xué)未來技術(shù)學(xué)院,北京 100049

3 中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049

1 發(fā)展歷程

虛單元方法(virtual element method,VEM)的基本驅(qū)動來源于對任意形狀多邊形的處理.傳統(tǒng)意義的有限元或有限差分需要將具備顯著幾何特征的物理實(shí)體離散化而求解,這在一定程度上喪失了對實(shí)體幾何信息的“宏觀”描述.而實(shí)際的工程需求中,越來越多的計(jì)算涉及到處理特定幾何體結(jié)構(gòu),如非凸性多邊形的變形以及復(fù)雜結(jié)構(gòu)的接觸等問題.2013年,由意大利 Milano-Bicocca 大學(xué)的 L.Beir?o da Veiga 教授組首先提出可適用于任意形狀多邊形的虛單元數(shù)值方法(Beir?o Da Veiga et al.2013a,2013b;Beir?o da Veiga et al.2014).由于該方法和有限元方法的高度兼容性,同時在處理懸掛節(jié)點(diǎn)(適用于含任意數(shù)目節(jié)點(diǎn)的單元)、接觸(均可轉(zhuǎn)化為點(diǎn)-點(diǎn)接觸)、晶體變形(對形狀、晶格匹配等無特殊要求)等問題方面具有特定優(yōu)勢,近年來得到計(jì)算力學(xué)領(lǐng)域的廣泛關(guān)注與發(fā)展.國際上不同的團(tuán)隊(duì)將該方法從理論上進(jìn)行了拓展,并應(yīng)用到不同的力學(xué)問題,如彈性問題 (Gain et al.2014,Artioli et al.2017a)、接觸問題 (Wriggers et al.2016,Wriggers and Rust 2019)、高階單元 (Beir?o da Veiga et al.2017a)、大變形 (Beir?o da Veiga et al.2015,Wriggers and Rust 2019)、斷裂問題 (Hussein et al.2018,Aldakheel et al.2019a),等工程科學(xué)/力學(xué)問題的計(jì)算.

對比有限元方法,虛單元法也需要對幾何空間進(jìn)行離散,通過形成并求解線性方程組對實(shí)際問題進(jìn)行近似.不同之處在于: (1)有限元中所采用的近似函數(shù)為顯式多項(xiàng)式函數(shù);虛單元法中,近似函數(shù)除了多項(xiàng)式函數(shù)外,還有在單元邊界為連續(xù)多項(xiàng)式且單元內(nèi)部滿足某些條件(如進(jìn)行拉普拉斯運(yùn)算后為多項(xiàng)式)的函數(shù).這些函數(shù)在單元域內(nèi)無顯式表達(dá),這也是虛單元方法中“虛”字的由來.(2)有限元方法中自由度均是近似函數(shù)的取值,虛單元方法通過定義合理的自由度,在形成剛度矩陣時避免計(jì)算單元內(nèi)部近似函數(shù)的取值;(3)有限元的剛度矩陣中僅含有一項(xiàng),虛單元法的剛度矩陣包含了協(xié)調(diào)矩陣和穩(wěn)定矩陣,以此來保證計(jì)算的收斂.

為了實(shí)現(xiàn)摘要中所提出的目標(biāo),本文將采取有別于傳統(tǒng)綜述論文的行文方式,大幅度地增加了關(guān)于虛單元法基本原理和最新理論的詳細(xì)介紹.這樣處理的目的是雙重性的: 一方面,讀者能夠了解該方法在一些典型問題中所體現(xiàn)的有別于有限元方法的能力;與此同時,對該方法感興趣的讀者將可以通過對本文的通讀和理解,全面掌握虛單元法的理論基礎(chǔ)和實(shí)現(xiàn)途徑,為利用該方法來發(fā)展適用于自身面臨的科學(xué)與工程計(jì)算問題提供一份全面的參考資料.基于這一設(shè)想,文章的基本框架為: 第 2 章簡略介紹有限元法與虛單元法的區(qū)別與聯(lián)系,以便具有力學(xué)背景的同行對虛單元法有個直觀了解;第 3 章詳細(xì)介紹虛單元法基本原理方面的發(fā)展,以熟知的泊松方程為對象,介紹了虛單元法的核心理論;第 4 章針對線彈性問題的求解,全面綜述了求解過程中涉及自由度、單元剛度矩陣構(gòu)建、相應(yīng)計(jì)算流程的具體方法和步驟,并給出了典型的計(jì)算實(shí)例來考察這一計(jì)算方法;第 5 章詳細(xì)介紹了當(dāng)前虛單元在非線性問題、斷裂問題、接觸問題、多場耦合等涉及復(fù)雜幾何邊界處理的典型應(yīng)用;第 6 章總結(jié)該方法的利弊并展望虛單元法的發(fā)展前景.

為便于介紹,需要在此做一些符號及其運(yùn)算的基本設(shè)定.以二維問題為例,下文中黑斜體表示張量（張量函數(shù)）或矩陣,比如位移場u,剛度矩陣K;斜體表示標(biāo)量函數(shù),比如泊松方程的解u;若帶有下標(biāo)則表示矢量的分量形式,比如u=uiei(i=1,2),下文為描述簡單省去單位基矢量ei.?()表示空間梯度運(yùn)算,比如?u=du/dxi為矢量;?·() 表示散度運(yùn)算,比如?·u=dui/dxi為標(biāo)量,這里重復(fù)下標(biāo)采用了愛因斯坦求和約定;Δ()=?·?() 表示拉普拉斯運(yùn)算,比如為標(biāo)量.

2 虛單元法與有限元法的異同概述

在詳細(xì)闡述虛單元法的基本理論前,需要對比有限元法并指出兩種方法的異同,以便工程力學(xué)背景的同行對該方法有個直觀了解.需要指出,該章意為對比基本概念,相關(guān)定義并不嚴(yán)格,更嚴(yán)格推導(dǎo)與解釋詳見下兩章.

考慮Ω域內(nèi)定義的標(biāo)量場u(x),x ∈Ω,為得到該場的近似解uh(x),兩種方法均需進(jìn)行空間離散Ω ≈∪iEi,其中Ei為單元網(wǎng)格.

網(wǎng)格要求的不同: 有限元需要網(wǎng)格為特定節(jié)點(diǎn)數(shù)目的凸單元,如三角形單元、四邊形單元等;虛單元法中單元可為任意節(jié)點(diǎn)數(shù)目,且對單元凸/凹性沒有限制.其原因?yàn)? 有限元中需要進(jìn)行單元映射,單元為凸才能保證映射矩陣的雅可比行列式為正,即單元體積恒正;而虛單元法無需進(jìn)行單元映射,可直接采用空間離散單元的幾何信息進(jìn)行近似.

在單元E上,近似函數(shù)uh(x) 均表示成形函數(shù)與待解自由度的乘積,即

其中,nd為待解自由度的數(shù)目,?i(x) 為節(jié)點(diǎn)i的形函數(shù),dofi(uh) 為節(jié)點(diǎn)i處的自由度.根據(jù)變分原理,剛度矩陣分量的一般形式為

式中推導(dǎo)采用了高斯公式,f,g,l為根據(jù)具體問題確定的函數(shù).有限元確定剛度矩陣主要依賴式(2)中第二個等號左邊直接進(jìn)行體積積分,而虛單元主要依賴式(2)中第二個等號右邊,需要在單元邊界上進(jìn)行線積分以及單元內(nèi)部進(jìn)行體積積分(通過自由度設(shè)定可避免).因而,形函數(shù)與自由度的設(shè)定均有不同.

形函數(shù)?i(x)形式的不同: 有限元中單元形函數(shù)為顯式的多項(xiàng)式函數(shù);虛單元法形函數(shù)除了多項(xiàng)式函數(shù)也可以為在單元邊界上?E為多項(xiàng)式、單元內(nèi)部滿足拉普拉斯運(yùn)算后 Δ?i為多項(xiàng)式的函數(shù),該函數(shù)可沒有顯式表達(dá).如此,虛單元方法將形函數(shù)的形式進(jìn)行了擴(kuò)充.

自由度設(shè)定的不同: 有限元中待解量為函數(shù)在節(jié)點(diǎn)處的取值,即:dofi(uh)dofi(uh)=uh(xi);虛單元法中,在單元邊界上的節(jié)點(diǎn)處的待解量為函數(shù)取值,但在單元內(nèi)部的待解量為形函數(shù)的拉普拉斯運(yùn)算與多項(xiàng)式的乘積.如此,可避免對沒有顯式表達(dá)的形函數(shù)進(jìn)行體積分.

具體實(shí)施起來,有限元法由于形函數(shù)形式確定,可直接對式(2)中第二個等號左邊進(jìn)行體積積分,剛度矩陣因而僅有一項(xiàng).需要指出,在特定單元與積分格式會造成“沙漏”,需要進(jìn)行積分格式的調(diào)整或者增加穩(wěn)定項(xiàng).而對于虛單元方法需進(jìn)行線積分以及體積分,由于形函數(shù)在邊界上為多項(xiàng)式,因而線積分容易求得;通過自由度設(shè)定可避免進(jìn)行單個形函數(shù)的體積分,但對于牽扯兩個形函數(shù)?i,?j的積分便需要引入映射操作符 Π 來使得 Π?i為多項(xiàng)式,如此

式中,為協(xié)調(diào)剛度矩陣分量,為穩(wěn)定剛度矩陣分量,其具體形式可根據(jù)收斂準(zhǔn)則給出.

剛度矩陣的不同: 有限元的剛度矩陣可直接通過體積積分獲得,僅含有一項(xiàng);虛單元法的剛度矩陣包含協(xié)調(diào)剛度矩陣和穩(wěn)定剛度矩陣,其中穩(wěn)定剛度矩陣通過收斂準(zhǔn)則確定.

至此,通過對比有限元中的概念,簡述了虛單元法的基本理論,第三章從數(shù)學(xué)角度探討虛單元法的完備性;第四章從應(yīng)用角度介紹虛單元方法在彈性問題中的應(yīng)用流程,若對數(shù)學(xué)完備不感興趣的讀者,可直接跳至第四章.

3 虛單元法的基本原理與其在泊松方程中的應(yīng)用

為解釋虛單元法的基本理論,特別是數(shù)學(xué)上的完備性,需要補(bǔ)充部分基本理論,涉及到函數(shù)空間、離散方程推導(dǎo)的基本流程、解的存在唯一性與收斂條件、以及多項(xiàng)式函數(shù)空間與虛單元法特定的近似函數(shù)形式等的討論.為便于理解且從信息的全面性考慮,這里作一些基本介紹.

3.1 函數(shù)空間

實(shí)數(shù)的集合稱為實(shí)數(shù)集,類似的,函數(shù)的集合稱為函數(shù)空間.比如函數(shù)u(x) 和v(x) 隸屬于同一個函數(shù)空間(為表述更簡潔,下文中省去自變量x),若在作用域Ω內(nèi)定義了兩者的內(nèi)積

則稱這樣的函數(shù)空間為內(nèi)積空間.對于內(nèi)積空間的任意函數(shù)u,進(jìn)一步定義一個實(shí)數(shù)‖u‖Lp(Ω)表示其大小,稱為函數(shù)的模

其中,p為整數(shù).若‖u‖Lp(Ω)＜∞,則稱 Lebesgue 空間(下標(biāo)L的來源).在本文討論中僅涉及L2空間(即p=2,并略去了作用域Ω的符號).進(jìn)一步,若函數(shù)及其導(dǎo)數(shù)滿足

其中,k為整數(shù),則稱為 Hilbert 空間(符號H的來源).在本文討論中僅涉及H1空間(即k=1 ).因?yàn)榭紤]了函數(shù)的導(dǎo)數(shù),需要將函數(shù)的模進(jìn)一步細(xì)化為半模(單豎線表示)與模(雙豎線表示),分別定義為

其中,下標(biāo)表示最高求導(dǎo)階次.

3.2 離散化方程組(泊松方程及其離散化)

無論何種數(shù)值方法最終都將變?yōu)榍蠼饩€性方程組,因而本小節(jié)將以簡單的泊松方程為例,簡述推導(dǎo)離散方程組的基本流程,以便后續(xù)的討論.

記空間域Ω的邊界為?Ω,考慮

將包含真實(shí)物理解的函數(shù)空間記為V,則方程(8)對應(yīng)的弱形式可通過分部積分并考慮邊界條件推得,即尋找u ∈V,使得

其中,a(u,v),(f,v) 分別稱為雙線性格式與線性格式,與式(4)相對應(yīng),分別定義為

將包含有限個元素的近似函數(shù)空間記為Vh ?V,待解問題變?yōu)? 尋找uh ∈Vh,使得

需要指出,由于Vh ?V,對于?uh,vh ∈Vh,除了在Vh空間內(nèi)定義雙線性格式ah(uh,vh),還可以在V空間內(nèi)定義相應(yīng)的a(uh,vh).為避免引入過多的數(shù)學(xué)定義,而不能專注于問題本身,在下文不引起混淆的情況下,把a(bǔ)h中的下標(biāo)h略去①h 一般代表空間離散尺寸,考慮到單元離散后必對應(yīng)近似函數(shù)空間,為避免引入過多符號,直接采用h符號表示近似空間,且在定義變量時是上標(biāo)、在定義運(yùn)算時是下標(biāo).,但在考慮解的收斂性時需要加以區(qū)分.

進(jìn)一步,將Vh中線性無關(guān)的基函數(shù)記為?i(i=1,2,···,n) (考慮到Vh的定義,n為有限值),則對于函數(shù)uh,vh ∈Vh,將有如下線性組合形式

其中,Uj,Vi為實(shí)數(shù).將式(12)代入到式(11)中,并考慮到vh的任意性,則有

由于雙線性格式的對稱性,記Kij=a(?i,?j),Fi=(f,?i),若采用矩陣寫法K=(Kij),U=(Uj),F=(Fi),則式(13)可寫為

再進(jìn)一步,將空間域Ω離散為有限個空間單元,Ω ≈Ωh ≡∪iEi,此時,作用于整個空間域Ω上的函數(shù)空間Vh將由作用在單元E上的函數(shù)空間Vh|E集合構(gòu)成,對應(yīng)著整體剛度矩陣與右端項(xiàng)將由單元內(nèi)積分(亦可以理解為全域積分,但在單元外部的函數(shù)取值為零)組裝得到,即

其中,∑為組裝操作,由于不易產(chǎn)生歧義,直接采用了相加符號.

至此,上述理論是有限元法與虛單元法共同的理論基礎(chǔ).

3.3 存在性、唯一性和收斂條件

對于數(shù)值求解,自然地需要考慮式(9)、式(11)是否存在解且解是否唯一;當(dāng)將計(jì)算域進(jìn)行空間離散后(記空間離散特征尺寸為h),還需要考慮收斂性,即h→0 時,‖u-uh‖→0 是否成立.分別討論如下:

(1)解的存在性與唯一性

對于式(9),真實(shí)解所在的空間為在邊界取值為零的一階 Hilbert 空間,即V=H01,下標(biāo) 0 表示函數(shù)在邊界處取值為零.如果存在C ＞0,使得

稱雙線性格式a(·,·) 為連續(xù)的,類似地,連續(xù)性需要|(f,v)|≤||f||·||v||m.進(jìn)一步如果存在α ＞0,使得

稱a(·,·) 為強(qiáng)制的(coercive) 或者Hm-橢圓的(elliptic).可以證明強(qiáng)制的雙線性格式對應(yīng)的方程存在且具有唯一解,即著名的 Lax-Milgram 定理.需要指出,式(16)、式(17)中不等式右端項(xiàng)采用模度量或半模度量是等價的.對于式(9),可以證明

因而,式(9)存在唯一解.上述理論與相關(guān)證明,可參閱 (Braess 2007,Brenner and Scott 2008).對于式(11),其解的存在性亦需要驗(yàn)證ah(uh,vh) 是否是強(qiáng)制的,即: 是否滿足式(16)和式(17),相關(guān)討論將在下一節(jié)展開.

(2)解的收斂性

判斷解是否收斂,需要驗(yàn)證協(xié)調(diào)性(consistency)和穩(wěn)定性(stability)(Ralston and Rabinowitz 2001).協(xié)調(diào)性指當(dāng)h→0 時,式(14)的解滿足微分方程與邊界條件,即趨向真實(shí)解;穩(wěn)定性指式(14)有解,即:K非奇異.因而可以通過給定一些簡單的真實(shí)解,來判斷方程的解是否滿足協(xié)調(diào)性與穩(wěn)定性,即進(jìn)行分片測試(patch test) (Taylor et al.1986).

3.4 多項(xiàng)式函數(shù)空間與虛單元法特定的近似函數(shù)形式

3.4.1 多項(xiàng)式函數(shù)空間

函數(shù)空間中,最直觀的為多項(xiàng)式函數(shù)構(gòu)成的空間.記Pk為最高階次不超過k的多項(xiàng)式函數(shù)構(gòu)成的空間,Pk為對應(yīng)的矢量函數(shù)空間.

(1) 在一維問題中,記ξ為無量綱坐標(biāo),則有

Pk共有k+1 個線性無關(guān)函數(shù),此外,記P-1=0;

(2) 在二維幾何空間的標(biāo)量問題中,記ξ,η為無量綱坐標(biāo),則有

Pk共有 (k+2)(k+1)/2 個線性無關(guān)的函數(shù);

(3) 在二維幾何空間的矢量問題中,則有

共有 (k+2)(k+1) 個線性無關(guān)的矢量函數(shù).對于小變形下的彈性問題,將P1變?yōu)?/p>

如此,并不改變函數(shù)空間內(nèi)線性無關(guān)多項(xiàng)式矢量函數(shù)的數(shù)目,但利用前三項(xiàng)對位移矢量u進(jìn)行近似時,可對應(yīng)著應(yīng)變?yōu)榱愕那闆r,即發(fā)生剛體位移.比如,若u=C(-η,ξ)T,其中,C為常數(shù),對應(yīng)的應(yīng)變?yōu)?=[?u+(?u)T]/2=0.如此調(diào)整,在生成剛度矩陣時,僅需少量操作,便可避免剛性位移,具體討論將在下一章給出.

3.4.2 虛單元法特定的近似函數(shù)形式

首先考慮如下問題,在二維多邊形E域內(nèi),標(biāo)量函數(shù)u ∈H1滿足

即函數(shù)u進(jìn)行拉普拉斯運(yùn)算后得到最高階次不超過k-2 的多項(xiàng)式,并且在多邊形邊界上u為連續(xù)的且不超過k次的多項(xiàng)式,那么函數(shù)u一定是多項(xiàng)式么？

通過反證法,顯然答案是否定的.下文將介紹虛單元法的函數(shù)空間由滿足上述條件的函數(shù)構(gòu)成(多項(xiàng)式函數(shù)以及非多項(xiàng)式函數(shù)),因?yàn)樗鶎?yīng)的非多項(xiàng)式?jīng)]有顯式形式(或者不易求出),而被稱為虛單元法.與之對應(yīng),有限元法的函數(shù)空間僅為多項(xiàng)式空間,因而,不同方法的核心區(qū)別是近似函數(shù)空間Vh的差異.

3.5 虛單元法在泊松方程中的應(yīng)用

由于泊松方程待解量為標(biāo)量函數(shù),相對簡單,本節(jié)將以泊松方程為例,介紹虛單元法的基本思路與流程.相比于下一節(jié)介紹的彈性方程,本節(jié)略側(cè)重于數(shù)學(xué)理論,以保證方法的數(shù)學(xué)完備.本節(jié)按照逆序的思路,首先給出式(11)存在唯一解以及收斂的條件,其次引出虛單元方法來實(shí)現(xiàn)上述條件,最后給出具體的實(shí)施細(xì)節(jié).

3.5.1 目標(biāo)

空間離散后,對于任意的Vh ?V,假定

即整體的雙線性格式可由單元內(nèi)的雙線性格式組裝得到,如此僅需對單元內(nèi)的雙線性格式進(jìn)行討論.對于任意的多邊形單元E,相應(yīng)的空間離散特征尺度h,以及定義在單元上的近似空間Vh|E,如果存在一個正整數(shù)k≥1,使得Pk(E)?Vh|E,且有

(1) 對于?p ∈Pk(E) 以及?vh ∈Vh|E,滿足

(2) 存在兩個不依賴于E以及h的正數(shù)α*和α*,對于?vh ∈Vh|E,滿足

此時,可以證明式(11)存在唯一解,且有足夠的精度,保證

條件(1)保證真實(shí)解為最高階次不超過k的多項(xiàng)式時,近似解可逼近真實(shí)解,即協(xié)調(diào)性條件;條件(2)保證ah的度量與a相同,此時能保證生成的剛度矩陣是滿秩的,稱為穩(wěn)定性條件.上述證明可參考 (Beir?o Da Veiga et al.2013a).

為滿足條件(1),首先定義函數(shù)空間VE,k為式(19)的解所構(gòu)成的空間.其次,定義映射操作符 ΠEk:VE,k→Pk(E)?VE,k,滿足: 對于?v ∈VE,k有

至此,對于?u,v ∈Vh|E,令aEh(u,v)=aE(ΠEk u,ΠEk v),則能滿足條件(1),即式(21).

為滿足條件(2),假定一個對稱正定的雙線性格式SE(u,v),對于?v ∈VE,k滿足 ΠEk v=0 的前提下,滿足

其中,正數(shù)c0,c1與單元和單元尺度無關(guān).根據(jù)式(25),對于?u ∈VE,k,則有:令v=u-ΠEk u,代入式(26),則對于?u ∈VE,k,有

至此,令

此時,考慮式(25),易證明式(28)的定義滿足條件(1),即式(21);考慮到式(26)以及aE(ΠEk v,ΠEk v)+aE(v-ΠEk v,v-ΠEk v)=aE(v,v),亦可證明式(28)滿足條件(2),即式(22),其中α*=max{1,c1},α*=min{1,c0}.

基于上述討論,虛單元法的總體思路可以簡述為: 在滿足式(19)的VE,k空間中,定義滿足式(24)的映射操作符 ΠEk:VE,k→Pk(E)?VE,k,并選用滿足式(27)的SE來定義雙線性格式(28).

3.5.2 實(shí)施方案

虛單元法的具體實(shí)施中,包括三個核心細(xì)節(jié): 根據(jù)VE,k設(shè)定自由度;構(gòu)造SE;以及計(jì)算線性格式(離散方程組的右端項(xiàng)).現(xiàn)分別予以討論.

(1)根據(jù) VE,k設(shè)定自由度

考慮多邊形單元E,記其形心為xE,邊數(shù)為n,特征尺寸為hE.

對于式(19)的第一個條件,對于多邊形的任意邊,需該邊上k+1 個點(diǎn)的函數(shù)取值則可確定最高階次為k的多項(xiàng)式函數(shù),這些點(diǎn)可由該邊的 2 個頂點(diǎn),以及邊內(nèi)k-1 個平均分布的點(diǎn)(或位于邊內(nèi)線積分點(diǎn)的 位置)構(gòu)成(當(dāng)k=1 時,每個多項(xiàng)式函數(shù)均為線性函數(shù),則無需邊內(nèi)節(jié)點(diǎn)).由于多邊形的頂點(diǎn)被 2 條邊共用,因而,共需要n+(k-1)n=nk個自由度即可確定?E上連續(xù)的多項(xiàng)式函數(shù).

對于式(19)的第二個條件,對于任意給定的f,則式(19)的解是唯一的.因而若滿足f ∈Pk-2(E),由于Pk-2共有k(k-1)/2 個線性無關(guān)的多項(xiàng)式函數(shù)(預(yù)備知識 3.4.1 中已予以討論),需要增加k(k-1)/2個自由度來保證式(19)的解的函數(shù)空間的維度,這些自由度設(shè)定在單元內(nèi)部點(diǎn)上.

因而,共需Ndof=dimVE,k=nk+k(k-1)/2 個自由度.圖1表示了五邊形(n=5 )單元在k=2時,自由度的設(shè)定位置.

圖1 k=2的單元自由度設(shè)定位置(E:單元； e:單元邊界)

在確定完自由度位置后,現(xiàn)在確定自由度的取值.將節(jié)點(diǎn)與單元邊界上的設(shè)定自由度的點(diǎn)的集合記為nf,單元內(nèi)部設(shè)定自由度的點(diǎn)的集合記為nb,并記xi處自由度取值為 dofi(vh).

若xi ∈nf,則該點(diǎn)處自由度取值為近似函數(shù)的取值,即: dofi(vh)=vh(xi);

若xi ∈nb,則該點(diǎn)處自由度取值需要考慮雙線性格式的定義,即,對于?uh,vh ∈VE,k

上式推導(dǎo)中采用了分部積分以及高斯定理,n為單元的外法向,且有:n·?vh=?vh/?n.由于VE,k為滿足式(19)的解構(gòu)成的空間,因而,在?E上vh,?uh/?n為多項(xiàng)式,可由xi ∈nf上設(shè)定的自由度確定(可設(shè)定邊內(nèi)自由度的點(diǎn)位于線積分點(diǎn)的位置,從而可以精確計(jì)算上式右端第一項(xiàng)的單元邊界積分);單元內(nèi)部積分僅包含右端第二項(xiàng),而考慮到式(19)中 Δuh ∈Pk-2(E),因而把xi ∈nb處的自由度取值設(shè)定為

其中,AE為單元面積.需要指出,一個xi ∈nb對應(yīng)一個p,Pk-2共有k(k-1)/2 項(xiàng),則共有k(k-1)/2個單元內(nèi)部點(diǎn)設(shè)定了自由度.如此,可以避免體積分來確定剛度矩陣,且避免了確定uh在單元內(nèi)的取值.因而,虛單元的一個優(yōu)勢為: 通過設(shè)定合理的自由度取值,在規(guī)避求解形函數(shù)在單元內(nèi)部的取值的同時,避免進(jìn)行體積分來確定剛度矩陣.

對自由度的取值進(jìn)一步討論如下.

① 式(30)中p為無量綱的多項(xiàng)式,因而式(30)的量綱與vh的量綱相同.特別的,當(dāng)p=1時,該自由度對應(yīng)著vh在單元內(nèi)部的平均值.在單元E內(nèi)的無量綱度量可以為

② 與有限元類似,每個自由度對應(yīng)著一個基函數(shù)?i(x),使得

此時,亦有:dofi(?j)=δij.由于vh(x)無顯式表達(dá)式,則?i(x)亦無顯式表達(dá)式.

③ 考慮到式(28),與上述自由度對應(yīng)的局部剛度矩陣為

(2)穩(wěn)定項(xiàng)SE的構(gòu)造

式(33)右端第一項(xiàng)生成的剛度矩陣并不滿秩,因而需要增加穩(wěn)定項(xiàng).考慮到穩(wěn)定性需要滿足式(27),最直觀的設(shè)定為

需要指出,在上述推導(dǎo)中,?i和 ΠEk ?i被表示成基函數(shù)?r的線性組合,即:對于 ΠEk ?i同理.在合理的hE的設(shè)定下,可以保證aE(?r,?r)=O(1),因而可以設(shè)定

來滿足式(27).需要指出穩(wěn)定項(xiàng)的設(shè)定有多種形式,更多的討論可參考 (Beir?o da Veiga et al.2017b).

(3)離散方程組右端項(xiàng)的構(gòu)造

考慮到單元離散,有

根據(jù)k的取值不同,fh的定義與 (fh,?i)E的構(gòu)造可分為三種情況

① 對于k=1,fh為f在單元E內(nèi)的平均值,即

此時

上述推導(dǎo)中,考慮了?i(xv)=δiv.

② 對于k=2,類似式(24),需要定義線性格式的映射操作符 Π0k:L2(E)→Pk(E),對于?w ∈L2(E),滿足

此時,fh=Π0kf,類似可定義 Π0k?i.如此,反復(fù)利用式(39),可構(gòu)造單元右端項(xiàng)為

其中,Π0k可以通過單元內(nèi)部自由度式(30)求出,更多實(shí)施細(xì)節(jié)可參考 4.2.4 節(jié).

③ 對于斜體>k＞2,類似的,定義映射操作符Π0k-2:L2(E)→Pk-2(E)fh=Π0k-2f,以及Π0k-2?i.構(gòu)造單元右端項(xiàng)為

其中,Π0k-2的計(jì)算可參考 4.2.4 節(jié).

如此,可證明(Beir?o Da Veiga et al.2013a)

至此,虛單元求解泊松問題時單元內(nèi)部的剛度矩陣為式(33),右端項(xiàng)根據(jù)k的取值在式 (38)、式(40)、式(41)三者中選擇合適的公式計(jì)算,進(jìn)行組裝后,便可進(jìn)行近似求解.

4 虛單元法處理線彈性問題

上一章以泊松方程為例,介紹了虛單元法的核心理論,側(cè)重解釋為什么引入一些處理,比如自由度的設(shè)定、多項(xiàng)式映射等,并強(qiáng)調(diào)數(shù)學(xué)完備.本章以彈性問題為例,通過對比有限元,具體介紹虛單元法的應(yīng)用過程.為獨(dú)立于上一章節(jié),本章直接從應(yīng)用角度理解虛單元,部分概念略有重復(fù),但更側(cè)重物理意義.

4.1 彈性問題描述

以二維連續(xù)彈性體為例,如圖2所示,物體Ω在準(zhǔn)靜態(tài)、小變形條件下的平衡條件為

圖2 二維彈性邊值問題

其中,σ為柯西應(yīng)力張量,f是為單位體積的體力.物體邊界記為?Ω,邊界條件為

其中,n為邊界法向,為指定應(yīng)力,且有?Ω=?Ωt ∪?Ωu以及?Ωt ∩?Ωu=?.本構(gòu)方程為σ=C?,其中C為四階材料彈性常數(shù),?為應(yīng)變張量,定義為位移梯度的對稱部分,即?=[?u+(?u)T]/2.

控制方程的弱形式為: 尋找u ∈V使得①計(jì)算數(shù)學(xué)的微元符號常采用 dx(如上章所述),本章更側(cè)重于力學(xué)解釋,因而體積微元為 dΩ,界面微元為d?Ω

其中,v為試矢量,V為包含真實(shí)解的矢量函數(shù)空間,其中每個分量隸屬于H1(Ω).通過分部積分,考慮應(yīng)力的對稱性以及代入邊界條件,式(45)可寫成如下形式

其中,a(u,v) 與L(v) 分別為彈性問題的雙線性與線性形式,定義為

同樣,記含有限個基函數(shù)的近似空間為Vh,相應(yīng)地,近似的雙線性格式ah(uh,vh) 和線性格式Lh(vh)形式與式(47)、式(48)類似,這里不再重復(fù).

4.2 離散方程

如前所述,虛單元法與有限元法類似,也需要將計(jì)算域近似為有限個單元構(gòu)成的區(qū)域,即:Ω ≈Ωh ≡∪iEi.不同于有限元中單元需要為凸多邊形來保證雅可比矩陣的行列式為正,在虛單元中單元Ei可以為任意形狀多邊形(凸、凹).本節(jié)將重復(fù)虛單元法的映射操作符與自由度設(shè)定,并給出剛度矩陣與右端項(xiàng)的具體矩陣表達(dá).

4.2.1 映射操作符

對于矢量函數(shù)空間,同樣定義映射操作符 ΠEk,將定義在單元E上的近似函數(shù)空間Vh(E) 映射到最高階次不超過k的多項(xiàng)式空間Pk(E),ΠEk:Vh(E)→Pk(E),下文中在不產(chǎn)生歧義的情況下,Π≡ΠEk.式(24)亦稱為正交條件,即:

由于雙線性形式在某種程度上表示了單元內(nèi)部的能量,因而正交條件表示了以能量為度量時,uh與 Πuh的誤差并不能通過不超過k次的多項(xiàng)式空間進(jìn)行度量.

如此,利用式(49),定義在單元上的雙線性格式為:uh,vh ∈Vh

右端第一項(xiàng)為映射操作符對uh,vh操作后在單元E內(nèi)的積分,第二項(xiàng)為進(jìn)行映射操作后沒能考慮的能量,又被稱為穩(wěn)定項(xiàng),需要能夠確保aE(uh,vh) 為強(qiáng)制的,保證存在唯一解.

4.2.2 自由度設(shè)定

在彈性問題中,近似函數(shù)空間中的分量亦需滿足式(19).對于nv邊形單元,且邊界上近似函數(shù)為最高階次不超過k的多項(xiàng)式時,需在

●nv個節(jié)點(diǎn)布置函數(shù)取值的自由度;

● 每條邊界上存在k-1 個邊界內(nèi)部節(jié)點(diǎn),在其上布置函數(shù)取值的自由度;

● 在單元內(nèi)部布置vh與最高階次不超過k-2 的多項(xiàng)式內(nèi)積的自由度(矩形式)

如第三章關(guān)于自由度的討論所述,對于標(biāo)量問題共需nd=nvk+k(k-1)/2 個自由度,因而對于二維矢量問題,共需 2nd個自由度.

在一個維度的函數(shù)空間內(nèi)定義操作符 dofi(vh),i=1,2,···,nd,表示vh在第i個自由度的取值,?i(x) 為第i個自由度對應(yīng)的基函數(shù),如前所述,此時有

當(dāng)設(shè)定自由度的點(diǎn)處均有兩個自由度(二維矢量)時,?表示基函數(shù)的向量表達(dá),有?1=[?1,0]T,?2=[0,?1]T,···,?2nd-1=[?nd,0]T,?2nd=[0,?nd]T,且有

需要指出,盡管虛單元法的空間離散格式(53)與有限元法相同,但由于節(jié)點(diǎn)形函數(shù)并無顯示表達(dá),單元內(nèi)部任意位置處的位移并不能直接獲取.但網(wǎng)格節(jié)點(diǎn)的位移是待解量,可通過節(jié)點(diǎn)位移后處理插值獲取位移場、應(yīng)力場等信息.

4.2.3 單元剛度矩陣

首先回顧傳統(tǒng)有限元中單元剛度矩陣的構(gòu)造.二維彈性問題的自由度為節(jié)點(diǎn)位移

其中,下標(biāo) 1,2 分別表示兩個維度,上標(biāo) 1,2,···,nd為單元節(jié)點(diǎn)編號.形函數(shù)矩陣為

應(yīng)變?yōu)??=BuE=?NuE,其中?在二維 Voigt 表示下,為

從式(47)出發(fā),并考慮vh的任意性,可以推得有限單元法中單元剛度矩陣為

其中,C是材料彈性常數(shù)矩陣.

虛單元法中的剛度矩陣構(gòu)造過程與有限元法類似,但如式(50)所示,需包含兩部分

其中,KcE為協(xié)調(diào)剛度矩陣(與映射操作符的協(xié)調(diào)性相關(guān)),KsE為穩(wěn)定剛度矩陣(與穩(wěn)定項(xiàng)相關(guān)).

(1)協(xié)調(diào)剛度矩陣

最高階次不超過k的多項(xiàng)式基矢量集為:{pα}α=1,2,···,nk,且pα ∈Pk(E),如 3.4.1 中所述,nk=(k+2)(k+1).映射操作符 Π 對式(53)中的矢量基函數(shù)?i的操作定義為多項(xiàng)式基的線性組合,即

式中,si,β為常數(shù).按照正交條件式(49),需要滿足

令α=1,2,···,nk,則可以得到如下矩陣形式

其中

以及

代入所有的矢量基函數(shù)?i,i=1,2,···,2nd,bi組裝成矩陣,映射系數(shù)si組裝成矩陣Π?,即

此時,由式(61),得

該條件可視為vh和 Πvh在pα(α=1,2,3) 度量下的平均值相等.如此,矩陣和矩陣G需要做如下修正 (Mengolini et al.2019)

對于矩陣G,前三行外的其他元素均可通過數(shù)值積分直接求得,從而可完全確定矩陣G的各個元素.

其中,為單元邊界?E的外法向量.

對于式(70)右端第一項(xiàng),如前所述,pα ∈Pk(E),則:?·(C ?(pα))∈Pk-2(E),進(jìn)一步考慮線性組合

其中,dα,β為有量綱的系數(shù).考慮到內(nèi)部自由度的定義式(51),可以看出,式(70)右端第一項(xiàng)與單元內(nèi)部自由度相關(guān),且有 dof2nvk+β(?i)=δ2nvk+β,i.

對于式(70)右端第二項(xiàng),單元邊界積分為各條邊上線積分的貢獻(xiàn)之和,每條邊上的線積分按如下數(shù)值積分計(jì)算

其中,ej為單元邊界?E的第j條邊,是該邊外法線矢量.可將邊界與頂點(diǎn)設(shè)定的自由度置于每條邊界上的數(shù)值積分點(diǎn)ξr(其相對應(yīng)的權(quán)重為wr),如此,可提高計(jì)算精度.

作為總結(jié),矩陣G的前三行按照式(69)求得,剩余行可直接數(shù)值積分得到;矩陣的前三行按照式(68)求得,剩余行通過式(70)～式(72)求得;矩陣按照式(66)求得.進(jìn)一步,將矩陣G作一變換,令其前三行所有元素設(shè)為 0,其他行所有元素不變,記此矩陣為.

如此,協(xié)調(diào)剛度矩陣的分量形式可表示為

相應(yīng)地,協(xié)調(diào)剛度矩陣為

(2)穩(wěn)定剛度矩陣

與有限元法不同,虛單元法中需要增加穩(wěn)定剛度矩陣.

首先,定義矩陣D2nd×nk為多項(xiàng)式基pα在自由度設(shè)定位置上的自由度取值,即

對應(yīng)的矩陣形式為

考慮到自由度取值可為函數(shù)本身取值,亦可為函數(shù)與多項(xiàng)式乘積的積分(矩形式).由此,可以確定矩陣D的各個元素:

① 對于邊界自由度: 1 ≤i≤2nvk,dofi(pα) 定義為基矢量pα在第i個自由度所屬節(jié)點(diǎn)的矢量分量值.

②對于單元內(nèi)部自由度: 2nvk+1 ≤i≤2nd,考慮到單元內(nèi)部自由度的定義式(51)

按照數(shù)值積分即可求得矩陣D的剩余元素.

進(jìn)一步,定義

與式(35)相對應(yīng),二維彈性問題的穩(wěn)定項(xiàng)可以取

將式(78)中Πij的定義,以及單位矩陣Iij=dofi(?j) 的定義代入式(79),則穩(wěn)定項(xiàng)的矩陣形式為

可以證明(Beir?o da Veiga et al.2014,2013b),若存在一個因子σE ＞0,滿足

其中,常數(shù)σ*,σ*均不依賴于單元特征尺寸h,則穩(wěn)定項(xiàng)式(80)乘以該因子σE依然能夠保證計(jì)算結(jié)果收斂.通常取σE=τh·tr(kcE)＞0,從而得到穩(wěn)定剛度矩陣

式中,τh為穩(wěn)定系數(shù),對于彈性問題,τh=0.5,更多選擇可參考 (Gain et al.2014);tr(·) 為矩陣的跡,引入此項(xiàng)可保證單元內(nèi)部的能量相對于單元尺寸和材料彈性常數(shù)為相對正確的比例(Artioli et al.2017a).

最后,需要指出矩陣G恰好可以表示成矩陣和矩陣D的乘積:G=.因此,如果在計(jì)算過程中首先求出了矩陣和矩陣D,就無需再通過式(63)進(jìn)行求解.當(dāng)然,按式(63)求出矩陣G,可幫助檢查兩種方式的求解結(jié)果是否相同.

4.2.4 右端項(xiàng)

如式(48),單元內(nèi)部的線性格式為

由于近似函數(shù)在單元邊界上為多項(xiàng)式,式中第二項(xiàng)無需特別處理,但對于第一項(xiàng),近似函數(shù)在單元內(nèi)部無顯式表達(dá),因此需要進(jìn)行近似處理.在不產(chǎn)生歧義的情況下,本節(jié)下述右端項(xiàng)均指與體積力f對應(yīng)的右端項(xiàng).

與泊松方程類似,對于彈性問題,右端項(xiàng)也需要進(jìn)行合理近似

根據(jù)k的取值不同,fh的定義與 (fh,?i)E的構(gòu)造可分為三種情況

(1) 對于k=1,fh為f在單元E內(nèi)的平均值,即

(2) 對于k=2,類似式(49)定義的雙線性格式的正交條件,首先定義線性格式的映射操作符Π0k:VE,k→Pk(E),對于??i ∈VE,k,滿足

其中

式中,si,β為常數(shù).下面將討論如何計(jì)算si,β.

將式(88)代入式(87),可得方程組

將其寫成矩陣形式

其中

以及

代入所有矢量基函數(shù)?i,i=1,2,···,2nd,將qi0組裝成矩陣Q0,將映射系數(shù)si0組裝成矩陣Πk0,即

方程組(90)可整理為

至此,系數(shù)矩陣Π0k可由式(95)求得.下面討論H0,Q0的求解.對于矩陣可通過已知的基矢量pα(α=1,2,···,nk) 完全確定.對于Q0的前nk-2行

其中

考慮到單元內(nèi)部自由度的定義式(51),可將基函數(shù)在單元內(nèi)部的自由度寫成

因此,矩陣Q中的元素均可以確定

對于Q0的第nk-2+1 行至第nk行(需要特別注意nk-2和nk分別是Pk-2和Pk空間的線性無關(guān)基函數(shù)的個數(shù),兩者并非相差 2),需要計(jì)算基函數(shù)?i和基矢量pα在單元內(nèi)部的積分,然而?i沒有顯式表達(dá),考慮到式(60)與式(89),由于已按照式(66)求得,可直接采用的第nk-2+1行至第nk行進(jìn)行替換,即

如此,矩陣Q0nk×2nd的計(jì)算公式如下

通過式(95)求出Π0k,反復(fù)利用式(87),可構(gòu)造右端項(xiàng)為

(3) 對于k ＞2,類似的,定義映射操作符 Π0k-2:VE,k→Pk-2(E),對于??i ∈VE,k,滿足

其中

式中,ri,β為常數(shù).將式(104)代入式(103),則有方程組

其矩陣形式為

代入所有矢量基函數(shù)?i,i=1,2,···,2nd,將qi組裝成矩陣Q,恰好為式(97),將映射系數(shù)ri組裝成矩陣Π0k-2,即

方程組(106)可以整理為

其中,矩陣Hnk-2×nk-2可通過已知的基矢量pα(α=1,2,···,nk-2) 完全確定;矩陣Qnk-2×2nd可通過式(97)、式(99)確定.

此時,通過式(109)求出操作符 Π0k-2的矩陣表達(dá)式Π0k-2,可構(gòu)造右端項(xiàng)

如上所述,可以在三種不同情況下求得基函數(shù)?i所對應(yīng)的離散方程右端項(xiàng).如此,單元內(nèi)與體積力f對應(yīng)的右端項(xiàng)為

4.3 計(jì)算流程

本小節(jié)給出單元矩陣與右端項(xiàng)的計(jì)算流程,如統(tǒng)一輸入量為: 內(nèi)插階數(shù)k,多邊形單元E,單元節(jié)點(diǎn)坐標(biāo)X=[x1,x2,···,xnv],積分點(diǎn)坐標(biāo)為ξ,權(quán)重為γ,材料彈性常數(shù)矩陣C,體積力f;輸出量為: 單元矩陣kE,單元右端項(xiàng)rE(包含與體積力f對應(yīng)的右端項(xiàng)rfE).

Π ←D ?ΠVh(E)(10) //映射到 空間的操作符式(78)(11) //計(jì)算剛度矩陣kc E ←?ΠT ?G ?Π(12) //協(xié)調(diào)剛度矩陣式(74)ks E ←τhtr(kc E)(I -Π)T(I -Π)(13) //穩(wěn)定剛度矩陣式(82)kE ←kc E +ks (14)(15) // 計(jì)算右端項(xiàng)rf E,rE ←02nd×1(16) //初始化k=1 E(17) if: then Xfrf ErE(18) 根據(jù) 和 計(jì)算,//右端項(xiàng)式(86)、式(83)(19) end(20) else if then H0αβ ←∫E pα·pβ dEH0nk×nk(21) // 矩陣式(92)Qαi ←∫E pα·?i dEQnk-2×2nd k=2(22) // 矩陣式(97)Q0 ←Q,H0,G, ?BQ0(23) // 矩陣式(101)Π0k ←(H0)-1 Q0nk×2ndΠ0k(24) // 的 矩陣式(95)XfΠ0krf ErE(25) 根據(jù),,計(jì)算,//右端項(xiàng)式(102)、式(83)(26) end(27) else Hαβ ←∫E pα·pβ dEHnk-2×nk-2 nk×2nd(28) // 矩陣式(105)Qαi ←∫E pα·?i dEQnk-2×2nd (29) // 矩陣式(97)Π0 k-2 ←H-1Qnk-2×2ndΠ0k-2(30) // 的 矩陣式(109)k-2rf ErE(31) 根據(jù),,計(jì)算,//右端項(xiàng)式(110)、式(83)Xf Π0(32) end

4.4 應(yīng)用實(shí)例

如前所述,虛單元法可處理任意節(jié)點(diǎn)數(shù)目的多邊形且對單元凸凹性沒有限制.為做出直觀展示,本小節(jié)以含有“力”字單元的懸臂梁模型為例,討論收斂性.本小節(jié)結(jié)果為階次為k=1 的情況.

考慮到懸臂梁的位移場存在理論解(Timoshenko and Goodier 1951).如圖3所示,在平面應(yīng)變狀態(tài)下單位厚度的懸臂梁左端,根據(jù)理論解施加位移邊界條件,右端施加沿 y 軸分布的剪應(yīng)力:

圖3 懸臂梁模型

此時,水平向與垂向的位移為

為展示虛單元的優(yōu)勢,采用含“力”字單元的網(wǎng)格,分別計(jì)算四種網(wǎng)格密度,如圖4(a) 所示.圖4(b) 為計(jì)算得到的水平向位移場,可以看出位移場連續(xù)且光滑.

圖4 網(wǎng)格與水平向位移場.(a)網(wǎng)格,(b)水平向位移

圖5 收斂曲線

圖6展示了不同縱向剖面(如圖3所示x=2.0 m,x=4.0 m,x=7.2 m)的水平向位移與理論解的比較,可以看出虛單元方法在不均勻的多邊形網(wǎng)格下,隨著網(wǎng)格的加密能夠趨于收斂.

圖6 不同縱向剖面的水平向位移與理論解比較

5 虛單元方法的深入應(yīng)用

上一章系統(tǒng)梳理了理解虛單元基本理論的預(yù)備知識;側(cè)重?cái)?shù)學(xué)推導(dǎo),闡述了虛單元在泊松方程的應(yīng)用;側(cè)重物理解釋,介紹了虛單元在二維彈性問題的應(yīng)用過程.如前言所述,虛單元方法的相關(guān)研究已從計(jì)算數(shù)學(xué)領(lǐng)域逐漸拓展到工程科學(xué)/力學(xué)領(lǐng)域.本章對虛單元方法在計(jì)算力學(xué)領(lǐng)域的研究進(jìn)展進(jìn)行歸納與總結(jié),包括材料與幾何非線性、接觸與界面問題、斷裂問題等.為避免內(nèi)容過度重合,每部分僅引用近 5年的代表性工作.

5.1 材料與幾何非線性

本節(jié)首先梳理小變形的相關(guān)工作,側(cè)重于單元曲直、階數(shù)的研究進(jìn)展,其次按照準(zhǔn)靜態(tài)問題與動力學(xué)問題進(jìn)行區(qū)分,討論材料與幾何非線性問題的研究進(jìn)展.需要指出,晶體塑性的相關(guān)工作在下一節(jié)接觸與界面問題中予以討論.

5.1.1 小變形問題

對于低階直線單元,Artioli et al.(2017a) 提出虛單元方法求解線彈性問題的標(biāo)準(zhǔn)化代碼實(shí)現(xiàn)流程,詳細(xì)介紹了協(xié)調(diào)項(xiàng)、穩(wěn)定項(xiàng)和載荷右端項(xiàng)的構(gòu)造過程,結(jié)果表明,同階的虛單元方法與有限元方法具有相近的準(zhǔn)確性和相同的收斂階數(shù),且對網(wǎng)格畸變幾乎不敏感.此后,Artioli et al.(2017b) 將其工作擴(kuò)展到非線性材料中,研究了黏彈性模型、含各向同性與運(yùn)動硬化的 Mises 塑性模型、以及記憶合金本構(gòu)模型.Beir?o da Veiga et al.(2015) 對小變形條件下,虛單元的非線性彈性本構(gòu)與 Mises 屈服準(zhǔn)則的彈塑性本構(gòu)的應(yīng)用過程進(jìn)行了總結(jié).Gain et al.(2014) 將虛單元法應(yīng)用到三維的小變形彈性問題中,并對比了不同類型網(wǎng)格的收斂性.

針對單元階數(shù),Beir?o da Veiga et al.(2017a) 提出了高階單元方法求解對流反應(yīng)方程,Beir?o da Veiga et al.(2016)類比高階 Serendipity 有限元方法,提出了高階 Serendipity 虛單元方法求解彈性問題.與高階有限元方法類似,高階虛單元法對網(wǎng)格畸變問題具有更高的魯棒性,并且由于增加了自由度,計(jì)算精度更高.為避免離散曲線邊界時需要“以直代曲”,針對單元曲直,Beir?o da Veiga et al.(2019) 將單元邊界的描述轉(zhuǎn)化為自然坐標(biāo),提出了曲邊單元,并對橢圓方程進(jìn)行了求解.Artioli et al.(2020) 將其方法應(yīng)用到二維線彈性問題中,但曲線形狀不能任意.Wriggers et al.(2020) 提出可將參考構(gòu)型的單元(直線邊)映射到初始構(gòu)型(曲線邊)中,實(shí)現(xiàn)任意形狀的單元離散,對于映射操作可采用等參映射或者 NURBS 映射.而后,Wriggers et al.(2021b) 將該方法擴(kuò)展到高階單元中.

5.1.2 準(zhǔn)靜態(tài)問題的幾何與材料非線性

對于彈性條件下的幾何非線性問題,Wriggers et al.(2017) 通過對應(yīng)變能的討論,提出了有限變形的穩(wěn)定項(xiàng)構(gòu)建方法,通過可壓縮和不可壓縮材料的模擬,驗(yàn)證了其方法的收斂性和魯棒性,如圖7所示.同一時期,Chi et al.(2017) 通過構(gòu)造局部位移空間準(zhǔn)確求得體積變化,根據(jù)切線模量的跡構(gòu)建穩(wěn)定項(xiàng),討論了混合單元(位移、壓力)格式與純位移格式的有限變形模型.圖8是其模型對彈性多孔材料的模擬結(jié)果.針對不可壓縮材料,Taylor-Hood 是混合單元格式,采用不同的階數(shù)插值位移與壓力,以保證 LBB 穩(wěn)定條件,Wriggers et al.(2021a) 將 Taylor-Hood 單元引入到虛單元法中,構(gòu)造了不可壓縮材料的有限變形格式.

圖7 有限變形下的沖壓問題.(a)模型,(b)可壓縮材料的沖壓變形,(c)不可壓縮材料沖壓變形(修改自 (Wriggers et al.2017))

圖8 彈性多孔材料的有限變形.(a)模型,(b)剪應(yīng)變?yōu)?1.345 下的最大應(yīng)變分布(修改自(Chi et al.2017))

若進(jìn)一步考慮材料非線性,Wriggers and Hudobivnik (2017) 采用低階單元從最小能量增量出發(fā),構(gòu)造了二維彈塑性本構(gòu)下有限變形的虛單元格式,并對精度和收斂性進(jìn)行了討論.Hudobivnik et al.(2019) 將該方法擴(kuò)展到三維問題中,并對比了不同收斂算法對結(jié)果的影響.圖9是對沖壓以及扭轉(zhuǎn)過程的模擬結(jié)果.

圖9 準(zhǔn)靜態(tài)大變形彈塑性問題的等效塑性應(yīng)變分布.(a)沖壓,(b)扭轉(zhuǎn) (修改自(Hudobivnik et al.2019))

需要指出,大部分虛單元方法在處理非線性問題時,均采用低階單元,De Bellis et al.(2019)成功將高階單元應(yīng)用到非線性彈性問題中.

5.1.3 動力學(xué)問題的幾何與材料非線性

對于動力學(xué)問題,Cihan et al.(2021b) 指出不同于剛度矩陣需要穩(wěn)定項(xiàng),虛單元方法的質(zhì)量矩陣只需要進(jìn)行映射操作即可,其中時間積分可按照隱式 Newmark 格式進(jìn)行.圖10是不考慮阻尼條件下,懸臂梁的振動模擬結(jié)果.

圖10 懸臂梁振動問題.(a)模型,(b)垂向位移時間演化(修改自 (Cihan et al.2021b))

對于彈塑性的動力學(xué)問題,Cihan et al.(2021a) 按照 Hu-Washizu 泛函求極值的基本思路,提出了混合單元格式,以避免彈塑性不可壓條件引入的體積自鎖問題,其中的時間積分亦采用隱式 Newmark 格式.圖11為泰勒桿與沖壓過程中的塑性應(yīng)變.若采用顯式時間積分,Park et al.(2020) 對彈性不可壓縮材料、Park et al.(2019) 對彈塑性材料分別進(jìn)行了討論,其中臨界時間步長都是通過矩陣最大的特征值確定.

圖11 彈塑性動力學(xué)問題的塑性應(yīng)變.(a)泰勒桿,(b)沖壓 (修改自 (Cihan et al.2021a))

上述討論的都是各向同性材料,對于橫觀各向同性材料（比如纖維增強(qiáng)）,虛單元法也被逐步應(yīng)用到小變形 (Wriggers et al.2018b,Reddy and van Huyssteen 2019) 與有限變形模擬中(Wriggers et al.2018a,van Huyssteen and Reddy 2021).

5.2 接觸與界面問題

由于虛單元法對網(wǎng)格形狀沒有要求,因而在網(wǎng)格配對處理上可得到極大簡化,在處理接觸與界面問題上具有先天優(yōu)勢.本節(jié)分別從接觸問題、非均勻材料的界面問題以及衍生的晶體塑性三個方面總結(jié)虛單元的研究進(jìn)展.

5.2.1 接觸問題

采用虛單元法處理接觸問題,是該方法一個非常重要的應(yīng)用.處理接觸問題需要判斷兩個物體間的相對位置,對于有限元網(wǎng)格(含節(jié)點(diǎn)、網(wǎng)格邊界線),如圖12(a)～圖12(c)所示,常用的方法有點(diǎn)-點(diǎn)、點(diǎn)-線、Mortar 型,其中點(diǎn)-點(diǎn)型最容易實(shí)現(xiàn),但需要網(wǎng)格匹配.由于虛單元對單元節(jié)點(diǎn)數(shù)目沒有要求,可任意增加節(jié)點(diǎn),如圖12(d)所示,可以將潛在接觸面的節(jié)點(diǎn)向目標(biāo)面做映射生成新的節(jié)點(diǎn),采用點(diǎn)-點(diǎn)格式進(jìn)行接觸判斷,如此可避免復(fù)雜的網(wǎng)格匹配、助于各體間獨(dú)立建模.需要指出,虛單元法中增加了單元節(jié)點(diǎn)后,增加節(jié)點(diǎn)的單元形式會發(fā)生改變,待解自由度增加,而其他單元沒有變化,這是與有限元節(jié)點(diǎn)插值最大的區(qū)別.

圖12 判斷兩物體間相對位置的常用方法.(a)點(diǎn)-點(diǎn),(b)點(diǎn)-線,(c)Mortar 型,(d)虛單元法插節(jié)點(diǎn)

對于兩體接觸,Wriggers et al.(2016) 對虛單元法在接觸問題中的應(yīng)用進(jìn)行了梳理,分別采用拉格朗日乘子法和罰函數(shù)法施加約束,其結(jié)果表明虛單元法利用插點(diǎn)格式處理接觸問題,精度高、操作簡單.圖13是采用不同方法對接觸問題進(jìn)行的分片測試,可以看出相對于有限元中的點(diǎn)-線格式,虛單元法求得的應(yīng)力均勻,精度更高.Wriggers and Rust (2019) 將該方法擴(kuò)展到大變形的摩擦接觸問題中,圖14為其模擬的赫茲接觸過程.對于接觸體具有曲線邊界的情況,Aldakheel et al.(2020) 引入了曲邊單元處理接觸,其中插入的節(jié)點(diǎn)需要按照曲邊單元進(jìn)行處理.

圖13 接觸問題的分片測試.(a)虛單元法,(b)有限元點(diǎn)-線格式(修改自(Wriggers et al.2016))

圖14 大變形條件下的赫茲接觸問題(修改自 (Wriggers and Rust,2019))

顆粒材料中顆粒間的相互作用涉及多體接觸問題,常采用離散元進(jìn)行描述,其核心假定為顆粒為剛體.虛單元法在網(wǎng)格上的任意性以及接觸處理的優(yōu)勢,使得快捷準(zhǔn)確地捕捉變形顆粒間的相互作用變?yōu)榭赡?Gay Neto et al.(2021) 首次做出了嘗試,采用完全的隱式時間積分、并且研究了剛度與質(zhì)量矩陣參數(shù)的敏感性.圖15模擬多面體顆粒集合的沙漏過程.該方法為顆粒材料研究提供一個新的思路與方向.

圖15 顆粒材料沙漏過程(修改自 (Gay Neto et al.2021))

5.2.2 不考慮斷裂的非均勻材料界面問題

對于非均勻材料,考慮兩方面的研究: 材料界面與均質(zhì)化.

關(guān)于材料界面問題,Lo Cascio et al.(2021) 采用虛單元法描述基材、邊界元法描述內(nèi)嵌材料,邊界元獲取的應(yīng)力-位移方程轉(zhuǎn)化為力-位移方程,組裝到虛單元法的方程組中,實(shí)現(xiàn)兩者的耦合.該方法有效地描述了非均勻材料的變形及損傷演化,圖16為典型算例的模擬結(jié)果.

圖16 虛單元法與邊界元法模擬非均勻材料.(a)網(wǎng)格,(b)單軸拉伸條件下的應(yīng)力分布(修改自 (Lo Cascio et al.2021))

多晶材料的均質(zhì)化是進(jìn)行多尺度建模的核心.對于有限元,增加晶粒的各向異性需增加大量的自由度,而虛單元法將任意形狀的晶粒視為一個網(wǎng)格,自由度數(shù)目自然大大降低.Marino et al.(2019)對線性和非線性的均質(zhì)化格式開展研究.若考慮電-磁-力耦合的多晶材料,B?hm et al.(2021b) 計(jì)算了不同晶格結(jié)構(gòu)和不同程度的各向異性材料,發(fā)現(xiàn)虛單元方法在低自由度條件下計(jì)算的均質(zhì)化性能都表現(xiàn)出很好的精度,如圖17所示.

圖17 多晶結(jié)構(gòu)均質(zhì)化.(a)虛單元法網(wǎng)格,(b)粗糙的四面體網(wǎng)格,(c)精細(xì)的四面體網(wǎng)格,(d) BaNiO3，中度各向異性，六方晶系晶胞(修改自 (B?hm et al.2021b))

5.2.3 晶體塑性

晶體塑性采用虛單元方法進(jìn)行模擬具有巨大潛力.虛單元法中單元的任意性可以完美地契合晶粒結(jié)構(gòu),無需進(jìn)行精細(xì)的網(wǎng)格剖分,但該部分工作尚屬起步階段.B?hm et al.(2021a) 采用虛單元模擬了 FCC 晶格結(jié)構(gòu)受單向拉伸的滑移過程,其中映射操作在晶粒的表面進(jìn)行以產(chǎn)生協(xié)調(diào)矩陣,穩(wěn)定項(xiàng)需要將晶格內(nèi)部劃分為多個四面體.圖18為單元分解過程與單滑的模擬結(jié)果.工業(yè)應(yīng)用方面,Behrens et al.(2020)將基于虛單元的晶體塑性模型應(yīng)用到軸承襯套的生產(chǎn)設(shè)計(jì)中,以描述多種材料連接處多晶材料的損傷與疲勞性質(zhì),圖19為模擬結(jié)果.

圖18 虛單元法模擬晶體塑性.(a)單個晶格網(wǎng)格分解為界面網(wǎng)格與內(nèi)部四面體網(wǎng)格,(b)單軸拉伸條件下 FCC 晶格結(jié)構(gòu)的剪應(yīng)變(修改自 (B?hm et al.2021a))

圖19 虛單元方法分析軸承襯套中不同材料連接區(qū)域的多晶塑性模型.(a)軸承襯套示意,(b)三種材料連接區(qū)域的代表性體積單元,(c) FCC 晶格在 (111) 滑移面上剪切應(yīng)變率 的分布,(d)等效von Mises 應(yīng)力分布(修改自 (Behrens et al.2020))

5.3 斷裂問題

虛單元法中網(wǎng)格的任意性簡化了空間離散的難度、并能隨意增加/減少節(jié)點(diǎn).虛單元法處理斷裂問題可以從直接單元切割(強(qiáng)間斷)以及與其他方法(相場、內(nèi)聚力、擴(kuò)展有限元等)的結(jié)合兩個角度進(jìn)行總結(jié).

5.3.1 直接單元切割

針對直接單元切割,Hussein et al.(2019)提出了沿裂紋擴(kuò)展方向增加新的單元節(jié)點(diǎn)并進(jìn)行單元分割的裂紋擴(kuò)展技術(shù),對比不同形狀的網(wǎng)格,驗(yàn)證其方法的有效性.圖20為單元切割算法與單軸拉伸下的裂紋擴(kuò)展.直接切割算法的優(yōu)勢是可直接確定裂紋面的位移間隔,在采用位移差求解應(yīng)力強(qiáng)度因子時大大簡化了程序設(shè)計(jì),并且由于可任意增加節(jié)點(diǎn),虛單元法在裂紋交叉、分叉、三維曲面裂紋的擴(kuò)展模擬中具有潛力.

圖20 虛單元法在斷裂問題中的應(yīng)用.(a)單元切割,(b)裂紋擴(kuò)展(修改自 (Hussein et al.2019))

5.3.2 與其他方法結(jié)合模擬斷裂

斷裂問題一直是計(jì)算力學(xué)的熱點(diǎn),在虛單元法提出以前,適用于不同場景的描述斷裂的方法已被廣泛發(fā)展.

與相場法耦合.考慮到裂紋擴(kuò)展方向可以容易從相場方法獲得,Hussein et al.(2020) 將相場法與單元切割融合,采用虛單元法中單元切割表征裂紋擴(kuò)展,未開裂部分的網(wǎng)格采用自適應(yīng)方案,來求解相場.該方法為裂紋擴(kuò)展模擬提供了高效、魯棒性強(qiáng)的模擬方案.圖21展示了該方法模擬的裂紋擴(kuò)展過程.

圖21 單元切割與自適性相場耦合模擬裂紋擴(kuò)展.(a)自適應(yīng)網(wǎng)格求解相場,(b)單元切割模擬裂紋擴(kuò)展(修改自 (Hussein et al.2020))

與內(nèi)聚力模型結(jié)合.Marfia et al.(2022)將虛單元法與內(nèi)聚力模型相結(jié)合,按照單元域內(nèi)平均的最大主應(yīng)力確定裂紋方向,通過內(nèi)聚力模型以及單元切割,模擬裂紋的起裂與擴(kuò)展過程,如圖22所示.進(jìn)一步考慮非均勻材料的界面引起的斷裂問題,Benedetto et al.(2018) 應(yīng)用零厚度的界面單元模擬應(yīng)力-裂紋張開過程,結(jié)果如圖23所示.

圖22 虛單元法與內(nèi)聚力模型模擬 L 形狀板的裂紋擴(kuò)展(修改自 (Marfia et al.2022))

圖23 虛單元法與內(nèi)聚力模型模擬耦合.(a)接觸脫離模擬,(b)非均勻材料三點(diǎn)彎曲梁模擬(修改自(Benedetto et al.2018))

與擴(kuò)展有限元法結(jié)合.擴(kuò)展有限元法是將近似函數(shù)空間進(jìn)行擴(kuò)充,采用間斷函數(shù)以及奇異函數(shù)描述物理場的間斷.Benvenuti et al.(2019) 將擴(kuò)展有限元法與虛單元進(jìn)行結(jié)合,對拉普拉斯方程的裂紋尖端進(jìn)行了描述,其中滿足拉普拉斯方程的間斷函數(shù)以及奇異函數(shù)被擴(kuò)充到基函數(shù)空間中,并推導(dǎo)了擴(kuò)充映射操作符,以便將擴(kuò)充的基函數(shù)空間映射到線性多項(xiàng)式與擴(kuò)充函數(shù)上.圖24所示為含裂紋的薄膜受 III 型載荷作用的模擬結(jié)果.Benvenuti 等 (2022)將擴(kuò)展有限元法與虛單元結(jié)合的方法進(jìn)一步應(yīng)用到線彈性斷裂問題中,著重討論了應(yīng)力強(qiáng)度因子的求解.

圖24 虛單元法與擴(kuò)展有限元結(jié)合模擬含裂紋的薄膜在 III 型載荷作用下變形.(a)100 個網(wǎng)格,(b)1600 個網(wǎng)格(修改自 (Benvenuti et al.2019))

5.4 多場耦合、拓?fù)鋬?yōu)化與地下裂隙網(wǎng)絡(luò)流體流動

本節(jié)總結(jié)了虛單元法在其他領(lǐng)域的應(yīng)用,包括熱-力耦合作用、拓?fù)鋬?yōu)化、以及地下裂隙網(wǎng)絡(luò)流體流動模擬等.

5.4.1 熱-力耦合作用

對于熱-力耦合問題,Dhanush and Natarajan (2019) 采用 Matlab 生成多邊形網(wǎng)格以及Abaqus 的輸入文件,將虛單元法成功應(yīng)用到 Abaqus 中來處理小變形下的熱-力耦合問題,并給出了詳細(xì)的數(shù)據(jù)格式.進(jìn)一步,Aldakheel et al.(2019b) 將虛單元法中有限變形下的彈塑性模型引入到熱-力耦合問題中,如此,自由能包含彈性能、彈-熱耦合部分、純熱部分、以及塑性勢能.圖25為三維鋼螺栓成型過程模擬.

圖25 三維鋼螺栓成型的熱力耦合過程.(a)～(f)等效塑性應(yīng)變,(g)～(l)絕對溫度場 T (Aldakheel et al.2019b)

5.4.2 拓?fù)鋬?yōu)化

拓?fù)鋬?yōu)化問題需要求解各種幾何形狀的邊值問題,由于虛單元法放松了對網(wǎng)格形狀的要求,簡化了復(fù)雜形狀的網(wǎng)格剖分難度,特別適用于優(yōu)化問題.Gain et al.(2015) 對比了虛單元法與有限元法的優(yōu)化結(jié)果,指出虛單元法只需要少量的單元,便可得到足夠精度.圖26展示了對懸臂梁的優(yōu)化結(jié)果,六面體網(wǎng)格需要 50 000 多個網(wǎng)格,而非規(guī)則多面體僅需 10 000 個網(wǎng)格.對于復(fù)合材料的大變形優(yōu)化問題,Zhang et al.(2020) 引入虛單元法提高計(jì)算效率,其中材料性能的確定是通過插值能量函數(shù)確定,而非插值材料參數(shù),且該函數(shù)可以解決粗糙網(wǎng)格面臨的單元畸變問題.進(jìn)一步,對于纖維增強(qiáng)的軟物質(zhì)的拓?fù)鋬?yōu)化問題,Zhang 等 (2021) 將基材與纖維參數(shù)作為兩組設(shè)計(jì)參數(shù),其中纖維的方向是預(yù)先設(shè)定的一組離散方向,并采用虛單元法求解有限變形的邊值問題.圖27為三點(diǎn)彎曲梁的優(yōu)化結(jié)果.

圖26 拓?fù)鋬?yōu)化問題.(a)懸臂梁模型,(b)六面體網(wǎng)格,(c)非規(guī)則多面體網(wǎng)格(修改自 (Gain et al.2015))

圖27 拓?fù)鋬?yōu)化問題.(a)目標(biāo)承載體,(b)優(yōu)化結(jié)果(修改自 (Zhang et al.2021))

5.4.3 地下裂隙網(wǎng)絡(luò)流體流動模擬

地下裂隙網(wǎng)絡(luò)中流體的流動是地質(zhì)力學(xué)中的重要問題.相對于裂隙流動,孔隙滲流過程一般忽略不計(jì),但三維空間中二維裂隙平面相互交叉、獨(dú)立,若采用傳統(tǒng)有限元協(xié)調(diào)網(wǎng)格,無法對每個裂隙單獨(dú)建模,必需整體考慮,如此網(wǎng)格生成變得十分困難.考慮到虛單元法對網(wǎng)格節(jié)點(diǎn)數(shù)目沒有要求,因而處理連接處的懸掛節(jié)點(diǎn)等十分自然,可以大大簡化網(wǎng)格生成的難度.Benedetto et al.(2016) 采用虛單元法對多種裂隙網(wǎng)絡(luò)下的流動進(jìn)行了嘗試,如圖28所示.此后,Benedetto et al.(2017) 將該方法擴(kuò)展到混合單元與高階單元.

圖28 地下裂隙流體流動問題.(a)網(wǎng)格生成,(b)裂隙網(wǎng)絡(luò)中水頭分布(修改自 (Benedetto et al.2016))

6 結(jié)論與展望

正如在該綜述的開頭所強(qiáng)調(diào)的,通過有別于傳統(tǒng)綜述論文的行文方式,詳述了虛單元方法的基本原理和最新理論進(jìn)展,同時也展示了目前該方法在一些典型問題中所體現(xiàn)的有別于有限元方法的潛力: 包括它對單元凸凹性的處理,以及因此延伸到如懸掛節(jié)點(diǎn)、接觸、裂紋擴(kuò)展、多晶體變形等問題方面的應(yīng)用.基于處理非凸性單元任意多邊形單元的考慮,虛單元方法引入了和有限元方法迥異的處理方法,最為顯著的是(1)虛單元法中的近似函數(shù)既包含了多項(xiàng)式函數(shù),同時也引入了在單元域內(nèi)無法顯式表達(dá)的函數(shù),這也是虛單元方法中“虛”字的由來;(2) 自由度取值在節(jié)點(diǎn)以及單元邊界上為近似函數(shù)的取值,在單元內(nèi)部為矩形式(近似函數(shù)與多項(xiàng)函數(shù)的乘積在單元內(nèi)的積分),避免計(jì)算非多項(xiàng)式在單元內(nèi)部的取值;(3)為保證計(jì)算的收斂性,虛單元法的剛度矩陣包含了協(xié)調(diào)矩陣和穩(wěn)定矩陣兩部分,對應(yīng)于多項(xiàng)式的精確解和多項(xiàng)式未被考慮的部分,比有限元方法更為復(fù)雜.

虛單元法尚屬于發(fā)展階段,在幾何與材料非線性問題、接觸與界面問題、斷裂問題等諸多方面具有發(fā)展?jié)摿?尤其需要關(guān)注它在以下幾個方面的應(yīng)用.

(1) 裂紋擴(kuò)展: 斷裂問題一直是計(jì)算力學(xué)的熱點(diǎn),在虛單元法提出以前,適用于不同場景的描述斷裂的方法已被廣泛發(fā)展.當(dāng)前虛單元法在這一類問題方面展現(xiàn)了非常好的兼容性.將相場法與單元切割融合,采用虛單元法中單元切割表征裂紋擴(kuò)展 (Hussein et al.2020),可以為裂紋擴(kuò)展模擬提供了高效、魯棒性強(qiáng)的模擬方案.將虛單元法與內(nèi)聚力模型相結(jié)合,通過內(nèi)聚力模型以及虛單元法中的單元切割 (Marfia et al.2022),可高效模擬裂紋的起裂與擴(kuò)展過程.將擴(kuò)展有限元法與虛單元進(jìn)行結(jié)合 (Benvenuti et al.2019,2022),可模擬諸多的線彈性斷裂問題.

(2) 接觸問題: 虛單元法對單元形狀沒有要求,可任意增加接觸面節(jié)點(diǎn),因而在處理接觸與界面問題時,可避免網(wǎng)格配對帶來的挑戰(zhàn),簡化計(jì)算與單元劃分問題.(Wriggers et al.2016) 對虛單元法在接觸問題中的應(yīng)用進(jìn)行了梳理,分別采用拉格朗日乘子法和罰函數(shù)法施加約束,其結(jié)果表明虛單元法利用插點(diǎn)格式處理接觸問題,精度高、操作簡單.對涉及大變形的摩擦接觸問題Wriggers and Rust (2019)和具有曲線接觸邊界的情況 (Aldakheel et al.2020),虛單元法也具有高效的處理方案.同樣對于涉及多體接觸的問題,虛單元法在網(wǎng)格上的任意性以及接觸處理方面,也具超越離散元的優(yōu)勢,可快捷準(zhǔn)確捕捉變形顆粒間的相互作用 (Gay Neto et al.2021).工程中涉及到類似的接觸問題不勝枚舉.

(3) 多晶體變形: 利用虛單元法來模擬晶體塑性具有巨大潛力.這是因?yàn)樘搯卧ㄖ袑Χ噙呅螁卧陌輰⒖梢院途Я＝Y(jié)構(gòu)完美地契合,一方面無需進(jìn)行精細(xì)地網(wǎng)格剖分,可以確保具有晶體取向的晶粒的整體性,同時又可以保持晶粒內(nèi)部不同區(qū)域之間的非局域作用,使得變形梯度、位錯塞積等晶體變形特征能夠得到妥善處理.盡管目前一些研究組已經(jīng)采用虛單元方法開展了晶格塑性模擬方面的應(yīng)用 (B?hm et al.2021a,Behrens et al.2020),該方面的工作還亟待開發(fā),以充分利用虛單元法對非規(guī)則單元處理方面的優(yōu)勢.

以上是筆者從自身的科研背景和需求出發(fā)所展望的虛單元計(jì)算方法的發(fā)展?jié)摿?因此它不可能是全面的.同時需要強(qiáng)調(diào)的是,虛單元方法在這些特定問題方面,具有超越傳統(tǒng)有限元方法的某些優(yōu)勢,因此為飽受這些問題困擾的科研工作者提供了一種新的解決方案,但這并不意味著它能全面地替代有限元方法的成熟性、便利性以及歷經(jīng)長時間發(fā)展帶來的高可靠性.它作為一個特定模塊,結(jié)合相應(yīng)的有限元方法,將豐富和提升計(jì)算力學(xué)的處理能力和便利性.通過對兩者的融會貫通,促進(jìn)發(fā)展出適應(yīng)我國力學(xué)計(jì)算需求的新型算法與高性能軟件.

致 謝作者們感謝國家自然科學(xué)基金委基礎(chǔ)科學(xué)中心項(xiàng)目“非線性力學(xué)的多尺度問題研究”(資助號 11988102)支持和面上項(xiàng)目支持(資助號 12172368,劉傳奇),劉傳奇同時感謝中科院百人計(jì)劃項(xiàng)目支持;作者們對意大利 Milano-Bicocca 大學(xué) L.Beir?o da Veiga 教授和美國加州大學(xué)N.Sukumar 教授關(guān)于論文算法富有成效的討論深表謝意.