经典悬链曲线常采用解析方法[1]求解，具有弹性的悬链或一般力作用下的悬链也可用解析方法描述，但总的来说解析解适用性不强，不能很好地面对结构特征沿链非均匀或载荷条件特殊的悬链，如带有负重的悬索桥。在工程中常使用数值方法来求解结构变形。其中，无论是有限元方法还是传递矩阵法，分析结构系统力学特性时总是将模型离散划分为若干简单单元，组集单元的质量、刚度矩阵及位移向量等，获得整体的变形平衡方程并进行求解。但这2种方法更适用于小变形和线弹性假设下的问题，并不适合求解悬链线这种基于大变形、强非线性假设的问题。例如，使用杆单元[2]或索单元[3-4]，前者面临强非线性的困难，后者则不适用于复杂载荷条件分析且不满足C1连续条件[5]。
本文提出一种较为通用的数值方法，可以求解实际悬链结构在构形复杂或受力复杂情况下的曲线构形。基于传递矩阵法的思想，构建悬链的力学模型，将悬链分为首尾相接的悬链单元，分析解得单元首尾状态参数间的函数关系方程组(以下称特征函数)，将函数顺序嵌套获得系统整体的特征函数方程组。结合边界条件求解该方程组，获得模型各单元状态参数及曲线构形。利用解析方法验证该数值方法的正确性。进一步分析表明，该方法对具有复杂结构特征和载荷特征的悬链线问题有很好的适用性，同时对求解其他特征类似结构的广义变形也适用。
1 悬链力学模型及求解方法 图 1为悬链力学模型。沿悬链方向将悬链划分成特性统一的若干单元，相邻单元有共同的节点，设该结构系统沿轴向划分为n个单元，n+1个节点，首尾两端为边界。其中，第i个单元左右2个节点的编号分别为i和i+1。
 |
| 图 1 悬链结构系统力学模型 Fig. 1 Mechanical model of cable chain structure system |
| 图选项 |
在载荷之下，悬链产生变形，节点发生位移，如图 2所示，各单元在节点位置处有力和位移的相互联系，而后传载到边界上，最终达到系统内力与载荷的平衡。模型中每个单元是规范简单的，在已知一端节点广义位移和单元所受广义外力的情况下，可以求解单元的变形与另一端的广义位移。
 |
| 图 2 单元受力分析 Fig. 2 Load analysis of element |
| 图选项 |
1.1 系统特征函数的建立 图 2为单元受力与变形示意图。第i单元的左端受到来自节点的力为Fi，右端受到来自节点的力为F′ i+1。连接2个单元的节点的受力必是平衡的，故Fi+1与F′i+1大小相等，方向相反。单元在特定载荷条件下受到载荷Qi。设节点拥有m个自由度，则上述力的向量都有m个分量。
依据理论力学或弹性力学的方法，可以建立描述单元受力、变形和节点几何联系的方程组[6]，如下：
 | (1) |
 | (1b) |
式中：yi, 1, …, yi, m, yi+1, 1, …, yi+1, m分别为悬链受载后第i单元左右节点m个自由度上的位移；Fi, 1, …, Fi, m, Fi+1, 1, …, Fi+1, m为对应节点对应自由度上的广义力。
式(1)融合了2类方程：第1类方程为单元本构-几何方程，描述了单元在载荷之下的变形和节点间的几何联系；第2类方程为广义力的平衡方程。每个单元的平衡方程组(1)包括m个(节点自由度个数)力平衡方程和m个本构-几何方程。
令第i个节点的状态参数向量为zi，该向量由广义位移和广义力组成，即
 | (2) |
在各参数为线性关系下，可以用矩阵表达已知单元一端的状态向量zi与单元另一端状态向量zi+1间的关系，这时的方法就是传递矩阵法[6-7]，但一般情况下应为复杂函数关系，也即式(1)。简写式(1)，得到
 | (3) |
式中：CFi为单元特征函数组。
这一过程是具体结构在具体载荷下发生特定变形的体现，取决于所研究的串列结构系统的变形特性。式(3)约束了每个单元的左右端受力和变形状态，又因为悬链系统是串联排布首尾相接的，所以系统首尾节点的状态参数也存在约束关系，这种关系是整个系统特性受载和变形的体现。式(4)所示为系统首尾状态向量间的联系：
 | (4) |
CF为整体特征函数组。式(4)为具有4m个变量、2m个方程组成的方程组，该方程组就是悬链系统整体特征函数方程组。只需要根据边界条件，找到2m个边界条件，方程的解便是系统首尾两端未知的状态参数。
1.2 边界条件确定 从弹性力学的角度理解，结构系统的边界条件决定微分方程在积分求解时的积分常数，进而影响整体的变形效果。这里边界条件设置与传递矩阵法[8]中的基本一致。
常见的边界条件主要有简支、自由、弹性支承等，其状态参数如表 1所示。表中：yx为节点未知挠度，θx为节点未知挠角，Fx为节点未知剪力，Kspring为弹性支承刚度。实际问题中，悬链的边界条件一般是简支，但理论上允许其他边界条件，并不局限于表 1所述的几种情况，应根据具体问题来分析获得。一般来说，边界节点的2m个状态参数中，一半参数是独立的，另一半为0或者可以用已知独立的参数来表示。
表 1 常见边界条件参数 Table 1 Common boundary condition parameters
| 边界条件 | 挠度 | 挠角 | 弯矩 | 剪力 |
| 简支 | 0 | θx | 0 | Fx |
| 自由 | yx | θx | 0 | 0 |
| 弹性支承 | yx | θx | 0 | Kspringyx |
表选项
1.3 方程组求解 方程组(4)一般是非线性的。求解非线性方程组的数值方法在数学上是比较完备的，如迭代法、同伦法[9]、人工智能法等[10]，其中后2种属于大范围收敛法。大范围收敛法对初值的选取几乎没有特别的要求，常可以求解方程全部的解。但是由于悬链线构形往往可以根据经验获得，这里选择较为实用的、以Newton迭代法为基础派生的离散Newton迭代法进行求解[11]，该方法的收敛速度较快，一般是超线性收敛[10]。
令方程组(4)中的2m个未知状态参数向量为zx，改写该方程的形式为
 | (5) |
设zx*为方程组精确解，zx(k)为第k个近似解。
由一阶Taylor公式，可得
 |
现用线性方程组
 | (6) |
近似代替非线性方程组(4)。但由于非线性函数组偏导数的计算比较复杂，G′(zx(k))难以获取，故以差商J(zx(k), h(k))代替偏导数。令
 |
 | (7) |
式中：h(k)=c‖G(zx(k)) ‖，c为绝对值较小的非零常数，能使差商更接近导数；ej为第j个2m维基本单位向量。
最终，系统变形求解的迭代公式为
 | (8) |
在经验解或预测解附近选取初始解zx(0)，而后迭代计算式(8)，可快速获得满足误差要求的精确解zx*。
对于某些线性变形求解问题，实际上只需要一步迭代就能解得变形结果，因为此时差商等于导数，所以一阶Taylor公式是严格成立的。
在获得系统左端状态向量z1后，利用式(3)，可以顺序求解每个节点的状态向量，进而可得每个单元与整个系统的变形和受力。
1.4 悬链线问题求解流程 根据上述数学求解方法，建立使用特征函数传递法求解悬链线构形的流程，如图 3所示。该求解过程主要分为5部分：
 |
| 图 3 结构变形求解流程 Fig. 3 Flowchart for solving structure deformation |
| 图选项 |
3) 代入边界条件和初始解，使用数值方法求解非线性方程组，获得首尾节点状态参数。
4) 对于某些非线性或大变形问题，变形结果会一定程度改变系统结构特征或载荷特征，故需要根据初算结果对模型作变形修正，迭代收敛后方可获得最终变形。
5) 计算系统各节点受力和位移，输出或绘图。
2 方法验证 使用经典悬链线解析方法对函数传递法及其求解流程进行验证。
经典悬链线具有以下3个假设特征：①悬链是理想柔性的，不能承受弯矩，只能承受沿索链的拉力；②索链的密度等特征参数沿链均匀；③索链在重力作用下自然变形，受力前后长度不变、截面特征不变。
2.1 经典悬链线特征函数 将悬链分为n个单元，第i个单元的长度和质量分别为li和mi，视单元为刚体杆件，单元质量集中在单元中间。单元间的连接方式视为铰接，只传递力而不传递力矩。单元受力如图 4所示。
 |
| 图 4 刚体杆件单元受力分析 Fig. 4 Load analysis of rigid rod element |
| 图选项 |
每个单元均视为刚体杆件，整个单元各处挠角相同，节点施加到杆件上的弯矩为0。这些可作为简化方程的依据。但应注意，每个节点两侧单元的挠角一般不相同。
单元节点具有3个自由度，列出力平衡和几何关系方程组(9)，该方程组由6个方程组成。
 | (9) |
式中：g为重力加速度。
对方程组(9)进行变化，消去θi，方程组由6项变4项，将i+1项全部放到等号左侧，i项全部放到等号右侧，为
 | (10) |
即
 |
上式为经典索链单元特征函数组。
由式(4)可进一步得到系统特征函数组。
2.2 结果验证 令悬链线两端支承等高，索链边界条件为简支，故位移参数x1、xn+1、y1、yn+1为已知条件，力参数T1, x、Tn+1, x、T1, y、Tn+1, y为未知条件。索链参数如表 2所示，重力加速度g取9.8 m/s2。
表 2 悬链算例参数 Table 2 Example parameters of catenary
| 参数 | 数值 |
| 跨度Lspan/m | 1.5 |
| 索链总长Lline/m | 2.0 |
| 横截面圆半径Rline/m | 1.0×10-3 |
| 材料密度ρ/(kg·m-3) | 7 800 |
表选项
将悬链划分为首尾相接的20个单元，各单元尺寸一致。按图 3所示流程编程计算得到悬链数值解，并与解析解[1]做对照，如图 5所示，关键数据结果如表 3所示。
 |
| 图 5 悬链线解析解与数值解曲线 Fig. 5 Analysis and numerical solution curves of catenary |
| 图选项 |
表 3 悬链算例结果 Table 3 Example results of catenary
| 解析/数值解 | 高低差Δy/m | 索端水平拉力Fx-end/N | 索端角θend/(°) |
| 解析解 | 0.588 6 | 0.0137 8 | 60.96 |
| 数值解 | 0.589 5 | 0.0135 9 | 59.72 |
| 相对误差/% | 0.15 | -1.4 | -2.0 |
表选项
2种方法的曲线结果几乎重叠。高低差的数值解相较解析解误差仅约0.15%，索端水平拉力的误差约为-1.4%，索端角略微大一点，达到-2.0%。索端角误差略大实际上是悬链划分单元太粗糙导致的，因为这时的数值索端角为解析悬链线的割线角而非切线角。
结果的一致性表明了特征函数传递方法的正确性。
实际上，该算例整体特征函数方程组可以获得2组解，如图 6所示，分为受拉解和受压解。受拉解正是图 5所示结果，该结果每个结构单元两端受拉力而平衡，受压解则是每个单元两端受压力而平衡的结果，此时悬链线的最高点在跨中位置。显然，实际索链并不存在受压平衡这种不稳定的结果，然而拱桥往往是这种受压构形。
 |
| 图 6 两组悬链线解 Fig. 6 Two groups of catenary solutions |
| 图选项 |
3 悬链线特性分析 基于数值方法的正确性，可以进一步做悬链线构形参数变化对力学特性影响的探讨。
图 7和图 8为两端支点等高情况下，不同跨度Lspan与长度Lline比值(以下简称跨长比)的悬链线的曲线造型及索端角θend和跨度与最大横向挠度比值(简称跨横比)的变化曲线。
 |
| 图 7 不同跨长比的悬链线 Fig. 7 Catenary with different span/length ratios |
| 图选项 |
 |
| 图 8 悬链线构形参数间的变化关系 Fig. 8 Relationship between configuration parameters of catenary |
| 图选项 |
由图 7和图 8可知，当跨长比从0逐渐增加到1的过程中，索端角由90°渐变为0°，最大横向挠度趋小，跨横比趋大。由于索链支点竖直方向的力Fy=mg/2保持不变，水平力为Fx=Fycot θend，跨长比接近1而索端角趋于0°时水平力Fx急剧增大，图 7中力三角形的底边随索端角变小而急剧变长说明了这一点。用公式表达即为
 |
因此理论上说，很难用力把悬链完全水平拉直，索端角只能趋近于0而不能为0。
4 方法适用性的进一步探讨 函数传递法不单适用于复杂条件悬链线问题，还可延伸到更多问题中去。此外，若在某种情况下，zi+1=CFi(zi)可以用线性代数方程组来表达，则方法就演变成经典的传递矩阵法[12-14]。
4.1 关于悬链线问题的延伸 实际上，函数传递法相较于解析法更具实用性。原因在于：该方法可以方便调整参数以适应更复杂的结构和载荷，如悬链质量分布不均匀且具有拉伸性，或某节点具有集中力、载荷与挠度有关联性的情况。以下简述上述多种情况的解决方法：
1) 若索链质量或其他结构特征沿链不均匀，只需修改各单元初始特征参数li和mi即可；若索链受到外力，可将之等效为附加质量，即m=F/g，加到相应单元即可。
图 9为悬链某处附加重物前后曲线造型对比。附加重物后的位置曲线曲率更大，水平位置更低。
 |
| 图 9 经典悬链线解与附加重物悬链线解 Fig. 9 Classical catenary solution and catenary with additional weight |
| 图选项 |
2) 如果各索链单元具有可拉伸性(更具实际意义)，设各单元轴向拉伸刚度Ki，则在纯刚性悬链单元的受力基础上，算出各单元的拉伸变形量而得到各自新的长度，返回修改初始单元长度，迭代求解，直至收敛即可。给出各单元拉长量为
 | (11) |
图 10为不可伸长的经典悬链线与带弹性悬链线对比，后者水平位置更低。
 |
| 图 10 经典悬链线解与带弹性悬链线解 Fig. 10 Classical catenary solution and catenary with elasticity |
| 图选项 |
3) 当载荷与悬索变形结果有关联性时，将质量力mig写成挠度的函数，而后再次解得单元特征函数即可。
例如，质量力与单元跨度成正比的情况(主要对应悬索桥[15])，改式(9)中重力项mig为ρ(xi+1-xi)即可，再解方程得到单元特征函数组；若质量力与挠度成正比的情况(对应离心力下跳绳的曲线问题[16-17])，需替换式(9)中重力项mig为离心力项-miΩ2(yi+1+yi)/2(Ω为跳绳进动角速度)。
另外一种处理载荷与变形结果存在关联性情况的手段是迭代。可先将问题视为经典悬链线条件进行初算，而后根据该初算变形结果计算载荷，求得等效质量，返回修改初始参数mi，迭代计算，直至收敛。图 11和图 12分别绘出了悬索桥曲线和跳绳曲线。由于这2种情况悬链跨中附近单位索长载荷集度较两端更大，跨中附近水平位置较经典悬链线解更低，曲线曲率更大。
 |
| 图 11 经典悬链线与悬索桥曲线 Fig. 11 Classical catenary and suspension bridge curve |
| 图选项 |
 |
| 图 12 经典悬链线与跳绳曲线 Fig. 12 Classical catenary and rope skipping curve |
| 图选项 |
4.2 方法与应用对象的演变和延伸 理论上说，本文所述的特征函数传递法对求解具有以下特征的结构的广义变形皆适用：结构单元首尾相接，结构各处状态参数存在可描述的力学或数学联系。例如，梁结构的变形便可使用该方法来求解：①在悬链线问题的基础上，在结构单元节点处加入扭簧[18]，注意此时结构单元仍旧是刚性的；②引入材料力学的梁单元，在小变形和线弹性的条件下，单元两端的结构特征参数(挠度、挠角、剪力、弯矩)具有线性关系，节点两侧的挠角是相等的，用矩阵表达式(3)为[8]
 | (12) |
式中：M为弯矩；Fτ为剪力；li为梁单元长度；Ei为材料弹性模量；Ii为梁截面惯性矩。
进而式(4)可用线性代数方程组来表达。求解该线性代数方程组是较为简单的，但这种线性表达只适用少数情况。
当式(3)可以用线性代数方程组表达时，函数传递法就演变成了传递矩阵法。函数传递法正是基于传递矩阵法思想提出的，相较传递矩阵法更为普适，它所描述的结构的状态参数间关系不限于线性关系。但这也决定了函数形式的复杂性，整体特征函数方程组用函数嵌套的形式来表达而非矩阵相乘，方程的求解略显繁琐。
5 结论 1) 提炼了悬链的力学模型，并给出单元特征函数及整体特征函数的建立方法，使用离散Newton迭代法对非线性整体特征函数方程组进行求解，进而获得索链受力和各处位移。经算例验证，该方法能够正确地求解悬链线问题。
2) 经典悬链线问题里，当跨长比趋向于1时，索端角趋于0，索端水平力趋于极大。
3) 函数传递法对复杂结构特征复杂载荷的悬链线问题具有很好的适用性，同时该方法对求解其他特征类似的结构的广义变形也适用。传递矩阵法实际上正是该方法在线性条件下的演变。
参考文献
| [1] | ROUTH E J. A treatise on analytical statics: With numerous examples[M]. Cambridge: Cambridge University Press, 2013. |
| [2] | 周强, 杨文兵, 杨新华. 斜拉桥索力调整在ANSYS中的实现[J]. 华中科技大学学报(城市科学版), 2005, 22(增刊): 81-83. ZHOU Q, YANG W B, YANG X H. Simulation of tensing stayed cables in ANSYS[J]. Journal of Huazhong University of Science and Technology (Urban Science Edition), 2005, 22(Supplement): 81-83. (in Chinese) |
| [3] | THAI H T, KIM S E. Nonlinear static and dynamic analysis of cable structures[J]. Finite Elements in Analysis and Design, 2011, 47(3): 237-246. DOI:10.1016/j.finel.2010.10.005 |
| [4] | JAYARAMAN H B, KNUDSON W C. A curved element for the analysis of cable structures[J]. Computers and Structures, 1981, 14(3-4): 325-333. DOI:10.1016/0045-7949(81)90016-X |
| [5] | 郑平芳. 有限元中索单元的研究及应用进展[J]. 山西建筑, 2011, 37(28): 42-43. ZHENG P F. Research on cable element in finite element and its application progress[J]. Shanxi Architecture, 2011, 37(28): 42-43. (in Chinese) |
| [6] | PESTEL E, LECKIE F. Matrix methods in elastomechanics[M]. New York: McGraw-Hill Book Company, 1963. |
| [7] | RUBIN S. Transmission matrices for vibration and their relation to admittance and impedance[J]. Journal of Manufacturing Science and Engineering, 1964, 86(1): 9-21. |
| [8] | 顾家柳, 丁奎元, 刘启洲, 等. 转子动力学[M]. 北京: 国防工业出版社, 1985. GU J L, DING K Y, LIU Q Z, et al. Rotordynamics[M]. Beijing: National Defense Industry Press, 1985. (in Chinese) |
| [9] | 夏林林, 吴开腾. 大范围求解非线性方程组的指数同伦法[J]. 计算数学, 2014, 36(2): 215-224. XIA L L, WU K T. An exponential homotopy method for nonlinear equations in large scope[J]. Mathematica Numerica Sinica, 2014, 36(2): 215-224. (in Chinese) |
| [10] | 夏林林, 户晗蕾, 吴开腾. 非线性方程组数值方法的研究进展[J]. 内江师范学院学报, 2013, 28(10): 12-17. XIA L L, HU H L, WU K T. Progress in the research of nonliear equations via numerical methods[J]. Journal of Neijiang Normal University, 2013, 28(10): 12-17. (in Chinese) |
| [11] | 颜庆津. 数值分析[M]. 3版. 北京: 北京航空航天大学出版社, 2006. YAN Q J. Numerical analysis[M]. 3rd ed. Beijing: Beihang University Press, 2006. (in Chinese) |
| [12] | TAIPING H. The transfer matrix impedance coupling method for the eigensolutions of multi-spool rotor systems[J]. Journal of Vibration and Acoustics, 1988, 110(4): 468-472. DOI:10.1115/1.3269552 |
| [13] | LEE J W, LEE J Y. An exact transfer matrix expression for bending vibration analysis of a rotating tapered beam[J]. Applied Mathematical Modelling, 2018, 53: 167-188. DOI:10.1016/j.apm.2017.08.022 |
| [14] | ROSEN A, GUR O. A transfer matrix model of large deformations of curved rods[J]. Computers and Structures, 2018, 87(7-8): 467-484. |
| [15] | 王建国, 逄焕平, 李雪峰. 悬索桥线形分析的悬链线单元法[J]. 应用力学学报, 2008, 25(4): 627-631. WANG J G, PANG H P, LI X F. Catenary cable element method for cable curve or suspension bridges[J]. Chinese Journal of Applied Mechanics, 2008, 25(4): 627-631. (in Chinese) |
| [16] | 何竞飞. 跳绳曲线的数值解及误差估计的间接判定法[J]. 中南工业大学学报, 1995, 26(4): 532-535. HE J F. The numerical solution of rope skipping curve and the indirect evaluation of its error[J]. Journal of Central South University of Technology, 1995, 26(4): 532-535. (in Chinese) |
| [17] | 王庸禄. 对跳绳式捻股成绳机转子弓外形设计及钢丝(股)张力控制的探讨[J]. 金属制品, 1985(6): 8-17. WANG Y L. Discussion on shape design of rotor bow and tension control of steel wire (strand) of rope skipping stranding machine[J]. Metal Products, 1985(6): 8-17. (in Chinese) |
| [18] | 罗雯瑛. 有限到无限: 基于离散静力平衡推导悬链线方程和梁位移方程[J]. 力学与实践, 2020, 42(1): 85-91. LUO W Y. From finiteness to infinity: Deriving the catenary equation and the beam displacement equation based on the discrete static equilibrium model[J]. Mechanics in Engineering, 2020, 42(1): 85-91. (in Chinese) |
