基于Total-FETI法的混凝土損傷分析子區域剖分研究

邱莉婷 馬福恒 沈心哲

摘要：在整體有限元撕裂對接法和隱式梯度損傷模型結合的基礎上，針對損傷子區域模擬的高精度要求和混凝土損傷擴展路徑的不確定性問題，通過引入基于非局部損傷應變的子區域自適應更新方法，結合子區域更新誤判修正，實現了混凝土損傷失效過程非線性子區域的高精度有限元網格自適應更新。同時，對L型混凝土試件進行不同子區域剖分數量和子區域分解形式的數值試驗對比分析。研究結果表明：模型不存在子區域分解敏感性;子區域剖分數量越多，或者子區域分解方式與損傷擴展路徑相適應時，高精度有限元網格子區域的自適應更新數量越少，模型計算規模減小明顯;此外，子區域界面節點的增加對計算規模削減的影響較小，模型整體計算效率提高明顯。研究結果可為混凝土損傷分析的并發多尺度數值計算提供參考。

關 鍵 詞：混凝土損傷; 隱式梯度損傷模型; 整體有限元撕裂對接法; 子區域自適應更新; 子區域分解形式; 計算效率

中圖法分類號： TU43

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2022.04.028

0 引 言

混凝土結構服役過程中在結構突變區域容易出現局部應力集中現象，其工程結構的失效往往源于局部構件的細觀缺陷和損傷局部化行為[1-2]。進行混凝土結構的穩定性和強度分析時不僅要把握結構整體力學響應，同時也應該重視混凝土的損傷萌生、擴展和失效過程[3]。損傷力學作為模擬連續介質逐步劣化過程的有效分析方法，可反映結構從完整連續介質結構體到局部損傷萌生直至整體失效的全過程[4]。隨著損傷的不斷積累，混凝土因其準脆性特性將出現應變局部化和應變軟化行為[5]，為保證計算精度，該區域有限元網格的精細化程度也需要不斷提高。混凝土作為典型的準脆性材料[6]，其細觀結構由骨料、硬化水泥砂漿和二者間的界面黏結帶組成[7-8]。混凝土的斷裂破壞是細觀層次上損傷積累與宏觀尺度上裂紋失穩擴展交織發展的復雜過程[9]。混凝土結構的損傷擴展路徑具有不確定性[10]，數值計算過程需要不斷地進行網格自適應更新，以滿足對損傷局部化區域力學特性的準確把握。

本文將整體有限元撕裂對接法和隱式梯度損傷模型相結合，在損傷子區域高精度有限元網格自適應更新的基礎上，開展混凝土損傷失效過程自適應、多尺度區域分解模擬，并開展子區域分解影響研究，為模型混凝土損傷區域的高精度模擬和損傷擴展路徑的不確定性問題研究提供解決途徑。

1 混凝土損傷分析的Total-FETI法

1.1 基于隱式梯度損傷模型的Total-FETI法

隱式梯度損傷模型屬于特殊的非局部損傷模型[11-12]，其將非局部模型中的權函數替換成格林（Green）函數[13]，該模型結合了空間相互作用與梯度公式的計算效率，操作更為簡單，也稱修正的Helmholtz’s方程：

εeq-clSymbolQC@2εeq=ε（1）

式中：非局部等效應變εeq不是由等效局部應變ε和它的導數顯式描述的，而是作為包含方程（1）的邊值問題和近似邊界條件的解[14];梯度參數為內部長度參數l的函數。

有限元撕裂對接法（finite element tearing and interconnecting method，FETI）屬于采用局部邊界條件的不重疊區域分解法[15-16]，其并行計算的收斂速度與子區域的數量相互獨立，可以將原求解問題劃分成更多的子區域進行求解，且可以保證較高的求解效率[17]。整體有限元撕裂對接法（Total-FETI）在有限元撕裂對接法基礎上將拉格朗日算子同時用于子區域界面連接和狄利克雷邊界條件的施加[18-19]，可簡化子區域剛度矩陣的求逆過程。非線性有限元撕裂對接法的線性方程系統和子區域Ω（s）間的位移連續條件可表述為

A（s）kδu（s）k+1+B（s）Tδλk+1=f（s）ext-B（s）Tλk-f（s）int，k（u（s）k）（2）

Nss=1B（s）u（s）k+Nss=1B（s）δu（s）k+1=0（3）

式中：A（s）是切線剛度矩陣;δu（s）和δλ（s）分別為位移增量和拉格朗日乘子增量;布爾矩陣B（s）由子區域界面的布爾矩陣和狄利克雷邊界條件對應的布爾矩陣按行串連構成;f（s）ext為外力向量;f（s）int為內力向量;Ns為子區域數量。

隱式梯度損傷模型的非局部等效應變εeq在Total-FETI中可通過對給定子區域界面的拉格朗日算子進行修正實現：

λd=λλεeq（4）

式中：λεeq為非局部等效應變對應的自由度。

再通過擴展布爾矩陣B（s）將λd裝配到相應位置，布爾矩陣Bd （s）由原來Total-FETI的子區域界面的布爾矩陣B（s）和非局部等效應變對應的布爾矩陣B（s）eq按行串連構成：

B（s）d=B（s）B（s）eq（5）

1.2 FETI的局部子區域自適應更新

本文混凝土結構的損傷分析采用隱式梯度損傷模型。對于各個非線性子區域的實時判別，可以根據特定加載步的各子區域的非局部等效應變最大值ε（s），maxeq與損傷發生時的應變閾值k0的相對關系進行。僅當子區域的非局部等效應變最大值ε（s），maxeq大于等于應變閾值k0時，損傷才會發生。這里，對子區域Ω（s）遍歷其所有節點n得到其非局部等效應變最大值ε（s），maxeq，t和非局部等效應變最小值ε（s），mineq，t：

ε（s），maxeq，t=max（εeqi），i=1，2，…，n（6）

ε（s），mineq，t=min（εeqi），i=1，2，…，n（7）

記加載步t和加載步t+1所對應的非局部等效應變差值的預測值為Δε（s）t，則加載步t+1的非局部等效應變最大值ε（s），peq，t+1的實時預測值為[20]

Δε（s）t=max（ε（s），maxeq，t-ε（s），mineq，t-1，ε（s），mineq，t-ε（s），maxeq，t-1）（8）

ε（s），peq，t+1=ε（s），maxeq，t+Δε（s）t（9）

引入基于層級多尺度的子區域自適應更新方法[21]進行非線性子區域實時更新。首先，通過邊值問題求解以獲取被替換子區域的最近變形;然后，進行整體區域的重新平衡迭代以消除由于子區域更新導致的殘余應力;進一步地，采用尺度間的線性多點約束（interscale linear multi-point constraints，ILMPCs）[21-22]進行子區域更新后的非協調界面連接，使得不同精細度有限元網格子區域可在同一個有限元計算模型中共存。

對于兩個網格精細度不一致的子區域，其尺度間的線性多點約束用矩陣形式可表示如下[23]：

Pu=[P（1） P（2）]u（1）u（2）=0 （10）

式中：u（s）為子區域Ω（s）的位移向量;矩陣P（s）由線性多點約束組成，建立了子區域更新后新增界面節點與鄰接子區域節點間的聯系，即高精度有限單元的節點自由度通過鄰接子區域交界面處單元的形函數插值求得。

這一系列尺度間的線性多點約束施加可以通過在Total-FETI中用拉格朗日乘子添加額外數量的方程完成。通過定義一個修正的布爾矩陣B（s）來實現。修正后布爾矩陣B（s）由B（s）d和約束矩陣P（s）按行串連組成：

B（1）B（2）=B（s）dB（s）dP（1）P（2）（11）

擴展后的拉格朗日乘子U包含鄰接子區域間節點的界面約束λd以及與更新子區域新增節點間的約束η：

U=λdη（12）

基于修正后的布爾矩陣B（s），非線性整體有限元撕裂對接法方程（2）和方程（3）可以改寫為

A（s）kδu（s）k+1+B（s）TδUk+1=f（s）ext-B（s）TUk-f（s）int，k（u（s）k）（13）

Nss=1B（s）u（s）k+Nss=1B（s）δu（s）k+1=0（14）

1.3 非線性子區域誤判修正

為防止出現某個非線性子區域的誤判，在加載步t+1迭代收斂后，每個子區域Ω（s）的非局部等效應變預測值ε（s），peq和非局部等效應變計算值ε（s），maxeq都將與損傷發生的非局部等效應變閾值k（s）0進行對比分析：

Δε（s），peq，t+1=ε（s），peq-k（s）0（15）

Δε（s）eq，t+1=ε（s），maxeq-k（s）0（16）

一旦出現Δε（s），peq，t+1<0且Δε（s）eq，t+1≥0的情況，則表明子區域Ω（s）的非線性出現誤判，此時應該將子區域Ω（s）判定為非線性子區域，再返回上一個加載步t進行重新迭代計算，以修正加載步t+1的求解。

2 L型混凝土試件損傷分析的數值試驗

2.1 Total-FETI模型及計算參數

為分析子區域分解方式對模型計算結果和計算效率的影響，對圖1所示的三維L型混凝土試件采用3種不同的子區域分解方式進行計算對比。試件幾何尺寸如圖1所示，邊界條件為在x=100 mm的R面上施加約束，ux=0 mm和uz=0 mm;在y=0 mm的B面上施加約束uy=0 mm和uz=0 mm;在y=100 mm的頂面施加ux=-0.6 mm和uz=-0.6 mm的位移約束。計算采用隱式梯度損傷模型，其參數說明如下：采用修正的von Mises模型計算等效應變，其混凝土的壓縮和拉伸強度比值η為10;采用Mazars模型[24]中的損傷演化法則，其控制軟化速率的β取值為50，控制殘余應力的α取值為0.999，損傷發生時的應變閾值k0取值為0.000 5;彈性模量E取值為40 000 MPa，泊松比υ為0.2，梯度參數c=1 mm2。將各個子區域的剛度方程和描述各子區域界面和邊界條件的布爾矩陣直接組裝成有限元支配方程的整體系數矩陣。線性方程組系統采用直接求解法進行求解，求解器為并行直接求解器Pardiso[25]。同時，采用1.3節所述非線性子區域誤判修正對每個加載步非線性判斷進行復核。

多尺度計算采用粗網格和細網格兩套有限元模型，其中細網格模型是在粗網格模型上進行細分得到，如圖2所示，均采用八節點六面體單元進行網格剖分。粗網格模型圖2（a）的節點總數為3 000，單元總數1 536;細網格模型圖2（b）節點總數117 912，單元總數98 304。

為對比三維情況下不同子區域剖分方式對子區域自適應更新、計算內存和計算時間的影響，將L型混凝土試件分別分解為如圖3所示的24個正方體子區域試件，48個長方體子區域試件以及192個正方體子區域試件。需要指出的是，圖3（b）所示試件缺口兩側的長方體子區域分別采用水平向布置和豎直向布置，這與其他區域均采用水平布置有所不同。

2.2 試驗結果分析

各試件在其y=100 mm試件頂面位置的荷載-ux曲線如圖4所示。可以發現，3種不同子區域分解方式下試件的荷載-ux曲線保持基本一致，子區域分解數量增加時，試件剛度呈現微小幅度地增大。接下來，通過選取曲線上的各試件的第一次子區域網格自適應更新的A點、峰值荷載所對應的B點和試件基本完全破壞所對應的C點的各子區域網格自適應更新情況進行對比分析。

如圖5所示，3個試件的第一次基于子區域的自適應網格更新數量不一樣。直觀來看在這個階段呈現的規律是子區域分解數量越多，則自適應更新的網格數量越少。定量分析表明：其中該階段L-24的網格總數為13 632，為全局細觀網格總數的13.87%，內存占用量2 286 MB;L-48的網格總數為7 584，為全局細觀網格總數的7.71%，內存占用量1 542 MB;L-192的網格總數為3 048，為全局細觀網格總數的3.10%，內存占用量1 785 MB。定量分析結論和直觀規律基本一致，但是內存占用量反而是L-192較L-48增大了15.8%，可見在此階段大幅增加子區域分解數量會使得界面自由度加大，削減了其減少自適應網格更新數量帶來的效率。

如圖6所示為3個不同子區域剖分形式試件在峰值荷載階段的自適應網格更新數量。直觀來看仍是隨子區域分解數量加大自適應更新的網格數量減少。定量分析表明：該階段L-24的網格總數為45 888，為全局細觀網格總數的46.68%，內存占用量7 358 MB;L-48的網格總數為33 792，為全局細觀網格總數的34.38%，內存占用量5 728 MB;L-192的網格總數為20 184，為全局細觀網格總數的20.53%，內存占用量4 607 MB。定量分析結論和直觀規律一致，同時內存占用量也呈現隨子區域剖分數量增加而遞減的情況，說明子區域界面自由度增加所導致的計算效率下降影響在子區域數量更新中期已經不再明顯。

如圖7所示，3個不同子區域分解形式試件在基本完全損傷階段的自適應網格更新總數量與前述階段的規律一致，即均隨子區域數量遞增而減少。定量分析表明：L-24的網格總數為57 984，為全局細觀網格總數的58.98%，內存占用量9 271 MB;L-48的網格總數為43 872，為全局細觀網格總數的44.63%，內存占用量7 338 MB;L-192的網格總數為34 800，為全局細觀網格總數的35.40%，內存占用量7 014 MB。定量分析結論和直觀規律一致，子區域分解數量的增加可明顯減小計算網格模型規模以及內存使用量，子區域界面自由度增加對計算效率降低的影響較小。同時，圖7（b）可解釋為何將L-48試件缺口兩側子區域分別設置為水平向和豎直向布置，因為本文自適應網格更新是基于子區域網格更新進行的，根據損傷可能的擴展形式將兩側子區域設置為不同方向可使得損傷影響的子區域盡量少。再從計算時間來看，L-24計算時間為13.79 h，L-48計算時間為10.68 h，L-192計算時間為8.72 h。可知L-48和L-192的計算時間相差較小，進一步說明了合理劃分子區域形式也是提高計算效率的有效手段。

綜上可知，不同子區域分解方式的損傷擴展模式以及荷載-ux曲線均保持良好的一致性，子區域的不同形式以及不同數量的分解對計算結果影響不大，模型不存在子區域分解敏感性。對比分析不同子區域分解數量的網格更新結果和計算內存需求可知，子區域網格分解數量越多所需要的計算網格數量越小，計算所需內存也呈現同樣規律，且此優勢在子區域自適應更新的中后期更為明顯。因為從網格增長速率來看，L-24和L-48網格在后期增長較慢，說明自適應更新中期精細網格已經基本上被替換，整體計算效率要自然比L-192網格更新隨損傷發展均勻替換要低。

3 結 論

本文采用隱式梯度損傷模型作為混凝土損傷分析的本構模型，并采用Total-FETI方法進行混凝土損傷失效過程的自適應并發多尺度區域分解模擬。針對混凝土損傷子區域的高精度模擬要求和損傷擴展路徑的不確定性問題，通過L型混凝土試件數值試驗開展了子區域分解對模型規模及計算效率的影響研究，主要得出以下結論。

（1） 對混凝土試件采用不同的子區域分解形式和分解數量進行對比分析，計算所得損傷模式一致，可知模型無區域分解敏感性。

（2） 由于模型各子區域均采用粗網格，僅局部非線性子區域進行高精度細觀網格替換，界面節點數量增加對計算規模和計算效率影響較小。

（3） 子區域數量分解越多，自適應更新的高精度子區域越少，有限元模型網格數量減少，從而計算效率提高。

（4） 子區域分解方式與損傷可能擴展路徑相適應，損傷區域涉及的高精度子區域越少，模型耗費計算資源削減，計算效率改進。

綜上可知，子區域分解方式和子區域分解數量對自適應網格更新總數及計算所需內存影響顯著，合理的子區域剖分方式和適當的子區域分解數量可以有效降低計算規模和計算內存需求。

參考文獻：

[1] 吳佰建，李兆霞，湯可可.大型土木結構多尺度模擬與損傷分析：從材料多尺度力學到結構多尺度力學[J].力學進展，2007，37（3）：321-336.

[2] 胡少偉，陸俊，范向前.混凝土損傷斷裂性能試驗研究進展[J].水利學報，2014，45（增1）：10-18.

[3] 林皋，劉軍，胡志強.混凝土損傷類本構關系研究現狀與進展[J].大連理工大學學報，2010，50（6）：1055-1064.

[4] 張立，李廣凱，馬懷發.混凝土類材料彈塑性損傷問題的全隱式迭代法[J].水利學報，2020，51（8）：947-956.

[5] 汪忠明，牛曉玉，楊伯源.基于隱式梯度公式的各向異性混凝土損傷研究[J].工程力學，2011，28（2）：159-164.

[6] 任宇東，陳建兵.基于一類非局部宏-微觀損傷模型的混凝土典型試件力學行為模擬[J].力學學報，2021，53（4）：1196-1211.

[7] 肖詩云，朱梁.混凝土初始損傷細觀結構特性數值試驗研究[J].大連理工大學學報，2017，57（1）：78-86.

[8] 徐青，周祥森，程志誠.基于Ansys的混凝土隨機骨料模型及細觀力學分析[J].武漢大學學報（工學版），2019，52（12）：1035-1040，1047.

[9] 王康，吳佰建，李兆霞.損傷跨尺度演化導致的混凝土強度尺寸效應[J].東南大學學報（自然科學版），2014，44（6）：1230-1234.

[10] 程道輝，王苗，卿龍邦，等.混凝土Ⅰ型裂縫隨機損傷斷裂分析研究[J].混凝土，2020（11）：40-42，56.

[11] PEERLINGS R R.Enhanced damage modelling for fracture and fatigue[D].Eindhoven：Technische Universiteit Eindhoven，1999.

[12] SIMONE A，WELLS G N，SLUYS L J.From continuous to discontinuous failure in a gradient-enhanced continuum damage model[J].Computer Methods in Applied Mechanics and Engineering，2003，192（41-42）：4581-4607.

[13] 汪忠明，牛曉玉，楊伯源.基于非局部隱式梯度模型的混凝土斷裂損傷[J].中國科學技術大學學報，2010，40（6）：651-656.

[14] 汪忠明.基于非局部增強梯度模型的混凝土斷裂損傷研究[D].合肥：合肥工業大學，2010.

[15] FARHAT C，ROUX F X.A method of finite element tearing and interconnecting and its parallel solution algorithm[J].International Journal for Numerical Methods in Engineering，1991，32（6）：1205-1227.

[16] 張云鵬，楊新生，邵定國，等.基于自適應區域分解的電磁場有限元求解方法研究[J/OL].高電壓技術.https：∥doi.org/10.133361/j.1003-6520.hev.20201807.

[17] 張曉虎，孫秦.結構非匹配網格區域分裂的并行數值計算研究[J].西北工業大學學報，2019，37（4）：650-655.

[18] DOSTL Z，HORK D，KUERA R.Total FETI—an easier implementable variant of the FETI method for numerical solution of elliptic PDE[J].International Journal for Numerical Methods in Biomedical Engineering，2006，22（12）：1155-1162.

[19] ERMK M，HAPLA V，HORK D，et al.Total-FETI domain decomposition method for solution of elasto-plastic problems[J].Advances in Engineering Software，2015，84.

[20] LLOBERAS-VALLS O，RIXEN D J，SIMONE A，et al.Domain decomposition techniques for the efficient modeling of brittle heterogeneous materials[J].Computer Methods in Applied Mechanics and Engineering，2011，200（13-16）：1577-1590.

[21] LLOBERAS-VALLS O，RIXEN D J，SIMONE A，et al.Multiscale domain decomposition analysis of quasi-brittle heterogeneous materials[J].International Journal for Numerical Methods in Engineering，2012，89（11）：1337-1366.

[22] HAIKAL G，HJELMSTAD K D.An enriched discontinuous Galerkin formulation for the coupling of non-conforming meshes[J].Finite Elements in Analysis & Design，2009，46（6）：496-503.

[23] LLOBERAS-VALLS O，RIXEN D J，SIMONE A，et al.On micro-to-macro connections in domain decomposition multiscale methods[J].Computer Methods in Applied Mechanics and Engineering，2012，225-228：177-196.

[24] 柯楊，姜勝濤，李佳，等.不同骨料體積率的輕混凝土損傷-斷裂分析[J].混凝土，2020（7）：15-19，24.

[25] 曹大志，強洪夫，任革學.基于OpenMP和Pardiso的柔性多體系統動力學并行計算[J].清華大學學報（自然科學版），2012，52（11）：1643-1649.

（編輯：鄭 毅）

Research on sub-domain meshing of concrete damage analysis based on Total-FETI method

QIU Liting，MA Fuheng，SHEN Xinzhe

（Nanjing Hydraulic Research Institute，Nanjing 210029，China）

Abstract：

Based on combination of total finite element tearing-interconnecting method and gradient-enhanced continuum damage model，aiming at the high-precision requirements of the damage sub-domain simulation and the uncertainty of concrete damage propagation path，through introducing the sub-domain adaptive updating method based on non-local equivalent strain，and incorporating the sub-domain update misjudgment correction，the high-precision finite element mesh adaptive update of nonlinear sub-domain during concrete damage and failure process simulation were realized.Moreover，numerical tests having different domain meshing numbers and domain meshing forms were carried out on L-shaped concrete specimen.The results show that there was no sub-domain meshing sensitivity for the model.When the number of sub-domain subdivision increased or the meshing mode of sub-regions was adapted to the damage propagation path，the number of adaptive updating of the sub-domain with high-precision finite element mesh decreased，and the calculation scale of the model was significantly reduced.In addition，the increment of sub-regional interface nodes had little effect on the reduction of calculation scale，and the overall calculation efficiency of the model was significantly improved.This study can provide reference for concurrent multi-scale numerical calculation of concrete damage analysis.

Key words：

concrete damage;gradient-enhanced continuum damage model;total finite element tearing-interconnecting method;sub-domain adaptive updating;sub-domain decomposition pattern;computational efficiency