郭小剛, 金 星, 周 濤, 宋曉東, 鄧旭輝
(1.湘潭大學(xué) 土木工程與力學(xué)學(xué)院,湘潭411105;2.長(zhǎng)沙礦冶研究院、深海礦產(chǎn)資源開發(fā)利用技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)沙410012)
張拉結(jié)構(gòu)的應(yīng)用日益廣泛,已經(jīng)覆蓋如大跨度空間結(jié)構(gòu)的索網(wǎng)、海洋工程中的軟管和柔索[1-4]等結(jié)構(gòu)。目前國(guó)內(nèi)外對(duì)張拉結(jié)構(gòu)的分析手段亦日趨復(fù)雜,采用兩節(jié)點(diǎn)直線單元[5-7]、三節(jié)點(diǎn)拋物線單元甚至五節(jié)點(diǎn)等參數(shù)索單元[8]的非線性有限元分析,特別在動(dòng)力時(shí)域分析時(shí)沿用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)斜率兩個(gè)未知數(shù)懸而未決。鄭麗鳳等[13]選取水平張力和左端邊界垂直力這兩個(gè)參數(shù)作為獨(dú)立未知量,建立了二元非線性牛頓迭代法方程組來求解懸鏈線的解,但左端邊界垂直力未必適宜作為獨(dú)立的參數(shù)。楊孟剛等[14]采用UL列式增量法,提出了兩節(jié)點(diǎn)懸鏈線索元非線性切線剛度矩陣,但未明確如何計(jì)算最重要的水平初始張力。周緒紅等[15]采用的懸索理論實(shí)質(zhì)上是懸鏈線理論的拋物線逼近,是工程上可行的方法。靳明君等[16]推導(dǎo)了懸鏈線索長(zhǎng)度與水平張力的函數(shù)關(guān)系,但將水平張力作為已知量似有因果倒置之嫌。張卓杰等[17]在懸鏈線的求解中選取左邊界的斜率即廣義傾角α的映射作為初始試探值,反復(fù)迭算從而求得懸鏈線的解,計(jì)算過程過于迂回。馮海暴[18]采用懸鏈線理論對(duì)35m作業(yè)水深鋪排管線進(jìn)行了力學(xué)分析,但采用的方程要求廣義傾角α=0,這一前提很難適用于流體動(dòng)力的變化。劉婷等[19]采用懸鏈線二維解作為跨接管安裝的柔性段三維分析的初始構(gòu)型,但在理論上沒有詳述如何求解二維懸鏈線解的水平張力這一關(guān)鍵問題。呂玉蘭等[20]在求解懸鏈線方程時(shí),以水平張力和懸鏈線最低點(diǎn)的x坐標(biāo)為未知數(shù),算法過于復(fù)雜,且懸鏈線最低點(diǎn)能否選為獨(dú)立變量也有待商榷。龐曉旭等[21]探討了水平懸鏈線在小垂度下彈性模量對(duì)力學(xué)特性的影響。周陽等[22]以給定上部頂張力和懸掛角作為邊界條件,研究了下部海床土非線性剛度對(duì)鋼懸鏈線立管的影響。吳衛(wèi)國(guó)等[23]采用懸鏈線靜態(tài)解研究了船只補(bǔ)給懸鏈線系統(tǒng)的動(dòng)力特性。
總之對(duì)于懸鏈線的數(shù)學(xué)解析解即雙曲余弦函數(shù),如何簡(jiǎn)潔地得到其中最重要的未知參數(shù)即水平張力,依然困擾著工程界,而懸鏈線的數(shù)值求解在張拉結(jié)構(gòu)中具有舉足輕重的作用。小垂度下懸鏈線的拋物線近似解目前在工程上得到了可行的解決途徑[1,15,21,23],但鮮有文獻(xiàn)探討懸鏈線在全局范圍任意垂度下的近似解,本文在此方面作出了有益的嘗試。
若結(jié)構(gòu)任意截面上的彎矩可假設(shè)為0,則結(jié)構(gòu)稱為完全柔性,若結(jié)構(gòu)變形前后可認(rèn)為長(zhǎng)度保持不變,則結(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ì)。
對(duì)于僅有垂直力作用的結(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表示懸鏈線左右兩端的相對(duì)水平距離,簡(jiǎn)稱水平距離;c=y(tǒng)b-ya表示懸鏈線左右兩端的相對(duì)垂直距離,簡(jiǎn)稱垂直距離。
懸鏈線在水平距離l變動(dòng)時(shí)存在一個(gè)不能逾越的極限水平距離d,簡(jiǎn)稱極限距離,
極限距離d對(duì)應(yīng)的線是一條直線,這條直線是懸鏈線的水平距離l無限逼近d時(shí)懸鏈線的漸近線,簡(jiǎn)稱極限漸近線。s為懸鏈線的索長(zhǎng)。
懸鏈線解(1)中高度綁定的表達(dá)式h/q具有深刻的內(nèi)在涵義,表示懸鏈線曲率半徑中的最小值。因此,引入曲率半徑r或曲率直徑D為
r=h/q,D =2r=2h/q (3)
曲率半徑r和懸鏈線形態(tài)本身僅取決于懸鏈線段的索長(zhǎng)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正比于單位均布載荷單位水平距離對(duì)應(yīng)的水平張力,可稱為相對(duì)水平張力,表示懸鏈線的張馳程度。
不失一般性,將坐標(biāo)系放置于左端點(diǎn)a處,利用式(1,3),從起點(diǎn)到任意位置x處的索長(zhǎng)度為
定義廣義傾角
θ=α+β (7)
根據(jù)兩點(diǎn)邊值約束條件,即c=y(tǒng)(l)=y(tǒng)b-ya,與不可拉伸弧長(zhǎng)約束條件,即s=s(l),且利用定義式(3,5)可得到索長(zhǎng)s和垂直距離c與廣義傾角θ和β的關(guān)系為
由式(8,9)可知廣義傾角θ的計(jì)算式為
將式(5,10)代入式(7),得到廣義傾角α的計(jì)算式為
從廣義傾角α的計(jì)算式(11)可知,α不是獨(dú)立的未知量,獨(dú)立的未知量只有一個(gè),即水平張力h。
利用式(2~5,8,9)運(yùn)算得到一組等價(jià)的隱含水平張力的超越方程為
一般地,懸鏈線水平張力h取決于三個(gè)參數(shù),即水平距離l、懸鏈線段的索長(zhǎng)s和垂直距離c。對(duì)方程(12)的無量綱參數(shù)b或n與無量綱參數(shù)β或z進(jìn)行仔細(xì)分析可知,懸鏈線相對(duì)水平張力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ù)等價(jià)的超越方程(12)可得到無量綱參數(shù)β與無量綱參數(shù)b之間相互關(guān)系的離散數(shù)值解,或可得到相對(duì)水平張力z與無量綱參數(shù)n之間的離散數(shù)值解,兩者都是超越方程的精確解或真實(shí)解。
在用C語言進(jìn)行計(jì)算機(jī)數(shù)據(jù)處理時(shí),8字節(jié)的double型實(shí)數(shù)不超過10的±309次方冪,因此計(jì)算中將接近0的double型的比極距離n=l/d的下限取為10-308。計(jì)算程序定義因變量0<n=l/d<1的計(jì)算域?yàn)椋?0-308,0.999999],該計(jì)算域是真小數(shù)區(qū)間(0,1)內(nèi)C語言仿真計(jì)算工程上可行的全域。
圖1將log10(l/d)作為橫坐標(biāo),顯示相對(duì)水平張力精確解的變化情況。n=l/d的計(jì)算域截取為[10-308,0.999999],log10(l/d)的值域?yàn)椋郏?08,0)。因變量n=l/d的計(jì)算域[10-308,10-7]內(nèi),相對(duì)水平張力在接近于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。相對(duì)水平張力精確解隨因變量n=l/d的變化情況如圖1所示。
對(duì)于β>0,數(shù)學(xué)上存在thβ<β,利用懸鏈線段靜力平衡關(guān)系式(13)有
從式(15)可確定相對(duì)水平張力上界z0:
將雙曲正弦函數(shù)泰勒級(jí)數(shù)展開,保留方程(12)的二次冪得到近似根β1為
保留二次冪得到的近似根β1大于真實(shí)的根,即β1>β,考慮到關(guān)系β1z1=zβ=1,因此z1<z;考慮到關(guān)系nb=1,可確定相對(duì)水平張力的下界之一為z1,如式(18)所示。相對(duì)水平張力z介于式(18)中上下界的變動(dòng)范圍:
分析發(fā)現(xiàn),式(16)中z0對(duì)應(yīng)的水平張力上界h0過高地偏離了真實(shí)的水平張力,引入相對(duì)水平張力z2:
計(jì)算表明實(shí)際相對(duì)水平張力z介于式(20)中相對(duì)水平張力下界z1和上界z2之間的變動(dòng)范圍:
以無量綱參數(shù)β=2l/qh和相對(duì)水平張力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%。
相對(duì)水平張力z1和z2是水平距離l趨于極限距離d的近似解。取下界z1和上界z2,兩者的平均值z(mì)m為
計(jì)算表明相對(duì)水平張力z介于相對(duì)水平張力下界z1和上界zm之間的變動(dòng)范圍:
即上界zm是比原上界z2更小的上界,在因變量n的計(jì)算域[0.4,0.999999]內(nèi),zm與精確解z 的誤差不超過10%。n≥0.62時(shí),zm與精確解z的誤差不超過5%;n≥0.68時(shí),zm與精確解z的誤差不超過4%;n≥0.75時(shí),zm與精確解z的誤差不超過3%;n≥0.91時(shí),zm與精確解z的誤差不超過1%;n≥0.991時(shí),zm與精確解z的誤差不超過0.1%。
水平距離l接近極限距離d時(shí)的近似解z1,zm,z2和精確解z在n=l/d介于[10-308,0.95]內(nèi)的變化趨勢(shì),如圖2所示。當(dāng)n=0.95時(shí),z1=1.779513,z=1.793394,zm=1.802627,z2=1.825742。
圖2 相對(duì)水平張力的近似解和精確解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)的變化趨勢(shì)。n=10-308時(shí),βzm=145.025;n=10-7時(shí),βzm=4.04;n=0.001時(shí),βzm=2.084;n=0.99時(shí),βzm=1.001。若工程誤差許可不超過5%,則當(dāng)水平距離l≥0.62d時(shí),水平張力可以由式(4,5,22)表示為
當(dāng)水平距離l接近極限距離d時(shí),水平張力變得充分大。式(24)清晰地表明了水平距離l趨于極限距離d時(shí),水平張力趨于無窮的變化趨勢(shì),公式簡(jiǎn)潔可行。從式(24)或圖2可以看出,保持d不變,l增大,水平張力隨之增大。保持l不變,d增大(s增大,c減?。?,水平張力隨之減小。式(24)在小垂度情形與文獻(xiàn)[1]拋物線近似解趨于一致。
探討l趨于0時(shí)的近似解。根據(jù)懸鏈線段平衡方程sh(ql/2h)=qd/2h可知,當(dāng)水平距離l趨于0時(shí),精確解的水平張力較l緩慢地趨于0。
鑒于C語言8字節(jié)浮點(diǎn)數(shù)的下界值大約為10-308量級(jí),下面計(jì)算上可行地找出無量綱參數(shù)或比極距離n=l/d在計(jì)算域[10-308,10-7]內(nèi)的近似解,該近似解視為l趨于0時(shí)的近似解。
通過多次試探計(jì)算,尋找到比極距離n=l/d在計(jì)算域[10-308,10-7]內(nèi)的近似解zl為
bzl=zl/z是近似解zl與精確解z之比,表示兩者的逼近程度。
水平距離l接近0時(shí)的近似解zl與精確解z之 比bzl,在n=l/d介于[10-20,0.009]內(nèi)的變化趨勢(shì)如圖4所示。其中n=10-308時(shí),b zl=1.000394;n=10-7時(shí),b zl=1.1565;n=0.009時(shí),bzl=1.296896。
圖4 n接近0的近似解zl和精確解之比Fig.4 Ratio of approximate solution zlto the exact one with nnear zero
從計(jì)算結(jié)果可知,近似解zl在[10-308,10-7]內(nèi)的誤差不超過16%,這是一個(gè)初始的近似解。
調(diào)整參數(shù)或選擇另外的函數(shù)等可以找到分段區(qū)域內(nèi)誤差更小、逼近程度更高的近似解。如函數(shù)z4=[1.0-0.5·sqrt(n)]/[1.0-log(n)]較好地模擬了n接近10-308時(shí)精確解z的變化趨勢(shì),將z4適當(dāng)減小約10%,得到近似解zi為
近似解zi與精確解z之比zi/z=bzi在n=l/d介于[10-308,0.25]內(nèi)的變化趨勢(shì)如圖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時(shí),bzi=0.89;n=10-7時(shí),bzi=1.029;n=0.001~0.009時(shí),bzi≈1.1;n=0.25時(shí),bzi=0.91。
圖6將log10l/d作為橫坐標(biāo),顯示近似解和精確解之比zi/z的變化情況,n=l/d的計(jì)算域截取為[10-100,0.999999],log10l/d 的值域?yàn)椋郏?00,0)??梢钥闯觯瑹o量綱水平距離n介于該范圍內(nèi)時(shí),bzi=zi/z的細(xì)微變化情況。
圖5 近似解zi和精確解之比Fig.5 Ratio of approximate solution zito the exact one
圖6 近似解和精確解之比zi/z的對(duì)數(shù)顯示Fig.6 Logarithmic displaying the ratio zi/zof approximate solution to the exact one
從圖3可以看出,當(dāng)n=10-308時(shí),zm/z=145.025278,即當(dāng)n接近于0時(shí),上界zm過高估算了精確解。因此將zm適當(dāng)減小,找到逼近精確解的近似解zn
關(guān)系式lgn=lgelnn=Mlnn表示常用對(duì)數(shù)與自然對(duì)數(shù)之間是線性映射關(guān)系。其中模數(shù)[25]M為M=lge=1/ln10≈0.434294481903…。
計(jì)算結(jié)果表明,zn與精確解z的比值bzn存在幾個(gè) 關(guān) 鍵 點(diǎn) 函 數(shù) 對(duì) 應(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時(shí)略微大于1,然后單調(diào)增加,在n=10-100附近時(shí)接近最大值1.086182,隨后單調(diào)減少,在n=0.04附近,bzn接近最小值0.87771,之后bzn隨n 單調(diào)增加,n趨于1時(shí)bzn趨于1。如圖7和圖8所示。
圖7 近似解zn和精確解之比Fig.7 Ratio of approximate solution znto the exact one
圖8 近似解和精確解之比zn/z的對(duì)數(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的計(jì)算域截取為[10-308,10-25],log10l/d 的值域?yàn)椋郏?08,-25],可以看出近似解和精確解之比bzn=zn/z的細(xì)微變化情況。圖8(b)中無量綱水平距離n計(jì)算域截取為[10-100,0.999999],log10l/d的值域?yàn)椋郏?00,0),可以看出bzn=zn/z的大致變化趨勢(shì)。
從計(jì)算結(jié)果分析可知,比值bzn在n的全計(jì)算域內(nèi)介于0.877<bzn(n)<1.0862;因此由公式(27)計(jì)算的相對(duì)水平張力zn是因變量n=l/d全域范圍內(nèi)誤差不超過13%的近似解。
首先,在兩點(diǎn)位移約束方程基礎(chǔ)上,自然地補(bǔ)充引進(jìn)弧長(zhǎng)不可拉伸約束方程,從而簡(jiǎn)單解決了求解經(jīng)典懸鏈線數(shù)學(xué)解的兩個(gè)未知參數(shù)即水平張力h和表示左邊界處斜率的廣義傾角α這一非線性數(shù)學(xué)問題。其次,推導(dǎo)出求解懸鏈線水平張力的超越方程。先指明h/q為懸鏈線的最小曲率半徑,后引進(jìn)無量綱參數(shù)z=2h/ql和n=l/d求解超越方程中的水平張力,使得水平張力形式上具有最簡(jiǎn)單的參數(shù)依賴關(guān)系。探討了廣義傾角β,θ和α的內(nèi)在含義以及與幾何參數(shù)的相互關(guān)系,創(chuàng)新地給出θ和α的計(jì)算公式。從廣義傾角α的計(jì)算公式可知,廣義傾角α不再是獨(dú)立的原始未知量,獨(dú)立的原始未知量只有一個(gè)即水平張力h。
隨之,通過對(duì)平衡方程的泰勒展開,經(jīng)過理論分析和反復(fù)計(jì)算得到水平張力的上下包絡(luò)界,為確定水平張力的精確數(shù)值解提供算法上的理論指導(dǎo)。
最后,提出了水平距離l分別趨于0與趨于極限距離d的近似解zl,zi和zm等以及在真小數(shù)計(jì)算范圍[10-308,0.999999]內(nèi)全局可行的近似解zn,分析計(jì)算了這些近似解關(guān)于精確解的誤差程度,近似解以函數(shù)的形式粗顯地刻畫了水平張力隨n=l/d的因果變化趨勢(shì),對(duì)于水平張力工程的初步快速估算具有重要的意義,也可對(duì)數(shù)值計(jì)算的發(fā)展扮演拋磚引玉的作用。這些近似解可作為初值加快精確數(shù)值解的非線性迭代計(jì)算。
綜上所述,本文關(guān)于懸鏈線提出的相關(guān)公式和近似解,在工程應(yīng)用上具有實(shí)際的價(jià)值。