循環荷載下混凝土開裂-閉合行為計算方法研究

李曉琴,張 田

(昆明理工大學 建筑工程學院，昆明 650500)

隨著有限元理論和計算技術的發展，有限元數值模擬方法已成為研究各種受力情況下RC結構響應的主要手段之一。ABAQUS中的CDP模型因具有較好的屈服準則、流動法則和滯回準則被廣泛用于各類RC結構的抗震分析[1]。相對主要用于單調荷載下分析的脆性開裂模型[2]和彌散裂縫模型[3]，CDP模型通過引入損傷參數(dt、dc)和剛度恢復因子(wt、wc)來分別描述混凝土在卸載時表現出的剛度退化和應力反轉時混凝土存在的剛度恢復現象[4]，該模型能較好地模擬循環荷載下混凝土的力學特性，如循環荷載下混凝土的開裂與壓碎、裂縫張開與閉合、拉壓應力反轉和損傷累積等特征。然而，對于有限元分析，定義合理的材料參數是模擬結果正確反映真實結構力學行為的關鍵所在。

目前，大量學者對CDP模型的材料參數進行了相關研究，如Syed等[5]對混凝土本構關系中的膨脹角Ψ與參數wt、wc的組合取值給模擬結果帶來的不確定性進行了研究；Yu等[6-7]對CDP模型采用的非關聯塑性流動勢和Drucker-Prager屈服準則進行了詳細討論；Alfarah等[8]基于混凝土破碎斷裂能Gc和開裂斷裂能GF闡述了混凝土單軸受壓、受拉材料參數的計算；聶建國等對CDP模型的單軸應力-應變關系和裂縫模型等做了詳細介紹。事實上，在CDP模型中，除了以上的材料參數需要明確以外，損傷參數(dt、dc)和剛度恢復因子(wt、wc)的確定也需要關注，對于損傷參數進行相關研究中，Alfarah等基于Gc和GF推導了新的拉壓損傷計算公式；Long等[9]采用累計耗散能與斷裂能的比值來建立損傷計算模型；方自虎等[10]基于軟化應力與峰值應力的比值建立了適合于高強混凝土模擬分析的損傷計算公式；Birtel等[11]基于試驗數據建立了與混凝土塑性應變有關的損傷計算公式；曹明[12]以能量等價原理為基礎，并結合GB 50010—2010《混凝土結構設計規范》中提供的混凝土應力-應變曲線建立了損傷參數計算公式。然而，在CDP模型眾多材料參數中，針對描述循環荷載下混凝土開裂-閉合行為的受壓剛度恢復因子wc的取值至今未有深入研究。目前一些關于wc的研究工作主要有：Syed等研究指出wc取值很難直接從試驗中直接測試得到，wc取值需要通過對比模擬與試驗結果來標定；Zhao等[13]著重考慮wc對洞室結構動力特性的影響，研究發現當wc取值增大時洞室的損傷面積呈線性減小；盡管Karimi等[14]取wc=0.3來模擬填充墻和拱形砌體墻在循環荷載下的力學響應，雖然兩個有限元模型的模擬值與試驗結果的均具有較好的吻合度，但對wc的取值原因尚未說明；Alfarah等[15]為了模擬RC結構中嚴重破損區域處裂紋閉合時剛度恢復較小的現象，采用較小wc值定義構件嚴重破損區域。顯然，對于wc的取值與應用問題，目前還沒有足夠的研究以一致的解決方法，這給CDP模型準確預測和模擬循環荷載下RC結構的力學響應帶來了困難，尤其是在模擬對象尚未完成試驗的情況下。

鑒于此，為解決上述問題，本研究針對wc的計算方法展開了理論研究和數值模擬，首先建立了考慮開裂斷裂能GF的指數型軟化應力-開裂位移曲線，并在現有的拉/壓損傷計算模型中比選出適用于鋼筋混凝土結構分析的損傷計算模型。在此基礎上，以受拉殘存斷裂能密度gFR與開裂斷裂能密度gf之比的θ次方(θ≥1)建立了wc計算模型，該模型同時考慮了單元特征長度leq和dt對wc的影響，并以單個單元為例說明了參數θ對wc的影響。進一步地，以循環荷載下嚴重破損RC節點為研究對象，基于現有混凝土裂縫特征與受拉損傷狀態的統計關系[16]，提出了一種基于拉伸損傷水平dt的模型簡化處理方法，以此解決ABQUAS中CDP模型wc只支持某一常數所存在的局限性。通過對多組循環荷載下的RC節點進行了數值模擬并與試驗數據比較，驗證了循環荷載下混凝土開裂-閉合行為計算方法的正確性以及循環荷載下有限元模型簡化處理方法的可行性。

1 CDP模型

圖1(a)顯示了混凝土單軸受拉時的力學行為，認為混凝土包括彈性階段(0～σt0)和軟化階段(下降段)。當混凝土超過破壞應力σt0時，混凝土將進入軟化和剛度退化階段。圖1(b)顯示了混凝土單軸受壓時的力學行為，認為混凝土受壓時分別經歷了彈性階段(0～σc0)、強化階段(σc0～σcu)和軟化階段。

(a) 受拉

現假定E0為材料的初始彈性剛度，則上述的應力-應變關系可用式(1a)～(1b)描述。

(1a)

(1b)

(2a)

(2b)

(3a)

(3b)

試驗表明，在單軸循環荷載下，荷載改變方向后材料的彈性剛度將得到部分恢復[20-21]。CDP模型引入剛度退化變量d描述混凝土的剛度退化行為，其中剛度退化變量d是關于應力狀態和拉壓損傷變量的函數。

E=(1-d)E0

(4)

(1-d)=(1-stdc)(1-scd)

(5)

式中：st和sc為考慮應力方向和剛度恢復效應的應力狀態函數，根據式(6a)～(6b)確定。

st=1-wtr*(σ11) 0≤wt≤1

(6a)

sc=1-wc(1-r*(σ11)) 0≤wc≤1

(6b)

式中：wt和wc分別為受拉、受壓剛度恢復因子，分別表示混凝土從受拉過渡到受壓和從受壓過渡到受拉時混凝土剛度的恢復程度；r*(σ11)為混凝土單軸應力因子，是主應力狀態的函數，由式(7)確定。

(7)

式中：σ11為混凝土的主應力狀態，σ11>0表示受拉，σ11<0表示受壓。

以單軸工況為例，CDP模型關于循環荷載下混凝土的應力-應變曲線如圖2所示。

圖2 單軸循環荷載作用下應力-應變關系(默認剛度恢復因子wt=0、wc=1)Fig.2 The stress-strain relationship under uniaxial cyclic loading (wt=0,wc=1)

CDP模型的塑性流動假定為非關聯流動，塑性勢G采用Drucker-Prager雙曲線函數描述。

(8)

式中：Ψ為高圍壓下在p-q面上的膨脹角(對混凝土取值在25°～35°區間，本文模型采用30°)；ξ為偏心率(默認值0.1)，表示該函數接近漸近線的速率(當偏心率趨于零時，流動勢趨于直線)。

(9)

式(9)中，各參數的計算公式如下

2 混凝土單軸受拉/抗壓本構關系

混凝土的單軸受壓力學行可通過三段函數定義，第一段為線彈性階段式(11a)；第二段為強化段[22]式(11b)；第三段為軟化段[23]式(11c)、式(11d)。

σc=E0εc

(11a)

(11b)

(11c)

(11d)

對于混凝土的單軸拉伸力學行為，本研究基于混凝土的開裂斷裂能GF提出了適用于混凝土非線性分析的指數型軟化段應力-開裂位移曲線(圖3)，其表示如下

圖3 應力-開裂位移曲線Fig.3 The stress-crack opening relationship

σt=ftme-cw

(12)

式中：w為開裂位移，單位mm；ftm為混凝土抗拉強度平均值，根據式(13a)計算，單位MPa；c為軟化參數，與混凝土的強度等級有關。經數據擬合后，參數c按式(13b)確定。

(13a)

(13b)

如圖3所示，軟化應力-開裂位移曲線是根據開裂斷裂能GF(曲線陰影面積Area)、開裂應力ftm和最大開裂位移w0這三參量確定的，能夠反映混凝土受拉應力隨著開裂位移的增加而逐漸退化(軟化行為)和裂縫開展過程中的能量耗散等特征。此外，該應力-開裂位移曲線適用于強度在12～120 MPa下的混凝土材料。

根據Bazant等[24]提出的裂縫帶理論，應力-開裂位移曲線可根據式(14)[25]轉換成應力-開裂應變曲線。

(14)

混凝土的開裂斷裂能GF根據式(15)計算，混凝土受壓破碎斷裂能Gc根據式(16)[26]確定。

GF=0.073(fcm)0.18

(15)

(16)

3 混凝土拉/壓損傷參數評價與選用

由式(1a)～(1b)和式(4)～(5)可知，描述混凝土卸載時剛度弱化響應的損傷參數包括受拉dt、受壓dc損傷參數。由于目前ABAQUS用戶手冊中沒有指明損傷變量的具體計算方法，因此現將已有的幾種損傷參數計算方法進行對比評價(表1)，以期選擇適用的損傷參數模型。

表1 損傷參數計算模型Tab.1 Calculation models of damage parameters

本研究以C30混凝土為例，采用上述5種計算模型確定的拉、壓損傷參數曲線如圖4所示。由圖4(a)可知，應變等效模型[27]計算的受壓損傷參數dc在所有曲線的最上方，損傷發展最為激烈。方自虎模型[10]下的受壓損傷參數dc從峰值應力開始計算，沒有考慮混凝土在強化階段的受壓損傷發展與累積。當非彈性應變小于0.005時，Birtel模型[11]與能量等效模型[28]計算的受壓損傷參數dc比較接近；當非彈性應變超過0.005時，Birtel模型[11]計算的受壓損傷參數dc比能量等效模型更大，說明損傷發展更快。從圖4(b)可知，應變等效模型[27]和Birtel模型[11]計算的受拉損傷參數dt較接近且在所有曲線的最上方，損傷發展最為激烈。方自虎模型[10]計算的受拉損傷參數dt在所有曲線的最下方，損傷發展較慢，不能較好描述混凝土受拉時表現出來的脆性特性。能量等效模型[28]計算的受拉損傷參數dt介于其他3個模型之間，損傷發展偏緩慢。

(a) θ=1

(a) 受壓損傷

以上的對比的四種損傷計算模型，主要是基于混凝土材料而建立的。然而，對于RC結構分析，有必要考慮受力鋼筋對混凝土損傷發展的影響，特別是對于受拉的情況。此外，ABAQUS用戶手冊還指出，對RC結構分析時應適當考慮鋼筋與混凝土協同工作的“拉伸硬化”效應。考慮到直接建立針對RC結構的損傷計算模型難度較大，因此本研究暫未深入研究損傷計算模型，而是通過對比現有的損傷計算模型的發展形態，選擇損傷隨非彈性應變和開裂應變發展速度處于中間水平的Sidoroff能量等效模型作為基準模型，以此來計算RC結構中混凝土的受壓dc、受拉dt損傷參數。

4 開裂-閉合受壓剛度恢復因子wc計算模型

4.1 wc計算模型的建立

CDP模型中通過引入受拉wt、受壓wc剛度恢復因子來描述混凝土在循環荷載下混凝土裂縫張開、閉合及裂縫之間的相互作用，以期更好模擬循環荷載下混凝土的反應。試驗研究表明當混凝土從受壓進入到受拉狀態，若進入受拉前已有受壓微裂縫形成時，應力反轉時混凝土的受拉剛度不能恢復；當混凝土從受拉進入受壓狀態時，隨著裂縫的閉合混凝土的受壓剛度可以恢復。因此，在本文研究中，不考慮混凝土從受壓進入受拉狀態時的剛度恢復效果(wt=0為默認值)，重點關注混凝土從受拉進入為受壓狀態時的剛度恢復程度(其中wc=1為默認值)。圖5給出了混凝土從受拉進入到受壓狀態時受壓剛度恢復因子wc的變化。

圖5 受壓剛度恢復因子wc的影響Fig.5 Effect of the compression stiffness recovery factor wc

關于圖5中wc的變化討論，本研究先假定混凝土在初始受拉時沒有受壓損傷(無受壓塑性應變)，有dc=0，代入式(5)和式(6b)得

(1-d)=(1-scdt)=

1-(1-wc(1-r*(σ11)))dt

(17)

由式(17)知，當混凝土受拉(σ11>0)時，有r*(σ11)=1，因此d=dt。當混凝土受壓(σ11<0)時，有r*(σ11)=0，因此d=(1-wc)dt。若wc=0，則d=dt，此時材料剛度沒有恢復；wc=1，有d=0，材料完全恢復壓縮剛度；0

為對比不同wc取值(0≤wc≤1)對混凝土剛度恢復程度和混凝土單元強度大小，本研究以幾何尺寸為150 mm、強度為C30的立方體混凝土單元為例，分別模擬了混凝土相同dt下wc為0、0.5和1時混凝土單元從受拉進入為受壓狀態時的力學行為(圖6)。

圖6 相同dt下不同wc對混凝土“由拉進入壓”后的力學行為影響Fig.6 Concrete compressive behaviour after cracking under the same dt with different wc values

由圖6可知，隨著wc取值增大，混凝土單元的彈性模量逐漸增大，單元的峰值強度也逐漸遞增，并且混凝土單元在相同變形處卸載時材料的卸載響應也更加緩慢。wc=1時對應的彈性模量為26 209 MPa，分別是wc=0.5(15 722 MPa)和wc=0(5 211 MPa)的1.67倍和5.03倍。wc=1時對應的單元峰值強度為20.1 MPa，分別是wc=0.5(12.1 MPa)和wc=0(4.0 MPa)的1.67倍和5.03倍。可見，wc的取值會直接影響到混凝土單元的彈性模量與峰值強度。

綜上可知，wc值影響混凝土的剛度恢復程度主要體現在彈性模量和峰值強度兩方面。因此，為了準確預測循環荷載作用下RC結構的宏觀承載力，有必要建立確定wc取值的計算模型，以期獲得合理的wc值描述循環荷載下混凝土裂縫開裂后閉合行為后的剛度恢復效應。如圖5所示，對于循環荷載下的混凝土，當受拉混凝土在軟化段B點卸載時，可認為混凝土的微裂縫將沿著卸載路徑B→C逐漸閉合。通常混凝土在拉→壓過渡區存在一個剛度恢復過程，即受拉應力卸載至C點時微裂縫并不會完全閉合。為簡便，本研究假設微裂紋在C點處完全閉合，并將多邊形BCEF的面積視為一種材料性質，即裂縫完全閉合點(C點)處對應的受拉殘存斷裂能密度gFR。進一步地，本研究基于混凝土受拉應力-應變關系(圖5)，以受拉殘存斷裂能密度gFR與開裂斷裂能密度gf比值的θ次方(θ為實數)建立了wc計算模型，其表示如下

(18)

式中：θ(θ≥1)為受壓剛度恢復修正參數；gf為混凝土開裂斷裂能密度(gf=GF/leq)，單位N/mm2；s1與s2分別為三角形BCD與多邊性BDEF區對應的面積(受拉殘存斷裂能密度)，單位N/mm2；s1與s2根據式(19)進行計算。

s1=

(19a)

(19b)

式中：wt為B點對應的開裂位移，單位mm；dt為描述卸載支路(B→C)上的受拉損傷參數；c為式(13)中的軟化參數。

4.2 wc計算模型的演化規律討論

以C30立方體混凝土單元為例，并選擇能量等效模型[28]計算混凝土的受拉損傷參數dt。圖7給出了wc計算模在不同θ取值下wc與dt的變化規律圖。

圖7 θ值對受拉損傷參數-受壓剛度恢復因子wc的影響Fig.7 Effects of θ to the compressive stiffness recovery factor wc

如圖7所示，wc與dt呈非線性遞減關系；在相同dt數值下，wc數值會隨著θ取值增大而減小。可見，θ會直接影響到wc-dt曲線的衰減速度。

同樣以C30混凝土為例，模擬并對比了混凝土在不同dt下wc分別取默認值(wcd1/2/3=1)和本研究提議的計算值下的滯回行為(圖8)。關于計算值選取如圖7所示，以θ=1、3下的wc-dt曲線為基礎，模擬了dt為0.58、0.75和0.83下混凝土卸載后再反向加載時剛度恢復效應。

4.3 wc計算模型在RC節點模型中的應用

以上工作中，已完成了循環荷載下混凝土開裂后閉合行為的wc計算模型的建立，并對其演化規律進行討論。由于ABAQUS中CDP模型下的wc只支持某一常數定義，不能描述混凝土在不同dt下具有不同的剛度恢復效應。為解決這一問題，本研究選擇循環荷載下嚴重破損的RC節點作為研究對象，基于現有混凝土裂縫特征與受拉損傷狀態的統計關系，提出了一種基于拉伸損傷水平dt的模型簡化處理方法。以0≤dt1≤0.75和0.75≤dt2≤dtm(dtm≤1)將RC節點簡單劃分為輕微損傷和損傷-破壞區兩個區域，并分別采用特征受壓剛度恢復因子wcz1和wcz2來定義兩個區域的材料屬性。輕微損傷和損傷-破壞區域范圍分別采用構件的非塑性鉸和塑性鉸區域進行度量，采用由Paulay等[29]提議的塑性鉸長度計算公式確定RC構件的塑性鉸區域(式20)。

Lp=0.08L+0.022fyd

(20)

式中：Lp為塑性鉸長度；L為構件剪跨；d受拉鋼筋直徑，單位mm;fy為受拉鋼筋強度，單位MPa。

受拉損傷參數dt的區間劃分和不同區域的受壓特征剛度恢復因子wcz1與wcz2的計算方法如圖9所示，受拉損傷分區分為0≤dt≤0.75與0.75≤dt≤dtm(dtm≤1)兩個區間，根據面積等效的原則計算每個區間對應的wcz1與wcz2，其中參數wcz1與wcz2分別用來定義圖11(a)中對應的黑色與灰色區域。

圖9 特征受壓剛度恢復因子Fig.9 Characteristic value of the compressive stiffness recovery factor

5 循環荷載下RC節點的數值模擬應用研究

5.1 有限元模型

以Yang等[30]研究中的CL1節點為研究對象，采用ABAQUS軟件進行三維有限元建模分析。混凝土采用實體單元C3D8R模擬，鋼筋采用桁架單元T3D2模擬，并采用Embedded命令將鋼筋內置于混凝土單元內。混凝土采用CDP模型，鋼筋選用由Clough提出的材料模型(圖10)模擬[31-32]，該模型通過削減鋼筋的卸載剛度來替代RC結構中因鋼筋與混凝土間的黏結滑移效而導致的強度退化和剛度退化。

圖10 鋼筋本構關系Fig.10 The reinforcement constitutive relationship

CL1節點的幾何尺寸和配筋詳細見文獻[30]，混凝土150 mm立方體抗壓強度fcu-150為35.5 MPa，縱向鋼筋和橫向鋼筋的屈服強度分為528 MPa和407 MPa。根據試驗的邊界條件，將梁兩端50mm區段采用剛體約束，并指定剛體的轉動中心為梁端面的形心位置處；柱子底部采用固定鉸約束；柱頂先施加豎直向下的軸壓荷載(1 420 kN)，然后在柱頂面施加平行于梁軸線方向的水平循環位移荷載，循環加載制度如圖11所示，每級加載位移循環2次。采用ABAQUS/Standard模塊中的Newton-Raphson算法求解，通過增量步施加荷載，經過若干次平衡迭代逐步獲得解答。有限元元模型邊界和網格如圖12所示。

圖11 CL1節點加載制度Fig.11 Loading system for CL1 joint

(a) 混凝土

5.2 受壓剛度恢復修正參數θ取值討論

如4.2節所述，對一個確定的有限元模型，受壓剛度恢復修正參數θ是一個確定的常量。從wc-dt曲線中可以看出θ的取值會直接影響到wc-dt曲線的衰減速度。因此，本研究將參數θ作為參數變量進行數值模擬研究，并分別取θ等于1、2和3時進行數值計算，通過對比有限元的模擬結果來標定θ的取值范圍，以期得到較好的模擬結果。CL1節點在不同θ取值下模擬結果見圖13和圖14。

圖13 不同θ取值下的模擬結果對比Fig.13 Comparison of the FE results with the different values of θ

(a) θ=2

由圖13可知，當θ增大時RC節點的峰值承載力降低，整個滯回環的面積減小，節點的耗能能力降低。相反地，θ太小可能會高估RC節點的承載力和耗能能力。當θ取2或3時，模擬結果在峰值承載力和滯回環面積方面相差較小，說明數值模擬結果趨于穩定。經有限元分析發現，當θ≥3時，混凝土wc-dt曲線衰減加快，從而導致wc1和wc2過小，這會使得模擬結果嚴重低估結構的承載力甚至導致有限元模型分析不收斂。可見，θ取值應該控制在1～3之間(1≤θ≤3)比較合理。

如圖14所示，通過對比θ=2(或3)下的模擬結果與試驗結果發現，當θ=3時的模擬結果更加逼近試驗結果。因此，對于循環荷載下RC節點，θ=3可作為確定wc計算模型的默認值。

5.3 循環荷載下不同RC節點數值模擬

基于以上循環荷載下混凝土開裂閉合行為相關演化算法，包括受壓剛度恢復因子wc演化規律的合理性和簡化有限元計算方法的可行性，進一步對其他類型框架十字節點CL2、邊節點E1[33]和門式閉合節點F0[34]作為算例進行有限元建模分析，有限元模型的單元類型與材料模型的選擇、算法控制與CL1節點模型保持一致，模擬結果見圖15。

(a) 十字節點CL2

通過骨架曲線，可以得到RC節點的極限荷載，有限元結果與試驗結果對比如表2。

表2 不同RC節點的極限荷載對比Tab.2 Comparison of the ultimate loads for the different type of the RC joints

由表2可知，有限元計算的極限承載力普遍小于試驗值，誤差在2.8%～11.7%。綜合考慮模擬骨架曲線與試驗骨架曲線的整體變化趨勢以及承載力誤差范圍，有限元模型均能較好的反映結構的實際力學響應。由此可見，本文建立的wc演化算法和循環荷載下簡化有限元計算方法的可行性，也驗證了θ=3時能獲得較好的模擬結果。

6 結 論

(1) 本研究基于混凝土開裂斷裂能GF提出了適用于混凝土非線性分析的指數型軟化段應力-開裂位移曲線，該曲線是根據GF、開裂應力ftm和最大開裂位移w0這三參量確定的，能夠反映混凝土受拉應力隨著開裂位移的增加而逐漸退化(軟化行為)和裂縫開展過程中的能量耗散等特征。此外，該應力-開裂位移曲線適用于強度在12～120 MPa下的混凝土材料。

(2) 通過對比評價現有的混凝土拉dt、壓損傷參數dc計算模型，選擇了適用于RC結構分析的能量等效dt和dc計算模型。

(3) 基于混凝土軟化應力-應變關系建立了受拉殘存斷裂能密度gFR與開裂斷裂能密度gf比值的θ次方(θ為實數)的wc計算模型，該計算模型同時考慮了單元特征長度leq和dt對變量wc的影響，同時本文還討論了wc計算模型的演化規律和受壓剛度恢復修正參數θ的取值。

(4) 基于現有混凝土裂縫特征與受拉損傷狀態的統計關系，提出了一種基于拉伸損傷水平dt的模型簡化處理方法。以0≤dt1≤0.75和0.75≤dt2≤dtm(dtm≤1)將循環荷載下嚴重破損的RC節點簡單劃分為輕微損傷(非塑性鉸區)和損傷-破壞(塑性鉸區)兩個區域，分別采用特征受壓剛度恢復因子wcz1和wcz2來定義兩個區域的材料屬性。采用ABAQUS對CL1節點進行建模分析，討論了受壓剛度恢復修正參數θ的取值(1≤θ≤3)。經過分析討論后以θ=3可作為確定wc計算模型的默認值。

(5) 通過對循環荷載下的十字節點、邊節點和門式框架RC節點進行數值模擬并與試驗數據比較，有限元計算的極限承載力與試驗值的誤差在2.8%～11.7%，綜合考慮模擬骨架曲線與試驗骨架曲線的整體變化趨勢以及承載力誤差范圍，有限元模型均能較好的反映結構的實際力學響應，驗證了循環荷載下混凝土開裂-閉合行為計算方法的正確性以及有限元模型簡化處理方法的可行性。