經(jīng)典懸鏈線理論精確解與近似解的非線性數(shù)值計算

郭小剛， 金 星， 周 濤， 宋曉東， 鄧旭輝

（1．湘潭大學(xué) 土木工程與力學(xué)學(xué)院，湘潭411105；2．長沙礦冶研究院、深海礦產(chǎn)資源開發(fā)利用技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，長沙410012）

1 引 言

張拉結(jié)構(gòu)的應(yīng)用日益廣泛，已經(jīng)覆蓋如大跨度空間結(jié)構(gòu)的索網(wǎng)、海洋工程中的軟管和柔索［1－4］等結(jié)構(gòu)。目前國內(nèi)外對張拉結(jié)構(gòu)的分析手段亦日趨復(fù)雜，采用兩節(jié)點(diǎn)直線單元［5－7］、三節(jié)點(diǎn)拋物線單元甚至五節(jié)點(diǎn)等參數(shù)索單元［8］的非線性有限元分析，特別在動力時域分析時沿用Wilson法Newmark法內(nèi)蘊(yùn) Newton－Raphson非線性迭代［9－12］，使得工程應(yīng)用復(fù)雜繁瑣，且處理的大多是小垂度的問題。

懸鏈線問題的初步數(shù)學(xué)解答已有300多年的歷史，關(guān)于采用懸鏈線解答求解工程問題的文獻(xiàn)層出不窮［1，12－24］。 沈 世 釗 等［1］介 紹 了 單 索 均 布 力 作用下的懸鏈線理論以及懸鏈線解的拋物線近似，但留下水平張力和左邊界端點(diǎn)斜率兩個未知數(shù)懸而未決。鄭麗鳳等［13］選取水平張力和左端邊界垂直力這兩個參數(shù)作為獨(dú)立未知量，建立了二元非線性牛頓迭代法方程組來求解懸鏈線的解，但左端邊界垂直力未必適宜作為獨(dú)立的參數(shù)。楊孟剛等［14］采用UL列式增量法，提出了兩節(jié)點(diǎn)懸鏈線索元非線性切線剛度矩陣，但未明確如何計算最重要的水平初始張力。周緒紅等［15］采用的懸索理論實(shí)質(zhì)上是懸鏈線理論的拋物線逼近，是工程上可行的方法。靳明君等［16］推導(dǎo)了懸鏈線索長度與水平張力的函數(shù)關(guān)系，但將水平張力作為已知量似有因果倒置之嫌。張卓杰等［17］在懸鏈線的求解中選取左邊界的斜率即廣義傾角α的映射作為初始試探值，反復(fù)迭算從而求得懸鏈線的解，計算過程過于迂回。馮海暴［18］采用懸鏈線理論對35m作業(yè)水深鋪排管線進(jìn)行了力學(xué)分析，但采用的方程要求廣義傾角α＝0，這一前提很難適用于流體動力的變化。劉婷等［19］采用懸鏈線二維解作為跨接管安裝的柔性段三維分析的初始構(gòu)型，但在理論上沒有詳述如何求解二維懸鏈線解的水平張力這一關(guān)鍵問題。呂玉蘭等［20］在求解懸鏈線方程時，以水平張力和懸鏈線最低點(diǎn)的x坐標(biāo)為未知數(shù)，算法過于復(fù)雜，且懸鏈線最低點(diǎn)能否選為獨(dú)立變量也有待商榷。龐曉旭等［21］探討了水平懸鏈線在小垂度下彈性模量對力學(xué)特性的影響。周陽等［22］以給定上部頂張力和懸掛角作為邊界條件，研究了下部海床土非線性剛度對鋼懸鏈線立管的影響。吳衛(wèi)國等［23］采用懸鏈線靜態(tài)解研究了船只補(bǔ)給懸鏈線系統(tǒng)的動力特性。

總之對于懸鏈線的數(shù)學(xué)解析解即雙曲余弦函數(shù)，如何簡潔地得到其中最重要的未知參數(shù)即水平張力，依然困擾著工程界，而懸鏈線的數(shù)值求解在張拉結(jié)構(gòu)中具有舉足輕重的作用。小垂度下懸鏈線的拋物線近似解目前在工程上得到了可行的解決途徑［1，15，21，23］，但鮮有文獻(xiàn)探討懸鏈線在全局范圍任意垂度下的近似解，本文在此方面作出了有益的嘗試。

2 懸鏈線超越方程與數(shù)值解

若結(jié)構(gòu)任意截面上的彎矩可假設(shè)為0，則結(jié)構(gòu)稱為完全柔性，若結(jié)構(gòu)變形前后可認(rèn)為長度保持不變，則結(jié)構(gòu)具有不可拉伸性。經(jīng)典懸鏈線理論采用了完全柔性和不可拉伸性假設(shè)，因其采用完全柔性假設(shè)，相應(yīng)的懸掛結(jié)構(gòu)稱為柔索結(jié)構(gòu)或柔性結(jié)構(gòu)，但因?yàn)椴豢衫煨约僭O(shè)，柔索結(jié)構(gòu)沿縱向具有內(nèi)在的剛性特質(zhì)。

對于僅有垂直力作用的結(jié)構(gòu)，其任一截面水平方向的合內(nèi)力處處相等。柔索在僅有均布垂直載荷q作用下經(jīng)典的懸鏈線解［1，4，24］為

式中ch為雙曲余弦函數(shù)；xa和ya分別為懸鏈線左端點(diǎn)a處的x方向坐標(biāo)和y方向坐標(biāo)；獨(dú)立的待求量h為沿懸鏈線不變的水平張力，h＞0；待求量α為左邊界a處斜率的廣義傾角，shα為邊界a處的斜率。線均布力q的單位是N／m，水平張力h的單位是N。

懸鏈線右端點(diǎn)b處的x方向坐標(biāo)為xb，y方向坐標(biāo)為yb。l＝xb－xa表示懸鏈線左右兩端的相對水平距離，簡稱水平距離；c＝y(tǒng)b－ya表示懸鏈線左右兩端的相對垂直距離，簡稱垂直距離。

懸鏈線在水平距離l變動時存在一個不能逾越的極限水平距離d，簡稱極限距離，

極限距離d對應(yīng)的線是一條直線，這條直線是懸鏈線的水平距離l無限逼近d時懸鏈線的漸近線，簡稱極限漸近線。s為懸鏈線的索長。

懸鏈線解（1）中高度綁定的表達(dá)式h／q具有深刻的內(nèi)在涵義，表示懸鏈線曲率半徑中的最小值。因此，引入曲率半徑r或曲率直徑D為

r＝h／q，D ＝2r＝2h／q （3）

曲率半徑r和懸鏈線形態(tài)本身僅取決于懸鏈線段的索長s、水平距離l和垂直距離c，而與作用其上的均布載荷q無關(guān)。定義無量綱參數(shù)b和n為

無量綱參數(shù)b和n為互逆關(guān)系，bn＝1。無量綱參數(shù)n＝l／d，為位于區(qū)間（0，1）的正小數(shù)，可稱為比極距離。無量綱參數(shù)n越大，曲率半徑r越大，n越接近于1，曲率半徑r趨于無窮大，懸鏈線越逼近懸鏈線的極限漸近線即越逼近直線；無量綱參數(shù)n越小，曲率半徑r越小，水平距離l越接近于0，曲率半徑r趨于0。

定義無量綱參數(shù)β和z為

無量綱參數(shù)β和z為互逆關(guān)系，βz＝1。無量綱參數(shù)z正比于單位均布載荷單位水平距離對應(yīng)的水平張力，可稱為相對水平張力，表示懸鏈線的張馳程度。

不失一般性，將坐標(biāo)系放置于左端點(diǎn)a處，利用式（1，3），從起點(diǎn)到任意位置x處的索長度為

定義廣義傾角

θ＝α＋β （7）

根據(jù)兩點(diǎn)邊值約束條件，即c＝y(tǒng)（l）＝y(tǒng)b－ya，與不可拉伸弧長約束條件，即s＝s（l），且利用定義式（3，5）可得到索長s和垂直距離c與廣義傾角θ和β的關(guān)系為

由式（8，9）可知廣義傾角θ的計算式為

將式（5，10）代入式（7），得到廣義傾角α的計算式為

從廣義傾角α的計算式（11）可知，α不是獨(dú)立的未知量，獨(dú)立的未知量只有一個，即水平張力h。

利用式（2～5，8，9）運(yùn)算得到一組等價的隱含水平張力的超越方程為

一般地，懸鏈線水平張力h取決于三個參數(shù)，即水平距離l、懸鏈線段的索長s和垂直距離c。對方程（12）的無量綱參數(shù)b或n與無量綱參數(shù)β或z進(jìn)行仔細(xì)分析可知，懸鏈線相對水平張力z＝2h／q l僅依賴于水平距離l與水平極限距離d之比n＝l／d，這是懸鏈線兩點(diǎn)邊值約束問題相當(dāng)明確且可靠的結(jié)論。

根據(jù)雙曲函數(shù)的性質(zhì)［25］和式（12），關(guān)于廣義傾角β存在關(guān)系為

式（13）可視為懸鏈線段兩點(diǎn)邊界約束的靜力平衡關(guān)系，其中

根據(jù)等價的超越方程（12）可得到無量綱參數(shù)β與無量綱參數(shù)b之間相互關(guān)系的離散數(shù)值解，或可得到相對水平張力z與無量綱參數(shù)n之間的離散數(shù)值解，兩者都是超越方程的精確解或真實(shí)解。

在用C語言進(jìn)行計算機(jī)數(shù)據(jù)處理時，8字節(jié)的double型實(shí)數(shù)不超過10的±309次方冪，因此計算中將接近0的double型的比極距離n＝l／d的下限取為10－308。計算程序定義因變量0＜n＝l／d＜1的計算域?yàn)椋?0－308，0．999999］，該計算域是真小數(shù)區(qū)間（0，1）內(nèi)C語言仿真計算工程上可行的全域。

圖1將log10（l／d）作為橫坐標(biāo)，顯示相對水平張力精確解的變化情況。n＝l／d的計算域截取為［10－308，0．999999］，log10（l／d）的值域?yàn)椋郏?08，0）。因變量n＝l／d的計算域［10－308，10－7］內(nèi)，相對水平張力在接近于0的附近緩慢增加。其中z（10－308）＝0．001408，z（10－200）＝0．002141，z（10－150）＝0．00284，z（10－100）＝0．00423，z（10－7）＝0．050513。相對水平張力精確解隨因變量n＝l／d的變化情況如圖1所示。

3 水平張力上下界的確定

對于β＞0，數(shù)學(xué)上存在thβ＜β，利用懸鏈線段靜力平衡關(guān)系式（13）有

從式（15）可確定相對水平張力上界z0：

將雙曲正弦函數(shù)泰勒級數(shù)展開，保留方程（12）的二次冪得到近似根β1為

保留二次冪得到的近似根β1大于真實(shí)的根，即β1＞β，考慮到關(guān)系β1z1＝zβ＝1，因此z1＜z；考慮到關(guān)系nb＝1，可確定相對水平張力的下界之一為z1，如式（18）所示。相對水平張力z介于式（18）中上下界的變動范圍：

分析發(fā)現(xiàn)，式（16）中z0對應(yīng)的水平張力上界h0過高地偏離了真實(shí)的水平張力，引入相對水平張力z2：

計算表明實(shí)際相對水平張力z介于式（20）中相對水平張力下界z1和上界z2之間的變動范圍：

4 懸鏈線超越方程的近似解

以無量綱參數(shù)β＝2l／qh和相對水平張力z＝D／l表示方程sh（l／D）＝d／D 或sh（β）＝bβ的精確解，βz＝1。若能找到無量綱函數(shù)zk（n），在0＜n1≤n≤n2＜1域內(nèi)滿足：

則zk（n）就是方程sh（l／D）＝d／D 或sh（β）＝bβ在n1≤n≤n2域內(nèi)的無量綱近似解或分段近似解。若近似解與精確解的比值bzk（n）介于0．8≤bzk（n）≤1．2，那么可以概略地說近似解的誤差不超過20%，若0．89≤bzk（n）≤1．11，則近似解的誤差大約不超過11%。

相對水平張力z1和z2是水平距離l趨于極限距離d的近似解。取下界z1和上界z2，兩者的平均值zm為

計算表明相對水平張力z介于相對水平張力下界z1和上界zm之間的變動范圍：

即上界zm是比原上界z2更小的上界，在因變量n的計算域［0．4，0．999999］內(nèi)，zm與精確解z 的誤差不超過10%。n≥0．62時，zm與精確解z的誤差不超過5%；n≥0．68時，zm與精確解z的誤差不超過4%；n≥0．75時，zm與精確解z的誤差不超過3%；n≥0．91時，zm與精確解z的誤差不超過1%；n≥0．991時，zm與精確解z的誤差不超過0．1%。

水平距離l接近極限距離d時的近似解z1，zm，z2和精確解z在n＝l／d介于［10－308，0．95］內(nèi)的變化趨勢，如圖2所示。當(dāng)n＝0．95時，z1＝1．779513，z＝1．793394，zm＝1．802627，z2＝1．825742。

圖2 相對水平張力的近似解和精確解Fig．2 Exact and approximate solutions of relative horizontal tension

圖3 近似解zm和精確解z之比Fig．3 Ratio of approximate solution zmto the exact solution z

圖3 顯示了近似解zm與精確解z之比zm／z＝βzm在n＝l／d 介于［0．001，0．99］內(nèi)的變化趨勢。n＝10－308時，βzm＝145．025；n＝10－7時，βzm＝4．04；n＝0．001時，βzm＝2．084；n＝0．99時，βzm＝1．001。若工程誤差許可不超過5%，則當(dāng)水平距離l≥0．62d時，水平張力可以由式（4，5，22）表示為

當(dāng)水平距離l接近極限距離d時，水平張力變得充分大。式（24）清晰地表明了水平距離l趨于極限距離d時，水平張力趨于無窮的變化趨勢，公式簡潔可行。從式（24）或圖2可以看出，保持d不變，l增大，水平張力隨之增大。保持l不變，d增大（s增大，c減?。?，水平張力隨之減小。式（24）在小垂度情形與文獻(xiàn)［1］拋物線近似解趨于一致。

探討l趨于0時的近似解。根據(jù)懸鏈線段平衡方程sh（ql／2h）＝qd／2h可知，當(dāng)水平距離l趨于0時，精確解的水平張力較l緩慢地趨于0。

鑒于C語言8字節(jié)浮點(diǎn)數(shù)的下界值大約為10－308量級，下面計算上可行地找出無量綱參數(shù)或比極距離n＝l／d在計算域［10－308，10－7］內(nèi)的近似解，該近似解視為l趨于0時的近似解。

通過多次試探計算，尋找到比極距離n＝l／d在計算域［10－308，10－7］內(nèi)的近似解zl為

bzl＝zl／z是近似解zl與精確解z之比，表示兩者的逼近程度。

水平距離l接近0時的近似解zl與精確解z之 比bzl，在n＝l／d介于［10－20，0．009］內(nèi)的變化趨勢如圖4所示。其中n＝10－308時，b zl＝1．000394；n＝10－7時，b zl＝1．1565；n＝0．009時，bzl＝1．296896。

圖4 n接近0的近似解zl和精確解之比Fig．4 Ratio of approximate solution zlto the exact one with nnear zero

從計算結(jié)果可知，近似解zl在［10－308，10－7］內(nèi)的誤差不超過16%，這是一個初始的近似解。

調(diào)整參數(shù)或選擇另外的函數(shù)等可以找到分段區(qū)域內(nèi)誤差更小、逼近程度更高的近似解。如函數(shù)z4＝［1．0－0．5·sqrt（n）］／［1．0－log（n）］較好地模擬了n接近10－308時精確解z的變化趨勢，將z4適當(dāng)減小約10%，得到近似解zi為

近似解zi與精確解z之比zi／z＝bzi在n＝l／d介于［10－308，0．25］內(nèi)的變化趨勢如圖5所示。因變量n在區(qū)域［10－308，0．25］內(nèi)，zi與精確解z 的比值bzi介于0．89＜bzi（n）＜1．10154，在此范圍內(nèi)zi與精確解的誤差不超過11%。其中n＝10－308時，bzi＝0．89；n＝10－7時，bzi＝1．029；n＝0．001～0．009時，bzi≈1．1；n＝0．25時，bzi＝0．91。

圖6將log10l／d作為橫坐標(biāo)，顯示近似解和精確解之比zi／z的變化情況，n＝l／d的計算域截取為［10－100，0．999999］，log10l／d 的值域?yàn)椋郏?00，0）。可以看出，無量綱水平距離n介于該范圍內(nèi)時，bzi＝zi／z的細(xì)微變化情況。

圖5 近似解zi和精確解之比Fig．5 Ratio of approximate solution zito the exact one

圖6 近似解和精確解之比zi／z的對數(shù)顯示Fig．6 Logarithmic displaying the ratio zi／zof approximate solution to the exact one

從圖3可以看出，當(dāng)n＝10－308時，zm／z＝145．025278，即當(dāng)n接近于0時，上界zm過高估算了精確解。因此將zm適當(dāng)減小，找到逼近精確解的近似解zn

關(guān)系式lgn＝lgelnn＝Mlnn表示常用對數(shù)與自然對數(shù)之間是線性映射關(guān)系。其中模數(shù)［25］M為M＝lge＝1／ln10≈0．434294481903…。

計算結(jié)果表明，zn與精確解z的比值bzn存在幾個 關(guān) 鍵 點(diǎn) 函 數(shù) 對 應(yīng) 關(guān) 系，bzn（10－308）＝1．076153，bzn（10－300）＝1．085261，bzn（10－100）＝1．086182，bzn（10－7）＝1．000548，bzn（0．04）＝0．87771，bzn（0．999990）＝0．999999。近似解zn與精確近似解zn與精確解z的比值b zn在n靠近0時略微大于1，然后單調(diào)增加，在n＝10－100附近時接近最大值1．086182，隨后單調(diào)減少，在n＝0．04附近，bzn接近最小值0．87771，之后bzn隨n 單調(diào)增加，n趨于1時bzn趨于1。如圖7和圖8所示。

圖7 近似解zn和精確解之比Fig．7 Ratio of approximate solution znto the exact one

圖8 近似解和精確解之比zn／z的對數(shù)顯示Fig．8 Logarithmic displaying the ratio zn／zof approximate solution to the exact one

圖8 將log10l／d作為橫坐標(biāo)，圖8（a）顯示近似解和精確解之比zn／z的變化。n＝l／d的計算域截取為［10－308，10－25］，log10l／d 的值域?yàn)椋郏?08，－25］，可以看出近似解和精確解之比bzn＝zn／z的細(xì)微變化情況。圖8（b）中無量綱水平距離n計算域截取為［10－100，0．999999］，log10l／d的值域?yàn)椋郏?00，0），可以看出bzn＝zn／z的大致變化趨勢。

從計算結(jié)果分析可知，比值bzn在n的全計算域內(nèi)介于0．877＜bzn（n）＜1．0862；因此由公式（27）計算的相對水平張力zn是因變量n＝l／d全域范圍內(nèi)誤差不超過13%的近似解。

5 結(jié) 語

首先，在兩點(diǎn)位移約束方程基礎(chǔ)上，自然地補(bǔ)充引進(jìn)弧長不可拉伸約束方程，從而簡單解決了求解經(jīng)典懸鏈線數(shù)學(xué)解的兩個未知參數(shù)即水平張力h和表示左邊界處斜率的廣義傾角α這一非線性數(shù)學(xué)問題。其次，推導(dǎo)出求解懸鏈線水平張力的超越方程。先指明h／q為懸鏈線的最小曲率半徑，后引進(jìn)無量綱參數(shù)z＝2h／ql和n＝l／d求解超越方程中的水平張力，使得水平張力形式上具有最簡單的參數(shù)依賴關(guān)系。探討了廣義傾角β，θ和α的內(nèi)在含義以及與幾何參數(shù)的相互關(guān)系，創(chuàng)新地給出θ和α的計算公式。從廣義傾角α的計算公式可知，廣義傾角α不再是獨(dú)立的原始未知量，獨(dú)立的原始未知量只有一個即水平張力h。

隨之，通過對平衡方程的泰勒展開，經(jīng)過理論分析和反復(fù)計算得到水平張力的上下包絡(luò)界，為確定水平張力的精確數(shù)值解提供算法上的理論指導(dǎo)。

最后，提出了水平距離l分別趨于0與趨于極限距離d的近似解zl，zi和zm等以及在真小數(shù)計算范圍［10－308，0．999999］內(nèi)全局可行的近似解zn，分析計算了這些近似解關(guān)于精確解的誤差程度，近似解以函數(shù)的形式粗顯地刻畫了水平張力隨n＝l／d的因果變化趨勢，對于水平張力工程的初步快速估算具有重要的意義，也可對數(shù)值計算的發(fā)展扮演拋磚引玉的作用。這些近似解可作為初值加快精確數(shù)值解的非線性迭代計算。

綜上所述，本文關(guān)于懸鏈線提出的相關(guān)公式和近似解，在工程應(yīng)用上具有實(shí)際的價值。