改进力密度法在悬索桥主缆找形中的应用

2022-03-18詹海鹏

同济大学学报（自然科学版） 2022年3期

刘超，詹海鹏，朱怡

（1.同济大学土木工程学院，上海 200092；2.无锡市政设计研究院有限公司，江苏无锡 214072）

悬索结构的工程应用历史悠久，是最常见的大跨度结构形式，应用范围覆盖到桥梁工程、房屋建筑等。悬索桥是最典型的悬索结构，通过其主缆的轴向力来抵抗外荷载作用，可以充分利用钢丝等材料的拉伸强度。当主缆采用高强度材料的时候，能够大幅度减轻自重比例，使得悬索桥能够实现较大的跨度。悬索桥主缆不同于一般的刚性结构，在无应力状态下的结构初始形状是不确定的，在使用状态下的主缆轴力提供几何刚度。因此，悬索桥主缆存在显著的非线性特征，这也使得主缆找形问题成为其结构分析的重点和难点。针对悬索桥结构，Ochsendorf等［1］指出，Euler提出的抛物线理论，即主缆在沿跨均布荷载下的形状为抛物线，其水平分力为恒定值，忽略了主缆重力沿着主缆曲线分布的特征，仅适用于跨度较小的桥梁。随着悬索桥跨度的不断增大，抛物线形状的主缆受力与实际受力不相符，误差较大，不能满足设计要求。因此，悬链线模型［2］、分段悬链线法［3］、整体力学分析法［4］等解析方法被先后提出，并广泛应用至今。近些年，有限元法［5-7］得到快速发展，成为分析工具之一。但非线性有限元方法在悬索结构找形过程中存在两个问题：一是多迭代子步带来的误差积累，二是收敛的不稳定性［8］。此外，宋旭明等［9］提出需要仔细考虑主梁的压缩对主缆无应力长度产生的影响。Kim等［10］考虑加劲梁和主塔的轴向压缩，基于弹性悬链线单元和牛顿迭代法，以无应力长度为未知数，提出了TUCD（target configuration under dead load）及其改进算法，但是该算法对初始值要求较高。我国悬索桥规范规定的缆形计算方法即为分段悬链线法［11］，此方法迭代参数较少，求解速度快，是目前应用最普遍的方法。但该方法对迭代初值敏感，迭代过程容易发散，且需要推导复杂的悬链线方程。

力密度法［12-17］由Schek于1974年提出，主要用于索网找形分析，因其简易的计算理念而得到快速发展。动力松弛法［18］自从1965年提出以来不断得到发展，其本质是利用特殊条件下的动力方法来解决静力问题。二者均广泛应用在索网和膜结构设计中。力密度法优势在于能够利用线性手段求解非线性问题，计算更加简便。索网结构找形与悬索桥主缆找形原理相似，本质上都是寻找适应边界条件的平衡状态。本文基于传统的力密度法，结合主缆无应力长度的概念，增加迭代求解过程的控制条件，提出了适合悬索桥主缆找形的改进力密度法。利用静力平衡和几何协调关系，采用力密度作为初始参数，形成力与节点坐标的线性关系，建立整体平衡方程，直接求解主缆节点坐标。该方法将非线性问题转化为线性方程的求解，步骤简便，收敛稳定，计算精度满足工程需求。

1 力密度法基本原理

力密度法［12］基本参数有结构的拓扑关系、边界约束条件、荷载以及平衡条件。悬索结构的拓扑关系规定了节点单元的排列规则和序列，采用单元-节点的关联作为其拓扑矩阵C t，定义为

式中：i，j为节点号，i，j=1～nt，节点分为两部分，nf表示固定节点，nd表示自由节点，二者的关系可以表示为nt=nd+nf；e为单元编号，e=1～m；拓扑矩阵C t的行数对应着单元号，列数对应着节点号，即为m行和nt列。根据固定节点和自由节点可将C t分块为C和Cf，即C t=[CCf]。

设结构中所有节点的三维坐标为x t、y t、z t，均为nt维矢量，坐标矢量同样可分为固定节点与自由节点相邻节点的坐标差矢量u，v，w为

设单元长度和轴力分别为m维矢量l和t，与之相对应的m维对角矩阵设为L和T，对应3个坐标方向的外荷载矢量分别为p x、p y、p z。结构处于平衡状态下，每个节点力的矢量和为零，因此可以得到以下平衡方程：

式中：U、V、W为相邻节点的坐标差矢量；u、v、w对应的m维对角矩阵。

令q=L-1t，其物理含义为单元轴力和单元长度的比值，即力密度的定义，可以得到

利用关系式

将式（2）和式（5）带入式（4）得

设D=CTQC，即为C的广义高斯变换，同理设Df=CTQCf，可得

式（7）为力密度法求解的基本方程，并且是一组能够求解自由节点坐标的线性方程。对于一个给定的荷载体系和固定节点，对应每一组力密度，就可以求得该组力密度下结构的唯一平衡形状。由于对角矩阵Q元素为力密度值，对于纯受拉的主缆体系没有交叉节点，且设置受拉为正，即满足每个单元力密度大于0，则D是正定矩阵。

对于给定的结构拓扑矩阵，在已知的荷载和边界条件下，不同的力密度对应着不同的平衡状态。力密度的变化与平衡状态的变化是相对应的，同时对于线性方程的求解效率较高。因此，力密度是悬索结构适合的描述参数。

2 改进力密度法

工程实践已经证明，悬索桥的主缆在单元长度合适的情况下，采用只受轴向力的杆单元模拟主缆索结构具有较高的精度，并且单个杆单元的垂度效应可以忽略。因此，主缆找形计算假定如下：

（1）主缆为完全柔性，只受轴向拉力作用，忽略主缆抗弯刚度的影响。

（2）单个杆单元的自重均布荷载根据有限元基本理论等分到单元两端节点形成等效节点力。

（3）杆单元张拉过程中截面的面积保持不变。

（4）主缆单元应力在弹性范围内满足胡克定律。

（5）主缆锚固点位置、吊杆位置、跨中垂点以及吊杆力为设计已知量。

传统的力密度法只是简单给定力密度，尚无法满足悬索桥主缆找形特定的限制条件，此时找形问题就转化为需满足附加控制条件的非线性问题，需采用特定的方式进行迭代求解。针对主缆体系，引入无应力长度的概念，基于不同的设计需求建立满足限定条件的迭代格式，从而求解出主缆找形问题中的成桥缆形和空缆缆形。

2.1 成桥缆形计算

在平衡方程中，两侧锚点和塔顶IP点为已知固定节点，其余节点可作为自由节点，吊杆力作为外荷载的一部分也是已知量。在成桥缆形求解问题中，主缆无应力长度和自由节点坐标为未知量，迭代过程中保持更新，其余已知条件保持不变。

由于主缆缆形未知，其自重也是未知的，需要设置初始缆形提供初始单元长度。初始单元轴力可设为一个较大值，易于收敛。依据二者数值得到初始力密度。初始缆形根据设计需求可设置为锚点、塔顶IP点与垂点之间的相连直线，以吊杆位置为界将主缆分成若干个单元和节点。初始缆形中每一段的主缆长度已知，可以用主缆的截面积乘上节段长度再乘上重度表示节段主缆的重力，便可得到成桥缆形计算的平衡方程。根据式（8）求解得到新的自由节点坐标。以节点z坐标为例，其平衡方程如式（10）所示。

式中：zf为固定节点z坐标向量，始终保持不变；qe为各单元力密度值，e是单元编号；z n为第n阶段待求

解的自由节点z坐标向量为平衡矩阵，由第n-1阶段的单元力密度qen-1确定；p z(l en-1)为结构z方向外荷载，由于吊杆力已知，所以外荷载由第n-1阶段的单元长度l en-1确定。

由于假设的直线缆形不可能保持平衡，因此需要进入下一迭代。以本阶段的缆形作为初始状态，重新代入求解平衡方程。主缆重力是对应上一阶段缆形的节段重力，在计算得到新的缆形之后，由于主缆重力发生变化，结构体系将再次不平衡，需要反复迭代计算，直到相邻节段的缆形坐标变化量达到收敛要求。在主缆无应力长度和坐标均未知的情况下，锚点和塔顶的边界条件不足以确定主缆的成桥线型。理论上，可以有无穷组力密度满足方程的要求，即有无穷组成桥缆型解，因此需要根据已知设计条件（垂点坐标和吊杆纵向坐标）来获得跟设计参数匹配的成桥缆形。在每次迭代过程中设定这些自由节点的设计坐标值，进入下一阶段的迭代计算，最终平衡方程求解获得的自由节点坐标会逐渐收敛至满足设计参数的缆形。

式中：x ni、yi n、zni为第n阶段i节点的三维坐标；x ni-1、为第n-1阶段i节点的三维坐标。

成桥缆形求解具体流程如图2所示。

图2 成桥缆形求解流程图Fig.2 Flowchart of cable shape finding in finished stage

2.2 空缆缆形求解

空缆缆形与成桥缆形的求解相似，初始缆形和初始轴力与成桥缆形分析中设置相同。空缆状态下，吊杆力为零，只受主缆重力作用，并且空缆缆形的垂度未知。由于从空缆到成桥状态主缆的锚点不变，因此全桥主缆的无应力长度不变。在成桥缆形确定后，每一个主缆单元轴力为已知，成桥状态的单元无应力长度可表示为

式中：s0e为单元无应力长度；se为成桥状态单元有应力长度；te为主缆单元轴力；A为主缆截面积；E为主缆弹性模量。

利用主缆在成桥状态和空缆状态时的无应力长度相等这一控制条件，即可求解空缆缆形。在迭代过程中保持固定节点不变，利用已知的单元无应力长度来控制迭代过程中单元轴力和长度。以节点z坐标为例，其平衡方程如式（13）所示。

式中：pz在外荷载迭代过程中保持不变；qen-1可以看作是单元无应力长度的函数；l e表示迭代过程中单元有应力长度由已知的无应力长度求得；te表示单元轴力。其余符号参考式（10）的说明，收敛要求见式（11）。

由于主缆在成桥和空缆状态下内力不同，因此需要设置索鞍预偏量，消除因受力不同导致的主缆对主塔的水平力。张拉过程中缆形会随着吊杆的逐步张拉，回到成桥状态的线形。对于索鞍的预偏量，可以在缆形计算的基础上，先假设预偏量为Ap，假设偏向边跨，预偏量调整值为Δ=Ap2。在此假设下，分别计算空缆状态和成桥状态下的主缆无应力长度。若前者大于后者，调整预偏量为Ap=Ap-Δ；若前者小于后者，调整预偏量为Ap=Ap+Δ，并设置下一阶段预偏量调整值为Δ=Δ2。如此反复计算，直到空缆缆形和成桥缆形的无应力长度之差满足精度要求。

3 工程应用

3.1 平面缆形悬索桥实例

如图4所示，京杭大运河悬索桥位于徐州市金山桥经济开发区境内，为双索面地锚式悬索结构。全桥长245.0 m，桥宽5.8 m，主跨为130.0 m的钢桁架梁。主塔采用钢桁架结构。主缆为1根预制平行钢丝索股，垂跨比1/10，垂度13 m，两主缆横桥向中心间距为5.8 m。全桥共有50根吊索，吊索间距为5 m，吊索索力在139～176 kN之间。该桥已于2020年11月施工完成。

利用本文方法计算该桥的成桥和空缆缆形。跨京杭大运河悬索桥的成桥状态如图5所示。在桥梁主缆架设和张拉期间对全桥缆形进测量，空缆和成桥实测缆形如图6所示。本桥主缆弹性模量为1.99×105MPa，截面积为0.01 m2。

图5 主缆成桥状态Fig.5 Main cable in finished satge

图6 空缆和成桥实测坐标Fig.6 Measured coordinates of main cable in unloaded stage and finished satge

根据流程图2、图3，基于Matlab软件对空缆缆形与成桥缆形进行编程计算。同时，通过Midas Civil软件建立有限元模型，如图7所示。比较两种算法之间的缆形误差，并得到本文程序所求得的空缆和成桥缆型与实测值的误差，如图8和图9所示。

图3 空缆缆形求解流程图Fig.3 Flowchart of cable shape finding in unloaded stage

图 4京杭大运河悬索桥总体布置图（单位：m）Fig.4 Layout of Beijing-Hangzhou Grand Canal suspension bridge(unit:m)

图7 有限元模型Fig.7 Finite element structural model

图8 空缆缆形比较Fig.8 Comparison of cable shapes in unloaded stage

从图8和图9可以看出，论文方法与有限元方法计算的空缆和成桥缆形误差最大分别为2.0 cm和2.3 cm。与实测值相比，空缆的缆形误差最大为5.6 cm。分析可知，产生此项误差的主要原因是空缆时荷载较小，导致主缆存在局部扭曲，因此测量存在一定误差。对于成桥缆形，论文方法计算的缆型理论值与实测缆形的最大误差为5.8 cm，主要原因是吊杆力理论值与实际值存在一定误差。总体而言，论文方法计算的缆形误差相对较小。

图9 成桥缆形比较Fig.9 Comparison of cable shapes in finished stage

3.2 空间缆形悬索桥实例

为验证论文方法对于空间缆形的适用性，选取韩国永宗大桥为工程背景，具体数据见参考文献［19］。永宗大桥为自锚式空间悬索桥，桥梁总体布置和具体参数见图10。吊索下锚点处桥面竖曲线半径为14 472 m，下锚点之间的横向距离为31.92 m；端吊索力竖向分量为3 905.30 kN，其余边跨吊索力竖向分量为3 142.86 kN，中跨吊索力竖向分量均为3 088.96 kN；边跨锚点和塔顶IP点高差为66.51 m。该桥主缆弹性模量为1.99×105MPa，截面积为0.135 5 m2。

图1 悬索桥示意图Fig.1 Schematic diagram of suspension bridge

图10 永宗桥总体布置图（单位：m）Fig.10 General layout of Yongzong Bridge(unit:m)

空间缆形相较于平面缆形增加了Y方向的自由度，同样选取锚点到IP点之间的连线作为找形程序的初始缆形，见图11中的虚线线形，边界条件和之前所述保持不变。已知吊索力竖向分量，但是吊索力横向分量未知，可以利用已知的吊索锚点与主缆节点坐标之间的位置关系，通过空间力矢量的分解，得到迭代过程中作用于主缆的吊索力横向分量。自锚式悬索桥梁体在巨大的索力作用下，会产生较大的压缩量，本文通过有限元计算位移，改变找形边界坐标，可以消除主梁压缩量对于缆形误差的影响。

通过本文程序求解的空间缆形如图11所示，由空缆的直线经过迭代收敛至成桥的缆形。图12和图13分别给出了详细的平面坐标、立面坐标以及线形迭代变化过程。