虛單元法

虛單元法（Virtual Element Method，VEM）是一種用於數值求解偏微分方程的新型方法，它的最大特點是可以使用多邊形/多面體格線，從而針對複雜幾何形狀以及斷裂、接觸等問題時具有優勢。它是有限元法（FEM）的推廣，旨在克服傳統有限元法在處理複雜幾何形狀和非標準單元格線時的某些局限性。

基本介紹

• 中文名：虛單元法(虛擬元法)
• 外文名：Virtual Element Method
• 所屬學科：數學、力學
• 其他譯名：虛擬元法

算法簡介

虛單元法（Virtual Element Method，VEM）是近年來提出的一種新的數值計算方法，主要用於解決複雜幾何形狀下的偏微分方程問題，尤其適用於計算力學、材料科學等領域。它可以被看作是有限元法（FEM）的擴展，但在處理複雜幾何形狀、斷裂力學、接觸力學、板殼等問題時更有優勢。

虛單元法在斷裂力學中的套用

基本原理

映射運算元

虛單元法的核心思想是通過虛擬（即沒有顯式格式）的方式定義多邊形或者多面體單元中的試函式。其方法是構造變數的映射,並通過映射對計算域進行插值離散。對於多邊形單元，其邊界單元（一維線單元）的基函式是已知的。通過自由度的定義以及正交條件

可計算出映射運算元。

虛單元法中的多邊形單元

通過映射運算元，可進一步構造離散雙線性格式。在計算映射運算元及離散雙線性格式的過程中，只需要沿著單元邊界進行積分，不需要多邊形內部的插值函式，因此可適用於任意多邊形單元。

離散格式

虛單元法將變數分解成映射項和餘項，因此可寫為。對於泊松方程，將代入至雙線性格式中，可得

因此，虛單元法所獲得的剛度矩陣分為兩部分：協調剛度矩陣和穩定剛度矩陣（穩定項）。這也是虛單元法和有限元法的另一個重要區別。最新的技術可以消去虛單元法的穩定項。

特點

多邊形單元

虛單元法的最大特點是自動適應任意形狀的多邊形單元（三維問題可適用多面體單元）。和眾多多邊形有限元法不同，虛單元法可以更靈活地處理非凸多邊形，例如懸掛節點單元，內凹節點單元等。這些特點使得虛單元法在處理斷裂、接觸等問題時更靈活。此外，也可以根據格線的靈活性，構造非匹配格線，從而更容易對複雜形狀進行格線劃分。格線的自適應也是虛單元法的一大重要優點。

虛單元法中的自適應格線

穩定項

虛單元法的另一個特點是剛度矩陣需要穩定項。由於使用映射運算元將變數投影到了多項式空間，因此會造成能量的損失。在剛度矩陣中的表現為協調剛度矩陣（直接通過映射形成）的秩不滿足條件，這將導致得到錯誤的計算結果。對於泊松方程，最簡單的穩定項構造方式為。在實際問題中，穩定項通常需要其他的參數。對於三角形單元，不需要額外的穩定項。穩定項的構造、穩定項參數的選擇等是比較棘手的問題。

虛單元法中不同空間的投影

C1連續

虛單元可以更靈活地定義自由度，因此在單元邊界定義變數的法嚮導數，從而進一步求解要求高階連續的方程（如板殼問題，高階相場模型等）

虛單元法求解Kirchhoff–Love殼問題

高階虛單元法

許多工程問題要求結果有較高的精度，因此需要高階虛單元法進行計算（一般是二階，）。

自由度

對於高階虛單元法，需要定義以下自由度：

• 單元頂點處的函式值
• 單元邊上個Gauss–Lobatto點的函式值
• 單元內部直到階的矩量

通過以上自由度的定義，可計算出虛單元法中的映射運算元，從而進一步計算出剛度矩陣。單元內部的額外自由度會造成計算量的上升。為此，可進一步構造serendipy單元，消去單元內部的額外自由度。

三階虛單元法中的自由度

投影運算元

考慮正交條件，其等號右端可進一步寫為

其中為單元邊界外法線。上式等號右端第一項不為0，可通過矩量自由度計算。右端第二項可通過高斯積分及邊界線單元的插值進行計算。因此映射運算元可解。

離散格式

高階虛單元法的離散格式和一階虛單元法的離散格式相同。對於泊松方程，其離散格式為

三維虛單元法

格線劃分

許多工程問題為三維問題。使用虛單元法求解三維問題時，需要將計算域劃分成多面體單元，從而計算映射運算元以及剛度矩陣。

三維問題中的格線劃分和多面體單元

投影運算元

對於三維多面體單元，其映射運算元的計算個二維相同：

通過分部積分及高斯散度定理，上式可寫為

上式等號右側第一項可通過自由度的定義計算。對於等號右側第二項，可以在單元邊界建立局部坐標系，將多面體單元表面視為二維多邊形單元（考慮增強空間），因此有

最終，可得

從而進一步求解出映射運算元。

多面體單元及多邊形表面

高階三維格式

大部分的工程問題為三維問題。為了較少計算中的“鎖定”現象，提高計算精度，往往需要用到二階單元進行計算。三維高階虛單元法的格式和程式都更複雜，僅有少數論文被發表。在力學問題中，可見文獻。

和高階二維格式相似，高階三維格式需要定義額外的自由度

單元頂點處的函式值

單元邊上個Gauss–Lobatto點的函式值

單元面內部直到階的矩量

單元內部直到階的矩量

二階三維單元的自由度定義

無穩定項虛單元法

傳統的虛單元法需要額外的穩定項，對於力學問題及非線性問題，穩定項的選擇對結果的影響較大。如何選擇合適的穩定項往往比較困擾研究者。因此，構造無穩定項的虛單元法(Stabilization-free virtual element method, SFVEM)顯得十分重要。

L2投影運算元

無穩定項虛單元法的基本思想是對變數的梯度構造高階的映射運算元，使得剛度矩陣的秩滿足要求，從而刪掉虛單元法中的穩定項。為了求解映射運算元，需要考慮如下正交條件

其中為階矢量多項式。等式的右端可進一步寫為

考慮到,則在當前的自由度定義下，映射運算元可求。根據單元的邊數選擇適當的多項式階數,可以得到無穩定項的虛單元格式。對於二維問題，多項式階數應滿足

定義映射運算元的矩陣形式為,則變數的梯度為

相關符號定義可見文獻

離散格式

無穩定項虛單元法（SFVEM）不需要額外的穩定項，其離散格式和傳統有限元法十分相似。對於泊松方程，其雙線性格式為

因此可得剛度矩陣為

可以看出單元剛度矩陣不需要任何穩定項。此外，對於其他有限元法可以求解的問題，基本上可以使用無穩定項虛單元法直接進行離散求解，如力學問題中的超彈性問題。

無穩定項虛單元法求解超彈性衝壓問題

三維格式

三維無穩定項虛單元法格式相關工作較少。到目前為止，主要工作由Bing-Bing Xu及Peter Wriggers完成，可見文獻。其基本原理與二維無穩定項格式相似，需要滿足以下正交條件

其中等號右側可繼續近似為

因此映射運算元可解。注意在此處的證明需要用到增強空間技術。

套用

泊松方程

傳熱問題、擴散-反應問題等控制方程為一系列泊松方程

虛單元法在此領域發表了大量論文，主要關注點在於算法格式的收斂性、穩定性等。

線彈性力學問題

虛單元法在力學問題中存在著廣泛的套用。線彈性力學問題的控制方程為

其中柯西應力滿足胡克定律,。其雙線性格式為

線彈性問題中，虛單元法的穩定項可選擇方法較多，較簡單的格式為

三維虛單元法求解複雜幾何模型的彈性變形

斷裂力學問題

在虛單元法中可以使用任意形狀的多邊形單元，因此非常適用於裂紋擴展問題。基於應力強度因子的計算和經典裂紋擴展算法，配合單元劃分算法，可以對裂紋擴展進行仿真計算。

虛單元法計算斷裂力學問題

有限變形彈性問題

虛單元法在求解超彈性問題時具有很大優勢。考慮Neo-Hookean超彈性本構，可以根據傳統連續介質力學理論構造其切線剛度矩陣。其中穩定項可選擇為

其中參數和當前狀態下的本構相關。與此同時，可以使用無穩定項虛單元法（SFVEM）求解超彈性問題，從而避免引入穩定項及相關參數。對於一階虛單元法，其計算過程中進需要用到點積分，因此計算效率高。在使用高階虛單元法時，會涉及到多邊形或者多面體單元上的域積分計算，會在一定程度上降低計算效率。

虛單元法求解超彈性桿的扭轉

動力學問題

動力學問題的方程為

其中為第一P-K應力張量。其相應弱形式為

可得及質量雙線性項為

其中穩定項可選擇為。配合時間差分格式，可計算出動力學回響。

虛單元法求解二維線彈性動力學問題

彈塑性力學問題

在彈塑性力學問題中，變形梯度被分解為彈性部分和塑性部分。左 Cauchy-Green張量為，其中塑性右 Cauchy-Green張量為。其勢能為

根據連續介質力學理論，可構造其切線剛度矩陣。對於無穩定項虛單元法，可根據傳統有限元的框架直接計算，見幾年Robert L. Taylor院士90歲生日的論文。

虛單元法和有限元法計算有限變形塑性問題

接觸力學問題

由於虛單元法可以使用多邊形單元，因此其在接觸力學（格線不匹配問題）中有較大的優勢。根據接觸面節點的插入算法，可以自動適應點-點接觸理論。此外，虛單元法還可以構造曲邊單元，模擬複雜情況下的接觸問題。由於虛單元法能夠承受較大的格線變形，因此使用第三介質接觸理論也會更有優勢。

旋轉體的接觸計算

相場斷裂模擬

虛單元法在求解相場斷裂問題時有著許多優勢。一階虛單元法具有比有限元法更高的效率，因此可以提高相場斷裂問題的計算速度。此外，虛單元法可以使用多邊形格線，因此在計算過程中，可以使用格線局部加密和自適應技術。

虛單元法計算相場斷裂問題

出版物

虛單元法及其工程問題的出版物較少，可見Peter Wriggers教授的相關工作：《Virtual Element Methods in Engineering Sciences》

虛單元法的相關出版物