幾何非線性平面梁考慮收縮徐變的算法研究*

鄧繼華，邵旭東，彭建新

(1．湖南大學 土木工程學院, 湖南 長沙 410082；2．長沙理工大學 土木與建筑工程學院, 湖南 長沙 410076)

隨著結構計算理論、高強材料及施工裝備的快速發展，混凝土斜拉橋等大跨柔性混凝土結構在跨度和高度上的記錄不斷被刷新,目前國內已建成的蘇通長江大橋主跨為1 088 m(主梁為鋼結構)，塔高達300 m(塔為混凝土結構).分析此類既有混凝土構件又比較柔性的結構，同時考慮幾何非線性和收縮徐變是很有必要的[1-2]，目前對幾何非線性平面梁單元已進行了大量研究[3-5]，理論已經很成熟，對混凝土結構的收縮徐變效應分析也有很多研究成果[6-8]，但對混凝土柔性結構同時考慮幾何非線性和收縮徐變效應的研究文獻非常少，筆者僅找到兩篇文獻[9-10]，文獻[9]介紹了幾何非線性結構進行徐變效應分析的原理和方法，提供了按施工階段進行幾何非線性結構徐變分析的增量形式和分析步驟．但該文采用按齡期調整的有效彈性模量法來求解結構的徐變問題，由于該方法既要形成結構的彈性剛度矩陣，又要形成徐變剛度矩陣，在每一迭代步內結構的平衡方程還要求解兩次，這在計算量本來就很大的非線性計算中是很不合適的．同時該文算法是基于增量平衡，在進行下一個增量步(施工階段)等效荷載列陣計算時未考慮上一個增量步(施工階段)的殘余力，更沒考慮由于徐變變形而導致上一個增量步(施工階段)已經平衡的構形被打破在下一個增量步(施工階段)會產生新的不平衡力，因此隨著施工階段數的增多誤差會越來越大．文獻[10]從虛功增量方程出發，建立了桿系結構的非線性與混凝土收縮徐變效應耦合分析的有限元方法．該文采用初應變法來計算分階段施工混凝土結構的徐變、基于全量平衡且考慮了由于徐變變形而導致上一個增量步(施工階段)已經平衡的構形被打破在下一個增量步(施工階段)會產生新的不平衡力問題．因此，從理論上講，相對于文獻[9]的算法而言，文獻[10]的算法無論在非線性計算效率還是精度方面均有提高．但仔細研究也會發現該文存在以下問題，該文建立的非線性平衡方程本質上是基于U.L列式的增量平衡方程，即計算單元切線剛度矩陣是基于每一迭代步的增量．但為了考慮由于徐變變形而導致的平衡構形被打破及提高計算精度對不平衡力計算又采用全量方法，顯然兩者不一致會導致各種狀態變量數增加從而使計算量增加的缺點．且該文是從張量分析出發，推導過程也較復雜，較難為工程技術人員理解和應用.鑒于此，本文在參考上述已有文獻的基礎上，首先基于微分法導出了隨轉坐標系下平面梁在大轉動小應變時計入初應變效應的幾何非線性平衡方程,結合用初應變法進行節段施工混凝土梁計算中收縮徐變等效節點力計算的有限元列式,再利用隨轉坐標系與結構坐標系下節點力之間和節點位移之間全量及增量的關系，最終獲得結構坐標系下平面梁單元幾何非線性分析中考慮收縮徐變效應影響的實用算法,給出了詳細的計算步驟，最后對文獻[10]的算例1進行了比較分析.

1 隨轉坐標系下梁單元的平衡方程

圖1(a)與(b)示意了在結構坐標系xy中初始時刻和計算t時刻混凝土梁元的幾何參數及位移的即時變量．設混凝土梁元的隨轉坐標系為XY，該坐標系是隨單元變形而轉動的，它始終以節點i為原點，以節點i到j的連線方向為X軸，由X軸逆時鐘轉90°為Y軸.

圖1 變形前后梁單元

由圖1(a)與(b)知：

(1)

混凝土梁單元在隨轉坐標系中的節點位移為：

(2)

(3)

在隨轉坐標系下單元的應變-位移關系只需考慮線性項[11]，有：

εe=B0dl

(4)

式中εe為節點位移引起的彈性應變;B0為線性應變矩陣.

根據廣義虎克定理，在單元內部具有初應變的情況下應力表達式為：

σ=D(εe-ε0)

(5)

式中ε0為混凝土收縮、徐變而引起的初應變.

基于虛功原理有：

(6)

將式(5)代入式(6)，并聯立式(4)得：

(7)

對于式(7)，顯然右邊第一項積分就是平面梁單元在局部坐標系下的線彈性剛度矩陣k0，為求簡便，可將式(7)改寫為：

f=fe+f0

(8)

對式(7)微分有：

δf=k0δdl

(9)

2 混凝土收縮、徐變的增量法計算

2.1 混凝土徐變計算

在常應力σ作用下，以徐變系數來表示混凝土徐變應變的表達式為：

(10)

式中εc為時刻τ始作用于混凝土的常應力σ至時刻t所產生的徐變應變;D為應力矩陣;φ(t,τ)為加載齡期為τ、計算齡期為t時的混凝土徐變系數.

由于大跨徑橋梁都是采用分節段多階段施工完成，是一個變結構變應力的問題，一般采用增量法的形式來進行徐變分析，將整個求解時間軸分為t0,t1,…,tn,tn+1等時刻，在各時刻作用有Δσ0,Δσ1,…,Δσn,Δσn+1等應力增量，下面以在tn時刻，在Δtn+1即tn→tn+1時段內為例，介紹基于初應變法的單元徐變等效節點力增量計算方法.

該時段內，在結構的應力和材料常數不隨時間而變化的假定下，由式(10)得徐變增量為：

(11)

式中Δσi為時間軸上ti時刻的應力增量.

(12)

Δσi=D(Δεe,i-Δεc,i)=DΔεe,i-DΔεc,i

(13)

式中Δεe,i如前述為節點位移增量引起的應變增量;Δεc,i為混凝土徐變(上一時步)初應變增量.

Δfe,i+Δfc,i=Δfi

(14)

將式(14)代入式(12)得：

(15)

(16)

式(16)中Ci,qi,(i=1,2,3)等6個待定系數的具體值可參見文獻[12].

(17)

將式(16)代入式(17)可得：

(18)

令

(19)

則式(18)可化為：

(20)

又有：

(21)

在不考慮幾何非線性的情況下，用初應變法計算多階段施工橋梁徐變效應的具體步驟如下．

首先是輸入基本參數，確定時間軸，然后進行時步循環，時步循環結束則徐變計算結束，其中在時步的每次循環中又包括7個主要步驟按順序依次為[14]:

2)計入上一時步的徐變等效結點荷載，同時將其反號疊加至Δfc,i中作為徐變初應變影響部分；

3)結構剛度矩陣計算；

4)求解總剛，獲得結點位移增量；

6)結點位移和桿端內力累加；

7)按式(15)計算下一時步的徐變等效結點荷載.

2.2 混凝土收縮計算

類似徐變過程分析，得到Δtn+1(tn-tn+1)時段內收縮應變增量為：

(22)

式中εs(∞)為收縮應變終極值;p為收縮應變增長速度系數;τ0為混凝土的硬化時間.其他符號意義同前，收縮應變增量引起的等效節點力增量為：

(23)

對式(23)的計算可采取與徐變計算相同的方法通過遞推得到，此處不贅述.

在程序中，對于收縮按初應變的考慮，計算方法與上述徐變計算步驟完全一樣.

3 結構坐標系下梁單元的平衡方程

聯立式(8),(14),(15)及(23)，顯然有：

fe=∑Δfe,i

(24)

F=tT·f

(25)

式中t為結構坐標系與隨轉坐標系之間的轉換矩陣，一般桿系有限元書上均可查到.

將式(25)兩邊微分可得：

δF=δtT·f+tT·δf

(26)

結合式(1)對式(26)中矩陣t微分，不難得到δtT用計算t時刻結構坐標系中的位移向量微分δd表達的形式，據此可將δtT·f寫成用δd表達的形式為：

δtT·f=kn·δd

(27)

δdl=T·δd

(28)

聯立式(27),(9)及(28)，可將式(26)改寫成：

δF=(kn+tTk0T)δd

(29)

因此，混凝土梁單元在結構坐標系下非線性切線剛度矩陣K為：

K=kn+tTk0T

(30)

從上述推導過程可得到僅考慮幾何非線性時計算步驟如下：

1)根據單元在結構坐標系下的總位移向量d由式(2)形成隨轉坐標系下的位移向量dl；

2)通過位移向量dl由式(8)形成單元在局部坐標系下的節點力向量f(此時因不計初應變有f0=0)，再由結構當前構形的幾何參數由式(27)形成kn；

3)kn結合式(25),(28)和k0由式(30)形成結構坐標系下單元的切線剛度矩陣K，計算涉及到的幾何參數都是結構當前構形下的值；

4)將單元在局部坐標系下的節點力向量f由式(25)轉換到結構坐標系下得到F；

5)對所有單元重復1)至4)的步驟，生成結構切線剛度矩陣[K]=∑K和節點合力[F]=∑F；

6)計算不平衡力ΔR=P-[F]，其中P為截止到計算施工階段施加的總外荷載所轉換成的等效節點力；

7)求解結構方程[K]·Δd=ΔR，得到節點位移增量Δd，并將其疊加到總位移向量d中；

8)收斂條件判斷，如收斂，則轉到下一施工階段計算，如不收斂，則返回1)，進行下一次迭代計算.

4 同時考慮幾何非線性與收縮徐變的計算步驟

從前面的整個推導過程及幾何非線性、收縮徐變效應單獨分析時的計算步驟，可知同時考慮幾何非線性和收縮徐變效應時的基本計算原理為：將不平衡力(由總外荷載減去節點位移產生的結構總抗力及收縮徐變等效節點力總量得到)作用下經過幾何非線性分析得到的變形作為每一工況初始瞬時彈性變形，徐變變形與之成線性關系，將上一階段的收縮徐變變形作為下一階段的初應變，計算得到的等效節點力增量總是作用在單元隨轉坐標系中，因此在各階段的累加計算過程中無須考慮坐標系的轉換，可直接累加形成收縮徐變等效節點力總量，將收縮徐變等效節點力總量由結構當前構形下的隨轉坐標轉換到結構坐標系下就可參與不平衡力的計算.

由以上分析，可知只要注意以下4項就可實現對僅考慮幾何非線性的程序稍加改造而獲得同時考慮幾何非線性和收縮徐變效應程序的目的.

1)在計算kn時涉及到的f須計入f0，而f0就是本階段前所有階段的收縮徐變等效節點力總量；

2)不平衡力由總外荷載減去節點位移產生的結構總抗力及收縮徐變等效節點力總量得到；

從上述計算步驟可看出：非線性方程的不平衡力是完全基于全量平衡計算的，能保證每個施工階段結構處于平衡狀態且不累積誤差.

5 算例分析

為驗證本文算法及計算程序，以單獨的幾何非線性考題、單獨的收縮徐變考題以及同時考慮幾何非線性與收縮徐變的考題進行了計算，篇幅所限，僅以文獻[10]中算例1的比較分析來說明.

算例 如圖2所示懸臂梁，L=300 m，彈性模量E=3.45×104MPa，面積A=54.906 m2， 慣性矩I=450.084 m4，自由端施加一水平集中力Px，豎直集中力Py，收縮徐變計算時間為1 000 d，加載齡期為10 d.分4 種工況計算：線性、線性+混凝土徐變、幾何非線性和幾何非線性+混凝土徐變.在水平集中力Px=1 200 kN，豎直集中力Py=50 000 kN作用下，按上述4種工況進行計算.在計算中，將時間軸上從10 d到1 000 d按文獻[14]所述采用對數函數來劃分為20個時步長，分別為10，13,16,20,25,32,40,50,63,79,100,126,159,200,251,316,398,501,631,794,1 000.

圖2 懸臂梁模型

計算結果如表1、圖3和圖4所示.

表1 梁自由端的位移與梁根部彎矩

從表1、圖3和圖4可看出，在不考慮幾何非線性的情況下，由于為一次形成的靜定結構，徐變不產生次內力，固計徐變與不計徐變的固定端彎矩相同，兩者位移值之比為3.27，該值實際上就是基于已知徐變參數基礎上的加載齡期為10 d,計算齡期為1 000 d的徐變系數值；在不考慮徐變情況下，基于幾何非線性和線性的水平位移比值為1.131，豎向位移比值為1.165，固定端彎矩比值為1.109，但考慮幾何非線性與徐變效應耦合后，無論是水平、豎向位移還是固定端彎矩比值，變化均很大，基于幾何非線性與徐變效應耦合和計入徐變效應的線性值在水平位移方面比值為1.472，在豎向位移方面比值為1.876，固定端彎矩比值為1.520；同時還可看出基于幾何非線性與徐變效應耦合的計算值并不是基于幾何非線性的計算值與徐變系數的簡單乘積，兩者差別較大，表明兩者之間存在明顯的耦合作用.

時間軸/d

時間軸/d

時間軸/d

還應指出的是，對于表2中最后一行的計算值，如按變形后受力平衡條件來計算梁的根部彎矩，其值應是527 636.68 kN·m，與表列的程序計算值547 048.6 kN·m相差3.68%，仔細研究后發現這是正常的非線性迭代計算誤差，如將時間軸上從10 d到1 000 d按1 d的間距來劃分成990個步長，相應的計算值為：懸臂梁端水平位移為3 455.7mm，豎向位移為51.0mm，梁的根部彎矩為532 862.22 kN·m，再按變形后受力平衡條件計算的梁根部彎矩為532 725.03kN·m，兩者相差0.026%，完全可以忽略不計.以上結論與國際材料和結構試驗室聯合會(PJLEM)提供的關于混凝土徐變行為的8個經典試驗中一細長偏心受壓混凝土柱徐變試驗結論[16]是基本一致的.

6 結 論

1)基于隨轉坐標系有關理論及收縮徐變分析的初應變法建立了同時考慮幾何非線性和收縮徐變效應的非線性平衡方程，為混凝土斜拉橋等大跨柔性混凝土結構分析打下了理論基礎.

2)本文算法采用外荷載、由于節點位移而產生的結構抗力、混凝土收縮徐變產生的等效節點力的全量來計算節點不平衡力，故能有效消除時步間誤差累積的問題．

3)算例表明，對于既有混凝土構件又比較柔性的結構進行幾何非線性和收縮徐變效應藕合分析是很有必要的，只考慮幾何非線性或考慮線性和徐變共同作用的計算結果與之差別很大.

4)按本文算法，可以方便地對現有平面桿系程序進行改造，使之具有同時考慮幾何非線性和收縮徐變效應的計算功能，使工程技術人員能更好地理解大跨柔性混凝土結構的受力行為.

5)由于同時考慮幾何非線性與收縮徐變效應共同作用的理論研究及結構試驗極少，而工程實際的發展已經要求開展這方面的理論研究及試驗驗證.

[1]MARU S,ASFAW M,SHARMA R K,etal. Effect of creep and shrinkage on RC frames with high beam stiffness[J].Journal of Structural Engineeirng, ASCE, 2003, 129(4):536-543.

[2]EDWARDS N P,BILLINGTON D P. FE analysis of tucker high school roof using nonlinear geometry and creep[J]. Journal of Structural Engineeirng,ASCE,1998, 124(9): 984-990.

[3]蔡松柏,沈蒲生.大轉動平面梁有限元分析的共旋坐標法[J].工程力學,2006,23(增I): 69-72.

CAI Song-bai, SHEN Pu-sheng. Co-rotational procedure for finite element analysis of plane beam element of large rotational displacement[J]. Engineering Mechanics, 2006, 23(z1): 69-72. (In Chinese)

[4]呂和祥, 朱菊芬, 馬莉穎. 大轉動梁的幾何非線性分析討論[J]. 計算結構力學及其應用, 1995,12(4):485-490.

LV He-xiang, ZHU Ju-fen, MA Li-ying. Discussion of analyzing of geometric non-linear beams with large rotations[J]. Chinese Journal of Computational Mechanics, 1995, 12(4): 485-490. (In Chinese)

[5]鄧繼華,邵旭東. 基于U.L列式的帶剛臂平面梁元非線性分析[J]. 湖南大學學報:自然科學版, 2012,39(5):8-12.

DENG Ji-hua, SHAO Xu-dong. Non-linear analysis of plane beam element with rigid arms based on U.L formulation[J]. Journal of Hunan University: Natural Sciences, 2012,39(5): 8-12. (In Chinese)

[6]TAYLOR R L，PISTER K S，GOUDREAU G L．Thermomechanical analysis of viscoeoelastic solids[J]．Num Meth Eng ，1970(2):45-60.

[7]方志，黃立葵，周樂農．砼結構徐變的灰色系統預測[J].湖南大學學報:自然科學版,1994，21(4):96-102.

FANG Zhi, HUANG Li-kui,ZHOU Le-long. Creep deformation prediction of concrete structure with grey system models[J]. Journal of Hunan University:Natural Sciences, 1994，21(4):96-102. (In Chinese)

[8]金成棣．混凝土徐變對超靜定結構變形及內力的影響——考慮分段加載齡期差異及延遲彈性影響[J].土木工程學報，1981，14(3):19-32.

JIN Cheng-di. Analysis of deformation and stress redistribution of continuous structure due to creep of concrete—taking account of loading at different age and delayed elasticity of concrete[J].China Civil Engineering Journal, 1981，14(3):19-32. (In Chinese)

[9]占玉林,向天宇,趙人達. 幾何非線性結構的徐變效應分析[J]. 工程力學,2006,23(7):45-48.

ZHAN Yu-lin,XIANG Tian-yu,ZHAO Ren-da. Creep effect analysis of geometric nonlinear structures[J]. Engineering Mechanics, 2006, 23(7): 45-48. (In Chinese)

[10]陳常松,顏東煌,李學文. 混凝土收縮徐變分析的虛功增量方程及應用[J]. 工程力學,2010,27(10):139-144.

CHEN Chang-song,YAN Dong-huang，LI Xue-wen.The incremental virtual work equation for concrete shrinkage and creep analysis and its applications[J]. Engineering Mechanics,2010,27(10): 139-144. (In Chinese)

[11]CRISFIELD M A. Nonlinear finite element analysis of solids and structures[M].Chichester: John Wiley & Sons Inc, 1997：26-45.

[12]李學文，姚康寧，顏東煌. 利用最小二乘法實現2004規范徐變系數的指數函數擬合[J]. 長沙交通學院學報，2006，22(3)：21-24.

LI Xue-wen,YAO Kang-ning,YAN Dong-huang. Using least square method fitting the creep coefficient functions of concrete listed in the bridge criterion (JTG D62—2004)with exponential function model[J]. Journal of Changsha Communications University, 2006，22(3)：21-24. (In Chinese)

[13]姚康寧.大跨度混凝土斜拉橋運營階段混凝土收縮徐變影響研究[D]. 長沙: 長沙理工大學土木與建筑學院, 2006:37-49.

YAO Kang-ning. The shrinkage and creep analysis of long-span concrete cable-stayed bridges during operations[D]. Changsha: School of Civil Engineering and Architecture,Changsha University of Science and Technology, 2006: 37-49. (In Chinese)

[14]顏東煌,田仲初,李學文,等. 混凝土橋梁收縮徐變計算的有限元方法與應用[J].中國公路學報,2004,17(2):55-58.

YAN Dong-huang,TIAN Zhong-chu,LI Xue-wen,etal. Finite element method and application for the shrinkage and creep of concrete bridge[J].China Journal of Highway and Transport, 2004,17(2):55-58. (In Chinese)

[15]蔡松柏，沈蒲生，胡柏學，等. 基于場一致性的2D四邊形單元的共旋坐標法[J].工程力學，2009，26(12)：31-34.

CAI Song-bai, SHEN Pu-sheng, HU Bai-xue,etal.A field consistency based co-rotational finite element procedure for 2D quadrilateral element[J]. Engineering Mechanics，2009，26(12)：31-34. (In Chinese)

[16]ESPION B Benchmark. Examples for creep and shrinkage analysis computer program[C]//Proceeding of the 5th International RILEM Symposium．1993:901-911.