0 引言

悬索桥的超强跨越能力是业界共识，成为跨越江河与山壑首要桥型选择^[1-2]。悬索桥建设高速发展趋势下，悬索桥体系设计阶段的技术参数备受桥梁建设者重视，其中主缆找形及索鞍设计位置是悬索桥设计阶段关键参数，关系到悬索桥设计初始构形^[3-4]。

悬索桥作为1种高次超静定柔性结构体系，设计阶段成桥状态主缆几何非线性效应显著^[5-11]。目前成桥状态主缆找形计算方法主要有分段悬链线法^[12-14]、有限元法^[15-16]，分段悬链线法计算过程复杂及收敛性受迭代初始值及迭代增量影响，对于迭代初始值及迭代增量取值不合适，会导致迭代计算不收敛，有限元法难以准确确定成桥状态设计索鞍设计位置，同时存在计算发散等问题。为解决主缆找形计算问题，罗凌峰等^[17]考虑主缆与索鞍摩擦效应，基于闭合条件推导出中、边跨缆索系统找形的目标函数，采用Nelder-Mead算法求解目标函数最优解, 实现主缆找形计算；王邵锐等^[18]提出基于弯矩平衡推导出末点和中间点标高与主缆左端水平和竖向分力的解析关系式，得到综合考虑末点和中间点标高误差影响的修正迭代计算公式；邓小康等^[19]提出基于悬链线理论的切线表达式，根据跨中点斜率最低为计算起点，建立主缆迭代计算方程，通过最小二乘法实现主缆找形计算。这3种算法的计算格式及计算理论较为复杂，设计人员难以掌握及运用其计算理论。国内研究人员展开关于成桥状态索鞍设计位置深入研究，李传习等^[20]从力学关系和几何关系出发，推导得到了复合圆曲线索鞍位置计算的解析方程组，通过牛顿拉斐森迭代算法求解解析方程组，实现索鞍设计位置分离算法，有效解决索鞍设计位置计算，但没有将散索鞍弯矩平衡纳入解析方程组中，散索鞍设计位置处主缆切点坐标与实际主缆找形下主缆在散索鞍切点坐标计算存在不闭合等问题。邓小康等^[21]基于悬链线斜率表达公式，根据主缆和索鞍几何关系建立非线性方程，采用二分法求解一元非线性方程，得出复合圆曲线散索鞍的设计位置，该法有效降低方程组的个数，从而提高计算效率，但同样存在不闭合问题。

以上研究推导悬索桥设计阶段理论发展，为克服国内外关于悬索桥主缆找形及索鞍设计位置计算研究的局限性，本研究基于悬链线方程理论及闭合条件建立主缆找形非线性方程组及索鞍设计位置非线性方程组，同时考虑散索鞍力矩平衡关系建立解析方程组，通过求解非线性方程组，实现闭合同步计算求解, 有效解决主缆找形迭代计算易发散、索鞍设计位置计算公式复杂及散索鞍设计位置与主缆找形不闭合等问题。

1 悬链线方程基本理论

缆索单元力学简图见图 1。其中q为悬索结构的自重荷载集度，s为自重荷载作用下的长度，T为索端张力，Tdy/ds为索端张力竖向分量，Tdx/ds为索端张力水平分量，取索段微分段根据力学平衡条件可得^[22]：

(1)

(2)

如图 1所示，索的水平分量表达式为H=Tdx/ds，则将T=Hds/dx代入式(2)可得：

(3)

式(3)为索段基本平衡微分方程，索段张力水平分量沿索长保持不变。该微分方程是建立在索段自重作用后的变形微分方程，表示索微元的二阶导数与竖向荷载计算成正比例关系。qds/dx可看作荷载集度q沿水平方向等效均布荷载q_y。

(4)

将式(4)代入式(3)可得到式(5)：

(5)

表达式(5)可进一步变换为：

(6)

将等式(6)左边进一步化简如式(7)所示：

(7)

将等式(7)进行积分可得式(8)。

(8)

式中，a与b为积分常数；c=-H/q。

如图 2所示，对于跨中缆索而言，以缆索的首端点B为坐标原点建立坐标轴，索段高程计算公式为^[7]：

(9)

式中，h_i为第i号索段高差，对于节点处作用节点力P_i，根据力学平衡条件推导以下计算公式：

(10)

式中，α, β分别为索段始端、末端缆索切线与水平方向夹角。

根据式(10)表达式，悬链线系数a_i+1计算公式如式(11)所示^[23]：

(11)

索段无应力长度S_i计算公式见式(12)：

(12)

2 非线性方程组的建立 2.1 中跨主缆非线性方程组

如图 2所示，节点B与C为中跨已知主缆IP点，N为中跨垂点。假定B处主缆张力为H₀及V₀，根据悬链线方程理论，则B处悬链线系数表达式如式(13)所示，各节点悬链线系数计算公式如式(14)所示。以中跨主缆各个节点为原点建立直角坐标系，根据式(9)求得各个节点间高差如式(15)所示，成桥状态主缆找形跨度l_i及节点力P_i为已知量，根据几何闭合条件，主缆BN及BC竖向高度表达式分别如式(16)~式(17)所示。

(13)

(14)

(15)

(16)

(17)

2.2 边跨及锚跨主缆找形非线性方程组

如图 3所示，A为左边跨已知主缆在锚跨处IP点，假定B处主缆张力为H₀及V_b0，根据悬链线方程理论，则节点B处悬链线系数如式(18)所示，各节点悬链线系数计算公式如式(19)所示。以左边跨主缆各个节点为原点建立直角坐标系，根据式(9)求得各个节点间高差如式(20)所示，成桥状态边跨主缆找形跨度l_bi及节点力Q_i为已知量，根据几何闭合条件，BA竖向高度表达式如式(21)所示。

(18)

(19)

(20)

(21)

根据右边跨竖向高度闭合条件建立解析方程式(22)。

(22)

式中，a_yi为左边跨悬链线系数，根据左锚跨主缆力与左边跨主缆力相等原则建立式(23)：

(23)

式中，H_m为锚跨索股水平缆力；V_m为锚跨索股竖向缆力；a_m为锚跨索股悬链线系数。

根据几何闭合条件，AM竖向高度表达式如式(24)所示：

(24)

右锚跨竖向高度相容解析方程如式(25)所示：

(25)

2.3 主索鞍非线性方程组建立

如图 4所示，主缆在主索鞍左、右切点到B点水平距离分别为l_c与l_e，切点坐标(x_c, y_c)及(x_e, y_e)，索鞍圆心C′的坐标为(x_f, y_f)。根据悬链线理论计算公式，切点坐标计算表达公式如式(26)~式(29)所示：

(26)

悬链线系数a_b0表达式见式(27)：

(27)

同理可得：

(28)

悬链线系数a₀可表示为：

(29)

假定主索鞍半径为R₁，根据几何关系建立切点到圆心距离平方和表达式如式(30)~ (31)所示：

(30)

(31)

根据主索鞍切点上的缆力与切点和圆心的连线相互垂直几何关系，建立方程组如式(32)~式(33)所示：

(32)

(33)

同理可建立右侧主索鞍非线性方程组。

2.4 散索鞍及锚跨非线性方程组建立

复合型圆曲线散索鞍示意如图 5所示。为便于问题描述，计算简图取复合圆曲线的段数为3段进行讨论。鞍座由半径为R₃，R₄，R₅的3段圆曲线组成，各段圆曲线对应的圆心角分别为θ₃，θ₄，θ₅，索鞍首段圆曲线起始点与该段圆心连线水平角为θ₆，该角度为索鞍转动定位角度。在直角坐标系下：假定主缆在散索鞍处A点的坐标为(x_g，y_g)，第1段圆弧S₁圆心及第3段圆弧S₃圆心坐标分别为(x_n1，y_n1), (x_n3，y_n3)。坐标(x_n1, x_n1)与坐标(x_n3, x_n3)间计算表达式如式(34)~式(35)所示：

(34)

(35)

散索鞍IP点为点A，主缆与主索鞍左切点与右切点到A点水平距离分别为l_a，l_b。对于计算过程，本研究以右切点位于第1段圆弧，左切点位于第3段圆弧为例进行说明。根据悬链线计算公式可得切点坐标(x_a, y_a)及(x_b, y_b)计算表达公式如式(36)~(38)所示。根据悬链线理论计算公式，切点坐标计算表达公式如式(36)~式(39)所示。

(36)

锚跨及边跨主缆在散索鞍切点处缆力和散索鞍自重与索鞍转动中心取矩为0：

(37)

式中，e₁为散索鞍重心到转动中心J的水平距离；e₂，e₃分别为边跨切点处缆力到转动中心竖向水平距离；e₄，e₅分别为锚跨切点处缆力到转动中心竖向水平距离，具体计算公式可参考文献[24]，本研究不在此赘述; V_b(n+1)为节点A处主缆竖向分力。根据边跨主缆闭合条件可得式(38)：

(38)

悬链线系数a_b计算表达式为：

(39)

根据几何关系建立切点到圆心距离平方和表达式如式(40)~式(41)所示：

(40)

(41)

根据散索鞍切点处缆力与切点和圆心的连线相互垂直几何关系，建立方程组如式(42)~ 式(43)所示：

(42)

(43)

同理可建立右侧散索鞍非线性方程组。

2.5 解析方程组求解

联立非线形方程组进行主缆找形及索鞍设计位置计算，26个独立方程组含有一共26个未知数，则方程组存在有解情况。根据MATLAB软件编制非线性方程组求解程序，采用GRG方法进行求解^[23]，依次可求解各跨主缆力学要素，索鞍设计位置及主缆在索鞍处切点坐标。实现悬索桥主缆找形及索鞍位置同步快速求解，本研究算法只需要求解非线性方程组，与分段悬链线理论比较，本研究算法计算收敛性对于迭代初始值及迭代增量没有依耐性。

3 工程运用及算法比较 3.1 工程背景

如图 6所示，某悬索桥为本研究工程研究背景，桥梁构造为主跨1 038 m单跨地锚式，主梁采用钢箱梁，矢跨比为1/9，全桥共63对吊索。

3.2 有限元模型

如图 7所示，采用Midas Civil建立的悬索桥整体模型，全桥离散为1 330个节点，1 280个单元。主缆及吊杆弹性模量为E=1.98×10⁵ MPa。运用软件自带基于节线法的建模助手进行初步平衡态计算，根据悬索桥精确分析模块进行精细化找形计算分析。

3.3 分段悬链线理论

分段悬链线理论广泛应用于主缆线形找形计算，为1种经典理论计算方法。以中跨为例对主缆找形方法进行说明，基于抛物线理论计算假定的初始段迭代初始值及基于力学平衡条件计算节点处竖向力计算公式见式(44)~式(46)，依次代入各索段分段悬链线隐式方程，将求得各索段高度依次求和，根据公式(47)计算跨中垂度收敛误差值及缆索始末端高度误差值，将计算误差值与收敛精度比较，不满足收敛条件，则根据牛顿拉斐森切线刚度矩阵求取下一次初始值迭代增量，进一步修正初始迭代值，直到计算差值满足收敛精度要求。因此在使用分段悬链线计算理论进行主缆找形计算时，初始迭代值及迭代增量的选取决定收敛速度及收敛性。

(44)

(45)

(46)

收敛误差计算公式如式(47)所示，h₁与h₂分别表示为缆索垂度、缆索始末端高差。

(47)

3.4 算法验证

应用GRG算法对本研究计算方法提出的非线性方程组进行求解，可快速求解非线性方程组，根据非线性方程组的解与悬链线基本理论可实现主缆找形计算及索鞍设计位置计算。运用有限元法进行主缆找形，需要建立复杂的软件计算模型及进行多次迭代才能达到设计精度要求，同时存在无法直接求解索鞍设计位置及计算发散等问题。分段悬链线计算理论收敛性与初始迭代值及迭代增量相关，索鞍设计位置计算需在主缆找形计算完成后进行。为验证本研究提出的计算方法的可靠性及适用性，进行3种计算方法对比分析。吊杆无应力长度对比如图 8所示。

本研究计算值与经典理论计算值及有限元值较吻合，其中中跨主缆线形最大差值分别为0.8 mm及0.9 mm。在成桥状态主缆无应力计算结果对比中，本研究计算值与经典理论计算值及有限元值较为吻合，其中中跨主缆无应力长度误差最大，误差值分别为2.4 cm及0.8 cm，误差率在规范允许范围内，误差率满足规范范围要求。

在悬索桥吊杆无应力下料长度计算结果对比中，本研究计算值与经典理论计算值及有限元值较为吻合，本研究计算值偏短，其中1号吊杆误差值最大，误差控制值在1.0 mm以内，吊索无应力长度计算值满足工程精度要求。

成桥状态主缆在索鞍切点坐标与索鞍设计位置计算结果对比见表 1，本研究计算值与经典理论计算值误差控制在2.0 cm以内。由于本研究考虑散索鞍力矩平衡的完备性闭合计算，主缆找形理论IP点对应与索鞍设计位置的理论IP点为同一位置，实现计算的闭合，而经典理论计算值则采用常规计算理论，本研究算法为成桥状态找形与索鞍设计位置闭合同步算法，力学概念简洁清晰，同时求解计算过程中只需要求解多维非线形方程组。运用MATLAB软件基于GRG法进行编制求解非线性方程组程序，可快速实现迭代收敛计算。

4 结论

(1) 本研究提出基于悬链线方程理论及几何相容条件建立主缆找形解析方程组，基于索鞍几何关系及变形协调关系建立索鞍设计位置解析方程组，通过求解联立的26维非线性方程组，基于非线性方程组的解与悬链线基本理论实现主缆找形计算及索鞍设计位置同步计算。本研究算法在提高计算效率的同时满足精度要求。

(2) 本研究算法纳入散索鞍弯矩平衡方程，使得主缆设计位置下的主缆在散索鞍切点坐标与实际主缆找形计算坐标一致，由此计算得到的索鞍IP点位置对应与主缆找形时的IP点位置一致，实现闭合计算。

(3) 本研究成桥状态主缆线形计算值与经典理论计算值及有限元值较吻合，中跨主缆线形误差值在1.0 mm以内；本研究成桥状态主缆无应力计算值与经典理论计算值及有限元值较为吻合，其中中跨主缆无应力长度误差最大，相对误差值分别为2.4 cm与0.8 cm，误差率在规范允许范围内；本研究索鞍切点坐标与索鞍设计位置计算值与经典理论计算值误差控制在2 cm以内，存在差异性原因：本研究考虑散索鞍力矩平衡的完备性闭合计算，而经典理论计算值则采用常规计算理论，验证本研究算法的高效性及可靠性。本研究算法可推广到超大跨度悬索桥设计计算理论。