張 立，李廣凱，馬懷發(fā)

（1.山東泰山抽水蓄能電站有限責(zé)任公司，山東 泰安 271000）

（2.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038）

（3.中國水利水電科學(xué)研究院 工程抗震研究中心，北京 100048）

1 研究背景

混凝土類材料一般具有初始微缺陷或微孔洞，在外載荷作用下混凝土類材料內(nèi)部的微孔洞、微裂紋等初始缺陷延長、擴(kuò)展、交匯，從而形成微裂縫。微裂縫引起材料剛度劣化和強(qiáng)度降低。同時，微裂縫界面的摩擦滑動引起類似于金屬材料的晶格位錯，產(chǎn)生不可恢復(fù)殘余（塑性）變形[1-3]。基于金屬晶體滑移或位錯發(fā)展起來的彈塑性力學(xué)，能夠較好地描述混凝土類材料的彈性變形、破壞條件以及不可恢復(fù)變形的演化發(fā)展。盡管混凝土類材料的破壞機(jī)理與金屬材料有所不同，但塑性理論在處理不可恢復(fù)變形具有顯著優(yōu)勢，且損傷力學(xué)可以很好地描述其內(nèi)部微缺陷或微裂縫的演化規(guī)律及其對其宏觀力學(xué)性能的影響，二者結(jié)合形成的彈塑性損傷模型可以全面地描述混凝土類材料的非線性變形特征[4-6]。彈塑性問題是典型的材料非線性問題，其應(yīng)力一般與變形或加載歷史有關(guān)，需要采用增量法求解。因此，增量應(yīng)力應(yīng)變本構(gòu)方程與迭代法相結(jié)合被廣泛應(yīng)用于求解彈塑性問題。牛頓-拉斐遜法包括N-R和mN-R 算法，是求解非線性方程組最常用的增量迭代法。人們?yōu)榱烁纳剖諗克俣忍岢隽烁鞣N加速迭代的措施，艾特肯加速法[7]和線性搜索法[8-11]與牛頓法結(jié)合使用，以減少計算迭代次數(shù)。但這些算法都需要剛度矩陣顯式表達(dá)形式，往往還受限于本構(gòu)方程的積分形式。與上述方法不同的是，彈塑性隱式全迭代法[12]能夠聯(lián)立求解彈塑性問題的平衡方程、屈服函數(shù)和塑性流動方程，不需要剛度矩陣的顯式形式，也不拘泥于變形的增量和全量形式，并且不受屈服函數(shù)和塑性流動方程具體形式的限制，可擴(kuò)展應(yīng)用于求解混凝土類材料的非線性問題。

在極限荷載作用下，混凝土類材料的非線性變形一般伴隨著損傷演化和塑性流動的過程，即為彈塑性損傷過程。基于應(yīng)變等效假定[13-15]，在有效應(yīng)力空間利用塑性力學(xué)方法，用有效應(yīng)力張量代替經(jīng)典塑性力學(xué)中的柯西（名義/實際）應(yīng)力張量，以考慮不可恢復(fù)變形。由于有效應(yīng)力隨著彈性應(yīng)變的增加而單調(diào)增加，屈服面處于膨脹狀態(tài)，不會出現(xiàn)柯西應(yīng)力空間材料軟化所導(dǎo)致的屈服面收縮情況，因此，只需考慮應(yīng)力強(qiáng)化，避開了處理應(yīng)變軟化的問題。

本文首先介紹彈塑性問題的全隱式迭代法的基本思想，再引入局部回映迭代，以改進(jìn)其收斂性。然后，基于混凝土類材料彈塑損傷變形特性及其在有效應(yīng)力空間應(yīng)力強(qiáng)化的特點，將混凝土類材料彈塑性損傷問題分解為彈塑性問題和損傷問題。在有效應(yīng)力空間上采用彈塑性全隱式迭代法求解其彈塑性變形，由損傷參數(shù)描述材料剛度的弱化，并將有效應(yīng)力轉(zhuǎn)換為物理空間的實際（名義）應(yīng)力。這種求解方法即為本文提出的混凝土類材料彈塑性損傷問題的全隱式迭代法。

2 全隱式迭代法的基本思想

彈塑性材料的本構(gòu)關(guān)系通常由屈服函數(shù)和流動規(guī)律以增量形式給出。本構(gòu)關(guān)系一般表示為一種顯式形式，Δσ=DepΔε，其中的Δσ、Δε分別為應(yīng)力增量和應(yīng)變增量，Dep為彈塑性矩陣。彈塑性矩陣與加載過程有關(guān)，將根據(jù)應(yīng)力應(yīng)變狀態(tài)進(jìn)行調(diào)整。本文提出的全隱式迭代法的基本思想是將含有塑性因子未知量的屈服條件和塑性流動方程作為基本方程，并與平衡方程聯(lián)立，通過隱式迭代求解變形及其相關(guān)物理量。假設(shè)εn、εp n、σn和Δλn分別是第n個加載步的應(yīng)變、塑性應(yīng)變、應(yīng)力和塑性因子。在n+1個加載步，求解描述彈塑性問題的方程組為：

式中：Y為屈服函數(shù)，是應(yīng)力和內(nèi)變量κ 的函數(shù)；G為塑性流動勢函數(shù)；F為外分布荷載。

在n+1 加載步的應(yīng)變增量為Δεn、塑性應(yīng)變增量Δεp n、內(nèi)變增量Δκn和應(yīng)力增量Δσn，則總應(yīng)變、塑性應(yīng)變、應(yīng)力以及內(nèi)變量分別表示為：

其中De為彈性矩陣。

令：

塑性因子的計算公式為：

將式（6）代入式（3），可以得到：

即：

其中彈塑性矩陣為：

如果Yn=0，應(yīng)力增量為：

式（1）中的第一個方程即平衡方程，可以表示為：

將式（3）代入式（11），可以得到：

將式（6）代入式（12），并消去Δλn，就得到求解彈塑性問題基本方程的積分弱解形式：

3 彈塑性問題的全隱式迭代法

3.1 求解彈塑性方程的隱式迭代格式式（13）是式（1）的等價弱解形式。在彈性加載或塑性卸載時，塑性因子Δλn小于零，則將式（13）就轉(zhuǎn)化為線性彈性方程：當(dāng)應(yīng)力狀態(tài)處于塑性區(qū)時，塑性因子Δλn大于零，則需要求解彈塑性方程。

在每一加載步內(nèi)，保持外部荷載為常量，第k迭代步的方程（13）有如下形式[12]：

傳統(tǒng)的彈塑性問題求解，一般按式（8）計算彈塑性矩陣，按式（9）求解應(yīng)力增量。與傳統(tǒng)方法不同，本文的隱式迭代法不需用顯式求出彈塑性矩陣，總應(yīng)力和應(yīng)力增量在迭代過程結(jié)束由式（15）自動給出。式（14）的迭代格式是全量的形式，這樣在解決實際問題時其邊界條件可以用全量的形式給出。

3.2 回映迭代法上述迭代法認(rèn)為，當(dāng)時可以通過平衡方程迭代式的右端項拉回到屈服面，即當(dāng)?shù)剑?4）收斂時，如果則且有即可得到方程的平衡形式：當(dāng)應(yīng)力進(jìn)入屈服狀態(tài)，迭代過程實質(zhì)際上是應(yīng)力回歸屈服面的過程，也是應(yīng)力重分布的過程。

如圖1所示的受集中荷載作用的簡支梁，梁跨度為2l，高度為2h，厚度為b，跨中受集中荷載P作用，若梁的材料為理想彈塑性，并服從Mises 屈服準(zhǔn)則，按照塑性理論可知，當(dāng)梁跨中橫截面彈性高度he =h，即跨中橫截面完全處于彈性狀態(tài)，彈性極限彎矩當(dāng)he =0，即跨中橫截面處于完全塑性狀態(tài)，已形成塑性鉸，此時塑性極限彎矩則有在塑性極限彎矩作用下形成塑性鉸時，有如圖1所示的屈服塑性區(qū)。其塑性區(qū)邊界為在梁的上、下邊緣離跨中距離處，he =h，塑性區(qū)在上下邊緣的寬度為梁長度的但按式（14）模擬得到如圖2的屈服塑性區(qū)，由于應(yīng)力過于集中，超極限應(yīng)力不再收斂到屈服面，此時的計算得到塑性區(qū)與理論解差別較大。

圖1 集中荷載作用下簡支梁的理論塑性區(qū)分布

圖2 集中荷載作用下數(shù)值迭代不收斂的塑性區(qū)

分析其原因，可能是由于式（14）是針對全區(qū)域的迭代，很難保證應(yīng)力過于集中的局部小范圍屈服函數(shù)得到精確滿足。為了使屈服應(yīng)力回歸到屈服面，改善的收斂性，下面將在迭代式（13）中增加了應(yīng)力在屈服點的回映迭代，即回映算法。

回映算法實質(zhì)上是一種彈性預(yù)測和塑性修正過程，需要局部迭代調(diào)整塑性因子，使預(yù)測應(yīng)力返回到屈服面[16]。這里的回映迭代式采用了半隱式向后歐拉方法[17]，其具體步驟如下所述。

（1）步驟1。給定n+1步的初始條件：計算令

圖3 周圍受壓的圓形隧筒的塑性區(qū)

3.3 彈塑性全隱式迭代法的驗證算例1。在隧洞開挖時巖體伴隨應(yīng)力重分布，隧洞周邊局部可能進(jìn)入塑性狀態(tài)，其受力變形狀態(tài)可以等同于厚壁圓筒，外部受均勻壓力，如圖3（a）所示。在圍壓p作用下，如巖體滿足Mohr-Coulomb 準(zhǔn)則屈服條件，隧洞半徑為R0的圓孔周圍塑性區(qū)半徑具有理論解：其中，c為巖體的黏聚強(qiáng)度，φ為巖體的內(nèi)摩擦角。本算例中，圓形孔口半徑R0為1 m，若黏聚強(qiáng)度c取1.0 MPa，摩擦系數(shù)f即tanφ取0.2。按彈塑性全隱式迭代法采用兩種方案求解，即分別采用和不采用局部回映迭代，計算得到圖3（b）所示的塑性半徑與壓力關(guān)系曲線。由該曲線可以看出，施加外部周邊壓力至6.5 MPa之前，兩種方案得到的塑性半徑與壓力關(guān)系曲線，與理論曲線基本一致。但在壓力值大于6.5 MPa 之后，采用局部回映迭代的數(shù)值解與理論解仍然吻合較好。圖3（c）給出了在外部周邊壓力達(dá)到7.0 MPa時，采用局部回映迭代得到的塑性區(qū)分布情況，此時塑性半徑Rp為3.0 m，與理論值基本一致。

算例2。仍以如圖1所示的簡支梁為例，設(shè)梁材料的彈性模量取1GPa，泊松比取0.25，并為理想彈塑性材料，服從Mises 屈服準(zhǔn)則，屈服應(yīng)力σs取2.7 MPa。在上、下邊緣正應(yīng)力剛達(dá)到屈服，此時的荷載定義為Ps，即有令b=1 m，h=0.1 m，l=0.6 m，則

用0.1Ps增量荷載進(jìn)行靜態(tài)數(shù)值模擬，加載至1.5Ps，即MpMe =1.5時，計算得到塑性區(qū)分布如圖4所示。與如圖1所示的理論塑性區(qū)分布相比，此時計算得到的塑性區(qū)分布與理論塑性區(qū)分布幾乎完全相同。與圖2的計算結(jié)果相比可以看出，采用局部回應(yīng)迭代后使迭代收斂性得到了大大的改善。

圖4 集中荷載作用下的塑性區(qū)數(shù)值模擬

以上數(shù)值計算是在配置為Intel Core i5-2400 CPU@3.10GHz，MEM 2.99GB 的PC 機(jī)上進(jìn)行的。在上面的算例中，平衡迭代誤差屈服函數(shù)控制在簡支梁有限元數(shù)值模型共有13 202個自由度，15個加載步，耗時不到3 min；平衡迭代最多需要15次迭代；局部迭代只需3或4次即可趨近于零。這說明本文提出的迭代法具有較高的計算效率，并且迭代穩(wěn)定性良好。

4 彈塑性損傷問題的全隱式迭代法

4.1 求解彈塑性損傷方程的隱式迭代式在有效應(yīng)力空間中，有效應(yīng)力取代實際（名義）應(yīng)力σ。按照應(yīng)變等效假設(shè)，在實際物理空間里的實際應(yīng)力σ可用損傷參數(shù)d和有效應(yīng)力表示為混凝土非線性變形的分解，及其與彈性塑性損傷的關(guān)系如圖5所示，E0為初始彈性模量，為極限彈性強(qiáng)度，ε0為與極限彈性強(qiáng)度對應(yīng)的極限應(yīng)變，即有為塑性應(yīng)變（殘余應(yīng)變），εel為彈性應(yīng)變。由于有效應(yīng)力隨彈塑性應(yīng)變的增加而單調(diào)增加，因此屈服面總是膨脹狀態(tài)，因此，實際物理空間里的材料變形軟化不會導(dǎo)致有效應(yīng)力空間的屈服面收縮。

圖5 混凝土非線性變形與彈性塑性損傷的關(guān)系

由于有效應(yīng)力σˉ和塑性變形滿足塑性理論的屈服準(zhǔn)則、強(qiáng)化法則、加載-卸載準(zhǔn)則和流動準(zhǔn)則，因此，第3節(jié)的彈塑性問題的全隱式迭代法可以移植到有效應(yīng)力空間。在每一加載步內(nèi)，保持外部荷載Fn+1為常量，在第k迭代步的迭代式（14）可以改寫為如下求解彈塑性損傷方程問題的隱式迭代式：

另外，式（17）的迭代式是以全量的形式給出，這樣應(yīng)用該方法分析高壩地震響應(yīng)時，可以方便地以全量的形式輸入地震動荷載。

4.2 彈塑性損傷全隱式迭代法的驗證屈服準(zhǔn)則采用Lee和Fenves[18-19]基于Lubliner 等[20]提出的混凝土塑性損傷模型，塑性流動勢函數(shù)采用Drucker-Prager 雙曲線函數(shù)。計算模型采用混凝土試件尺寸和網(wǎng)格剖分如圖6。取文獻(xiàn)[21]的軸拉試驗和軸壓試驗觀測數(shù)據(jù)：彈性模量為53.6 GPa，抗拉強(qiáng)度為1.84 MPa，抗壓強(qiáng)度為25.28 MPa，泊松比為0.2，密度為2400 kg/m3。由文獻(xiàn)[22-23]試驗觀測數(shù)據(jù)分析得到如圖7和圖8給出的混凝土材料的峰（失效）后應(yīng)力、損傷變量與非線性（開裂）變形關(guān)系曲線。

圖6 混凝土試件有限元數(shù)值模型及網(wǎng)格剖分（單位：mm）

圖7 壓縮非線性參數(shù)曲線

圖8 拉伸非線性參數(shù)曲線

為了保證在跨中區(qū)域發(fā)生斷裂破壞，將梁跨中下邊緣寬45 mm×高22.5 mm 的區(qū)域的彈性模量分別取0.67 GPa和5.36 GPa，即為實測彈性模量的1/80和1/10。對圖6所示的混凝土彎拉試件，施加豎向位移荷載P（w），位移增量取為10-3mm，進(jìn)行彎拉破壞全曲線數(shù)值模擬。為了消除計算結(jié)果的網(wǎng)格敏感性，在有限元計算時，由圖8給出的損傷變量和峰后應(yīng)力與開裂位移關(guān)系，計算單元拉伸塑性應(yīng)變單元的特征長度he取單元面積的平方根。

圖9給出了隨豎向位移增加，梁的彎拉應(yīng)力由彈性階段到達(dá)峰值，然后進(jìn)入軟化階段的過程。梁（P為荷載，L支座間距，b為梁厚度，h為梁高度）與梁上邊緣中點豎向位移相對應(yīng)。

圖9 混凝土試件彎拉應(yīng)力與變形全曲線數(shù)值計算與試驗結(jié)果

表1 混凝土靜態(tài)彎曲破壞試驗值與計算值

表1列出了由同一批次混凝土澆筑兩個混凝土試件[21]：編號為F110509A22和F110517A03 彎拉峰值上邊緣中點彎拉應(yīng)力荷載及其對應(yīng)的彎拉強(qiáng)度的試驗值和數(shù)值計算結(jié)果，其中數(shù)值模擬-1計算得到彎拉應(yīng)力峰值為2.90 MPa，與試件F110517A03 試驗值2.49 MPa 比較接近；數(shù)值模擬-2 計算得應(yīng)力峰值為3.40 MPa，與試件F110509A22的試驗值3.45 MPa接近。

由圖9給出的混凝土試件彎拉應(yīng)力與變形全曲線中的數(shù)值計算與試驗結(jié)果對比可以看出，在軟化曲線的中間段，數(shù)值模擬與試驗觀測曲線存在一定差異，但最后與實測試件的軟化變形曲線趨于一致。另一方面，從梁跨中局部區(qū)域彈性模量不同取值所得兩條計算全曲線分別與同一批次混凝土兩試件的試驗結(jié)果相吻合，這也說明實際試件在局部材料性能的不均勻可能是導(dǎo)致試驗觀測結(jié)果差異的誘因。

圖10給出了試件彎拉損傷破壞數(shù)值模擬結(jié)果。圖10（a）為加載至荷載峰值狀態(tài)下的最大主應(yīng)力分布，最大主拉應(yīng)力接近抗拉極限強(qiáng)度1.84 MPa，位于裂紋的擴(kuò)展的前沿，說明最大應(yīng)力控制在屈服面上；圖10（b）給出的損傷云圖是梁彎曲變形接近失穩(wěn)時的損傷云圖，由圖10（b）看出，梁從跨中下邊緣開始起裂，向上邊緣擴(kuò)展。由以上混凝土試件彎拉損傷破壞全過程的數(shù)值模擬分析可以看出，計算結(jié)果與材料試驗觀測結(jié)果吻合較好，從而驗證了本文提出的彈塑性損傷問題的隱式迭代法的正確性。

圖10 混凝土試件彎拉損傷破壞數(shù)值模擬結(jié)果

5 結(jié)論

本文的全隱式迭代法將屈服函數(shù)和塑性流動方程作為基本方程，聯(lián)立隱式求解由平衡方程、屈服條件和塑性流動方程組成的方程組，無需顯式計算彈塑性矩陣。該算法在引入局部回映迭代法后，通過很少幾次回映迭代就可以獲得理想的收斂精度，并且具有很好的數(shù)值穩(wěn)定性。

由于在有效應(yīng)力空間里混凝土類材料的塑性損傷遵循塑性變形規(guī)律，本文基于等效應(yīng)變假設(shè)，將混凝土類材料的彈塑性損傷問題分解為彈塑性問題和損傷問題。應(yīng)用彈塑損傷問題全隱式迭代法在有效應(yīng)力（虛擬）空間計算塑性變形，再按應(yīng)變等價的原則，考慮材料的損傷程度，即借助于損傷參數(shù)將有效應(yīng)力轉(zhuǎn)換為物理空間的實際應(yīng)力。本文混凝土試件彎拉損傷數(shù)值計算表明，在有效應(yīng)力空間里，隨有效應(yīng)力的增加，塑性（殘余變形）屈服面處于膨脹狀態(tài)，混凝土彈塑性損傷問題的全隱式迭代是穩(wěn)定的。由于損傷參數(shù)為塑性應(yīng)變的單調(diào)增函數(shù)，因此損傷參數(shù)的增大只會影響實際應(yīng)力的大小，而不會影響塑性變形迭代求解的穩(wěn)定性。

最后強(qiáng)調(diào)指出，本文的隱式迭代格式是以全量的形式給出，這樣應(yīng)用該方法分析高混凝土壩非線性地震響應(yīng)時可以很方便地以全量形式實現(xiàn)地震動荷載輸入。另外，本文提出的彈塑性損傷問題全隱式迭代法不僅適應(yīng)于混凝土壩體，也適應(yīng)于壩基巖體的彈塑性損傷問題的求解。