懸鏈線拱橋主拱圈變截面高度的數值計算方法

由變截面拱圈的截面變化方式引入，依據微積分弧微分公式計算曲線長度的原理，推導出懸鏈線拱軸線長度的計算公式無法用初等函數表示，由此引入數值積分的方法求解軸線長度，進而深入了解自適應辛普森積分法，采用Excel中的Visual Basic語言自編程序驗證理論充分可行。

懸鏈線； 變截面； 合理拱軸線； 自適應辛普森積分； Visual Basic

U448.22+5 A

［定稿日期］2022-03-23

［作者簡介］白高（1995—），男，本科，助理工程師，從事市政、公路橋梁設計工作。

拱橋具有獨特的山區地形適應能力、強大的跨越能力以及可靠的使用性能，已經成為我國公路、鐵路廣泛使用的一種橋型。拱橋主拱圈作為主要受力構件，起著傳遞和承受荷載的重要作用，拱橋主拱圈線形又對主拱圈的受力合理性、建設經濟性起著至關重要的作用，拱軸線常用類型有圓弧線、拋物線、懸鏈線。圓弧線適用于承受徑向均布荷載的結構，圓弧線拱在水工結構中多見，在橋梁結構中為施工方便，圓弧線多用于跨徑20 m以下的小跨徑拱橋；拋物線形拱軸線適用于拱圈承受均布荷載，理論上拋物線適用空腹式拱；懸鏈線形拱軸線適合承受拱頂至拱腳連續變化的荷載，是實腹式拱橋的合理拱軸線。雖然空腹式拱橋采用懸鏈線形拱軸線在恒載作用下壓力線與拱軸線存在偏離，但計算表明該偏離通常對拱圈控制截面受力是有利的，并且，僅需使拱軸線與恒載壓力線在拱頂、1/4跨徑點和拱腳五點重合即可實現拱內彎矩最小化，因此，懸鏈線是目前大、中跨徑拱橋采用最普遍的拱軸線形［1］。

對大跨徑拱橋，拱圈截面慣性矩采用拱頂至拱腳變化的方式，拱圈慣性矩變化通常是為了適應拱內內力變化，慣性矩的變化通常表現為拱圈高度的變化（圖1），拱圈高度的合理變化有利于充分發揮拱橋的每個截面的材料強度［2］。對于無鉸拱，通常根據李特公式（Ritter 公式［1］）來改變截面慣性矩，在設計工作中，直接應用李特公式逐個計算截面高度不便操作，原因是工程中拱圈截面形狀復雜、不同的材料在截面中分布的規律性與截面高度之間的關系無法通過簡單的數學模型描述（板拱除外），因此對復雜截面李特公式不能直接快速得到截面高度；為簡單實用起見，設計大跨徑拱橋時工程師常采用沿拱軸線長度（弧線）方向，拱圈截面高度采用拋物線漸變的方式，即拱圈截面高度是拱軸線弧長的拋物線函數，這就是要解決的問題：如何計算變截面中各變化點的截面高度，得到拱圈設計坐標。

將自適應辛普森算法的原理通過Excel中VBA語言編制成程序，將程序計算結果與實際工程中的設計值對比驗證，論證方法的可行性和正確性。

1 理論計算原理

一元函數積分計算原理是牛頓-萊布尼茲原理，即找到被積函數的原函數，將積分上下限代入原函數作差，這是微積分的基本原理。牛頓法計算定積分的前提是被積函數的原函數能用初等函數表示，下面討論懸鏈線方程（圖2）。

由《橋梁工程》教材［1］，懸鏈線的方程為式（1）。

y=fm-1（coshkξ-1）（1）

其中：

k=cosh-1m=ln（m+m2+1）（2）

其中：

ξ=xl1=xl/2=2xl（3）

由《高等數學》［3］教材得懸鏈線的弧微分為：

ds=1+y′2dx（4）

將式（1）、式（3）代入式（4），借助數學專業程序Wolfram Mathematica 12.1計算得到式（4），結果為式（5）。

ds=1+4f2k2Sin（2kx/l）2l2（m-1）2dx（5）

于是弧長計算的原函數用Wolfram Mathematica 12.1程序計算為：

∫1+4f2k2Sin2kxl2l2（m-1）2dx=-lEllipticE2ikxl，4f2k2l2（m-1）22k+C（6）

式（6）中i為虛數符號，EllipticEφ，m為第二類橢圓積分，第二類橢圓積分是無法通過初等函數給出解析式的，牛頓法解決不了這個問題。

考慮用數值積分來解決問題，一元函數積分幾何意義是函數y=f（x）的圖形，在積分上下限區間上與x軸所圍成的曲邊梯形的面積（圖3）。數值積分的目的還是求這個曲邊梯形的面積，并且需要滿足一定的精度要求。

［a，b］是一個較小的區間，見圖4，由于［a，b］區間很小，［a，b］中曲邊梯形的面積可以近似采用所示的梯形陰影面積來表示，即∫baf（x）dx≈b-a2［f（a）+f（b）］=（b-a）12f（a）+12f（b），這就是數值計算梯形公式的原理。

上式中a、b稱積分點，f（a）、f（b）是積分點對應的函數值，2個“12”是積分點處函數值所占的權重，顯然梯形公式的核心是求2個積分點的函數值的加權平均值；加權平均值越接近［a，b］內f（x）的真實平均值，計算精度也就越高。由此想到可以在［a，b］間等距插入有限個節點x0，x1…xn，然后采用合適的權值A0，A1，…，An，然后令（b-a）∑nk=0Akf（xk）更加精確地表達陰影部分的面積，觀察上式Ak f（xk）可利用n次拉格朗日型插值多項式構造［4］，得到牛頓-柯特斯公式見式（7）。

I（f）=（b-a）∑nk=0C（n）kf（xk）（7）

C（n）k為柯特斯系數，如n=1時，C（1）0=12，C（1）1=12，這就是梯形公式，當n=2時，C（2）0=16，C（2）1=46，C（2）2=16，這時得到的公式稱為辛普森公式［5］見式（8）。

I（f）=（b-a）6f（a）+4fa+b2+f（b）（8）

牛頓-柯特斯公式計算積分的精度不但與區間b-a的大小有關，與區間b-a內所劃分的節點數n也有關（n值越大越精確，但n≥8時，牛頓-柯特斯公式計算誤差不穩定）［6］，現今對n=2時的辛普森公式研究已經十分成熟，便于編程應用。

設拱軸線的曲線長度積分區間為［x1，x2］，（x1，x2∈0，l/2，x1

如何控制積分區間的劃分是一個重要的問題，令積分區間［a，b］中點為c=a+b2，用S（a，c）表示采用辛普森公式計算區間［a，c］的積分值，S（a，b）、S（c，b）同理表示用辛普森公式計算的不同區間的積分值，由數值分析的知識可以得到S（a，b）的絕對誤差近似計算為［4］式（9）。

eps=|S（a，b）-（S（a，c）+S（c，b））|15（9）

將eps作為程序循環計算的精度判斷條件，可自行控制計算精度，由此理論完備。

2 程序實現

用Excel中的Visual Basic語言編制程序，曲線長度計算的核心代碼如圖5所示，代碼能實現給定積分區間，返回該區間的懸鏈線弧長，代碼中的yfsenal變量即上述eps，至此引言問題得解。

3 數據驗證

筆者找到了廣西南寧市某大橋拱肋線形設計圖，該橋拱肋高度按照2.5次拋物線變化，由作者自編程序計算得到的拱圈坐標，與設計文件所給坐標對比完全吻合，選取有代表性的7個點對比（表1）。

4 結論

通過數學推導，得出懸鏈線的曲線長度是不可通過初等函數表達的結論，利用數值分析中求一元函數積分的辛普森公式，在辛普森公式中又引入了復化和自適應的算法，結合Excel中的Visual Basic語言編制成實用的計算程序，程序通過與實橋的設計文件對比驗證，論證了理論的正確性，并且達到了足夠的精度。橋梁設計是一門綜合了數學、力學、計算機科學等多自然學科的綜合技術，在橋梁設計工作中，跨學科的綜合應用有很多，工程師不可能熟練掌握各個學科的各種理論，因此，工程師們在遇到復雜問題時，需要冷靜分析，大膽尋求解決問題的方法，小心求證方法的可行性和正確性，準確有效地解決橋梁設計中的疑難問題，保證橋梁結構安全可靠的使用性能，是廣大橋梁工程師的使命。

參考文獻

［1］ 顧安邦，向中富.橋梁工程（下冊）［M］.3版.北京：人民交通出版社股份有限公司，2017（7）：83，135-137.

［2］ 宋功譚.700m級鋼管混凝土拱橋主拱圈構造研究［D］.重慶：重慶交通大學，2021.

［3］ 同濟大學數學系.高等數學（上冊）［M］.6版.北京：高等教育出版社，2007.4：169-170.

［4］ （美）薩奧爾（Sauer， T.）著；裴玉茹，馬賡宇譯.數值分析（原書第2版）［M］.北京：機械工業出版社，2014.10：124-126，240.

［5］ 劉佳霖，臧婷.兩種定積分數值計算方法的比較［J］.電子技術，2021，50（7）：252-253.

［6］ 林承初.牛頓—柯特斯求積公式穩定性分析［J］.社科與經濟信息，2001（9）：172-173.

［7］ 唐珺，楊道杰.自適應積分的遞歸實現研究［J］.中國高新技術企業，2016（27）：13-16.