高慶飛,洪能達,郭斌強,劉 洋,馬其魯
(1.哈爾濱工業大學 交通科學與工程學院,哈爾濱 150090;2. 浙江無限元組合結構橋梁設計有限公司,杭州 310000;3.浙江省交通規劃設計研究院,杭州 310000)
懸索橋由于具有受力性能合理、抗震性能較好、大跨輕盈的特點,是特大跨徑橋梁的首選橋型. 當橋梁跨徑超過1 200 m時,懸索橋被認為是最具有競爭力的橋型方案. 隨著跨徑的增大,懸索橋的優勢就越明顯[1]. 主纜是懸索橋至關重要的受力構件. 懸索橋施工質量的指標很大程度上要考慮橋梁施工完成后的實際主纜線形與設計主纜線形是否一致[2-3]. 對于特大跨徑懸索橋,主纜的安全系數要求不小于2.5,這必然要對懸索橋進行精確的求解. 目前的懸索橋計算分析方法已經比較成熟,但由于計算模型在假設時精度達不到要求,所以在處理一些局部受力問題時,懸索橋的計算內力和線形與實際受力與變形有一定偏差[4]. 懸索橋的索鞍起到了轉向或分散主纜的作用,并限制了主纜的變形. 對于懸索橋的成橋階段,主纜在進出鞍座處與鞍座是處于相切狀態的,兩者在切點間是緊密接觸的,主纜在鞍座內不能有任何的相對滑移. 而在施工階段,考慮到施工的實際情況,主纜與鞍座接觸的切點是不斷變化位置的,主纜與鞍座間存在著接觸非線性[5]. 在進行懸索橋計算時,如果不能準確考慮主纜與鞍座間復雜的接觸關系,主纜和吊索的下料長度將會減小,主纜架設控制標高的精度也會相應降低.
基于有限單元法的懸索橋的施工分析計算中,大多數方法是通過利用成橋理論點來建模分析. 但這不能確保主纜與鞍座處于相切狀態,從而造成了主纜與鞍座出現相交或分離[5]. 為了充分考慮鞍座對主纜的約束作用,文獻[6]采用桿單元來模擬. 橋塔頂部采用兩個桿單元模擬,通過主從約束的方式與鞍座桿單元連接. 通過修改桿單元的節點坐標,來模擬鞍座的頂推作用,改變鞍座桿單元的無應力長度,從而使得主纜在脫離點的斜率等于鞍槽曲線上此點處的斜率,進而得到切點的位置. 文獻[7]則采用4個梁單元來模擬鞍座及其頂推作用,用一個直桿代替散索鞍,該方法與文獻[6]的方法實質是相同的. 上述方法在修正懸空段主纜的無應力長度時,未考慮切點相對于鞍座位置的變化,使得懸索橋求解精度較低. 另外,跨內主纜線形由假定的切點坐標迭代求出后,還必須驗證主纜與鞍座是否相切,否則需重新迭代. 因此,該方法計算效率也不高.
有形狀的索鞍在其表面上包裹并支撐主纜,在施工過程中,主纜與索鞍間的切點位置在一直變化著. 在懸索橋計算分析中,由理論頂點假設推導出的主纜線形而得到的考慮索鞍曲線修正后實際的主纜形狀是比較繁瑣的,其計算關鍵在于找到主纜與索鞍間的相切點. 目前常用的解析算法是求解八元非線性方程組. 即借助懸鏈線公式得到4個幾何關系式、利用力平衡條件得到兩個平衡方程、依據切點和索鞍圓心的距離關系得到兩個公式,再加上一些工程經驗性的約束條件后通過數值分析或Mathcad等軟件進行計算. 多位學者對此均有論述[8-10],但該方法計算效率并不高. 另外,有學者采用通用有限元軟件建立主纜與索鞍的接觸關系進行數值仿真分析[11-15],也有部分學者考慮了主纜與鞍座之間的摩擦來提高計算精度[16-19]. 上述方法在精度方面雖然能滿足工程需求,但計算比解析法更困難,效率更低.
目前的索鞍曲線修正計算分析理論推導繁瑣、難以操作、效率較低. 所以在實際工程中,索鞍曲線修正的影響基本只在懸索橋的成橋狀態下考慮. 但是為了滿足人們建造更大跨徑懸索橋的需求,僅在成橋狀態進行索鞍修正無法滿足施工對線形的控制要求[20]. 因此,推導一種便于計算、效率更高的索鞍曲線修正算法是十分有必要的. 不僅可以實時進行索鞍曲線修正,還可為大跨徑懸索橋的施工控制提供一定的參考. 為此,本文從幾何關系角度對索鞍曲線修正算法進行推導,并借助牛頓-拉斐森迭代法,提供了一種可行的數值求解方法.
現代懸索橋主纜一般由鋼絲集束形成,且其線形平順,在轉折處曲率半徑很大,因此相比其軸向抗拉剛度,抗彎剛度很小,具有很強的柔性索結構特征. 本文在推導懸索橋主纜的算法時,對主纜作3項前提假設:1)主纜是理想柔性的索,只承受軸向拉力,不承受軸向壓力和彎矩;2)在設計懸索橋時通常會給主纜拉應力設置較大的安全系數,主纜在施工和使用過程中,材料均不會達到塑性階段,因此其符合胡克定律,是線彈性的;3)主纜橫截面積在外荷載作用下變化量微小,不予考慮.
將一段只受自重的纜索分割出來,如圖1所示,H和V分別表示主纜的水平分力和豎向分力,q為纜索沿弧長均布的自重荷載,l為主纜在x軸的投影長度,c為主纜在y軸的投影長度.

圖1 懸鏈線纜索示意圖
顯然,可以得到主纜在自重作用下的曲線形狀為懸鏈線[8],具體表達式為
(1)
其中,各參數表達式為
(2)
對于原點處H、V已知,則有
(3)
索鞍處的主纜線形如圖2所示,x為水平向,y為豎向. 整體分析時為方便計算,會假設一個虛擬的主纜交點,即圖中的理論頂點;本文將局部坐標系原點建立在理論頂點上,其坐標為(0, 0),左側主纜與索鞍的切點標為切點1,其局部坐標為(x1,y1),右側主纜與索鞍的切點標為切點2,其局部坐標為(x2,y2),索鞍與主纜接觸面為圓弧,索鞍圓心局部坐標為(x3,y3),半徑為R. 另外,x1、x2與y1、y2的關系滿足懸鏈線式(1),為推導過程簡潔,用y=f(x)來表示.

圖2 索鞍與主纜關系示意圖
根據左側主纜與索鞍相切可以得到,切點1與索鞍圓心的連線斜率k1表達式為
(4)
根據右側主纜與索鞍相切可以得到,切點2與索鞍圓心的連線斜率k2表達式為
(5)
其中k為理論切點位置索(索鞍)的切線斜率.
切點1與切點2同在以索鞍圓心為圓心,半徑為R的圓上,根據幾何關系,可得
(6)
左側主纜方程為
(7)
右側主纜方程為
(8)
其中,式(7)、(8)中的參數表達式為
(9)
另外,可以得到
(10)
將f1′(x1)、f2′(x2)、y1、y2代入式(6)后得到求解主纜與索鞍切點的二元非線性方程組為
(11)
通過牛頓-拉斐森迭代法[21]可對本文所述二元非線性方程組(11)進行數值求解,令
(12)
再令
其中Ai為雅克比矩陣,矩陣各元素表達式為
(13)
根據牛頓-拉斐森迭代,則有
Xi+1=Xi-Ai-1Gi.
(14)
取切點1、切點2的初始x坐標代入牛頓-拉斐森迭代,得到x1、x2的數值解,y1、y2代入前文提到的懸鏈線方程即可得到. 值得注意的是,傳統的八元非線性方程組的初始值選取原則較為繁瑣,而本文算法對初始值的選取基本沒有要求.
根據前述理論推導,索鞍曲線修正算法由傳統的八元非線性方程組(傳統方法)簡化為二元非線性方程組(本文方法),可較為方便地將其嵌入有二次開發功能的有限元分析軟件,如:TDV、Ansys、Abaqus等,也可通過MATLAB等編程語言方便地進行求解.
本文分別選取主索鞍與散索鞍兩個算例,對其正確性與計算效率進行驗證. 需要說明的是,表中左切點坐標(x1,y1)、右切點坐標(x2,y2)、索鞍圓心坐標(x3,y3)均為整體坐標系坐標.
某懸索橋的其中一個主索鞍理論頂點坐標為(230.000 m, 131.425 m),主纜面積A1=A2=0.408 973 m2,主纜自重集度為q1=q2=33 kN/m,彈性模量E=198 GPa. 索鞍半徑R=6 m,索鞍左右兩側的主纜索力水平分量H1=H2=189 500 kN,V1=90 622.7 kN,V2=73 504.1 kN.
通過編程,對傳統方法與本文方法(第2節的理論推導)進行對比,初始取值見表1,計算過程見表2,計算結果見表3.

表1 主索鞍曲線修正初始取值

表2 主索鞍曲線修正計算過程

表3 主索鞍位置及主纜曲線計算結果
某懸索橋的其中一個散索鞍理論頂點坐標為(0 m, 54 m),主纜面積A1=A2=0.408 973 m2,主纜自重q1=q2=33 kN/m,彈性模量E=198 GPa. 索鞍半徑R=6 m,索鞍左右兩側的主纜索力水平分量H1=H2=189 500 kN,V1=137 557 kN,V2=-41 804.3 kN.
通過編程,對傳統方法與本文方法(第2節的理論推導)進行對比,初始取值見表4,計算過程見表5,結果見表6.

表4 散索鞍曲線修正初始取值

表5 散索鞍曲線修正計算過程
從上述兩個算例的計算結果和計算效率對比可知,本文提出的二元非線性方程組求解方法與目前普遍采用的八元非線性方程組求解方法計算結果在0.01 mm的數量級上可以保持一致性,對于工程應用可以認為兩種方法的計算結果是沒有差別的. 計算效率方面,本文提出的方法迭代次數減少了50%,計算耗時僅為目前普遍采用方法的10%. 更重要的是,在設置初始參數時,本文提出的方法要簡便很多,只需要設置x1和x2兩個參數即可,且沒有特殊要求. 一般情況下,為便于編寫程序可以直接設置初始x1=x3-R/2,x2=x3+R/2,即可很快迭代收斂. 然而,參考文獻[8]中的推導過程表明,目前普遍采用方法,需要設置8個初始參數,這8個參數中只有x1、x2較易設置,其余6個參數均需有一定的經驗或者經過數次試算方可獲得較好的計算結果,這給該算法實現程序自動進行索鞍曲線修正計算帶來很大困難. 由此說明,本文所提出方法更為方便,大大提高了計算效率,且精度可滿足要求.

表6 散索鞍位置及主纜曲線計算結果
1)根據主纜與索鞍的幾何關系,推導了索鞍處主纜曲線的修正計算方法. 相比于傳統算法,減少了6個方程與6個初始輸入參數,表達形式更明確,物理含義更清晰.
2)利用牛頓-拉菲森迭代法,對所得二元非線性方程組進行求解. 其中,關于輸入參數的初始值設置沒有嚴格要求,具有很強的可操作性.
3)針對常見的主索鞍與散索鞍結構,相比于傳統算法,迭代次數減少50%,計算時間不足10%,大大提高了計算效率,且精度可滿足要求.