橋梁徐變效應的有限元算法研究及程序驗證

王雨權,廖立堅,李林安,霍學晉,李 黎

(1.中國鐵路設計集團有限公司,天津 300308； 2.天津大學力學系,天津 300350； 3.中鐵大橋勘測設計院集團有限公司,武漢 430050； 4.中鐵第六勘察設計院集團有限公司,天津 300308)

在應力狀態下，徐變是混凝土的固有特性。混凝土長期徐變會引起橋梁結構發生主要受力方向的位移，因此，徐變計算是橋梁，尤其是大跨度橋梁結構的控制性因素[1-5]。設計及研究人員對于徐變引起的各種不利影響進行了深入研究，如徐變引起CRTSⅡ型縱連板式軌道和橋面之間的脫空現象[6]，徐變對剪力滯效應影響[7]等。因此，鑒于徐變的重要性，很多學者也對如何更好地模擬徐變效應進行了深入研究[8-12]。

人們可通過徐變系數來描述徐變現象，查閱相關文獻可知，徐變系數的影響因素很多，其表達式的模型也很多[13-15]，國內公路規范[16]和鐵路規范[17]代表了兩種表達形式。

如何利用規范上的徐變系數進行精確合理地計算混凝土橋梁的徐變效應，是一個較難確定的復雜問題。由于商業軟件Midas Civil操作簡便性，在一般設計中，設計人員很少去深究徐變到底是如何計算的。而對于Midas Civil來說，徐變計算模塊是其核心功能模塊，運算機理一直處于保密狀態。因此，設計人員對于該模塊的計算過程其實是處于一個黑盒子狀態[18-20]。這對于準確把握設計精髓，進一步提高橋梁設計水平有阻礙作用。因此，進行徐變算法的相關研究工作十分必要。

從徐變系數出發，以Midas Civil幫助文件中涉及到的算法為參照，深入研究相關計算方法，并進行相應的公式推導。在確定算法路徑后，詳細給出了遞推格式的計算公式和步驟，并進行了相關案例測試。

1 徐變系數

影響混凝土徐變的因素很多，但主要因素有混凝土強度、受荷時的混凝土齡期、周邊環境相對濕度、溫度、混凝土體積、作用于混凝土結構的持續應力大小等。為計算徐變引起的變形，引入徐變系數φ(t，τ)的概念。

徐變系數φ(t，τ)可表示為從時刻τ開始作用一個應力σ(τ)后，混凝土在時刻t時的單軸應變與加載齡期τ時的應變之比。即

(1)

式中，εcr(t，τ)為齡期為τ的混凝土在時刻t時應變；εe(τ)為加載時刻的彈性應變。

目前，國內外對徐變效應考慮的因素不盡相同，采用了不同的數學表達方式，但歸納起來有以下兩大類。

一是將徐變系數表達為一系列系數的乘積，每一個系數表示一個影響因子。如BS5400英國規范(1984年版)、美國ACI209標準(1982年版)、CEB-FIP標準規范(1990年版)、JTG 3362—2018《公路鋼筋混凝土及預應力混凝土橋涵設計規范》。

二是將徐變系數表達為若干個性質互異的分項系數之和。如CEB-FIP標準規范(1978年版)、JTJ 023—85《公路鋼筋混凝土及預應力混凝土橋涵設計規范》(1985年版)、TB 10092—2017《鐵路橋涵混凝土結構設計規范》等。

目前，國內有關橋梁徐變系數的相關規定，主要采用公路和鐵路兩個行業規范。TB 10092—2017《鐵路橋涵混凝土結構設計規范》(以下簡稱“《2017年鐵路橋規》”)附錄A給出了徐變系數的相關規定。該規范沿用了歷年橋梁規范規定，徐變系數的計算表達式中，存在很多圖表。JTG 3362—2018《公路鋼筋混凝土及預應力混凝土橋涵設計規范》(以下簡稱“《2018年公路規范》”)附錄C也給出了徐變系數的相關規定，該規范在上一版規范的基礎上，增加了添加粉煤灰的混凝土徐變修正系數。對于一般混凝土材料，仍然維持上一版規定。下面對這兩本規范的內容進行簡要介紹。

1.1 公路規范相關規定

《2018年公路規范》的徐變系數采用乘積形式，具體表達式如下。

(2)

式中，t0為加載時的混凝土齡期，d；t為計算考慮時刻的混凝土齡期，d；φ(t，t0)表示加載齡期為t0，計算考慮齡期為t時的混凝土徐變系數；φ0為名義徐變系數；βc(t-t0)為加載后徐變隨時間發展的系數；βH為與年平均相對濕度相關的系數；fcm為強度等級C20～C50混凝土在28 d齡期時的平均立方體抗壓強度，MPa；fcm0取值10 MPa；RH為環境年平均相對濕度，%；RH0取值100%；h為構件理論厚度，mm，h=2A/μ，A為構件截面面積，μ為構件與大氣接觸的周邊長度；h0取100 mm；t1取1 d。

1.2 鐵路規范相關規定

《2017年鐵路規范》的徐變系數采用分項系數相加的形式，具體表達式如下。

(3)

式中，φ(t，τ)為徐變系數；βd(t-τ)為隨時間而增長的滯后彈性應變，取值查規范附錄圖A.0.2-1；fτ/f∞表示加載時混凝土的強度與最終強度R∞之比，查規范附錄圖A.0.2-2；φf為流塑系數，φf=φf1×φf2；φf1為依周圍環境而定的系數；φf2為依理論厚度h而定的系數，查規范附錄圖A.0.2-3。

兩本規范的表達式相差很大，公路規范的表達式明確，而鐵路規范的相關系數取值依賴于圖表的查詢，這容易產生徐變系數取值誤差。

2 徐變算法對比研究

混凝土徐變的計算方法大致有Dischinger微分求解法、Tr?st-Bazant代數求解法和有限元逐步增量法。對于結構形式復雜、施工步驟較多的橋梁結構，大多采用有限元逐步增量法。

目前，在設計過程中，由于商業軟件較為發達，以Midas Civil為代表的橋梁結構有限元軟件，已成為設計人員從事橋梁設計的必備工具。實踐工程中，Midas Civil設計的橋梁在最近十幾年時間內，也經受住了考驗。因此，在研究徐變算法時，以與Midas Civil的計算對比結果為基準，但Midas Civil并未給出詳細具體的計算算法，只能從其僅有的幫助文件或相關手冊內容進行研究，來論證可能的徐變計算方法。

下面分別對Midas Civil給出的徐變計算模型和常見的基于Tr?st-Bazant方法給出的工程實用分析基本方程的邏輯關系進行論證。并基于論證結果，給出具體的徐變計算算法。

2.1 Midas Civil相關算法研究

在Midas Civil的幫助文件及相關技術手冊中，該軟件給出了混凝土的應變方程，即

ε(t)=εe(t)+εc(t)

(4)

式中，ε(t)為混凝土總應變；εe(t)為混凝土彈性應變；εc(t)為混凝土徐變應變；t為徐變計算時刻。

當Midas Civil采用徐變系數計算方法時，其混凝土的徐變應變表達式為

(5)

式中，σe(t)為混凝土加載的初始應力；Δσ為每一個加載持續時間段內的應力增量；E為混凝土彈性模量；φ(t，τ0)為徐變系數；τ0為混凝土的加載齡期。

式(4)、式(5)是Midas Civil中混凝土材料徐變計算中應力應變增量關系，以這兩個表達式為基礎，將其轉換成全量形式，過程如下。

由式(4)可寫出對應的應力表達式

σ(t)=σe(t)+Δσ

(6)

式中，σ(t)為混凝土計算時刻的總應力；Δσ為每一個加載持續時間段內的應力增量。

聯立方程(4)～(6)可得

[1+φ(t，τ0)/2]

(7)

對式(7)進行展開，由此可得

(8)

由于σe(t)=E×εe(t)，對式(8)進行整理可得

(9)

對式(9)兩邊同乘E并移項，可得

(10)

該表達式即為從Midas Civil的幫助文件給出的混凝土應變方程出發，得到的混凝土徐變應力應變全量關系表達式。

將式(10)代入伽遼金平衡弱解方程，即可利用相關的弱解方程求解工具進行求解。

伽遼金平衡弱解方程的通式為

(σ;δε)=(f;δu)

(11)

將式(10)代入后可得

(12)

2.2 基于Tr?st-Bazant的工程實用分析算法研究

在國內外可查詢的相關資料中，混凝土徐變的應力-應變方程一般表述為

式中，ε(t)、σ(t)分別為混凝土計算時刻的總應變和總應力；E(τ0)為混凝土彈性模量；φ(t，τ0)為徐變系數；εe(τ0)為彈性應變；εs(t)為徐變應變；τ0為加載齡期；t為徐變計算時刻。

采用Tr?st-Bazant方法給出的工程實用分析基本方程為

(14)

式中，εs(t)為徐變應變；σs(t)為徐變應力；σ(τ0)為彈性應力；E為彈性模量；φ(t，τ0)為徐變系數；ρ(t，τ0)為老化系數。

下面以式(13)、式(14)為基礎，推導出對應的混凝土徐變應力應變全量關系表達式。

全量形式的徐變應變和徐變應力可表示為

(15)

將式(15)代入Tr?st-Bazant基本方程(14)可得

[1+ρ(t，τ0)φ(t，τ0)]

(16)

對式(16)兩邊進行整理可得

(17)

對式(17)兩邊同乘E并移項，可得

[1+ρ(t，τ0)φ(t，τ0)]σ(t)=Eε(t)-

σ(τ0)φ(t，τ0)[1-ρ(t，τ0)]

(18)

同理將式(18)代入伽遼金平衡弱解方程的通式(σ;δε)=(f;δu)，可得

([1+ρ(t，τ0)φ(t，τ0)]-1Eε(t);δε)=(f;δu)+

-(σ(τ0)φ(t，τ0)[1-ρ(t，τ0)]×

[1+ρ(t，τ0)φ(t，τ0)]-1;δε)

(19)

對比式(12)和式(19)可知，這兩個表達式只要滿足ρ(t，τ0)=1/2時，其方程是等價的。

因此，徐變計算要取得與Midas Civil一致的結果，只需按照基于Tr?st-Bazant基本方程推導而來的計算方法即可。

3 徐變的增量遞推算法

基于Tr?st-Bazant基本方程的徐變算法，國內相關學者[21]已進行較多研究，下面根據相關文獻總結出橋梁結構基于空間梁單元的徐變計算增量算法。

3.1 各子步徐變位移表達式

圖1 徐變計算應力與時間的關系

根據徐變的定義，時刻tn的徐變應變可表示為

(20)

時刻tn-1的徐變應變為

(21)

利用積分中值定理，可得第n個階段即tn-1→tn的徐變應變為

式中相關變量的表達式為

(23)

式中，γ(ti，ti-1)=1/[1+ρ(ti，ti-1)φ(ti，ti-1)]。

另外，根據虛功原理和Castigliano第一定理，可得到如下表達式

{F}=γ(t，τ0)[K]{δ}-γ(t，τ0)φ(t，τ0){F0}=

γ(t，τ0)[K]{δ}-γ(t，τ0)φ(t，τ0)[K]{δ0}=

[Kφ][{δ}-φ(t，τ0){δ0}]=

(24)

將上述平衡方程與式(20)結合，可得第n階段(tn-1→tn時段)的徐變位移為

(25)

式中，彈性位移{δ0i}由初始階段的位移法分析得到；總徐變位移{δi}由結構總體平衡方程解出。

3.2 徐變計算的遞推公式

直接利用式(25)來計算徐變位移并不方便，因為每一步都需存儲前一階段的中間參數{δ0i}、{δi}等參數，因此，需結合徐變系數構造一個遞推格式，以便快速計算徐變效應。

設徐變系數的通用公式為

(26)

式中，βa(τ)、βd(0)為徐變系數擬合時與初始時刻0、加載時刻τ有關的擬合系數；Ci(τ)為Dirichlet級數展開時與加載時刻τ相關的徐變系數形狀系數；qi為Dirichlet級數展開時與隨時間變化的徐變系數曲線形狀相關的系數。

彈性位移{δ0i}由初始階段位移法分析得到，總徐變位移{δi}由結構總體平衡方程解出。

將式(26)代入上述各表達式，并定義如下變量

(27)

則第n階段(tn-1→tn時段)的徐變位移表達式可表示為

{δ0n-1}[βa(tn-1)+βd(0)]

(28)

γ(tn-tn-1)Ci(tn-1/2)e-qi(tn-tn-1/2)

(29)

因此，徐變計算表達式即可用式(28)、式(29)對應的遞推公式表示。

接下來，只需將本文所述的公路規范和鐵路規范的徐變系數擬合成式(26)所示的通用表達格式，即可自動處理不同規范的徐變計算問題。

對于自定義的徐變系數，可通過如下加權最小二乘法進行擬合，加權最小二乘的格式如下

(30)

式中，w(xi)為加權值；s(xi)、f(xi)分別為擬合點處的擬合值及規范計算值。

3.3 徐變遞推計算基本步驟

在進行徐變效應分析時，可按以下3個步驟分別進行分析。

第一步：對徐變系數進行擬合。

式(2)、式(3)分別展示了公路和鐵路規范的徐變系數，但程序計算時，不能用該表達式進行計算，均需將其轉化為式(26)的公式。因此，本階段的主要任務即為用加權二乘法對徐變系數進行處理。

第二步：對n=1階段進行處理，主要按如下過程進行處理。

(1)計算初始階段的彈性變形。

(5)組集總剛和荷載列陣，列出總體平衡方程式，處理邊界條件后解出{δ1}。

(6)計算徐變引起的總桿端力[F1]和約束反力。

第三步：對n=2以后的階段進行遞推計算處理，主要按如下過程進行處理。

(1)計算現階段的彈性變形{δ01}。

(5)組集總剛和荷載列陣，列出總體平衡方程式，處理邊界條件后解出{δ2}。

(6)計算徐變引起的總桿端力[F2]和約束反力。

(7)繼續循環進入下一階段計算。

按照上述3個步驟，即可編制有限元程序，進行徐變位移計算。

最終的徐變位移需將每一施工階段計算得到的位移增量進行疊加。計算中，遇到體系轉化，只需將由此導致的二次力疊加入下一個階段即可。

筆者在研制自主有限元軟件“臥龍”時，按照上述計算方法和過程進行了徐變計算的探索，并對其結果的準確性進行了案例測試。下面給出基于上述算法的“臥龍”與Midas Civil的驗證案例。

4 徐變位移計算程序驗證

某三跨連續梁，梁跨組合為94.5 m+166.4 m+107.9 m，總梁長368.8 m。同時模擬樁基礎、承臺、墩身的施工過程，其對應混凝土為C40、C30混凝土，這部分混凝土不考慮收縮徐變效應。梁體采用C60混凝土并考慮收縮徐變效應，梁體截面為單箱單室截面。鋼束采用Strand1860，鋼束特性值有3組，鋼束數量有238組，施工階段共分40個階段。靜力荷載包括自重、掛籃、鋼束張拉力、二期恒載等，其計算模型如圖2所示。

圖2 三跨連續梁測試模型

測試過程，徐變系數采用《2018年公路規范》提供的徐變系數表達式。

采用“臥龍”與Midas Civil背靠背對比計算模式。計算過程中，考慮了鋼束的二次力、鋼束引起的截面特性變化影響，也考慮了梁體變截面的影響。由于施工階段比較多，限于篇幅，直接選取成橋后10年的數據進行對比。

分別提取兩個軟件的徐變位移，數據如表1所示，由于空間梁單元的每個節點有6個自由度，但對于連續梁來說，其主要受力為自重及鋼束，因此，可選取豎向和水平兩個方向的徐變位移作為對比自由度。由于全橋節點較多，下面僅選取第29～51號單元進行比較，這些單元代表了中跨靠近第一個橋墩這半個跨度的橋梁結構，具有典型性。

表1給出了“臥龍”與Midas Civil在豎向和水平向兩個方向的徐變位移對比結果。表中，Dzm、Dxm分別為Midas Civil計算的豎向和水平向徐變位移；Dzw、Dxw分別為臥龍計算的豎向和水平向徐變位移。從表1中可以看出，豎向徐變位移的最大誤差為3.3%，該節點號為29號，位于第一個主墩節點上。該節點的徐變位移為上拱狀態，與其余節點的徐變位移方向相反，相對來說，該點徐變絕對值較小。中跨跨中位置節點號為51號，其相對誤差為1.19%。所有節點的水平向徐變位移，其誤差都在1%以內。

表1 自編程序臥龍與Midas的徐變位移對比

將表中兩個方向的徐變計算值分別繪制成對比圖，如圖3所示。

圖3 徐變位移值對比

從圖3中也可以看出，不論是豎向徐變位移值還是水平向徐變位移值，“臥龍”和Midas Civil的結果都是非常相近的。

“臥龍”與Midas存在誤差的原因為：采用遞推格式時，需對徐變系數進行擬合，其擬合數值和規范表達式計算出的精準數值是存在誤差的；Midas的徐變系數計算值與規范精準值之間也會存在誤差。但總體來看，“臥龍”與Midas Civil之間的誤差是可以接受的，這也表明，上述基于遞推形式的徐變計算方法是可行的。

5 結論

通過系統梳理徐變系數的相關規定，給出了現行公路和鐵路規范的相關規定。針對Midas和基于Tr?st-Bazant基本方程的兩種徐變算法，對其進行全量形式的公式推導和對比分析。隨后，給出遞推形式的徐變計算方法及計算步驟，并通過一個較為復雜的連續梁案例進行測試，結論如下。

(1)現行鐵路規范和公路規范的徐變系數表達式區別很大，可通過加權最小二乘法將其轉化成徐變系數設定的通用格式，從而統一徐變計算過程。

(2)Midas幫助文件給出的徐變算法和基于Tr?st-Bazant基本方程的徐變算法本質上是一致的。當滿足ρ(t，τ0)=1/2時，兩者之間的方程是等價的。

(3)將徐變系數寫成設定的通用格式，可將基于空間梁單元的徐變算法改造成遞推形式，遞推形式的徐變計算過程將大大減少運算量和數據存儲量，加快運算速度。

(4)實例測試結果表明，在充分考慮鋼束二次力、結構體系轉化影響、變截面性質及鋼束引起的截面特性變化情況下，本文所述徐變算法計算得到的徐變位移和Midas計算得到的徐變位移結果非常接近。

(5)測試結果表明，自編程序“臥龍”的徐變計算精度，在成橋后10年，豎向徐變位移的最大誤差為3.3%，且徐變絕對值較小，水平向徐變位移誤差均在1%以內，計算精度良好。