经典悬链线理论精确解与近似解的非线性数值计算

郭小刚， 金 星， 周 涛， 宋晓东， 邓旭辉

（1．湘潭大学 土木工程与力学学院，湘潭411105；2．长沙矿冶研究院、深海矿产资源开发利用技术国家重点实验室，长沙410012）

1 引 言

张拉结构的应用日益广泛，已经覆盖如大跨度空间结构的索网、海洋工程中的软管和柔索［1－4］等结构。目前国内外对张拉结构的分析手段亦日趋复杂，采用两节点直线单元［5－7］、三节点抛物线单元甚至五节点等参数索单元［8］的非线性有限元分析，特别在动力时域分析时沿用Wilson法Newmark法内蕴 Newton－Raphson非线性迭代［9－12］，使得工程应用复杂繁琐，且处理的大多是小垂度的问题。

悬链线问题的初步数学解答已有300多年的历史，关于采用悬链线解答求解工程问题的文献层出不穷［1，12－24］。 沈 世 钊 等［1］介 绍 了 单 索 均 布 力 作用下的悬链线理论以及悬链线解的抛物线近似，但留下水平张力和左边界端点斜率两个未知数悬而未决。郑丽凤等［13］选取水平张力和左端边界垂直力这两个参数作为独立未知量，建立了二元非线性牛顿迭代法方程组来求解悬链线的解，但左端边界垂直力未必适宜作为独立的参数。杨孟刚等［14］采用UL列式增量法，提出了两节点悬链线索元非线性切线刚度矩阵，但未明确如何计算最重要的水平初始张力。周绪红等［15］采用的悬索理论实质上是悬链线理论的抛物线逼近，是工程上可行的方法。靳明君等［16］推导了悬链线索长度与水平张力的函数关系，但将水平张力作为已知量似有因果倒置之嫌。张卓杰等［17］在悬链线的求解中选取左边界的斜率即广义倾角α的映射作为初始试探值，反复迭算从而求得悬链线的解，计算过程过于迂回。冯海暴［18］采用悬链线理论对35m作业水深铺排管线进行了力学分析，但采用的方程要求广义倾角α＝0，这一前提很难适用于流体动力的变化。刘婷等［19］采用悬链线二维解作为跨接管安装的柔性段三维分析的初始构型，但在理论上没有详述如何求解二维悬链线解的水平张力这一关键问题。吕玉兰等［20］在求解悬链线方程时，以水平张力和悬链线最低点的x坐标为未知数，算法过于复杂，且悬链线最低点能否选为独立变量也有待商榷。庞晓旭等［21］探讨了水平悬链线在小垂度下弹性模量对力学特性的影响。周阳等［22］以给定上部顶张力和悬挂角作为边界条件，研究了下部海床土非线性刚度对钢悬链线立管的影响。吴卫国等［23］采用悬链线静态解研究了船只补给悬链线系统的动力特性。

总之对于悬链线的数学解析解即双曲余弦函数，如何简洁地得到其中最重要的未知参数即水平张力，依然困扰着工程界，而悬链线的数值求解在张拉结构中具有举足轻重的作用。小垂度下悬链线的抛物线近似解目前在工程上得到了可行的解决途径［1，15，21，23］，但鲜有文献探讨悬链线在全局范围任意垂度下的近似解，本文在此方面作出了有益的尝试。

2 悬链线超越方程与数值解

若结构任意截面上的弯矩可假设为0，则结构称为完全柔性，若结构变形前后可认为长度保持不变，则结构具有不可拉伸性。经典悬链线理论采用了完全柔性和不可拉伸性假设，因其采用完全柔性假设，相应的悬挂结构称为柔索结构或柔性结构，但因为不可拉伸性假设，柔索结构沿纵向具有内在的刚性特质。

对于仅有垂直力作用的结构，其任一截面水平方向的合内力处处相等。柔索在仅有均布垂直载荷q作用下经典的悬链线解［1，4，24］为

式中ch为双曲余弦函数；xa和ya分别为悬链线左端点a处的x方向坐标和y方向坐标；独立的待求量h为沿悬链线不变的水平张力，h＞0；待求量α为左边界a处斜率的广义倾角，shα为边界a处的斜率。线均布力q的单位是N／m，水平张力h的单位是N。

悬链线右端点b处的x方向坐标为xb，y方向坐标为yb。l＝xb－xa表示悬链线左右两端的相对水平距离，简称水平距离；c＝yb－ya表示悬链线左右两端的相对垂直距离，简称垂直距离。

悬链线在水平距离l变动时存在一个不能逾越的极限水平距离d，简称极限距离，

极限距离d对应的线是一条直线，这条直线是悬链线的水平距离l无限逼近d时悬链线的渐近线，简称极限渐近线。s为悬链线的索长。

悬链线解（1）中高度绑定的表达式h／q具有深刻的内在涵义，表示悬链线曲率半径中的最小值。因此，引入曲率半径r或曲率直径D为

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

曲率半径r和悬链线形态本身仅取决于悬链线段的索长s、水平距离l和垂直距离c，而与作用其上的均布载荷q无关。定义无量纲参数b和n为

无量纲参数b和n为互逆关系，bn＝1。无量纲参数n＝l／d，为位于区间（0，1）的正小数，可称为比极距离。无量纲参数n越大，曲率半径r越大，n越接近于1，曲率半径r趋于无穷大，悬链线越逼近悬链线的极限渐近线即越逼近直线；无量纲参数n越小，曲率半径r越小，水平距离l越接近于0，曲率半径r趋于0。

定义无量纲参数β和z为

无量纲参数β和z为互逆关系，βz＝1。无量纲参数z正比于单位均布载荷单位水平距离对应的水平张力，可称为相对水平张力，表示悬链线的张驰程度。

不失一般性，将坐标系放置于左端点a处，利用式（1，3），从起点到任意位置x处的索长度为

定义广义倾角

θ＝α＋β （7）

根据两点边值约束条件，即c＝y（l）＝yb－ya，与不可拉伸弧长约束条件，即s＝s（l），且利用定义式（3，5）可得到索长s和垂直距离c与广义倾角θ和β的关系为

由式（8，9）可知广义倾角θ的计算式为

将式（5，10）代入式（7），得到广义倾角α的计算式为

从广义倾角α的计算式（11）可知，α不是独立的未知量，独立的未知量只有一个，即水平张力h。

利用式（2～5，8，9）运算得到一组等价的隐含水平张力的超越方程为

一般地，悬链线水平张力h取决于三个参数，即水平距离l、悬链线段的索长s和垂直距离c。对方程（12）的无量纲参数b或n与无量纲参数β或z进行仔细分析可知，悬链线相对水平张力z＝2h／q l仅依赖于水平距离l与水平极限距离d之比n＝l／d，这是悬链线两点边值约束问题相当明确且可靠的结论。

根据双曲函数的性质［25］和式（12），关于广义倾角β存在关系为

式（13）可视为悬链线段两点边界约束的静力平衡关系，其中

根据等价的超越方程（12）可得到无量纲参数β与无量纲参数b之间相互关系的离散数值解，或可得到相对水平张力z与无量纲参数n之间的离散数值解，两者都是超越方程的精确解或真实解。

在用C语言进行计算机数据处理时，8字节的double型实数不超过10的±309次方幂，因此计算中将接近0的double型的比极距离n＝l／d的下限取为10－308。计算程序定义因变量0＜n＝l／d＜1的计算域为［10－308，0．999999］，该计算域是真小数区间（0，1）内C语言仿真计算工程上可行的全域。

图1将log10（l／d）作为横坐标，显示相对水平张力精确解的变化情况。n＝l／d的计算域截取为［10－308，0．999999］，log10（l／d）的值域为［－308，0）。因变量n＝l／d的计算域［10－308，10－7］内，相对水平张力在接近于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，数学上存在thβ＜β，利用悬链线段静力平衡关系式（13）有

从式（15）可确定相对水平张力上界z0：

将双曲正弦函数泰勒级数展开，保留方程（12）的二次幂得到近似根β1为

保留二次幂得到的近似根β1大于真实的根，即β1＞β，考虑到关系β1z1＝zβ＝1，因此z1＜z；考虑到关系nb＝1，可确定相对水平张力的下界之一为z1，如式（18）所示。相对水平张力z介于式（18）中上下界的变动范围：

分析发现，式（16）中z0对应的水平张力上界h0过高地偏离了真实的水平张力，引入相对水平张力z2：

计算表明实际相对水平张力z介于式（20）中相对水平张力下界z1和上界z2之间的变动范围：

4 悬链线超越方程的近似解

以无量纲参数β＝2l／qh和相对水平张力z＝D／l表示方程sh（l／D）＝d／D 或sh（β）＝bβ的精确解，βz＝1。若能找到无量纲函数zk（n），在0＜n1≤n≤n2＜1域内满足：

则zk（n）就是方程sh（l／D）＝d／D 或sh（β）＝bβ在n1≤n≤n2域内的无量纲近似解或分段近似解。若近似解与精确解的比值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］内，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］内的变化趋势，如图2所示。当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＝10－308时，βzm＝145．025；n＝10－7时，βzm＝4．04；n＝0．001时，βzm＝2．084；n＝0．99时，βzm＝1．001。若工程误差许可不超过5%，则当水平距离l≥0．62d时，水平张力可以由式（4，5，22）表示为

当水平距离l接近极限距离d时，水平张力变得充分大。式（24）清晰地表明了水平距离l趋于极限距离d时，水平张力趋于无穷的变化趋势，公式简洁可行。从式（24）或图2可以看出，保持d不变，l增大，水平张力随之增大。保持l不变，d增大（s增大，c减小），水平张力随之减小。式（24）在小垂度情形与文献［1］抛物线近似解趋于一致。

探讨l趋于0时的近似解。根据悬链线段平衡方程sh（ql／2h）＝qd／2h可知，当水平距离l趋于0时，精确解的水平张力较l缓慢地趋于0。

鉴于C语言8字节浮点数的下界值大约为10－308量级，下面计算上可行地找出无量纲参数或比极距离n＝l／d在计算域［10－308，10－7］内的近似解，该近似解视为l趋于0时的近似解。

通过多次试探计算，寻找到比极距离n＝l／d在计算域［10－308，10－7］内的近似解zl为

bzl＝zl／z是近似解zl与精确解z之比，表示两者的逼近程度。

水平距离l接近0时的近似解zl与精确解z之 比bzl，在n＝l／d介于［10－20，0．009］内的变化趋势如图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

从计算结果可知，近似解zl在［10－308，10－7］内的误差不超过16%，这是一个初始的近似解。

调整参数或选择另外的函数等可以找到分段区域内误差更小、逼近程度更高的近似解。如函数z4＝［1．0－0．5·sqrt（n）］／［1．0－log（n）］较好地模拟了n接近10－308时精确解z的变化趋势，将z4适当减小约10%，得到近似解zi为

近似解zi与精确解z之比zi／z＝bzi在n＝l／d介于［10－308，0．25］内的变化趋势如图5所示。因变量n在区域［10－308，0．25］内，zi与精确解z 的比值bzi介于0．89＜bzi（n）＜1．10154，在此范围内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作为横坐标，显示近似解和精确解之比zi／z的变化情况，n＝l／d的计算域截取为［10－100，0．999999］，log10l／d 的值域为［－100，0）。可以看出，无量纲水平距离n介于该范围内时，bzi＝zi／z的细微变化情况。

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

图6 近似解和精确解之比zi／z的对数显示Fig．6 Logarithmic displaying the ratio zi／zof approximate solution to the exact one

从图3可以看出，当n＝10－308时，zm／z＝145．025278，即当n接近于0时，上界zm过高估算了精确解。因此将zm适当减小，找到逼近精确解的近似解zn

关系式lgn＝lgelnn＝Mlnn表示常用对数与自然对数之间是线性映射关系。其中模数［25］M为M＝lge＝1／ln10≈0．434294481903…。

计算结果表明，zn与精确解z的比值bzn存在几个 关 键 点 函 数 对 应 关 系，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，然后单调增加，在n＝10－100附近时接近最大值1．086182，随后单调减少，在n＝0．04附近，bzn接近最小值0．87771，之后bzn随n 单调增加，n趋于1时bzn趋于1。如图7和图8所示。

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

图8 近似解和精确解之比zn／z的对数显示Fig．8 Logarithmic displaying the ratio zn／zof approximate solution to the exact one

图8 将log10l／d作为横坐标，图8（a）显示近似解和精确解之比zn／z的变化。n＝l／d的计算域截取为［10－308，10－25］，log10l／d 的值域为［－308，－25］，可以看出近似解和精确解之比bzn＝zn／z的细微变化情况。图8（b）中无量纲水平距离n计算域截取为［10－100，0．999999］，log10l／d的值域为［－100，0），可以看出bzn＝zn／z的大致变化趋势。

从计算结果分析可知，比值bzn在n的全计算域内介于0．877＜bzn（n）＜1．0862；因此由公式（27）计算的相对水平张力zn是因变量n＝l／d全域范围内误差不超过13%的近似解。

5 结 语

首先，在两点位移约束方程基础上，自然地补充引进弧长不可拉伸约束方程，从而简单解决了求解经典悬链线数学解的两个未知参数即水平张力h和表示左边界处斜率的广义倾角α这一非线性数学问题。其次，推导出求解悬链线水平张力的超越方程。先指明h／q为悬链线的最小曲率半径，后引进无量纲参数z＝2h／ql和n＝l／d求解超越方程中的水平张力，使得水平张力形式上具有最简单的参数依赖关系。探讨了广义倾角β，θ和α的内在含义以及与几何参数的相互关系，创新地给出θ和α的计算公式。从广义倾角α的计算公式可知，广义倾角α不再是独立的原始未知量，独立的原始未知量只有一个即水平张力h。

随之，通过对平衡方程的泰勒展开，经过理论分析和反复计算得到水平张力的上下包络界，为确定水平张力的精确数值解提供算法上的理论指导。

最后，提出了水平距离l分别趋于0与趋于极限距离d的近似解zl，zi和zm等以及在真小数计算范围［10－308，0．999999］内全局可行的近似解zn，分析计算了这些近似解关于精确解的误差程度，近似解以函数的形式粗显地刻画了水平张力随n＝l／d的因果变化趋势，对于水平张力工程的初步快速估算具有重要的意义，也可对数值计算的发展扮演抛砖引玉的作用。这些近似解可作为初值加快精确数值解的非线性迭代计算。

综上所述，本文关于悬链线提出的相关公式和近似解，在工程应用上具有实际的价值。