于 翔1，曹明輝2，譚建平1

（1、長沙理工大學土木工程學院 長沙410114；2、中國建筑第五工程局有限公司 長沙410004）

0 概述

結(jié)構(gòu)的變形是一個動態(tài)的平衡過程，結(jié)構(gòu)實際發(fā)生的大轉(zhuǎn)動、大變形效應，使得結(jié)構(gòu)變形不能單純依靠線性方程進行求解，需在幾何方程中添加位移二次項式。一部分學者通過齡期調(diào)整有效模量法將混凝土的徐變效應與混凝土彈性模量相結(jié)合，以此考慮結(jié)構(gòu)幾何非線性與徐變效應［1］。該法在短期徐變分析中精度較高，考慮結(jié)構(gòu)荷載長期作用下徐變效應時就存在較為明顯的誤差，而且?guī)缀畏蔷€性的變形與徐變變形應當是互相影響，不斷迭代，最后達到動態(tài)平衡的過程，如果將其作為一個整體來分析計算，則會產(chǎn)生相應偏差［2］。另一部分學者通過推導考慮幾何非線性的結(jié)構(gòu)變形平衡方程，自編程序?qū)ζ溥M行分析，在計算精度與效率方面較為符合工程需要［3］。

1 平面梁單元幾何非線性理論

桿單元在不同時刻的變形可以用圖1 表示，其中XY為全局結(jié)構(gòu)坐標系，指向始終一致。引入共旋坐標系作為局部的單元坐標系，該坐標系隨時刻發(fā)生的變形而改變，并且始終以單元i節(jié)點指向j節(jié)點為x軸，y軸為x軸逆時針旋轉(zhuǎn)90°而成［4-7］。

初始時刻單元在全局坐標系中坐標分別為（xi，yi）和（xj，yj），計算時刻t梁單元總位移在全局坐標系中表示為d=［ui vi θ i uj vj θ j］T，共旋坐標系下的單元坐標系中表示為dl=［ui'vi'θi'uj'vj'θj'］T。其中u、v、θ分別表示節(jié)點水平位移、豎向位移、轉(zhuǎn)角位移。單元變形前后情況如圖1所示。

圖1 不同時刻梁單元Fig.1 Beam Element at Different Time

共旋法的局部坐標下平面梁單元各節(jié)點位移可簡化為：

式中：0l為初始時刻梁單元長度；tl為計算時刻t梁單元長度。其計算公式為：

式中：EA為抗壓剛度；EI為抗彎剛度。

平面梁單元在大位移、小應變的幾何非線性下的切線剛度矩陣Kt為：

式中具體表達式由文獻［7］推導給出。

2 初應變法計算徐變效應

共旋坐標系下單元的應變位移關系無需考慮非線性項，即：

式中：εe為節(jié)點位移產(chǎn)生的彈性應變矩陣；B0為線性應變矩陣。

考慮初應變的共旋單元坐標系下應力可表示為：

式中：ε0為混凝土徐變引起的初應變；D為應力矩陣。混凝土線性徐變應變εc表達式為：

幾何非線性下徐變效應的分析是在疊加原理條件下在每一個時刻中進行迭代，直到其迭代精度達到要求后才會執(zhí)行下一時刻的計算。鄧繼華等人［8］推導得出tn?tn+1時段內(nèi)，應變增量△εn+1c與應力增量△σ i表達式為：

由計算原理可明顯得出：

經(jīng)推導后的遞推公式可得出：

當n=0時：

當n=1時：

常規(guī)方法和遞推公式都可得出如下公式：

3 算例

3.1 有限元模型對比

本文采用工程實例為湖南省南縣的某大橋，以此建模并展開相應計算分析。該大橋為主跨368 m的大跨度自錨中承式鋼管混凝土系桿拱橋，主橋跨徑組合為（80+368+80）m。主跨線形為懸鏈線拱軸線，m=1.543，計算跨徑L=356 m，計算矢高71.2 m，矢跨比為1∕5。利用ANSYS 中BEAM3 單元將三維立體的模型壓縮成平面梁單元的茅草街模型，且數(shù)據(jù)與BEAM188 基本無異。然后通過格式化輸出模型數(shù)據(jù)，導入以Fortran匯編了可分析幾何非線性與徐變共同作用的程序。

整個模型共分為94 個節(jié)點，307 個單元。將空間模型壓縮成平面模型，BEAM3 單元不支持蠕變功能，因此以BEAM3 單元自重作用下是否考慮幾何非線性的跨中撓度為基準進行核對，以保證整個計算的說服性。具體如圖2、表1所示。

圖2 不同單元計算的撓度Fig.2 The Maximum Deflection of Midspan Calculated by Different Software，Element Type and Calculation Method （cm）

表1 不同軟件、單元類型、計算方式所計算的跨中最大撓度Tab.1 The Maximum Deflection of Midspan Calculated by Different Software，Element Type and Calculation Method （cm）

由圖2 和表1 可知，ANSYS 與Fortran 平面梁單元計算結(jié)果差值1%以內(nèi)，說明本文采用共旋坐標法的計算精度足以支撐整個計算。以Fortran 平面梁與ANSYS 空間梁BEAM188 的計算結(jié)果計算誤差5%左右，后續(xù)計算徐變效應分析時將以BEAM188 計算數(shù)據(jù)作為參為依據(jù)，以加強數(shù)據(jù)的說服性。

3.2 誤差對比

采用文獻［10］中時步劃分方法，以對數(shù)形式將1000d計算齡期劃分為20個時步。初始加載齡期為10 d。

圖3 和表2 顯示2 種計算不同軟件、不同單元類型、不同徐變分析模式計算下的差值從自重靜力分析下的4.45%最終趨于0.61%的差值，此間差異是越來越小的。整個誤差是正常的迭代誤差，從數(shù)據(jù)層面上是可以說明Fortran 自編程徐變效應計算程序具有說服性。

圖3 ANSYS與Fortran不同齡期跨中拱頂豎向撓度徐變比對Fig.3 Comparison of Vertical Deflection and Creep of Mid Span Vault in Different Ages between ANSYS and Fortran （cm）

表2 ANSYS與Fortran不同齡期跨中拱頂豎向撓度徐變比對Tab.2 Comparison of Vertical Deflection and Creep of Mid Span Vault in Different Ages between ANSYS and Fortran （cm）

3.3 差異分析

在原有是否考慮幾何非線性的自重作用下計算分析基礎上，增加了是否考慮幾何非線性與徐變效應的雙重非線性數(shù)值計算。并比較四者的豎向撓度，尋找并發(fā)現(xiàn)差異所在（見圖4、表3）。

圖4 各齡期下跨中拱頂豎向撓度Fig.4 Vertical Deflection of Middle Span Vault at Different Ages

表3 4種計算方式下跨中拱頂豎向撓度比Tab.3 Vertical Deflection Ratio of Mid Span Vault under Four Calculation Methods （cm）

不考慮徐變的撓度與軸力不隨齡期的發(fā)生而變化。隨著齡期的增長，徐變速率由最初的急劇增長之后趨于平緩，同時鋼管應力隨著齡期的增長而增長，同時混凝土應力在徐變效應的影響下開始下降。

4 結(jié)論

通過基于共旋法推導了可考慮大變形、轉(zhuǎn)動的平面梁單元的切線剛度矩陣，并在此基礎上將徐變效應以初應變法作為等效節(jié)點荷載施加于梁上。以某大橋為實例建模，在ANSYA 建?；A上壓縮為平面模型，并進行格式化輸出，并將Fortran自編程序與ANSYS對比，以保證數(shù)據(jù)準確性與說服性。結(jié)論如下：

⑴將ANSYS 建成的BEAM188 空間模型轉(zhuǎn)換成BEAM3 平面模型，并將關鍵數(shù)據(jù)格式化輸出，導入Fortran自編程序。并與Fortran計算結(jié)果相比較，保證了計算精確性。

⑵以對數(shù)形式將1 000 d 計算齡期劃分為20 個時步進行考慮，對自編程序進行了是否考慮幾何非線性的徐變效應分析。結(jié)果發(fā)現(xiàn)，考慮幾何非線性較之線性下的徐變分析而言，其數(shù)值相對于自重作用下的數(shù)值差異為5%左右。當計算荷載更為龐大、跨徑更大等情況下，結(jié)構(gòu)非線性更為明顯，其雙重非線性考慮則更有必要。