邱莉婷，馬福恒，沈振中，張 湛，霍吉祥

(1.南京水利科學(xué)研究院，江蘇 南京 210029； 2.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098； 3.河南省西霞院水利樞紐輸水及灌區(qū)工程建設(shè)管理局，河南 鄭州 450008)

混凝土作為由硬化水泥砂漿、骨料及二者界面過渡區(qū)組成的非均質(zhì)三相復(fù)合材料，在澆筑過程中由于溫度應(yīng)力和收縮問題會形成大量的孔隙和原始微裂紋[1]。其斷裂是在外界荷載作用下，局部微缺陷發(fā)展成核、擴展匯聚，宏觀裂紋形成，裂紋失穩(wěn)擴展的交織發(fā)展過程。整個斷裂過程無大比例塑性流動產(chǎn)生，僅吸收了裂縫表面形成過程中的剩余能量，是典型的準(zhǔn)脆性斷裂[2]。準(zhǔn)脆性混凝土材料由于細觀結(jié)構(gòu)的多相性和不均勻性，其宏觀力學(xué)響應(yīng)十分復(fù)雜，表現(xiàn)為非線性、應(yīng)變軟化和應(yīng)變局部化等材料特性。其非線性力學(xué)特性主要有彈塑性、黏性、損傷和斷裂等幾種，損傷占有重要地位。損傷力學(xué)是研究材料的微、細觀缺陷及裂紋萌生和發(fā)展對其宏觀力學(xué)特性的影響，以及材料和結(jié)構(gòu)損傷演化過程及機理的破壞力學(xué)。當(dāng)有限元計算采用局部損傷模型時會存在網(wǎng)格敏感性和零能耗問題[3]。隱式梯度損傷模型[4]屬于非局部損傷模型，以變形的梯度形式反映微結(jié)構(gòu)變化和相互作用，把空間相互作用與梯度公式的計算效率相結(jié)合，適用于同非線性有限元法結(jié)合開展混凝土損傷失效的數(shù)值模擬。

同時，傳統(tǒng)單一網(wǎng)格的有限元模型無法對混凝土線彈性區(qū)域和損傷集中區(qū)域進行區(qū)別處理，不能將有限的計算資源用于需要重點分析的局部易損區(qū)域。而采用過渡單元進行線彈性區(qū)域粗網(wǎng)格和非線性區(qū)域細網(wǎng)格的連接，當(dāng)網(wǎng)格尺寸跨度較大時，過渡單元的網(wǎng)格剖分質(zhì)量難以保證。如進行自適應(yīng)有限元分析[5]往往需要進行多次自適應(yīng)網(wǎng)格剖分才能達到要求的計算精度。對此，并發(fā)多尺度方法[6]可以使高效率的粗網(wǎng)格和高精度的細網(wǎng)格在同一有限元模型中共存，且對混凝土的應(yīng)變軟化行為不存在層級多尺度方法[7]所面臨的尺寸依賴性問題，可用于混凝土的損傷失效全過程分析。但計算過程中直接包含細網(wǎng)格子區(qū)域仍會使計算資源消耗較大。整體有限元撕裂對接法（Total Finite element tearing and interconnecting method,Total-FETI）[8]屬于局部邊界條件不重合的區(qū)域分解法，基于分而治之思想，具有高度的并行性和并行效率，適合于大規(guī)模數(shù)值計算。

本文采用隱式梯度損傷模型描述混凝土材料的非線性本構(gòu)關(guān)系，將并發(fā)多尺度方法和整體有限元撕裂對接法結(jié)合，在兼顧計算精度和計算效率的基礎(chǔ)上，構(gòu)建混凝土損傷失效分析的并發(fā)多尺度區(qū)域分解模型，并將該模型應(yīng)用于混凝土預(yù)制切口四點彎曲梁試驗的數(shù)值模擬。

1 混凝土單邊切口梁四點彎曲試驗數(shù)值模型

1.1 隱式梯度損傷模型在Total-FETI中的實現(xiàn)

有限元撕裂對接法[9]將所研究區(qū)域分解為多個不重疊的獨立子區(qū)域，通過拉格朗日乘子去除各子區(qū)域界面間的連續(xù)性約束。其求解由兩部分組成：首先是子區(qū)域界面問題求解，即采用預(yù)條件共軛梯度算法（Preconditioned conjugate gradiant，PCG）求解僅含界面未知量的并由狄利克雷邊界條件約束的縮減系統(tǒng)；然后再進行由紐曼邊界條件約束的各個獨立子區(qū)域問題求解。對于界面問題，在PCG迭代求解過程將誤差進行全局化處理，從而加快計算收斂速度。同時，相比傳統(tǒng)的子結(jié)構(gòu)法，F(xiàn)ETI的子區(qū)域劃分不需要人為干預(yù)，也不要求滿足子結(jié)構(gòu)內(nèi)部節(jié)點先編號、界面節(jié)點后編號的編號原則，無需再進行節(jié)點編號排序優(yōu)化，且子區(qū)域間通信規(guī)模也較小。FETI法的內(nèi)在機制使得計算時的收斂快慢和子域劃分的數(shù)量多少無直接聯(lián)系，可以根據(jù)計算需求靈活地進行子區(qū)域劃定。進一步地，整體有限元撕裂對接法（Total-FETI）[8]在有限元撕裂對接法基礎(chǔ)上，將拉格朗日算子不僅用于子區(qū)域交界面的約束，同時也用于邊界子區(qū)域的狄利克雷邊界條件施加，從而達到簡化子區(qū)域剛度矩陣的求逆過程，這也是本文多尺度區(qū)域分解算法構(gòu)建的基礎(chǔ)。

本文主要研究混凝土損傷失效的非線性演變過程，因此，這里給出非線性Total-FETI離散后的牛頓拉弗森迭代（Newton-raphson）求解公式。假定在某加載步的第k+1牛頓拉弗森迭代步，其位移向量和拉格朗日乘子的更新如下：

式中：u為位移向量；拉格朗日算子λ同時考慮了子區(qū)域界面位移自由度和狄利克雷邊界條件的施加。

非線性有限元撕裂對接法的線性方程系統(tǒng)和子區(qū)域Ω(s)間的位移連續(xù)條件可表述為：

式中：A(s)為切線剛度矩陣；δu(s)和δλ分別為位移增量和拉格朗日乘子增量；布爾矩陣B(s)由子區(qū)域界面的布爾矩陣和狄利克雷邊界條件對應(yīng)的布爾矩陣按行串連組成；fext(s)為外力向量；fint(s)為內(nèi)力向量；Ns為子區(qū)域數(shù)量。

隱式梯度損傷模型屬于特殊的非局部損傷模型，其求解非局部等效應(yīng)變的隱式梯度方程將非局部模型中的權(quán)函數(shù)替換成格林（Green）函數(shù)，也稱修正的Helmholtz’s方程[10]。

式中：非局部等效應(yīng)變?eq不是由等效局部應(yīng)變?及其導(dǎo)數(shù)顯式描述，而是作為包含方程（4）的邊值問題和近似邊界條件的解；梯度參數(shù)為內(nèi)部長度參數(shù)l的函數(shù)，其取值為c(l)=l2/2。

Total-FETI由于考慮了子區(qū)域界面相容條件，其控制方程的變分形式以混合變分形式表述。Total-FETI與隱式梯度損傷模型結(jié)合后，由于子區(qū)域界面的非局部等效應(yīng)變?eq的存在，需要對給定子區(qū)域界面的拉格朗日算子進行如下修正：

式中： λεeq考慮了非局部等效應(yīng)變的自由度。

通過擴展布爾矩陣B(s)將λd裝配到合適位置，布爾矩陣Bd(s)由原來Total-FETI的子區(qū)域界面的布爾矩陣B(s)和非局部等效應(yīng)變對應(yīng)的布爾矩陣Beq(s)按行串連組成：

1.2 并發(fā)多尺度方法與Total-FETI的結(jié)合

細網(wǎng)格子區(qū)域和粗網(wǎng)格子區(qū)域在其連接界面需滿足位移連續(xù)條件。細網(wǎng)格子區(qū)域界面節(jié)點位移可以通過粗網(wǎng)格子區(qū)域界面單元的形函數(shù)和界面節(jié)點位移進行插值求解。這里將粗網(wǎng)格子區(qū)域界面節(jié)點稱為主節(jié)點，將細網(wǎng)格子區(qū)域界面未能與主節(jié)點匹配的其余節(jié)點稱為從節(jié)點。從節(jié)點的自由度由線性多點約束以不同的插值方式實現(xiàn)與主節(jié)點的耦合。其中主節(jié)點間的約束為線性多點約束，從節(jié)點和主節(jié)點之間的約束稱為尺度間的線性多點約束。對于兩個有限元網(wǎng)格尺寸不一致的子區(qū)域界面，其尺度間的線性多點約束用矩陣形式[11]可表示如下：

式中：u(s)為子區(qū)域 Ω(s)的位移向量；矩陣P(s)由線性多點約束組成，建立了從節(jié)點和主節(jié)點間的聯(lián)系，即細網(wǎng)格節(jié)點自由度通過粗網(wǎng)格子區(qū)域交界面處單元的形函數(shù)插值求得。

尺度間的線性多點約束可通過拉格朗日乘子法實現(xiàn)[12]。采用拉格朗日乘子在Total-FETI中添加額外數(shù)量方程，通過定義一個修正的布爾矩陣來實現(xiàn)，修正后布爾矩陣由Bd(s)和約束矩陣P(s)按行串連組成：

擴展后的拉格朗日乘子U包含鄰域間主節(jié)點的界面約束λd及主節(jié)點與從節(jié)點間的約束η：

在引入隱式梯度損傷模型后，Total-FETI的界面問題求解需要重新修正。隱式梯度損傷模型在每個節(jié)點需要額外的自由度來描述其非局部等效應(yīng)變，所以子區(qū)域的剛度矩陣是非對稱矩陣。從而柔度矩陣和子區(qū)域界面問題的拉格朗日矩陣均為非對稱矩陣，對其進行迭代求解時只有雙共軛梯度穩(wěn)定解法（BiCGStab）和廣義極小剩余法（GMRES）兩種迭代求解器可滿足計算要求?？紤]到子區(qū)域界面方程的迭代求解難以構(gòu)建合適預(yù)處理算子，在整體有限元撕裂對接法求解框架下，不妨回到有限元撕裂對接法的最初形式：

式（12）可以通過串聯(lián)不同子區(qū)域的矩陣和向量改寫成熟悉的有限元支配方程形式Au=f。方程組采用分塊矩陣的形式表述如下：

將各個子區(qū)域的剛度方程和描述各子區(qū)域界面和邊界條件的布爾矩陣直接組裝成有限元支配方程的整體系數(shù)矩陣，稱為雙重組裝法[13]。再采用非對稱稀疏線性方程組的并行直接求解器Pardiso[14]進行大型線性方程組的求解。采用雙重組裝法的整體有限元撕裂對接法不再需要進行子區(qū)域的界面問題和剛體位移求解。

1.3 并發(fā)多尺度區(qū)域分解算法建立流程

混凝土損傷分析的并發(fā)多尺度算法建立流程主要包括子區(qū)域分解及多尺度有限元網(wǎng)格剖分，子區(qū)域界面信息和狄利克雷邊界信息提取，布爾矩陣及其轉(zhuǎn)置矩陣建立及各加載步的牛頓拉弗森迭代求解。并發(fā)多尺度有限元網(wǎng)格剖分，首先按照構(gòu)件幾何特征和重點分析區(qū)域進行子區(qū)域分解，其子區(qū)域的幾何形狀可以是不規(guī)則的；再對全部子區(qū)域進行粗網(wǎng)格剖分，粗網(wǎng)格的剖分要求能較好地反映構(gòu)件幾何及子區(qū)域邊界形狀，也要能準(zhǔn)確捕捉構(gòu)件線彈性階段位移場的空間分布和變形梯度變化；重點分析區(qū)域在保留粗網(wǎng)格子區(qū)域邊界節(jié)點的基礎(chǔ)上進行指定精度的細網(wǎng)格剖分，細網(wǎng)格子區(qū)域的內(nèi)部可以按照計算精度要求進行不同類型的有限元網(wǎng)格剖分，不受原粗網(wǎng)格子區(qū)域內(nèi)部節(jié)點限制。

2 計算分析

為驗證模型的合理性，在前人研究[15-16]的基礎(chǔ)上，對圖1所示的混凝土單邊切口梁四點彎曲試驗進行數(shù)值模擬驗證?；炷猎嚰叽鐬?00 mm×100 mm（長×高），凈跨為450 mm，底部預(yù)制切口高度d為10 mm。模型中預(yù)制切口采用去除節(jié)點水平向約束的零厚度切口模擬。試件底部右側(cè)支撐處和左側(cè)支撐處節(jié)點均施加豎直向和水平向約束。數(shù)值計算過程采用位移加載方式控制，在梁頂部兩個三分點處的節(jié)點分別施加強制位移荷載。為與文獻[16]作對比分析，采用試件底部預(yù)制切口偏左側(cè)7.5 mm處節(jié)點的豎向位移作為梁的撓度。

圖1 混凝土單邊切口梁四點彎曲試驗的幾何及邊界條件示意（單位：mm）Fig.1 Geometry and boundary conditions for four-point bending test (unit: mm)

首先，進行并發(fā)多尺度區(qū)域分解法建模。如圖2(a)所示，將試件分解為10個獨立的子區(qū)域，其中試件中部的子區(qū)域數(shù)量較多，主要是有限元網(wǎng)格精度區(qū)分需要，以及驗證子區(qū)域分解是否影響模型計算結(jié)果。各子區(qū)域根據(jù)計算精度要求進行不同尺寸的有限元網(wǎng)格剖分，剖分情況如下：左、右兩側(cè)的子區(qū)域Ω(1)、Ω(10)均采用h=10 mm的有限元網(wǎng)格；子區(qū)域 Ω(2)、Ω(5)、Ω(6)、Ω(9)均采用h=5 mm 的有限元網(wǎng)格；梁中部可能發(fā)生應(yīng)變局部化的4個方形子區(qū)域Ω(3)、Ω(4)、Ω(7)、Ω(8)的網(wǎng)格尺寸一致，為驗證模型是否存在網(wǎng)格敏感性問題，將這4個子區(qū)域分別剖分為5.000、1.250 和0.625 mm共3種不同精度的有限元網(wǎng)格進行對比分析。支撐點和加載點均采用寬為20 mm，頂部高出梁10 mm的剛體進行模擬，其彈性模量取值比混凝土梁的要大一個數(shù)量級，支撐點處的網(wǎng)格進行適當(dāng)局部加密處理。本算例假設(shè)為平面應(yīng)力問題進行計算，各子區(qū)域有限元計算網(wǎng)格均采用四節(jié)點的四邊形單元。為直觀地說明有限元網(wǎng)格剖分情況，在圖2(b)中給出了梁中部子區(qū)域網(wǎng)格為1.250 mm時的有限元網(wǎng)格剖分示意。

圖2 混凝土單邊切口梁子區(qū)域分解及有限元網(wǎng)格剖分示意（單位：mm）Fig.2 Domain decompostion and finite element of four-point bending beam (unit: mm)

隱式梯度損傷模型計算所采用的參數(shù)說明[16]如下：采用修正的Von Mises模型計算等效應(yīng)變，其混凝土的壓縮和拉伸強度比值η為10；采用Mazars模型中的損傷演化法則，其控制軟化速率的β取值為300，控制殘余應(yīng)力的α取值為0.92，損傷發(fā)生時的應(yīng)變閾值k0取值為0.000 075；彈性模量E取值為40 000 MPa，泊松比?為0.2，梯度參數(shù)c=4 mm2。雙重組裝后線性方程組系統(tǒng)的求解采用直接求解法進行求解，求解器為Pardiso。

計算得知，各加載步的牛頓拉弗森迭代步數(shù)均在4個加載步以內(nèi)計算收斂，可知模型魯棒性好。梁中部3種不同精度有限元網(wǎng)格剖分模型所得的損傷分布差別不明顯，這里僅將中等精度網(wǎng)格h=1.250 mm的損傷分布和發(fā)展過程在圖3中展示。由圖3可知，損傷萌生于預(yù)制切口的尖端，并迅速沿梁中軸線向梁頂部擴展。由于試件結(jié)構(gòu)和加載方式的對稱性，損傷場關(guān)于梁中軸線對稱分布，損傷模式與Simone等[16]的模擬結(jié)果（見圖4，圖例為損傷變量，0為未損傷，1為完全損傷）基本一致。同時，損傷分布在子區(qū)域界面連續(xù)性良好，未表現(xiàn)出子區(qū)域分解敏感性。

圖3 混凝土單邊切口梁四點彎曲試驗的損傷擴展過程（h=1.250 mm）Fig.3 Damage evolution of four-point bending beam (h=1.250 mm)

圖4 Simone文獻混凝土單邊切口梁四點彎曲試驗的損傷擴展過程（h=5.000 mm）[16]Fig.4 Damage evolution of four-point bending test in the paper of Simone (h=5.000 mm)[16]

圖5所示非局部等效應(yīng)變分布和發(fā)展過程表明，應(yīng)變集中區(qū)域與損傷集中區(qū)域一致，均位于梁純彎段的底部預(yù)制切口尖端區(qū)域。且非局部應(yīng)變的分布在子區(qū)域界面保持完好連續(xù)性，受子區(qū)域分解的影響較小。

圖5 混凝土單邊切口梁四點彎曲試驗的非局部等效應(yīng)變發(fā)展（h=1.250 mm）Fig.5 Nolacal equivalent strain evolution of four-point bending beam (h=1.250 mm)

梁中部選用不同精度網(wǎng)格模型所得到的荷載-豎向位移曲線分別如圖6中的實線所示。峰值荷載前的線性增長階段，試驗結(jié)果與數(shù)值模擬結(jié)果的一致性較好；對于峰值荷載，其中試驗峰值荷載為3.68 kN；文獻數(shù)值模擬峰值荷載為3.96 kN；本文h=5.000 mm峰值荷載為3.81 kN，h=1.250 mm峰值荷載為3.60 kN，h=0.625 mm峰值荷載為3.50 kN?？芍逯岛奢d隨著網(wǎng)格精細程度的增加有小幅減小，但與試驗相差2.2%～4.9%；軟化階段h=5.000 mm與參考文獻的數(shù)值模擬結(jié)果相近，h=1.250 mm和h=0.625 mm的軟化曲線幾乎重合，試驗軟化曲線介于h=5.000 mm和h=0.625 mm的軟化曲線之間。綜上可知，基于隱式梯度損傷模型的并發(fā)多尺度區(qū)域分解法不存在網(wǎng)格敏感性，且計算結(jié)果受子區(qū)域分解影響有限，這表明模型是合理的。

圖6 混凝土單邊切口梁四點彎曲試驗的荷載-豎向位移曲線Fig.6 Load-deflection curves of four-point bending beam

3 結(jié) 語

（1）采用隱式梯度損傷模型描述混凝土材料的非線性本構(gòu)關(guān)系，并在易損傷的局部子區(qū)域預(yù)設(shè)高精度有限元網(wǎng)格，同時，采用尺度間線性多點約束進行不同網(wǎng)格精度子區(qū)域的界面協(xié)調(diào)連接，構(gòu)建了混凝土損傷分析的并發(fā)多尺度區(qū)域分解模型。

（2）將模型用于混凝土單邊切口梁的四點彎曲試驗?zāi)M，對可能損傷的局部子區(qū)域分別采用了3種不同精度的有限元網(wǎng)格計算，結(jié)果表明模型不具網(wǎng)格敏感性。同時，損傷及非局部等效應(yīng)變分布受子區(qū)域分解影響不顯著，在子區(qū)域界面的連續(xù)性良好，荷載-豎向位移曲線的計算結(jié)果與數(shù)值模擬結(jié)果的一致性較好。計算結(jié)果表明模型是合理的，可用于混凝土損傷分析的并發(fā)多尺度區(qū)域分解計算。

（3）本文建立的并發(fā)多尺度區(qū)域分解法基于局部易損子區(qū)域的高精度有限元網(wǎng)格預(yù)設(shè)。考慮到混凝土材料細觀結(jié)構(gòu)的高度非均質(zhì)性，其裂縫擴展路徑具有不確定性，數(shù)值計算過程中需要不斷地進行網(wǎng)格自適應(yīng)更新，以滿足損傷局部化區(qū)域的高精度網(wǎng)格剖分要求。因此，建立子區(qū)域的自適應(yīng)更新方法，實現(xiàn)損傷子區(qū)域高精度有限元網(wǎng)格的實時自適應(yīng)更新是下一步研究工作的重點。