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

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

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

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

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

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

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

1 徐變系數(shù)

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

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

(1)

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

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

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

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

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

1.1 公路規(guī)范相關規(guī)定

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

(2)

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

1.2 鐵路規(guī)范相關規(guī)定

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

(3)

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

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

2 徐變算法對比研究

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

目前，在設計過程中，由于商業(yè)軟件較為發(fā)達，以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采用徐變系數(shù)計算方法時，其混凝土的徐變應變表達式為

(5)

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

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

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

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

(6)

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

聯(lián)立方程(4)～(6)可得

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

(7)

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

(8)

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

(9)

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

(10)

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

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

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

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

(11)

將式(10)代入后可得

(12)

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

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

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

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

(14)

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

下面以式(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的徐變應變?yōu)?/p>

(21)

利用積分中值定理，可得第n個階段即tn-1→tn的徐變應變?yōu)?/p>

式中相關變量的表達式為

(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)來計算徐變位移并不方便，因為每一步都需存儲前一階段的中間參數(shù){δ0i}、{δi}等參數(shù)，因此，需結合徐變系數(shù)構造一個遞推格式，以便快速計算徐變效應。

設徐變系數(shù)的通用公式為

(26)

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

彈性位移{δ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)對應的遞推公式表示。

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

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

(30)

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

3.3 徐變遞推計算基本步驟

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

第一步：對徐變系數(shù)進行擬合。

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

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

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

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

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

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

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

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

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

(7)繼續(xù)循環(huán)進入下一階段計算。

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

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

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

4 徐變位移計算程序驗證

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

圖2 三跨連續(xù)梁測試模型

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

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

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

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

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

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

圖3 徐變位移值對比

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

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

5 結論

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

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

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

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

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

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