基于有限元列式的懸索橋主纜找形算法

孫 遠(yuǎn)， 羅文孝， 劉 梅

(1. 華中科技大學(xué) 土木工程與力學(xué)學(xué)院， 湖北 武漢 430074；2. 貴州高速公路集團(tuán)有限公司， 貴州 貴陽(yáng) 550009)

懸索橋主纜屬于柔性結(jié)構(gòu)，在無(wú)應(yīng)力狀態(tài)下沒有確定的初始位形，其幾何剛度完全來(lái)源于初始內(nèi)力，因此，確定主纜平衡狀態(tài)的方法不同于框架結(jié)構(gòu)：一般會(huì)在設(shè)計(jì)初始單獨(dú)對(duì)主纜進(jìn)行分析以確定目標(biāo)線形，即主纜找形。盡管過去的研究成果表明,一般懸索橋的主纜找形并不是一個(gè)非常困難的問題[1,2],但隨著近些年復(fù)雜纜索體系市政橋梁的快速發(fā)展,傳統(tǒng)方法已體現(xiàn)出若干劣勢(shì)，追求整體性更好、效率更高的找形方法仍然是研究者們較為感興趣的課題[3]。

主纜找形方法主要分為解析法、非線性有限元法和數(shù)值迭代法，其中后兩者較為常用。這是因?yàn)槠胀ń馕鲱惙椒ㄍǔ＜俣ǖ跛髁鶆蚍植己椭骼|為拋物線，與實(shí)際情況不能完全符合；日本的 Ohtsuki 博士針對(duì)三維彎曲的主纜體系提出了節(jié)線法[4]，通過建立節(jié)點(diǎn)力的平衡方程組求解，但在無(wú)應(yīng)力長(zhǎng)度未知時(shí)只能作為近似分析，僅適用于初步設(shè)計(jì)。非線性有限元法的代表是借鑒于各類索網(wǎng)結(jié)構(gòu)找形[5～7]的“預(yù)加應(yīng)力法”，能夠整體地分析主纜，但有兩個(gè)問題：一是多迭代子步帶來(lái)的誤差積累，這是由非線性計(jì)算本身的特點(diǎn)所決定的；二是收斂的不穩(wěn)定性。該方法基本能夠處理加載路徑單一，成橋纜形易被估計(jì)的平面纜形，但有時(shí)仍會(huì)遇到迭代更新位移過大導(dǎo)致不收斂的情況，于是有學(xué)者通過分類討論和引入優(yōu)化算法等形式提高收斂效率[8]，取得了一定的成效。對(duì)此，Kim等[4]在永宗大橋的主纜體系中利用懸鏈線柔度矩陣對(duì)節(jié)線法的單元無(wú)應(yīng)力長(zhǎng)度進(jìn)行修正，并根據(jù)非線性位移迭代計(jì)算取得平衡。該方法一度被編入軟件，但被用戶發(fā)現(xiàn)在處理各種不同類型的纜索時(shí)極易發(fā)散，這是因?yàn)榈兄骼|水平分力保持不變卻并不準(zhǔn)確時(shí)，滿足設(shè)計(jì)條件的平衡狀態(tài)往往無(wú)法到達(dá)?；谙嗤脑恚琄im和Lee[9]，以及Myung等[10]以懸鏈線柔度矩陣中的單元無(wú)應(yīng)力長(zhǎng)度為未知數(shù)提出了TCUD (Target Configuration under Dead load)及其改進(jìn)算法，該方法計(jì)算結(jié)果準(zhǔn)確但對(duì)初值的要求較高。

當(dāng)前數(shù)值迭代法的優(yōu)勢(shì)是計(jì)算精確，其代表是分段懸鏈線法。唐茂林[11]提出了二維的分段懸鏈線法，羅喜恒等[12，13]則將該方法從二維擴(kuò)展至三維。該方法也存在兩個(gè)問題：一是限定較多，假定吊索必須與順橋向垂直，并且各跨主纜需要單獨(dú)計(jì)算，對(duì)于多向傾斜吊桿以及多跨懸索橋有所不便；二是迭代計(jì)算中需要不斷求解懸索的非線性狀態(tài)方程，容易遇到方程無(wú)解而導(dǎo)致不能收斂的情形。對(duì)此，一些學(xué)者采用了不同的數(shù)學(xué)修正模型以改善特定狀態(tài)下的收斂。李傳習(xí)[14]則直接指出了索段狀態(tài)方程無(wú)解時(shí)的若干情形，繼而采用成橋時(shí)具有近似線性剛度的假設(shè)以及懲罰因子來(lái)建立修正規(guī)則，提出了改進(jìn)的分段懸鏈線法，該方法收斂穩(wěn)定，但步驟略為復(fù)雜。

針對(duì)當(dāng)前主纜找形方法存在的問題，本文提出一種新的數(shù)值迭代方法即坐標(biāo)迭代算法：將主、邊跨主纜看作一個(gè)整體，建立直接求解節(jié)點(diǎn)平衡坐標(biāo)的迭代表達(dá)式并統(tǒng)一處理兩類主纜找形問題。該方法完全基于結(jié)構(gòu)力學(xué)基本原理，將線性有限元求解思路巧妙的引入到對(duì)柔性結(jié)構(gòu)的計(jì)算中，迭代過程不需要借助于特別的數(shù)學(xué)方法，具有良好的整體性和穩(wěn)定的收斂特性，能夠部分地克服當(dāng)前計(jì)算方法中所存在的問題。

1 基本原理

新數(shù)值迭代算法的核心思想來(lái)源于合理拱軸線的力學(xué)特征。實(shí)腹式拱橋的合理拱軸線與單段懸索在自重下具有相同的懸鏈線微分方程[15]，區(qū)別僅在于恒載壓力線和拉力線，這表明合理拱軸線可以用來(lái)解釋懸索橋主纜的力學(xué)平衡。同時(shí)，主纜無(wú)論是否作用有吊索力，其平衡線形與恒載拉力線都是自然重合的；一方面說明各種平衡狀態(tài)下的主纜線形都與合理拱軸線有相似之處；另一方面也說明求解恒載拉力線的過程就是求解主纜線形的過程。計(jì)算中，類似于在拱上作用有若干集中力的空腹式拱橋只能考慮在部分點(diǎn)達(dá)到兩線重合[15]，吊索力作用下的主纜也只需以有限的單元和節(jié)點(diǎn)為對(duì)象進(jìn)行考慮，如圖1所示。因此，可基于合理拱軸線的彎矩平衡特征構(gòu)建直接求解主纜坐標(biāo)的平衡迭代方程。本節(jié)內(nèi)容包含坐標(biāo)迭代算法的計(jì)算假定以及平衡迭代方程的建立。

圖1 拱橋合理拱軸與懸索橋主纜的力學(xué)相似性

1.1 計(jì)算假定

實(shí)際懸索主纜既具有索的柔度，又具有梁的抗彎剛度，無(wú)論是單純的懸鏈線單元，修改了剛度的梁?jiǎn)卧?，抑或是其他特殊索單元都不能完全精確地模擬懸索的實(shí)際特性，但在設(shè)計(jì)上均能滿足要求。在單元長(zhǎng)度合理時(shí)，桁架單元在處理索結(jié)構(gòu)問題的工程實(shí)踐中已被廣泛證明具有較高的精度[6]。當(dāng)線性有限元被引入到基于桁架索單元的坐標(biāo)迭代算法時(shí)，單元?jiǎng)澐指屿`活，使得單個(gè)桁架索單元的垂度效應(yīng)完全可以被忽略。本文坐標(biāo)迭代算法的假定主要包含以下幾點(diǎn)：

(1)不考慮垂度效應(yīng)對(duì)主纜剛度的影響。相鄰吊索之間的主纜可以用一個(gè)或多個(gè)桁架單元模擬。

(2)單個(gè)桁架索單元的自重均布荷載可按照有限元的基本理論[16]轉(zhuǎn)化為等效節(jié)點(diǎn)力。

(3)不考慮張拉效應(yīng)對(duì)于主纜剛度的影響。

(4)不考慮鞍座形狀對(duì)于主纜目標(biāo)線形的影響，以IP點(diǎn)(鞍座兩側(cè)主纜切線交點(diǎn))計(jì)入體系。

1.2 平衡迭代方程

當(dāng)懸索橋主纜滿足力學(xué)平衡時(shí)，應(yīng)具有類似合理拱軸線的彎矩平衡特征，即外荷載對(duì)軸線上每點(diǎn)的彎矩為0(計(jì)算時(shí)只需考慮對(duì)每個(gè)節(jié)點(diǎn))，單個(gè)單元不能發(fā)生轉(zhuǎn)動(dòng)，可看作單元平衡；此外，所有單元在節(jié)點(diǎn)處還要滿足節(jié)點(diǎn)力平衡。因此，主纜平衡狀態(tài)應(yīng)包含相當(dāng)于充要條件的兩點(diǎn)：

(1)每個(gè)桁架索單元的軸線水平和豎向長(zhǎng)度與相應(yīng)桿端力成比例，滿足單元平衡。

(2)每個(gè)節(jié)點(diǎn)各坐標(biāo)方向上集中力的代數(shù)和為0，滿足節(jié)點(diǎn)力平衡。

這兩點(diǎn)是建立坐標(biāo)迭代格式的基礎(chǔ)。對(duì)于單個(gè)桁架索單元而言，是基于第n-1迭代階段主纜線形近似解求出的單元桿長(zhǎng)和軸力建立第n階段的單元關(guān)系矩陣，從而形成第n迭代階段的單元平衡迭代方程。由于該方程在形式上類似于有限元的單元?jiǎng)偠确匠?，可借用有限元法的總綱矩陣對(duì)號(hào)入座規(guī)則[15]建立總體平衡迭代方程進(jìn)行求解，獲得第n階段主纜線形的近似解。重復(fù)上述過程，當(dāng)主纜節(jié)點(diǎn)坐標(biāo)經(jīng)過若干次迭代達(dá)到收斂時(shí)，即為精確解。

圖2為包含迭代信息的單個(gè)桁架索單元AB的受力圖，上標(biāo)括號(hào)內(nèi)表示迭代的階段。

圖2 包含迭代信息的桁架索單元

根據(jù)上述基本原理，已知桿長(zhǎng)lAB(n-1)和軸力TAB(n-1)，若要保持單元平衡，單元各力素應(yīng)符合一定的三角函數(shù)關(guān)系；此外，在正確的收斂趨勢(shì)下，各力素符號(hào)的相對(duì)關(guān)系也是一致的。不妨令軸力與單元A端桿端力的符號(hào)一致，在實(shí)際中即受拉為負(fù)，可獲得桿端力F與坐標(biāo)的關(guān)系式

(1)

(2)

式中：θ為單元與x軸的夾角。式(1)，(2)經(jīng)過整理后可得到

(3)

即

(4)

其中，

(5)

方程(4)即為以坐標(biāo)為未知數(shù)的單元平衡迭代方程。式(5)考慮了符號(hào)一致性，在一般情況下成立，但為了使初始線形在迭代中收斂至目標(biāo)線形，需要一些根據(jù)設(shè)計(jì)要求和階段來(lái)確定的附加控制條件；基于不同的控制條件可以建立總體平衡迭代方程不同的計(jì)算模式。

為了使得整個(gè)結(jié)構(gòu)達(dá)到平衡，還需使各單元在節(jié)點(diǎn)處滿足節(jié)點(diǎn)力平衡。因?yàn)閱卧胶獾匠?式(4))從形式上類似于有限元的單元?jiǎng)偠确匠?，只是未知?shù)從位移變成了坐標(biāo)，所以滿足節(jié)點(diǎn)力平衡的過程實(shí)際上就是借用有限元理論形成總體平衡迭代方程的過程。在不同設(shè)計(jì)階段組集單元平衡迭代方程可獲得不同未知數(shù)和控制條件下的總體平衡迭代方程,其類型主要分為設(shè)計(jì)模式和計(jì)算模式。

(1)設(shè)計(jì)模式

該模式下單元無(wú)應(yīng)力長(zhǎng)度未知，控制條件為已知設(shè)計(jì)參數(shù)，涵蓋了成橋主纜的兩類找形問題。處理方法：除固定邊界外，節(jié)點(diǎn)坐標(biāo)和單元無(wú)應(yīng)力長(zhǎng)度均為未知數(shù)，在迭代過程中保持更新；已知設(shè)計(jì)參數(shù)(主控制點(diǎn)設(shè)計(jì)坐標(biāo)、跨中垂點(diǎn)坐標(biāo)或主纜水平分力等)在迭代中保持不變，總體平衡迭代方程如式(6)所示。

(6)

式中：X為總體節(jié)點(diǎn)坐標(biāo)向量；S0為總體無(wú)應(yīng)力長(zhǎng)度向量；K為總體關(guān)系矩陣，由上一迭代階段求得的坐標(biāo)與桿端力確定；F為總體桿端力向量，由節(jié)點(diǎn)坐標(biāo)和式(4)反算確定，可視為X的函數(shù)；P為總體外荷載向量，由單元自重和吊索力確定，可視為單元無(wú)應(yīng)力長(zhǎng)度的函數(shù)；Hα為已知參數(shù)向量(元素為力或坐標(biāo))。

(2)計(jì)算模式

該模式下單元無(wú)應(yīng)力長(zhǎng)度已知，控制條件為“單元無(wú)應(yīng)力長(zhǎng)度不變”，可用于求解空纜線形和鞍座預(yù)偏。處理方法：除固定邊界外，節(jié)點(diǎn)坐標(biāo)為未知數(shù)，在迭代過程中保持更新；基于已知單元無(wú)應(yīng)力長(zhǎng)度來(lái)控制迭代中的桿長(zhǎng)與軸力的關(guān)系保持不變，總體平衡迭代方程如式(7)所示。

(7)

式中：l為單元長(zhǎng)度向量，由總體桿端力向量F和無(wú)應(yīng)力長(zhǎng)度控制條件確定，可看作F的函數(shù)；li,s0i分別表示第i號(hào)單元的桿長(zhǎng)與無(wú)應(yīng)力長(zhǎng)度；FiAx,FiAy分別表示第i號(hào)單元x，y方向的桿端力；E,A為單元的彈性模量與截面積。

2 平衡方程求解

2.1 求解條件

總體平衡迭代方程在求解時(shí)需要設(shè)置迭代初始線形、初始軸力和收斂容差，在設(shè)計(jì)模式時(shí)還需要設(shè)置迭代初始無(wú)應(yīng)力長(zhǎng)度；作為以坐標(biāo)為未知數(shù)的迭代方程，在處理邊界條件時(shí)也與常規(guī)有限元方法有所區(qū)別。在單元節(jié)點(diǎn)編號(hào)時(shí)可按照順時(shí)針方向順序編號(hào)以便于利用有限元的高效求解方法。坐標(biāo)迭代算法對(duì)初值是不敏感的，這將在后續(xù)算例中進(jìn)行驗(yàn)證。

迭代初始線形的設(shè)置是比較自由的。一般說來(lái)，存在唯一解時(shí)迭代初始線形越接近于目標(biāo)線形，迭代次數(shù)越少，但坐標(biāo)迭代算法還有一個(gè)特點(diǎn)，即在控制條件下的解不唯一時(shí)，會(huì)就近收斂至滿足現(xiàn)有條件的一個(gè)解，這也為設(shè)計(jì)提供了靈活和便利，因此，即使在設(shè)計(jì)初始對(duì)于唯一解的控制條件不確定時(shí)，也可以先根據(jù)需要設(shè)置一個(gè)與目標(biāo)接近的迭代初始線形，基于計(jì)算結(jié)果調(diào)整設(shè)計(jì)參數(shù)。主纜在自重荷載下呈懸鏈線，在均布荷載作用下呈拋物線，如果吊索力比較均勻且遠(yuǎn)大于主纜自重，成橋線形更接近于拋物線；反之，最終線形將介于拋物線與懸鏈線之間。實(shí)際上，拋物線也是懸鏈線方程解的一種特殊情況[15]，因此無(wú)論是對(duì)于成橋還是空纜，均可將主跨和邊跨的迭代初始線形按照合理拱軸懸鏈線方程解[15]的形式進(jìn)行設(shè)置，如圖3所示，懸鏈線參數(shù)對(duì)收斂效果的影響將通過算例進(jìn)行驗(yàn)證。

圖3 迭代收斂示意

同迭代初始線形一樣，初始軸力的設(shè)置也是比較自由的，既可以按照拋物線或者懸鏈線的公式計(jì)算，也可以為了方便將所有單元的軸力統(tǒng)一設(shè)置為一個(gè)值。但仍然需要注意，為了保證單向收斂，迭代初始軸力應(yīng)盡量取一個(gè)較大值：迭代初始軸力變化對(duì)于收斂效果的影響將在算例中進(jìn)行驗(yàn)證。

在設(shè)計(jì)模式下，還需要設(shè)置迭代初始的單元無(wú)應(yīng)力長(zhǎng)度。為了方便，可以直接設(shè)置為IP點(diǎn)連線后吊索間的直線段長(zhǎng)度，如圖4所示。

圖4 迭代初始無(wú)應(yīng)力長(zhǎng)度

(8)

將方程中的對(duì)角元中充一大數(shù)Am(如1020)，并對(duì)荷載項(xiàng)進(jìn)行修改，上式變?yōu)?/p>

(9)

(10)

需要注意的是，達(dá)到平衡狀態(tài)時(shí)，塔頂IP點(diǎn)的縱向是不能有約束的。

主纜節(jié)點(diǎn)坐標(biāo)和單元無(wú)應(yīng)力長(zhǎng)度的收斂容差，可按式(11)設(shè)置

(11)

式中：eps，epx分別為單元無(wú)應(yīng)力長(zhǎng)度和節(jié)點(diǎn)坐標(biāo)的收斂容差；s0i(n)為迭代第n階段第i號(hào)單元的無(wú)應(yīng)力長(zhǎng)度；xi(n)，yi(n)，zi(n)為迭代第n階段的第i號(hào)節(jié)點(diǎn)坐標(biāo)。

2.2 計(jì)算流程

根據(jù)前述基本原理和求解條件，懸索橋成橋和空纜線形的求解流程分別如圖5a，5b所示。

圖5 坐標(biāo)迭代算法流程

圖5a中，節(jié)點(diǎn)荷載信息包含主纜自重等效節(jié)點(diǎn)力以及吊索力；在形成總體關(guān)系矩陣時(shí)，針對(duì)順序編號(hào)的主纜可采用一維變帶寬存儲(chǔ)上三角陣[17]的策略縮小計(jì)算帶寬，進(jìn)一步提高計(jì)算效率。后續(xù)算例驗(yàn)證時(shí)只考慮主纜水平分力未給定的情況。

圖5b中，首先需要明確空纜線形的概念：空纜線形實(shí)際上應(yīng)該包含空纜自重和空纜無(wú)應(yīng)力兩種線形?？绽|自重線形表示未施加吊索力之前的主纜在自重狀態(tài)下的線形，這是實(shí)際存在的，可用于計(jì)算掛索長(zhǎng)度。但對(duì)后者來(lái)說，真正無(wú)應(yīng)力的線形是沒有意義的，故“無(wú)應(yīng)力”實(shí)際上是“小應(yīng)力”，小到可以忽略主纜的彈性伸長(zhǎng)，但保留了荷載比例；空纜無(wú)應(yīng)力線形是作為起點(diǎn)為非線性計(jì)算服務(wù)的：以成橋荷載縮放(本文取103倍)獲得的無(wú)應(yīng)力線形可以配合非線性有限元增大荷載系數(shù)正裝計(jì)算到成橋狀態(tài)，理論上荷載全部卸去(包括主纜自重)之后的鞍座預(yù)偏量應(yīng)該基于該線形計(jì)算；同理，以空纜自重荷載縮放成的無(wú)應(yīng)力線形則可作為起點(diǎn)計(jì)算到空纜自重狀態(tài)。因?yàn)槿N空纜線形的計(jì)算原理相同，本文將在后續(xù)算例中著重驗(yàn)證按成橋荷載比例的空纜無(wú)應(yīng)力線形。圖5中后處理1～3的具體步驟分別見表1～3。

表1 后處理1

表2 后處理2

表3 后處理3

3 算例分析

3.1 成橋與空纜線形計(jì)算

以大跨度地錨式懸索橋北盤江大橋?yàn)槔捎米鴺?biāo)迭代算法計(jì)算其成橋和空纜線形，并用梁?jiǎn)卧蜅U單元的非線性有限元方法進(jìn)行驗(yàn)證。

該橋跨度為192+636+192 m，主梁為鋼桁架，邊跨無(wú)吊索，如圖6a所示。為便于不同程序的比較，按照對(duì)稱性取一半結(jié)構(gòu)，對(duì)主纜統(tǒng)一劃分如圖6b所示，共72個(gè)單元和73個(gè)節(jié)點(diǎn)。主纜相關(guān)參數(shù)與吊索力分別見表4，5，跨中點(diǎn)吊索力取1/2實(shí)際大小。本文算法通過fortran95編制程序?qū)崿F(xiàn)。

圖6 北盤江大橋算例/m

程序單元容重/(kN/m3)彈性模量E/Kpa剪切模量G/Kpa軸向面積A/m2截面慣性矩Iy,Iz/m4梁格梁78.52.05×1087.90×1071.68×10-15.61×10-4本文桿78.52.05×108—1.68×10-1—ANSYS桿78.52.05×108—1.68×10-1—

表5 吊索力參數(shù) kN

非線性梁?jiǎn)卧绦颍翰捎每臻g6自由度直梁?jiǎn)卧?；為避免造成剛度矩陣奇異以及盡量符合實(shí)際主纜具有一定彎曲剛度的情況，將截面慣性矩取為真實(shí)值的1/4。經(jīng)過計(jì)算可驗(yàn)證截面慣性矩取其它更小值對(duì)線形結(jié)果影響很小，而剪切模量則對(duì)結(jié)果幾乎無(wú)影響。為了獲得成橋主纜的真實(shí)內(nèi)力，需要從迭代初始線形開始通過逐步增大荷載系數(shù)的理想非線性施工過程來(lái)獲得成橋線形，再根據(jù)與設(shè)計(jì)控制坐標(biāo)的差值對(duì)迭代初始線形進(jìn)行修正，反復(fù)計(jì)算直到收斂；鞍座采用能夠縱向滑動(dòng)的橡膠單元模擬。該程序提供了一個(gè)主纜線形修改工具，以保證迭代初始線形在指定修正坐標(biāo)不變的情況下滿足荷載比例的力學(xué)平衡，使非線性計(jì)算順利進(jìn)行。一共經(jīng)過5次迭代修正后達(dá)到收斂，計(jì)算后同時(shí)獲得了空纜無(wú)應(yīng)力線形。如圖7所示,線形計(jì)算在全橋中進(jìn)行。

圖7 非線性梁?jiǎn)卧绦蛴?jì)算成橋與空纜線形

ANSYS：采用通用有限元程序ANSYS[12]的目的是基于同樣的桿單元來(lái)驗(yàn)證本文算法計(jì)算的空纜無(wú)應(yīng)力線形與成橋線形之間的力學(xué)相關(guān)性。采用空間桿單元LINK10模擬一半主纜結(jié)構(gòu)，所有節(jié)點(diǎn)施加面外位移約束，塔頂縱向自由。計(jì)算過程：從坐標(biāo)迭代算法提供的空纜無(wú)應(yīng)力線形(按自重荷載比例)開始，各單元設(shè)置一個(gè)較小的初始應(yīng)變(10-6)，打開大變形和應(yīng)力剛度開關(guān)，迭代子步數(shù)設(shè)為20；第一荷載步考慮自重，第二荷載步考慮吊索力，計(jì)算完畢后獲得成橋線形，模型如圖8所示。

圖8 基于桿單元的ANSYS計(jì)算模型

以A、B、C分別表示非線性梁?jiǎn)卧绦?、本文程序和ANSYS中非線性桿單元的計(jì)算結(jié)果，則關(guān)于成橋和空纜線形的比較見表6，關(guān)于成橋單元內(nèi)力的比較見表7。限于篇幅，只規(guī)律性地列出部分單元和節(jié)點(diǎn)的計(jì)算結(jié)果，其中第28號(hào)節(jié)點(diǎn)表示塔頂IP點(diǎn)，73號(hào)為跨中節(jié)點(diǎn)。

從表6，7中可知，成橋線形中梁?jiǎn)卧c本文算法豎向坐標(biāo)最大差值21 mm，而軸力最大差值不超過0.10%；空纜線形時(shí)x和y坐標(biāo)的最大誤差分別為20，26 mm；兩種計(jì)算方法得到的鞍座預(yù)偏量則基本相同(348，349 mm)。對(duì)于主跨636 m的大型懸索橋來(lái)說，這種誤差是非常小的，可認(rèn)為本文方法與梁?jiǎn)卧M的結(jié)果相近；從誤差趨勢(shì)來(lái)看，考慮了抗彎剛度的梁?jiǎn)卧獣?huì)具有相對(duì)較大的整體豎向剛度；此外，誤差還包括主纜單元初始無(wú)應(yīng)力長(zhǎng)度的誤差和非線性有限元迭代過程中的誤差積累。而ANSYS的成橋結(jié)果則與本文結(jié)果非常接近，最大豎向坐標(biāo)差只有2 mm，而軸力最大誤差則不足0.01%，不僅驗(yàn)證了本文算法的精確性，也進(jìn)一步驗(yàn)證了本文算法從空纜無(wú)應(yīng)力到成橋線形的力學(xué)關(guān)系是正確的。由于使用的單元特性相同并具有同樣的單元無(wú)應(yīng)力長(zhǎng)度，誤差只是來(lái)源于ANSYS中非線性有限元迭代計(jì)算的誤差積累，本文算法的結(jié)果反而更精確。

表6 坐標(biāo)迭代算法與非線性方法的線形比較

表7 坐標(biāo)迭代算法與非線性方法的內(nèi)力比較

注：梁?jiǎn)卧卧蠖溯S力

綜上所述，應(yīng)用非線性有限元方法步驟較為復(fù)雜，而坐標(biāo)迭代算法則簡(jiǎn)單直接；上述對(duì)比則說明了坐標(biāo)迭代算法具有較高的精度，并且在同等初始無(wú)應(yīng)力長(zhǎng)度和不考慮抗彎剛度影響的情況下，可認(rèn)為該方法的計(jì)算結(jié)果是精確的。

3.2 收斂特性分析

以北盤江大橋?yàn)槔?、邊跨均引用合理拱軸線方程的解[15]作為迭代初始線形,重點(diǎn)考察中跨迭代初始線形的參數(shù)變化對(duì)收斂結(jié)果的影響。如果單元軸力的迭代初始值Ti不變，則拱軸系數(shù)m和矢高f*成為影響收斂速度的兩個(gè)參數(shù)；同理，在m和f*不變的情況下，改變Ti同樣可以影響收斂結(jié)果。若成橋設(shè)計(jì)矢高為f，則隨著m，f*，Ti的改變，坐標(biāo)迭代算法計(jì)算成橋線形的總迭代次數(shù)如表8，9所示。其中，表8中Ti按拋物線公式近似取值；表9中取m=2時(shí)，f*=1/2f。

表9 迭代次數(shù)隨Ti的變化

從表8，9中可知，當(dāng)Ti和m一定，f*的設(shè)置并不影響迭代次數(shù)，故f*的的設(shè)置是比較自由的；當(dāng)m發(fā)生變化時(shí)，也只在兩種情況下會(huì)改變迭代次數(shù)：m=1和m≠1，且m=1時(shí)的迭代次數(shù)更少。根據(jù)合理拱軸方程：當(dāng)m=1時(shí)，線形退化為拋物線；m≠1，線形為普通懸鏈線。也就是說，在迭代初始軸力一定的情況下，影響迭代次數(shù)的因素僅跟線型有關(guān)，拋物線因?yàn)楦咏诔蓸蚓€形，所以迭代次數(shù)更少。盡管將迭代初始線形設(shè)置成普通懸鏈線迭代次數(shù)會(huì)略多，但本文程序是基于桿單元的線性有限元窄帶寬計(jì)算，在收斂速度上也僅僅只是相差1～2 s而已。

當(dāng)m和f*一定時(shí)，采用拋物線法計(jì)算的迭代初始軸力與統(tǒng)一設(shè)置的軸力對(duì)迭代次數(shù)的影響不大；只要Ti大于成橋目標(biāo)軸力(為保證單向收斂)，都能夠順利收斂且其取值對(duì)迭代結(jié)果影響較小。

綜上所述，基于結(jié)構(gòu)力學(xué)原理的坐標(biāo)迭代算法，對(duì)于初值是不敏感的，大部分情況下都能快速收斂到正確解。只需注意將初始軸力取為一個(gè)大于成橋目標(biāo)軸力的值，這是容易做到的。根據(jù)以上分析，力學(xué)概念明確和收斂容易也是坐標(biāo)迭代算法的優(yōu)勢(shì)所在。

4 結(jié) 論

本文提出了一種新的數(shù)值迭代方法坐標(biāo)迭代算法用于解決懸索橋的主纜找形問題。該方法基于合理拱軸線的彎矩平衡特征構(gòu)建直接求解主纜坐標(biāo)的平衡迭代方程，其核心是結(jié)構(gòu)力學(xué)的基本原理，具有清晰的力學(xué)概念。該算法將線性有限元的計(jì)算方法巧妙地引入到索結(jié)構(gòu)的計(jì)算中，從而可以整體求解主纜的成橋和空纜線形，在求解步驟的方便性上優(yōu)于現(xiàn)有非線性方法和分段懸鏈線法。通過算例，驗(yàn)證了本文算法具有較高的精度；收斂性分析則證明了該算法對(duì)于初值是不敏感的，即不依靠特定數(shù)學(xué)方法處理仍能穩(wěn)定、可控地收斂，因而在計(jì)算效率上優(yōu)于分段懸鏈線法。

因此，計(jì)算準(zhǔn)確、方便和收斂容易是坐標(biāo)迭代算法的特點(diǎn)。本文方法可為懸索橋的設(shè)計(jì)和施工提供參考；同時(shí)，稍加擴(kuò)展便可用于考慮鞍座形狀影響的纜形計(jì)算以及空間索纜懸索橋的計(jì)算。

[1] 羅喜恒, 肖汝誠(chéng), 項(xiàng)海帆. 懸索橋理想恒載狀態(tài)的計(jì)算方法研究[J]. 橋梁建設(shè), 2008, (4): 31-35.

[2] 張其林. 索和膜結(jié)構(gòu)[M]. 上海: 同濟(jì)大學(xué)出版社, 2002.

[3] 檀永剛. 自錨式懸索橋的力學(xué)特性與施工控制研究[D]. 大連: 大連理工大學(xué), 2007.

[4] Kim H K, Lee M J, Chang S P. Non-linear shape-finding analysis of a self-anchored suspension bridge[J]. Engneering Structures, 2002, 24(12), 1547-1559.

[5] Kong W, Xiong X, Xue Q L. Study of the Initial Prestress for Guyed Transmission Tower[C]//Proceedings of 2009 International Workshop on Information Security and Application. 2009: 405-408.

[6] 袁行飛, 董石麟. 二節(jié)點(diǎn)曲線索單元非線性分析[J]. 工程力學(xué), 1999, 16(4): 59-64.

[7] Knudson W C. Static and Dynamic Analysis of Cable Net Structures[D]. California: University of California, Berkeley, 1971.

[8] 景天虎, 李青寧. 懸索橋主纜成橋線形確定的有限元新算法[J]. 世界橋梁, 2012, 40(1): 42-46.

[9] Kim K S, Lee H S. Analysis of target configurations under dead loads for cable-supported bridges[J]. Computers & Structures, 2001, 79(29-30), 2681-2692.

[10] Kim H K, Kim M Y. Efficient combination of a TCUD method and an initial force method for determining initial shapes of cable-supported bridges[J]. International Journal of Steel Structures, 2012, 12(2), 157-174.

[11] 唐茂林. 大跨度懸索橋空間幾何非線性分析與軟件開發(fā)[D]. 成都: 西南交通大學(xué), 2003.

[12] 羅喜恒, 肖汝誠(chéng), 項(xiàng)海帆. 空間纜索懸索橋的主纜線形分析[J]. 同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版), 2004, 32(10): 1349-1354.

[13] 王曉明, 郝憲武, 段瑞芳. 基于Steffens-Newton法的空間索形懸索橋初始平衡狀態(tài)分析[J]. 計(jì)算力學(xué)學(xué)報(bào), 2011, 28(5): 717-722.

[14] 李傳習(xí), 劉光棟, 柯紅軍. 懸索橋主纜系統(tǒng)數(shù)值解析法計(jì)算的一種收斂算法[J]. 工程力學(xué), 2008, 25(7): 66-73.

[15] 姚玲森. 橋梁工程[M]. 北京: 人民交通出版社, 2005.

[16] 朱伯芳. 有限單元法原理與應(yīng)用(第3版)[M]. 北京: 中國(guó)水利水電出版社, 2009.

[17] 于二青. 快速有限元分析程序的系統(tǒng)設(shè)計(jì)與算法優(yōu)化[D]. 上海: 上海交通大學(xué), 2009.