基于Kelvin 鏈率型模型的混凝土徐變有限元計算

楊 超，陳夢成

（華東交通大學土木建筑學院，江西 南昌330013）

混凝土徐變是指混凝土結構在恒載作用下隨時間變化的變形。 混凝土徐變計算的理論基礎主要包括老化理論[1]、多孔介質理論[2]、黏彈性理論[3]、繼效流動理論[4]、固結理論[5]等。 基于上述徐變理論，演繹出了各種混凝土結構的徐變分析方法如：初應變法[6-7]、位移增量法[8]、Galerkin 加權殘值法、等效模量法[9]、齡期調整有效模量法[10-13]、逐步積分法[14]、全量遞推方法[15]、率型疊加算法[16]等。 2001 年Bazant[17]將等效模量法、按齡期調整有效模量法、徐變率法、流動率法及Levi’s & Arutyunian’s 方法等6 種徐變計算方法的計算結果與徐變試驗曲線進行比較，結果發現， 按齡期調整有效模量法優于其它所有方法。 混凝土材料在長期服役過程中，徐變性能因受眾多因素影響，目前尚未形成統一、有效的徐變計算方法。

采用Kelvin 鏈模型將應變Volterra 積分方程轉化為率型徐變積分微分方程[18]，在此基礎上，引入中間變量離散連續時間， 結合應力應變增量彈性本構關系， 最終建立含初應變的彈性結構有限元分析模型。 最后使用該有限元分析模型，對大型商用有限元軟件Ansys 進行二次開發，實現混凝土徐變計算。

1 基于Kelvin 鏈率型模型的徐變理論

根據微觀力學理論[19]，混凝土在工作應力范圍內，單軸作用下徐變模型可依據Stieltjes 積分公式，用以下兩種Volterra 積分方程形式表示

式中：t 為混凝土開始澆筑至徐變計算時間；t′為加載齡期；σ 為混凝土應力；ε 為混凝土應變；ε0（t）為非彈性初始應變（如溫度引起的熱膨脹和收縮等除作用力以外產生的應變）；J（t，t′）為徐變函數（也稱為徐變柔度函數）， 定義為單位應力作用下，t-t′時間內引起的應變；ER（t，t′）為松弛函數（也稱為松弛模量）， 定義為單位應變作用下，t-t′時間內引起的應力；1/J（t，t′）=ER（t，t′）=E（t）表示為瞬時彈性模量。

徐變函數通過Dirichlet 序列變化[20]可以表示為

式中：τi為黏滯時間（事先給定值），其中i=1，2，…，M；E0為瞬時彈性模量；Di（t′）為齡期相關模量，可以由確切的柔度函數通過最小二乘法擬合得到，也可以由連續黏滯譜[21]確定，在后文給出。

將式（2）代入式（1）中，可得

其中

式（3）可以理解為總應變由M+1 個元件的應變累加而成；式（4）理解為編號為“0”、彈性剛度為E0的彈簧元件；式（5）理解為流變模型中的其他M 個Kelvin 黏滯元件。

將式（5）對t 進行微分，得到結合公式（5），可以得到以下積分微分方程

1） 考慮非老化情況，滿足Di（t′）=Di=常數，式（7）改寫為微分方程（率型方程）

由式（8）可得，總的應力可以分為兩部分。 第一部分為由剛度為Di彈性元件引起的應力；第二部分為ηi=τiDi由線性阻尼為黏滯元件引起的應力，該阻尼元件的應變率為ε˙i。 此時對應圖1 中的（a）非老化Kelvin 鏈模型。

圖1 Kelvin 鏈模型Fig.1 Kelvin chain model

2） 考慮老化情況。通常情況需要考慮材料的老化。 為了消除式（7）中的積分，再次對該式對t 進行微分，得到微分方程（率型方程）

同樣，可以采用式（9）來表征具有時變性能的Kelvin 元件（包含其中的彈性元件和阻尼元件）。 參考文獻[18]，考慮老化時，Ei（t）=Di（t）-τi（t），ηi=τi-Di。 此時對應圖1 中的（b）老化Kelvin 鏈模型。

式（8）和式（9）即為不考慮材料老化與考慮材料老化情況下的徐變計算微分方程形式（也稱為率型）。 在考慮材料老化情況下，為了解上述二階微分方程，文獻[18]提出中間變量使計算便捷有效，引入中間變量

將式（10）代入式（9），即可得到滿足初始條件γi（0）=0 的一階微分方程

在復雜應力狀態下，一階微分方程（11）可改寫為

一般來說，微分方程（11）式的解需要對每個時間步Δt（從時刻tn到tn+1）進行積分計算，分析任意步計算時，中間變量γi（tn）=γi（n）的值總是依賴于上一步（如果是第一步，則依賴于初始條件）。 根據文獻[18]分析，時間步Δt 的選取對計算效率影響很大。 當Δt≥τ1時，采用積分算法求解常微分方程所得到的結果是不穩定且精度不高，為此，需要對時間步Δt 尋找其它選擇方法。Zienkiewicz 等[22]在分析非老化材料時對時間步Δt 的選取， 曾提出了指數法，后來這種時間步Δt 選取的指數法被Bazant[23]延伸到了老化材料中。

指數算法的關鍵在于當Δt≥τ1時， 若微分方程（11）的右邊在任意時間步保持為常數，則微分方程（11）可以有解析解。 為簡便計算，我們先考慮單向應力問題，然后再拓展至復雜應力狀態問題。令Δt=tn+1-tn，Δσ=σ（tn+1）-σ（tn），Di（n+1/2）=Di（tn+Δt1/2），則式（11）右邊可以寫成Δσ/ΔtDi（n+1/2），即式（11）可以表示為

求解常微分方程（13），得到

其中，βi和λi為

對式（10）積分求出應變增量，根據Kelvin 鏈模型，對應變增量進行求和可以得到總的應變增量計算公式

式（16）右邊第一項理解為施加應力時引起的彈性應變增量，該彈性應變增量等于彈性應力增量乘以等效彈性柔量。 等效彈性模量表示為

式（16）右邊第二項理解為恒載引起的徐變應變增量

為方便有限元子程序編寫，式（16）可以改寫為

將式（19）代入公式（14），中間變量可以由下式更新計算

以上算法適用于單向應力情況，但是也很容易通過乘以各項同性材料矩陣CV應用至復雜應力情況（即三維情況）。

基于Kelvin 鏈模型的計算連續黏滯譜有[21]

式中：C（k）（kτi）為徐變函數C（t，tn-1/2）k 階導數；C（t，t′）=J（t，t′）-1/E0；當k=3 的時候，計算可滿足精度要求[21]。 將連續的黏滯譜離散得到離散的黏滯譜

2 Ansys 有限元軟件中Usermat 子程序的算法設計

基于Kelvin鏈率型的混凝土徐變計算方法式（19），采用64 位Ansys15.0+Vsiual Studio 2010+Inter Fortran Composer XE 2013 SP1 對用戶材料子程序Usermat 進行二次開發，算法流程圖如圖2 所示。具體計算步驟如下。

1） 說明Ansys 有限元軟件中Usermat 子程序所需要的常數、變量、變量數組。

2） 從Ansys 主輸入數據中讀取材料常數，包括彈性模量、泊松比、加載齡期等。

3） 在三維復雜應力情況下， 定義CV和雅克比矩陣D=E0CV。

4） 離散時間tn，n＝1，2，…；Δt=tn-tn-1；tn-1/2=t0+[（tn-t0）（tn-1-t0）]1/2； 當n=1 時，t1/2=（t0-t1）/2，t0為加載齡期，開始時間步循環步。

5） 當計算時間小于加載齡期t0時，按式（6）計算瞬時加載彈性變形，否則（即當t 超過加載齡期t0時）按式（7）計算持荷過程中徐變變形。

6） 初始化應力增量，按彈性結構計算初應變增量Δε。

7） 按照率型徐變擬線彈性應力應變關系計算任意時刻徐變變形：

8） 循環計算下一分析步，當增量步達到設置的最大值時結束計算。

圖2 基于率型徐變模型的有限元徐變分析算法流程圖Fig.2 Flowchart of finite element creep analysis algorithm based on rate-type creep model

3 標準算例

算例1 計算恒定軸壓荷載作用下素混凝土方塊徐變[24]。 素混凝土28 d 彈性模量E28=30 MPa。 混凝土齡期為7 d 時， 方塊頂面施加1 MPa 恒定應力，隨后保持荷載恒定30 年左右。 本節有限元模型數值解與及現有結果[24]對比情況如圖3（b）所示。 比較結果表明，三者結果幾乎一致。

圖3 素混凝土方塊有限元模型及計算結果對比圖Fig.3 Finite element model and calculation results comparison of plain concrete brick

圖4 素混凝土長柱有限元模型及計算驗證Fig.4 The finite element model and calculation validations of plain concrete long column

算例2 計算逐步施加軸壓荷載的素混凝土長柱徐變[24]。素混凝土強度fc28=30 MPa，相對環境濕度h=70%。 混凝土齡期為7 d 時，柱端面進行第一次加載，應力值為1 MPa；隨后每隔7 d 對柱端面疊加1 MPa 應力；疊加次數共12 次；疊加完成后，保持荷載恒定約30 年。 有限元計算模型如圖4（a）所示，本節有限元模型數值解與及現有結果[24]對比情況如圖4（b）所示。 比較結果表明，三者結果幾乎一致。

4 結論

1） 從微觀力學角度出發， 建立Kelvin 鏈率型徐變模型。 基于該模型對Ansys 有限元商用軟件實施二次開發， 修正了Ansys 軟件中材料子程序Usermat，最后通過經典算例驗證了本文徐變計算模型的正確性和有效性， 彌補了Ansys 軟件中徐變模型的不足。

2） 提出的徐變模型也適用于其它有限元商用軟件， 如：Abaqus，Msc.Marc 和ADINA 等的二次開發。