虛單元計(jì)算方法的最新理論與應(yīng)用進(jìn)展

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

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

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

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

1 發(fā)展歷程

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

對(duì)比有限元方法,虛單元法也需要對(duì)幾何空間進(jìn)行離散,通過(guò)形成并求解線性方程組對(duì)實(shí)際問(wèn)題進(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)無(wú)顯式表達(dá),這也是虛單元方法中“虛”字的由來(lái).(2)有限元方法中自由度均是近似函數(shù)的取值,虛單元方法通過(guò)定義合理的自由度,在形成剛度矩陣時(shí)避免計(jì)算單元內(nèi)部近似函數(shù)的取值;(3)有限元的剛度矩陣中僅含有一項(xiàng),虛單元法的剛度矩陣包含了協(xié)調(diào)矩陣和穩(wěn)定矩陣,以此來(lái)保證計(jì)算的收斂.

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

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

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

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

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

網(wǎng)格要求的不同: 有限元需要網(wǎng)格為特定節(jié)點(diǎn)數(shù)目的凸單元,如三角形單元、四邊形單元等;虛單元法中單元可為任意節(jié)點(diǎn)數(shù)目,且對(duì)單元凸/凹性沒有限制.其原因?yàn)? 有限元中需要進(jìn)行單元映射,單元為凸才能保證映射矩陣的雅可比行列式為正,即單元體積恒正;而虛單元法無(wú)需進(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ù)具體問(wèn)題確定的函數(shù).有限元確定剛度矩陣主要依賴式(2)中第二個(gè)等號(hào)左邊直接進(jìn)行體積積分,而虛單元主要依賴式(2)中第二個(gè)等號(hào)右邊,需要在單元邊界上進(jìn)行線積分以及單元內(nèi)部進(jìn)行體積積分(通過(guò)自由度設(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)式的乘積.如此,可避免對(duì)沒有顯式表達(dá)的形函數(shù)進(jìn)行體積分.

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

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

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

至此,通過(guò)對(duì)比有限元中的概念,簡(jiǎn)述了虛單元法的基本理論,第三章從數(shù)學(xué)角度探討虛單元法的完備性;第四章從應(yīng)用角度介紹虛單元方法在彈性問(wèn)題中的應(yīng)用流程,若對(duì)數(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) 隸屬于同一個(gè)函數(shù)空間(為表述更簡(jiǎn)潔,下文中省去自變量x),若在作用域Ω內(nèi)定義了兩者的內(nèi)積

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)于式(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)制的雙線性格式對(duì)應(yīng)的方程存在且具有唯一解,即著名的 Lax-Milgram 定理.需要指出,式(16)、式(17)中不等式右端項(xiàng)采用模度量或半模度量是等價(jià)的.對(duì)于式(9),可以證明

因而,式(9)存在唯一解.上述理論與相關(guān)證明,可參閱 (Braess 2007,Brenner and Scott 2008).對(duì)于式(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 時(shí),式(14)的解滿足微分方程與邊界條件,即趨向真實(shí)解;穩(wěn)定性指式(14)有解,即:K非奇異.因而可以通過(guò)給定一些簡(jiǎn)單的真實(shí)解,來(lái)判斷方程的解是否滿足協(xié)調(diào)性與穩(wěn)定性,即進(jìn)行分片測(cè)試(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為最高階次不超過(guò)k的多項(xiàng)式函數(shù)構(gòu)成的空間,Pk為對(duì)應(yīng)的矢量函數(shù)空間.

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

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

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

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

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

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

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

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

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

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

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

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

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

3.5.1 目標(biāo)

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

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

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

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

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

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

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

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

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

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

至此,令

此時(shí),考慮式(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}.

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

3.5.2 實(shí)施方案

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

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

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

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

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

因而,共需Ndof=dimVE,k=nk+k(k-1)/2 個(gè)自由度.圖1表示了五邊形(n=5 )單元在k=2時(shí),自由度的設(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)處自由度取值需要考慮雙線性格式的定義,即,對(duì)于?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為單元面積.需要指出,一個(gè)xi ∈nb對(duì)應(yīng)一個(gè)p,Pk-2共有k(k-1)/2 項(xiàng),則共有k(k-1)/2個(gè)單元內(nèi)部點(diǎn)設(shè)定了自由度.如此,可以避免體積分來(lái)確定剛度矩陣,且避免了確定uh在單元內(nèi)的取值.因而,虛單元的一個(gè)優(yōu)勢(shì)為: 通過(guò)設(shè)定合理的自由度取值,在規(guī)避求解形函數(shù)在單元內(nèi)部的取值的同時(shí),避免進(jìn)行體積分來(lái)確定剛度矩陣.

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

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

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

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

③ 考慮到式(28),與上述自由度對(duì)應(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的線性組合,即:對(duì)于 ΠEk ?i同理.在合理的hE的設(shè)定下,可以保證aE(?r,?r)=O(1),因而可以設(shè)定

來(lái)滿足式(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)造可分為三種情況

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

此時(shí)

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

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

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

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

③ 對(duì)于斜體>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)

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

4 虛單元法處理線彈性問(wèn)題

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

4.1 彈性問(wèn)題描述

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

圖2 二維彈性邊值問(wèn)題

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

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

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

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

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

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

4.2 離散方程

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

4.2.1 映射操作符

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

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

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

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

4.2.2 自由度設(shè)定

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

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

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

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

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

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

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

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

4.2.3 單元?jiǎng)偠染仃?/p>

首先回顧傳統(tǒng)有限元中單元?jiǎng)偠染仃嚨臉?gòu)造.二維彈性問(wèn)題的自由度為節(jié)點(diǎn)位移

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

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

從式(47)出發(fā),并考慮vh的任意性,可以推得有限單元法中單元?jiǎng)偠染仃嚍?/p>

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

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

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

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

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

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

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

其中

以及

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

此時(shí),由式(61),得

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

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

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

對(duì)于式(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.

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

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

作為總結(jié),矩陣G的前三行按照式(69)求得,剩余行可直接數(shù)值積分得到;矩陣的前三行按照式(68)求得,剩余行通過(guò)式(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è)定位置上的自由度取值,即

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

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

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

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

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

進(jìn)一步,定義

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

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

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

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

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

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

4.2.4 右端項(xiàng)

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

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

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

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

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

(2) 對(duì)于k=2,類似式(49)定義的雙線性格式的正交條件,首先定義線性格式的映射操作符Π0k:VE,k→Pk(E),對(duì)于??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的求解.對(duì)于矩陣可通過(guò)已知的基矢量pα(α=1,2,···,nk) 完全確定.對(duì)于Q0的前nk-2行

其中

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

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

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

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

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

(3) 對(duì)于k ＞2,類似的,定義映射操作符 Π0k-2:VE,k→Pk-2(E),對(duì)于??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可通過(guò)已知的基矢量pα(α=1,2,···,nk-2) 完全確定;矩陣Qnk-2×2nd可通過(guò)式(97)、式(99)確定.

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

如上所述,可以在三種不同情況下求得基函數(shù)?i所對(duì)應(yīng)的離散方程右端項(xiàng).如此,單元內(nèi)與體積力f對(duì)應(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對(duì)應(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ù)目的多邊形且對(duì)單元凸凹性沒有限制.為做出直觀展示,本小節(jié)以含有“力”字單元的懸臂梁模型為例,討論收斂性.本小節(jié)結(jié)果為階次為k=1 的情況.

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

圖3 懸臂梁模型

此時(shí),水平向與垂向的位移為

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

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

5.1 材料與幾何非線性

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

5.1.1 小變形問(wèn)題

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

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

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

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

圖7 有限變形下的沖壓?jiǎn)栴}.(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)下有限變形的虛單元格式,并對(duì)精度和收斂性進(jìn)行了討論.Hudobivnik et al.(2019) 將該方法擴(kuò)展到三維問(wèn)題中,并對(duì)比了不同收斂算法對(duì)結(jié)果的影響.圖9是對(duì)沖壓以及扭轉(zhuǎn)過(guò)程的模擬結(jié)果.

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

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

5.1.3 動(dòng)力學(xué)問(wèn)題的幾何與材料非線性

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

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

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

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

上述討論的都是各向同性材料,對(duì)于橫觀各向同性材料（比如纖維增強(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èn)題

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

5.2.1 接觸問(wèn)題

采用虛單元法處理接觸問(wèn)題,是該方法一個(gè)非常重要的應(yīng)用.處理接觸問(wèn)題需要判斷兩個(gè)物體間的相對(duì)位置,對(duì)于有限元網(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)格匹配.由于虛單元對(duì)單元節(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)的單元形式會(huì)發(fā)生改變,待解自由度增加,而其他單元沒有變化,這是與有限元節(jié)點(diǎn)插值最大的區(qū)別.

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

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

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

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

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

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

5.2.2 不考慮斷裂的非均勻材料界面問(wèn)題

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

關(guān)于材料界面問(wè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)行多尺度建模的核心.對(duì)于有限元,增加晶粒的各向異性需增加大量的自由度,而虛單元法將任意形狀的晶粒視為一個(gè)網(wǎng)格,自由度數(shù)目自然大大降低.Marino et al.(2019)對(duì)線性和非線性的均質(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),無(wú)需進(jìn)行精細(xì)的網(wǎng)格剖分,但該部分工作尚屬起步階段.B?hm et al.(2021a) 采用虛單元模擬了 FCC 晶格結(jié)構(gòu)受單向拉伸的滑移過(guò)程,其中映射操作在晶粒的表面進(jìn)行以產(chǎn)生協(xié)調(diào)矩陣,穩(wěn)定項(xiàng)需要將晶格內(nèi)部劃分為多個(gè)四面體.圖18為單元分解過(guò)程與單滑的模擬結(jié)果.工業(yè)應(yīng)用方面,Behrens et al.(2020)將基于虛單元的晶體塑性模型應(yīng)用到軸承襯套的生產(chǎn)設(shè)計(jì)中,以描述多種材料連接處多晶材料的損傷與疲勞性質(zhì),圖19為模擬結(jié)果.

圖18 虛單元法模擬晶體塑性.(a)單個(gè)晶格網(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èn)題

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

5.3.1 直接單元切割

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

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

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

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

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

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

與內(nèi)聚力模型結(jié)合.Marfia et al.(2022)將虛單元法與內(nèi)聚力模型相結(jié)合,按照單元域內(nèi)平均的最大主應(yīng)力確定裂紋方向,通過(guò)內(nèi)聚力模型以及單元切割,模擬裂紋的起裂與擴(kuò)展過(guò)程,如圖22所示.進(jìn)一步考慮非均勻材料的界面引起的斷裂問(wèn)題,Benedetto et al.(2018) 應(yīng)用零厚度的界面單元模擬應(yīng)力-裂紋張開過(guò)程,結(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ù)描述物理場(chǎng)的間斷.Benvenuti et al.(2019) 將擴(kuò)展有限元法與虛單元進(jìn)行結(jié)合,對(duì)拉普拉斯方程的裂紋尖端進(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)用到線彈性斷裂問(wèn)題中,著重討論了應(yīng)力強(qiáng)度因子的求解.

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

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

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

5.4.1 熱-力耦合作用

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

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

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

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

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

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

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

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

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

6 結(jié)論與展望

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

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

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

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

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

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

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