改進(jìn)力密度法在懸索橋主纜找形中的應(yīng)用

劉 超，詹海鵬，朱 怡

（1.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2.無(wú)錫市政設(shè)計(jì)研究院有限公司，江蘇 無(wú)錫 214072）

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

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

1 力密度法基本原理

力密度法［12］基本參數(shù)有結(jié)構(gòu)的拓?fù)潢P(guān)系、邊界約束條件、荷載以及平衡條件。懸索結(jié)構(gòu)的拓?fù)潢P(guān)系規(guī)定了節(jié)點(diǎn)單元的排列規(guī)則和序列，采用單元-節(jié)點(diǎn)的關(guān)聯(lián)作為其拓?fù)渚仃嘋 t，定義為

式中：i，j為節(jié)點(diǎn)號(hào)，i，j=1～nt，節(jié)點(diǎn)分為兩部分，nf表示固定節(jié)點(diǎn)，nd表示自由節(jié)點(diǎn)，二者的關(guān)系可以表示為nt=nd+nf；e為單元編號(hào)，e=1～m；拓?fù)渚仃嘋 t的行數(shù)對(duì)應(yīng)著單元號(hào)，列數(shù)對(duì)應(yīng)著節(jié)點(diǎn)號(hào)，即為m行和nt列。根據(jù)固定節(jié)點(diǎn)和自由節(jié)點(diǎn)可將C t分塊為C和Cf，即C t=[CCf]。

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

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

式中：U、V、W為相鄰節(jié)點(diǎn)的坐標(biāo)差矢量；u、v、w對(duì)應(yīng)的m維對(duì)角矩陣。

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

利用關(guān)系式

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

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

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

對(duì)于給定的結(jié)構(gòu)拓?fù)渚仃嚕谝阎暮奢d和邊界條件下，不同的力密度對(duì)應(yīng)著不同的平衡狀態(tài)。力密度的變化與平衡狀態(tài)的變化是相對(duì)應(yīng)的，同時(shí)對(duì)于線性方程的求解效率較高。因此，力密度是懸索結(jié)構(gòu)適合的描述參數(shù)。

2 改進(jìn)力密度法

工程實(shí)踐已經(jīng)證明，懸索橋的主纜在單元長(zhǎng)度合適的情況下，采用只受軸向力的桿單元模擬主纜索結(jié)構(gòu)具有較高的精度，并且單個(gè)桿單元的垂度效應(yīng)可以忽略。因此，主纜找形計(jì)算假定如下：

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

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

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

（4）主纜單元應(yīng)力在彈性范圍內(nèi)滿足胡克定律。

（5）主纜錨固點(diǎn)位置、吊桿位置、跨中垂點(diǎn)以及吊桿力為設(shè)計(jì)已知量。

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

2.1 成橋纜形計(jì)算

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

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

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

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

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

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

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

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

2.2 空纜纜形求解

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

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

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

式中：pz在外荷載迭代過(guò)程中保持不變；qen-1可以看作是單元無(wú)應(yīng)力長(zhǎng)度的函數(shù)；l e表示迭代過(guò)程中單元有應(yīng)力長(zhǎng)度由已知的無(wú)應(yīng)力長(zhǎng)度求得；te表示單元軸力。其余符號(hào)參考式（10）的說(shuō)明，收斂要求見(jiàn)式（11）。

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

3 工程應(yīng)用

3.1 平面纜形懸索橋?qū)嵗?/h3>

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

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

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

圖6 空纜和成橋?qū)崪y(cè)坐標(biāo)Fig.6 Measured coordinates of main cable in unloaded stage and finished satge

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

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

圖 4京杭大運(yùn)河懸索橋總體布置圖（單位：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可以看出，論文方法與有限元方法計(jì)算的空纜和成橋纜形誤差最大分別為2.0 cm和2.3 cm。與實(shí)測(cè)值相比，空纜的纜形誤差最大為5.6 cm。分析可知，產(chǎn)生此項(xiàng)誤差的主要原因是空纜時(shí)荷載較小，導(dǎo)致主纜存在局部扭曲，因此測(cè)量存在一定誤差。對(duì)于成橋纜形，論文方法計(jì)算的纜型理論值與實(shí)測(cè)纜形的最大誤差為5.8 cm，主要原因是吊桿力理論值與實(shí)際值存在一定誤差?？傮w而言，論文方法計(jì)算的纜形誤差相對(duì)較小。

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

3.2 空間纜形懸索橋?qū)嵗?/h3>

為驗(yàn)證論文方法對(duì)于空間纜形的適用性，選取韓國(guó)永宗大橋?yàn)楣こ瘫尘?，具體數(shù)據(jù)見(jiàn)參考文獻(xiàn)［19］。永宗大橋?yàn)樽藻^式空間懸索橋，橋梁總體布置和具體參數(shù)見(jiàn)圖10。吊索下錨點(diǎn)處橋面豎曲線半徑為14 472 m，下錨點(diǎn)之間的橫向距離為31.92 m；端吊索力豎向分量為3 905.30 kN，其余邊跨吊索力豎向分量為3 142.86 kN，中跨吊索力豎向分量均為3 088.96 kN；邊跨錨點(diǎn)和塔頂IP點(diǎn)高差為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方向的自由度，同樣選取錨點(diǎn)到IP點(diǎn)之間的連線作為找形程序的初始纜形，見(jiàn)圖11中的虛線線形，邊界條件和之前所述保持不變。已知吊索力豎向分量，但是吊索力橫向分量未知，可以利用已知的吊索錨點(diǎn)與主纜節(jié)點(diǎn)坐標(biāo)之間的位置關(guān)系，通過(guò)空間力矢量的分解，得到迭代過(guò)程中作用于主纜的吊索力橫向分量。自錨式懸索橋梁體在巨大的索力作用下，會(huì)產(chǎn)生較大的壓縮量，本文通過(guò)有限元計(jì)算位移，改變找形邊界坐標(biāo)，可以消除主梁壓縮量對(duì)于纜形誤差的影響。

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

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

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

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

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

4 結(jié)論

（1）論文提出的改進(jìn)力密度法用于解決懸索橋主纜找形問(wèn)題。將該算法應(yīng)用于地錨式懸索橋?qū)嶋H工程的施工監(jiān)控，并與有限元結(jié)果進(jìn)行對(duì)比驗(yàn)證。結(jié)果顯示，論文方法計(jì)算的纜形與實(shí)測(cè)纜形較為接近。

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

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

作者貢獻(xiàn)聲明：

劉 超：提出改進(jìn)的力密度法，指導(dǎo)算法推導(dǎo)。

詹海鵬：算法的編程和算例的計(jì)算，論文撰寫(xiě)。

朱 怡：論文工程背景的設(shè)計(jì)。