劉 超，詹海鵬，朱 怡

（1.同濟大學土木工程學院，上海 200092；2.無錫市政設計研究院有限公司，江蘇 無錫 214072）

懸索結構的工程應用歷史悠久，是最常見的大跨度結構形式，應用范圍覆蓋到橋梁工程、房屋建筑等。懸索橋是最典型的懸索結構，通過其主纜的軸向力來抵抗外荷載作用，可以充分利用鋼絲等材料的拉伸強度。當主纜采用高強度材料的時候，能夠大幅度減輕自重比例，使得懸索橋能夠實現(xiàn)較大的跨度。懸索橋主纜不同于一般的剛性結構，在無應力狀態(tài)下的結構初始形狀是不確定的，在使用狀態(tài)下的主纜軸力提供幾何剛度。因此，懸索橋主纜存在顯著的非線性特征，這也使得主纜找形問題成為其結構分析的重點和難點。針對懸索橋結構，Ochsendorf等［1］指出，Euler提出的拋物線理論，即主纜在沿跨均布荷載下的形狀為拋物線，其水平分力為恒定值，忽略了主纜重力沿著主纜曲線分布的特征，僅適用于跨度較小的橋梁。隨著懸索橋跨度的不斷增大，拋物線形狀的主纜受力與實際受力不相符，誤差較大，不能滿足設計要求。因此，懸鏈線模型［2］、分段懸鏈線法［3］、整體力學分析法［4］等解析方法被先后提出，并廣泛應用至今。近些年，有限元法［5-7］得到快速發(fā)展，成為分析工具之一。但非線性有限元方法在懸索結構找形過程中存在兩個問題：一是多迭代子步帶來的誤差積累，二是收斂的不穩(wěn)定性［8］。此外，宋旭明等［9］提出需要仔細考慮主梁的壓縮對主纜無應力長度產生的影響。Kim等［10］考慮加勁梁和主塔的軸向壓縮，基于彈性懸鏈線單元和牛頓迭代法，以無應力長度為未知數，提出了TUCD（target configuration under dead load）及其改進算法，但是該算法對初始值要求較高。我國懸索橋規(guī)范規(guī)定的纜形計算方法即為分段懸鏈線法［11］，此方法迭代參數較少，求解速度快，是目前應用最普遍的方法。但該方法對迭代初值敏感，迭代過程容易發(fā)散，且需要推導復雜的懸鏈線方程。

力密度法［12-17］由Schek于1974年提出，主要用于索網找形分析，因其簡易的計算理念而得到快速發(fā)展。動力松弛法［18］自從1965年提出以來不斷得到發(fā)展，其本質是利用特殊條件下的動力方法來解決靜力問題。二者均廣泛應用在索網和膜結構設計中。力密度法優(yōu)勢在于能夠利用線性手段求解非線性問題，計算更加簡便。索網結構找形與懸索橋主纜找形原理相似，本質上都是尋找適應邊界條件的平衡狀態(tài)。本文基于傳統(tǒng)的力密度法，結合主纜無應力長度的概念，增加迭代求解過程的控制條件，提出了適合懸索橋主纜找形的改進力密度法。利用靜力平衡和幾何協(xié)調關系，采用力密度作為初始參數，形成力與節(jié)點坐標的線性關系，建立整體平衡方程，直接求解主纜節(jié)點坐標。該方法將非線性問題轉化為線性方程的求解，步驟簡便，收斂穩(wěn)定，計算精度滿足工程需求。

1 力密度法基本原理

力密度法［12］基本參數有結構的拓撲關系、邊界約束條件、荷載以及平衡條件。懸索結構的拓撲關系規(guī)定了節(jié)點單元的排列規(guī)則和序列，采用單元-節(jié)點的關聯(lián)作為其拓撲矩陣C t，定義為

式中：i，j為節(jié)點號，i，j=1～nt，節(jié)點分為兩部分，nf表示固定節(jié)點，nd表示自由節(jié)點，二者的關系可以表示為nt=nd+nf；e為單元編號，e=1～m；拓撲矩陣C t的行數對應著單元號，列數對應著節(jié)點號，即為m行和nt列。根據固定節(jié)點和自由節(jié)點可將C t分塊為C和Cf，即C t=[CCf]。

設結構中所有節(jié)點的三維坐標為x t、y t、z t，均為nt維矢量，坐標矢量同樣可分為固定節(jié)點與自由節(jié)點相鄰節(jié)點的坐標差矢量u，v，w為

設單元長度和軸力分別為m維矢量l和t，與之相對應的m維對角矩陣設為L和T，對應3個坐標方向的外荷載矢量分別為p x、p y、p z。結構處于平衡狀態(tài)下，每個節(jié)點力的矢量和為零，因此可以得到以下平衡方程：

式中：U、V、W為相鄰節(jié)點的坐標差矢量；u、v、w對應的m維對角矩陣。

令q=L-1t，其物理含義為單元軸力和單元長度的比值，即力密度的定義，可以得到

利用關系式

將式（2）和式（5）帶入式（4）得

設D=CTQC，即為C的廣義高斯變換，同理設Df=CTQCf，可得

式（7）為力密度法求解的基本方程，并且是一組能夠求解自由節(jié)點坐標的線性方程。對于一個給定的荷載體系和固定節(jié)點，對應每一組力密度，就可以求得該組力密度下結構的唯一平衡形狀。由于對角矩陣Q元素為力密度值，對于純受拉的主纜體系沒有交叉節(jié)點，且設置受拉為正，即滿足每個單元力密度大于0，則D是正定矩陣。

對于給定的結構拓撲矩陣，在已知的荷載和邊界條件下，不同的力密度對應著不同的平衡狀態(tài)。力密度的變化與平衡狀態(tài)的變化是相對應的，同時對于線性方程的求解效率較高。因此，力密度是懸索結構適合的描述參數。

2 改進力密度法

工程實踐已經證明，懸索橋的主纜在單元長度合適的情況下，采用只受軸向力的桿單元模擬主纜索結構具有較高的精度，并且單個桿單元的垂度效應可以忽略。因此，主纜找形計算假定如下：

（1）主纜為完全柔性，只受軸向拉力作用，忽略主纜抗彎剛度的影響。

（2）單個桿單元的自重均布荷載根據有限元基本理論等分到單元兩端節(jié)點形成等效節(jié)點力。

（3）桿單元張拉過程中截面的面積保持不變。

（4）主纜單元應力在彈性范圍內滿足胡克定律。

（5）主纜錨固點位置、吊桿位置、跨中垂點以及吊桿力為設計已知量。

傳統(tǒng)的力密度法只是簡單給定力密度，尚無法滿足懸索橋主纜找形特定的限制條件，此時找形問題就轉化為需滿足附加控制條件的非線性問題，需采用特定的方式進行迭代求解。針對主纜體系，引入無應力長度的概念，基于不同的設計需求建立滿足限定條件的迭代格式，從而求解出主纜找形問題中的成橋纜形和空纜纜形。

2.1 成橋纜形計算

在平衡方程中，兩側錨點和塔頂IP點為已知固定節(jié)點，其余節(jié)點可作為自由節(jié)點，吊桿力作為外荷載的一部分也是已知量。在成橋纜形求解問題中，主纜無應力長度和自由節(jié)點坐標為未知量，迭代過程中保持更新，其余已知條件保持不變。

由于主纜纜形未知，其自重也是未知的，需要設置初始纜形提供初始單元長度。初始單元軸力可設為一個較大值，易于收斂。依據二者數值得到初始力密度。初始纜形根據設計需求可設置為錨點、塔頂IP點與垂點之間的相連直線，以吊桿位置為界將主纜分成若干個單元和節(jié)點。初始纜形中每一段的主纜長度已知，可以用主纜的截面積乘上節(jié)段長度再乘上重度表示節(jié)段主纜的重力，便可得到成橋纜形計算的平衡方程。根據式（8）求解得到新的自由節(jié)點坐標。以節(jié)點z坐標為例，其平衡方程如式（10）所示。

式中：zf為固定節(jié)點z坐標向量，始終保持不變；qe為各單元力密度值，e是單元編號；z n為第n階段待求

解的自由節(jié)點z坐標向量為平衡矩陣，由第n-1階段的單元力密度qen-1確定；p z(l en-1)為結構z方向外荷載，由于吊桿力已知，所以外荷載由第n-1階段的單元長度l en-1確定。

由于假設的直線纜形不可能保持平衡，因此需要進入下一迭代。以本階段的纜形作為初始狀態(tài)，重新代入求解平衡方程。主纜重力是對應上一階段纜形的節(jié)段重力，在計算得到新的纜形之后，由于主纜重力發(fā)生變化，結構體系將再次不平衡，需要反復迭代計算，直到相鄰節(jié)段的纜形坐標變化量達到收斂要求。在主纜無應力長度和坐標均未知的情況下，錨點和塔頂的邊界條件不足以確定主纜的成橋線型。理論上，可以有無窮組力密度滿足方程的要求，即有無窮組成橋纜型解，因此需要根據已知設計條件（垂點坐標和吊桿縱向坐標）來獲得跟設計參數匹配的成橋纜形。在每次迭代過程中設定這些自由節(jié)點的設計坐標值，進入下一階段的迭代計算，最終平衡方程求解獲得的自由節(jié)點坐標會逐漸收斂至滿足設計參數的纜形。

式中：x ni、yi n、zni為第n階段i節(jié)點的三維坐標；x ni-1、為第n-1階段i節(jié)點的三維坐標。

成橋纜形求解具體流程如圖2所示。

圖2 成橋纜形求解流程圖Fig.2 Flowchart of cable shape finding in finished stage

2.2 空纜纜形求解

空纜纜形與成橋纜形的求解相似，初始纜形和初始軸力與成橋纜形分析中設置相同?？绽|狀態(tài)下，吊桿力為零，只受主纜重力作用，并且空纜纜形的垂度未知。由于從空纜到成橋狀態(tài)主纜的錨點不變，因此全橋主纜的無應力長度不變。在成橋纜形確定后，每一個主纜單元軸力為已知，成橋狀態(tài)的單元無應力長度可表示為

式中：s0e為單元無應力長度；se為成橋狀態(tài)單元有應力長度；te為主纜單元軸力；A為主纜截面積；E為主纜彈性模量。

利用主纜在成橋狀態(tài)和空纜狀態(tài)時的無應力長度相等這一控制條件，即可求解空纜纜形。在迭代過程中保持固定節(jié)點不變，利用已知的單元無應力長度來控制迭代過程中單元軸力和長度。以節(jié)點z坐標為例，其平衡方程如式（13）所示。

式中：pz在外荷載迭代過程中保持不變；qen-1可以看作是單元無應力長度的函數；l e表示迭代過程中單元有應力長度由已知的無應力長度求得；te表示單元軸力。其余符號參考式（10）的說明，收斂要求見式（11）。

由于主纜在成橋和空纜狀態(tài)下內力不同，因此需要設置索鞍預偏量，消除因受力不同導致的主纜對主塔的水平力。張拉過程中纜形會隨著吊桿的逐步張拉，回到成橋狀態(tài)的線形。對于索鞍的預偏量，可以在纜形計算的基礎上，先假設預偏量為Ap，假設偏向邊跨，預偏量調整值為Δ=Ap2。在此假設下，分別計算空纜狀態(tài)和成橋狀態(tài)下的主纜無應力長度。若前者大于后者，調整預偏量為Ap=Ap-Δ；若前者小于后者，調整預偏量為Ap=Ap+Δ，并設置下一階段預偏量調整值為Δ=Δ2。如此反復計算，直到空纜纜形和成橋纜形的無應力長度之差滿足精度要求。

3 工程應用

3.1 平面纜形懸索橋實例

如圖4所示，京杭大運河懸索橋位于徐州市金山橋經濟開發(fā)區(qū)境內，為雙索面地錨式懸索結構。全橋長245.0 m，橋寬5.8 m，主跨為130.0 m的鋼桁架梁。主塔采用鋼桁架結構。主纜為1根預制平行鋼絲索股，垂跨比1/10，垂度13 m，兩主纜橫橋向中心間距為5.8 m。全橋共有50根吊索，吊索間距為5 m，吊索索力在139～176 kN之間。該橋已于2020年11月施工完成。

利用本文方法計算該橋的成橋和空纜纜形?？缇┖即筮\河懸索橋的成橋狀態(tài)如圖5所示。在橋梁主纜架設和張拉期間對全橋纜形進測量，空纜和成橋實測纜形如圖6所示。本橋主纜彈性模量為1.99×105MPa，截面積為0.01 m2。

圖5 主纜成橋狀態(tài)Fig.5 Main cable in finished satge

圖6 空纜和成橋實測坐標Fig.6 Measured coordinates of main cable in unloaded stage and finished satge

根據流程圖2、圖3，基于Matlab軟件對空纜纜形與成橋纜形進行編程計算。同時，通過Midas Civil軟件建立有限元模型，如圖7所示。比較兩種算法之間的纜形誤差，并得到本文程序所求得的空纜和成橋纜型與實測值的誤差，如圖8和圖9所示。

圖3 空纜纜形求解流程圖Fig.3 Flowchart of cable shape finding in unloaded stage

圖 4京杭大運河懸索橋總體布置圖（單位：m）Fig.4 Layout of Beijing-Hangzhou Grand Canal suspension bridge(unit:m)

圖7 有限元模型Fig.7 Finite element structural model

圖8 空纜纜形比較Fig.8 Comparison of cable shapes in unloaded stage

從圖8和圖9可以看出，論文方法與有限元方法計算的空纜和成橋纜形誤差最大分別為2.0 cm和2.3 cm。與實測值相比，空纜的纜形誤差最大為5.6 cm。分析可知，產生此項誤差的主要原因是空纜時荷載較小，導致主纜存在局部扭曲，因此測量存在一定誤差。對于成橋纜形，論文方法計算的纜型理論值與實測纜形的最大誤差為5.8 cm，主要原因是吊桿力理論值與實際值存在一定誤差?？傮w而言，論文方法計算的纜形誤差相對較小。

圖9 成橋纜形比較Fig.9 Comparison of cable shapes in finished stage

3.2 空間纜形懸索橋實例

為驗證論文方法對于空間纜形的適用性，選取韓國永宗大橋為工程背景，具體數據見參考文獻［19］。永宗大橋為自錨式空間懸索橋，橋梁總體布置和具體參數見圖10。吊索下錨點處橋面豎曲線半徑為14 472 m，下錨點之間的橫向距離為31.92 m；端吊索力豎向分量為3 905.30 kN，其余邊跨吊索力豎向分量為3 142.86 kN，中跨吊索力豎向分量均為3 088.96 kN；邊跨錨點和塔頂IP點高差為66.51 m。該橋主纜彈性模量為1.99×105MPa，截面積為0.135 5 m2。

圖1 懸索橋示意圖Fig.1 Schematic diagram of suspension bridge

圖10 永宗橋總體布置圖（單位：m）Fig.10 General layout of Yongzong Bridge(unit:m)

空間纜形相較于平面纜形增加了Y方向的自由度，同樣選取錨點到IP點之間的連線作為找形程序的初始纜形，見圖11中的虛線線形，邊界條件和之前所述保持不變。已知吊索力豎向分量，但是吊索力橫向分量未知，可以利用已知的吊索錨點與主纜節(jié)點坐標之間的位置關系，通過空間力矢量的分解，得到迭代過程中作用于主纜的吊索力橫向分量。自錨式懸索橋梁體在巨大的索力作用下，會產生較大的壓縮量，本文通過有限元計算位移，改變找形邊界坐標，可以消除主梁壓縮量對于纜形誤差的影響。

通過本文程序求解的空間纜形如圖11所示，由空纜的直線經過迭代收斂至成橋的纜形。圖12和圖13分別給出了詳細的平面坐標、立面坐標以及線形迭代變化過程。

圖11 空間纜形求解過程圖Fig.11 Solving process diagram of spatial cable shape

圖12 主纜平面線形迭代過程Fig.12 Iterative process of main cable plane alignment

圖13 主纜立面線形迭代過程Fig.13 Iterative process of main cable elevation

圖14給出了本文程序和文獻［19］計算值的比較結果，Y坐標誤差最大為3.3 cm，Z坐標誤差最大為3.0 cm?？梢?，本文方法在解決空間纜形計算問題中同樣具有較高的精度。

4 結論

（1）論文提出的改進力密度法用于解決懸索橋主纜找形問題。將該算法應用于地錨式懸索橋實際工程的施工監(jiān)控，并與有限元結果進行對比驗證。結果顯示，論文方法計算的纜形與實測纜形較為接近。

（2）采用論文提出的方法求解自錨式懸索橋的空間纜形，并與文獻方法進行比較，驗證了該算法的適用性。

（3）論文提出的方法收斂穩(wěn)定，初始參數簡單，計算速度較快，適用于平面和空間懸索橋的主纜找形，值得進一步推廣應用。

作者貢獻聲明：

劉 超：提出改進的力密度法，指導算法推導。

詹海鵬：算法的編程和算例的計算，論文撰寫。

朱 怡：論文工程背景的設計。