謝雪峰，羅喜恒

(1.同濟(jì)大學(xué)橋梁工程系，上海 200092;2.同濟(jì)大學(xué)建筑設(shè)計(jì)研究院(集團(tuán))有限公司，上海 200092)

1 前言

懸索橋分析程序可以分為專用程序和通用程序。專用程序如平面纜索懸索橋數(shù)值分析程序SBP、CAP(construction analysis program for suspension bridge)［1］和懸索橋分析系列程序［2］以及懸索橋主纜設(shè)計(jì)與施工計(jì)算專用軟件SGK2000、橋梁結(jié)構(gòu)空間幾何非線性靜動(dòng)力分析軟件BCAS2000［3］等。通用程序主要有橋梁博士、Midas/Civil、ANSYS、SAP2000等。專用程序和通用程序在懸索橋分析方面各有優(yōu)缺點(diǎn)，如專用程序具有能考慮懸索橋的特殊構(gòu)造、加載方便等優(yōu)點(diǎn)，計(jì)算速度較快，計(jì)算精度較高，但前后處理一般配合不好，有些程序還有些缺陷;而通用程序前后處理方便，但一般沒有專門的鞍座單元、頂推單元，并需要人為給出主纜和吊索的無應(yīng)力長度。

筆者在ANSYS平臺基礎(chǔ)上進(jìn)行二次開發(fā)，結(jié)合懸索橋的特點(diǎn)構(gòu)造了鞍座單元和模擬鞍座頂推的算法，編寫了主纜和吊索無應(yīng)力長度的迭代求解APDL程序，使得ANSYS在懸索橋分析方面具有了專用程序的結(jié)構(gòu)模擬精度和計(jì)算精度，并且具有良好的前后處理功能。

2 分析方法

2.1 鞍座的模擬

鞍座是使主纜平順轉(zhuǎn)向的構(gòu)件，計(jì)算時(shí)一般認(rèn)為主纜應(yīng)總是與鞍座相切的，然而在施工階段和活載作用過程中，塔頂會(huì)有比較大的位移，主纜與鞍座的切點(diǎn)位置也有比較大的變動(dòng)，因此，正確模擬鞍座及切點(diǎn)的變化，是很有必要的。

以塔頂鞍座為例，如圖1所示，鞍座單元的構(gòu)造如下:A為鞍座標(biāo)識線(成橋時(shí)與主塔中心線重合)上的主纜點(diǎn)，為主纜與鞍座的固定點(diǎn);B為鞍座標(biāo)識點(diǎn)，成橋時(shí)與塔頂中心點(diǎn)重合;C、D為主纜與鞍座的切點(diǎn);E、F為主纜上的節(jié)點(diǎn)。單元AB、AC、AD為剛臂單元，采用mpc184單元。CE、DF為索單元，采用link10單元。BHG為頂推剛臂單元。

圖1 鞍座單元構(gòu)造形式Fig.1 Saddle element form

在鞍座頂推或活載作用下，鞍座兩側(cè)圓弧段的圓心相對于AB是不變的，也就是說，可以根據(jù)AB的轉(zhuǎn)角和索單元CE、DF在C、D點(diǎn)的切線角更新切點(diǎn)C、D的位置，直至前后兩次索單元CE、DF在C、D點(diǎn)的切線角相差小于給定的限值，則認(rèn)為相切。一般只需更新兩到三次就能達(dá)到較高精度。

以下以D點(diǎn)為例給出其原理。假設(shè)圓弧AD的圓心為O(Ox，Oy)，半徑為R，設(shè)AB長為l，且逆時(shí)針擺動(dòng)α角，主纜在D點(diǎn)的切線與水平向夾角此時(shí)為β，則更新D(xd，yd)點(diǎn)位置如下:

散索鞍的構(gòu)造同上，只是擺動(dòng)式散索鞍沒有頂推剛臂BHG。

2.2 頂推的實(shí)現(xiàn)

在施工過程分析時(shí)，需要模擬頂推鞍座這一懸索橋特有的工況。如圖2和圖3所示，頂推過程的模擬采用以下單元構(gòu)造形式:BGI為mpc184滑塊單元，B為滑塊;HGI為mpc184滑塊單元，H為滑塊;BH、GH、B’I、B’G 為 mpc184 剛臂單元，B’J為塔頂梁單元。賦予單元GH一個(gè)合適的線膨脹系數(shù)，通過改變GH的溫度即可實(shí)現(xiàn)鞍座頂推的模擬。如單元GH的長為2 m，線膨脹系數(shù)設(shè)為1，以當(dāng)前工況為基準(zhǔn)，降溫該單元0.5℃則頂推1 m。在施工階段分析時(shí)可以通過殺死GH單元來實(shí)現(xiàn)鞍座能隨施工過程自由滑動(dòng)的情況，從而獲得滑動(dòng)量和塔頂水平力隨施工過程的變化曲線，進(jìn)而合理設(shè)計(jì)頂推次數(shù)和頂推量。

圖2 頂推單元構(gòu)造形式Fig.2 Incremental launching element form schematic diagram

2.3 主纜和吊索無應(yīng)力長度的確定

主纜和吊索無應(yīng)力長度求解，一般采用文獻(xiàn)［1］中公式(2.1)和(2.2)進(jìn)行迭代求解。然而，在利用ANSYS進(jìn)行懸索橋分析時(shí)候，如果可以用AP-DL(ANSYS parametric design language)語言進(jìn)行迭代求解出無應(yīng)力長度，則模型一氣呵成，修改起來也很方便，提高了工作效率。以下是迭代過程的基本原理和流程圖。

圖3 頂推單元示意圖Fig.3 Incremental launching element

2.3.1 基本原理［4］

1)原理1:非線性方程求解的Newton-Raphson方法。采用Newton-Raphson方法，每步迭代都要計(jì)算導(dǎo)數(shù)值。在該問題中，是無法計(jì)算出導(dǎo)數(shù)值的，因此采用兩點(diǎn)割線法，以差商代替導(dǎo)數(shù)，見圖4。迭代公式如下:

圖4 牛頓－辛普森法Fig.4 Newton-Raphson method

2)原理2:非線性方程組求解的迭代法。F:D?Rn→ Rn是 n元實(shí)向量函數(shù)，把非線性方程組F(x)=0改寫成不動(dòng)點(diǎn)迭代x=G(x)，則有如下定理:

設(shè)G:D?Rn→Rn在區(qū)域D0?D上滿足G把D0映入它自身，即G(D0)?D0;G在D0上是壓縮映射，壓縮因子為L＜1，即對任意x，y∈D0有

則下列結(jié)論成立:a.G在D0上存在唯一的不動(dòng)點(diǎn)x*;b.對于任意的x(0)∈D0，不動(dòng)點(diǎn)迭代x(k+1)=G(x(k))產(chǎn)生的迭代序列 { x(k)}?D0且收斂于唯一的不動(dòng)點(diǎn)x*。

2.3.2 基本原理的應(yīng)用

2.3.2.1 原理1 的應(yīng)用

以中間跨為例，定義函數(shù)

即跨中節(jié)點(diǎn)在恒載作用下豎向位置與設(shè)計(jì)豎向位置的差值為跨中節(jié)點(diǎn)豎向坐標(biāo)的函數(shù)。式中，nym為跨中節(jié)點(diǎn)的豎向坐標(biāo);uym為跨中節(jié)點(diǎn)在恒載作用下的豎向位移;f0為跨中節(jié)點(diǎn)的設(shè)計(jì)豎向位置。

以邊跨為例，定義函數(shù)

即邊跨主纜在恒載作用下水平力與中間跨的成橋狀態(tài)水平力的差值為邊跨主纜單元初應(yīng)變的函數(shù)。錨跨和邊跨相同。式中，H為邊跨主纜在自重和索夾重等作用下的水平力;H0為中間跨的成橋狀態(tài)水平力;istn為邊跨主纜單元的初應(yīng)變。

2.3.2.2 原理2 的應(yīng)用

一般情況下，按照設(shè)計(jì)的成橋位置建立纜索系統(tǒng)的初始模型。在恒載的作用下，上下吊點(diǎn)會(huì)發(fā)生水平位移和豎向位移。為了使成橋時(shí)各吊點(diǎn)位于設(shè)計(jì)位置，在迭代過程當(dāng)中，需要修改上吊點(diǎn)的水平坐標(biāo)和下吊點(diǎn)的豎向坐標(biāo)。

以上吊點(diǎn)水平坐標(biāo)為例，設(shè)其個(gè)數(shù)為n組成向量 x(x1，x2，…，xm，…，xn) 。向量 nx0表示吊桿上節(jié)點(diǎn)的成橋狀態(tài)在恒載力作用下的水平位置，即目標(biāo)位置。構(gòu)造向量nx=nx0－ux，ux表示上吊點(diǎn)的水平位移。nx滿足原理2的條件，在迭代過程當(dāng)中，逐次更新吊桿上節(jié)點(diǎn)的水平坐標(biāo)為nx，直至收斂。

2.3.3 迭代流程

以控制點(diǎn)為依據(jù)建立成橋狀態(tài)初始模型，先假定各個(gè)切點(diǎn)的位置，中間跨以兩切點(diǎn)和跨中點(diǎn)為依據(jù)，按照拋物線來確定其他節(jié)點(diǎn)坐標(biāo);錨跨和邊跨以兩切點(diǎn)間連直線段為依據(jù)，確定其他節(jié)點(diǎn)坐標(biāo)。初始的纜索單元初應(yīng)變設(shè)為1.0E－8。以中間跨迭代求無應(yīng)力長度為例，迭代流程如圖5所示，其中Y0、Y1為跨中節(jié)點(diǎn)豎向坐標(biāo)，函數(shù) f即2.3.2節(jié)中的函數(shù)f1。根據(jù)2.3.2可知，纜索系統(tǒng)的未知量:纜索單元的無應(yīng)力長度，其個(gè)數(shù)與迭代方程個(gè)數(shù)相等［1，3］，可迭代求解出纜索系統(tǒng)的線形和無應(yīng)力長度。

圖5 中間跨迭代流程Fig.5 Main span iterative process

2.4 無應(yīng)力長度保持不變的方法

按照上述方法迭代出纜索系統(tǒng)的無應(yīng)力長度之后，在鞍座頂推或活載作用下，更新切點(diǎn)位置時(shí)，只要保證鞍座頂與緊鄰的吊桿之間主纜無應(yīng)力長度不變，則整個(gè)主纜的無應(yīng)力長度不變。以圖1所示鞍座為例，圓弧AC、AD上主纜的無應(yīng)力長度按照文獻(xiàn)［1］中鞍座內(nèi)主纜內(nèi)力分布的模式二進(jìn)行計(jì)算，即文獻(xiàn)［1］中式(2.14)、式(2.15)。施工計(jì)算或活載加載過程當(dāng)中，進(jìn)行迭代求解，逐次更新切點(diǎn)坐標(biāo)，并計(jì)算新的圓弧段上的主纜無應(yīng)力長度，然后修改主纜單元CE、DF的初應(yīng)變，以保證ACE、ADF這兩段的無應(yīng)力長度不變，直至鞍座與兩邊的主纜重新相切。在初應(yīng)變較小的情況下，這種做法精度是很高的。

圖6 索塔結(jié)構(gòu)的計(jì)算圖示Fig.6 Cable－tower structure schematic diagram

3 驗(yàn)證性算例［1］

以潤揚(yáng)大橋?yàn)槔?，按照上述方法建立ANSYS計(jì)算模型，原始數(shù)據(jù)文獻(xiàn)［1］中算例Z5一致。不考慮主塔自重的影響，計(jì)算時(shí)未考慮樁的影響?？紤]以下3個(gè)工況，比較CAP與ANSYS的計(jì)算結(jié)果。

1)工況1:塔頂鞍座的相對偏移量取1.787 5 m;

2)工況2:在工況1的基礎(chǔ)上頂推鞍座，使塔頂鞍座的相對偏移量為1.0 m;

3)工況3:相當(dāng)于成橋狀態(tài)，在工況2基礎(chǔ)上頂推鞍座，使塔頂鞍座的相對偏移量為0，并將全部荷載作用在結(jié)構(gòu)上。

表1、表2、表3給出了二者控制點(diǎn)位置、主纜水平力和塔頂內(nèi)力及無應(yīng)力長度的比較。從結(jié)果來看，兩者非常吻合，印證了上述方法的正確性。從塔頂內(nèi)力的情況來看，兩者剪力隨塔頂鞍座的相對偏移量的加大而差異變大，是因?yàn)锳NSYS中梁單元考慮了剪切剛度的影響，所以比 CAP中梁單元柔一些。

表1 控制點(diǎn)位置對比Table 1 Comparison of the location of control points

表2 主塔和主纜內(nèi)力對比Table 2 Comparison of the tower and main cable’s internal force

表3 無應(yīng)力長度對比Table 3 Comparison of unstressed length

4 結(jié)語

1)該方法精度較高，能充分滿足工程及科研精度要求。

2)該方法不需要借助于其他程序迭代出纜索無應(yīng)力長或者切線角，整體性強(qiáng)，能夠提高工作效率。

3)該方法可復(fù)制性強(qiáng)，參數(shù)修改也很方便，對于確定懸索橋合適的結(jié)構(gòu)參數(shù)十分有利。

［1］羅喜恒.復(fù)雜懸索橋施工過程精細(xì)化分析研究［D］.上海:同濟(jì)大學(xué)博士論文，2004.

［2］潘永仁.懸索橋的幾何非線性靜力分析及工程控制［D］.上海:同濟(jì)大學(xué)博士論文，1996.

［3］唐茂林.大跨度懸索橋空間幾何非線性分析與軟件開發(fā)［D］.成都:西南交通大學(xué)博士論文，2003.

［4］同濟(jì)大學(xué)計(jì)算數(shù)學(xué)教研室.現(xiàn)代數(shù)值數(shù)學(xué)和計(jì)算［M］.上海:同濟(jì)大學(xué)出版社，2004.