在超彈材料的結構有限元仿真中，單元的鎖死現象（Locking）時常困擾著工程與研發人員，是什么樣的機理導致了這種數值失真，我們如何來通過科學的方法解決，從而得到正確的計算結果，是我們一直在探索的。而F-Bar數值方法似乎就可以很好的解決這個問題，今天我們來了解一下什么是F-Bar單元，它是如何完美解決超彈大變形的鎖死問題的。

1.  什么是體積鎖死現象

圖1 平面應變模型

圖1所示的由兩個三角形單元組成的平面應變問題中，如果變形體材料是不可壓縮的，那么我們不管在加載點施加多大的力，從有限元法計算得到的所有節點的所有位移都為零。因為任意位移都會導致兩個三角形中的一個的體積變化。這種有限元單元精度不足以表現所需體積應變的問題就是體積鎖死現象。

圖2 受內壓的線彈性厚壁圓筒；理論解（點線）和有限元計算結果的比較[1]

圖2給出了另外一個例子。這里顯示了受內壓的厚壁圓筒在彈性圓筒的泊松比為0.3時，有限元計算結果逼近理論解；而當泊松比逼近0.5時，有限元計算結果與理論解相差極大（實際上當泊松比逼近0.5時，有限元計算位移逼近于零）。

1.1  為什么會發生鎖死現象

這里給出一個非常粗糙的解釋: 我們通?？梢园炎冃误w的變形功W分解為幾個部分，比如說體積變形部分和非體積變形部分

W = W1 + W2

由于數值計算的誤差，當不同部分變形功的差相差很大時，將會帶來很大的舍入誤差，使計算結果變得很不可靠。具體來說，由于

    1） 單元尺寸帶來的誤差: 殼，梁單元是其典型。由于殼單元面內的尺寸和殼厚相差很大，使得面外應變與面內應變的計算精度相差很大，帶來剪切鎖死現象。由于細分后網格的單元邊長和殼厚將會減小, 剪切鎖死現象可以通過細分網格來減輕。

   2）由于材質帶來的誤差: 不可或近似不可壓縮材料如泊松比趨于0.5的彈性材料，塑性材料，這時的體積應變對應于巨大的變形能，因而帶來體積鎖死。體積鎖死現象不能靠細分網格來回避。

   3）由于物理現象不同帶來的誤差。比如說在做熱-變形耦合分析時，溫度（在三角形單元中線性分布）和應變（在三角形單元中均勻分布）計算精度帶來的誤差。

   對于該問題的進一步解釋建議參看文獻[2]，在哪里給出了通俗和準確的解釋。

1.2  如何減輕鎖死

可以通過單元技術來減輕鎖死的影響。如Selective reduce integration, Incompatible，Assumed strain，Mixed formulation，MITC等等。對于近似不可壓縮材料的體積鎖死，最常見的對應方法是采用B-bar單元，這種單元的新式變種是F-bar單元，這是本文下面要講的內容。

2. B-bar單元[3,4]

一般的，單元內任意一點的應變ε表為單元節點i的位移u的函數。在B-bar單元中，將應變分解，其中為單元的平均體積應變或單元中心的體積應變值為偏差應變值。相應地，單元內任意一點的應變用下式計算

這里將體積應變做平均計算的方法緩和了體積鎖死。這種方法實現起來非常簡單，計算費用不高，效果也不錯，因此得到非常廣泛的應用。但是這種方法視乎不能滿足LBB條件[5]。

3. F-bar單元[6,7]

相應于B-bar單元對應變的加法分解，單元中的任意一點的變形梯度（deformation gradient）可以分解為體積變形梯度和偏差變形梯度兩個部分 F=Fs Fv。 理論上來說這種分解方式是嚴密的，而B-bar單元中的應變加法分解是一種近似，只是在應變很小時與變形梯度的乘法分解相近。因此把B-bar單元推廣到 F-bar單元是非常自然的，在變形較大時采用F-bar單元精度更好。

在F-bar單元中，單元中的任意一點的變形梯度表述為，J=detF. 其中

可以是單元中Jocobian J的平均值 ，也可以指單元中心的Jocobian值。

3.1  Total Lagrange算式展開

在Total Lagrange計算中，單元中變形能由Green應變E和第2 Piola-Kirchhoff 應力S的乘積在初始構型V中的積分計算得出，其中中的采用上述變形梯度表達式。如此，變形能變分為，單元的consistency stiff matrix可以由變形能變分的微分展開得出。

3.2 Updated Lagrange算式展開

在Updated Lagrange計算中，單元中變形能由Almansi應變A和Cauchy應力T的乘積在現在構型v中的積分計算得出。

, 其中。把這些變換帶入3.1的相應算式即可得出Updated Lagrange法的內力，單元的consistency stiff matrix計算公式。

3.3 微小變形計算的算式展開

如前所述，F-bar單元用于大變形計算。在微小變形計算中應用B-bar單元就可以了。不過由于應變的加法分解可以理解為變形梯度乘法分解的近似。那么可以從F-bar單元公式推出B-bar公式來。雖然從實用的角度來看這樣做似乎意義不大。參見文獻[9]。

3.4 F-Bar-patch 單元[8]

F-Bar單元不能應用于三角形一次和四面體一次單元，因為這些單元內的應變本身就是均勻分布的。為了將這一方法擴展到三角形一次和四面體一次單元，文獻[8]提出了將單元分解為相鄰的patch的方法，稱為F-Bar-patch單元。

4. 算例: Cook's membrane [10,11]

Cook's membrane benchmark見上圖，鍥形版材質設定為近似不可壓縮的Mooney Rivlin超彈性材: 其·彈性能函數, 給定參數C10=1.0，C01=0.925，D1=0.000245。 采用單元為六面體F-Bar單元，B-Bar單元和一般的等參六面體單元。計算結果如下：

（a）等參單元：受力方向最大位移值=0.2706

（b）B-bar單元：受力方向最大位移值=3.898

（c）F-Bar單元：受力方向最大位移值=4.918

  

通過數值的實例計算，可以明顯的看出，B-Bar單元在超彈材料的計算上要明顯優于等參數單元，而F-Bar單元，在變形較大時是優于B-Bar單元的。值得一提的，由中國人主導研發的大型通用有限元軟件 WELSIM 已經支持了F-Bar單元，這為我們在工程分析中得到更加精確的結果提供了有力工具。

---

參考文獻

1）Allan F. Bower: Applied Mechanics of Solids, CRC, 2010, p535

2） Gangan Prathap: FINITE ELEMENT ANALYSIS AS COMPUTATION， http://www.cmmacs.ernet.in/cmmacs/pdf/ch08.pdf

3）TJR Hughes: Centralization of selective integration procedure to anisotropic and nonlinear media, Intl. J. Numer. Methods Engng., 15(1980), 1413-1418

4) JC Simo, KS Pister, RA Taylor: Variational and projection methods for the volume constraint in finite deformation elastoplasticity, Comp. Methods Appl. Mech. Eng., 51(1985), 177-208

5) F Brezzi, M Fortin: Mixed and Hybrid Finite Element Method, Springer, 1991

6)  B Moran, M Oritz, CF Shin: Formulation of implicit finite-element methods for multiplicative finite deformation plasticity, Intl. J. Numer. Methods Engng., 29(1990), 483-514

7)  EA de Souza Neto, D Peric, M Dutko, DRJ Owen: F-bar based linear triangles and tetrahedral for finite element strain analysis of nearly incompressible solids. Part I: Formulation and benchmarking Intl. J. Numer. Methods Engng., 33(1996), 3277-3296

8)  EA de Souza Neto, D Peric, M Dutko, DRJ Owen: Design of simple low order finite element for large strain analysis of nearly incompressible solids, Intl. J. Numer. Methods Engng., 62(2005), 353-383

9）CD Foster, T Mohammad Nejad，Trilinwear Hexahedra with integral-averaged volumes for nearly incompressible nonlinear deformation, Engineering, 7(2015), 765-788

10)  Cook, R. D.; Improved two-dimensional finite element; ASCE J. Struct. Div., ST9, 1851-1863 (1974).

11) Brink,  U., and E. Stein, “On some Mixed Finite Element Methods for Incompressible and Nearly Incompressible Finite Elasticity,” Computational Mechanics, vol. 19, pp. 105–119, 1996.