悬索桥有限元计算的三节点空间鞍座单元

齐东春， 沈锐利， 刘章军， 谈云志

（1.三峡大学土木与建筑学院，湖北宜昌443002；2.西南交通大学土木工程学院，四川成都610031）

悬索桥有限元计算的三节点空间鞍座单元

齐东春1，2， 沈锐利2， 刘章军1， 谈云志1

（1.三峡大学土木与建筑学院，湖北宜昌443002；2.西南交通大学土木工程学院，四川成都610031）

为了进行空间缆索悬索桥主缆与塔顶鞍座间接触非线性的计算，开发了包括塔顶鞍座及两侧主缆在内的三节点空间鞍座单元.基于空间悬链线理论及主缆与鞍座的几何关系，对单元进行状态求解，得到主缆与鞍座的切点位置及切点索力，根据静力平衡条件计算单元精确的节点力；由增量代替微分，根据切线刚度矩阵的定义计算单元刚度矩阵的元素.空间鞍座单元自动满足主缆与鞍座相切，通过修改一个参数可实现鞍座顶推的计算.计算表明：计算结果与数值解析解结果完全相同，收敛速度较快，在每次状态求解时，迭代次数在12次以内.

悬索桥；空间鞍座单元；有限元法；接触非线性；切线刚度矩阵

主缆与鞍座的接触非线性问题一直是悬索桥计算的难题［1］.为解决这一问题，国内外学者采用解析法和有限元法两种不同方法对平面缆索体系悬索桥的主缆与鞍座接触非线性问题进行了探讨.文献［2-3］采用解析法，根据主缆与塔顶鞍座的几何关系和静力平衡条件建立了多元非线性方程组，求解方程组获得主缆与鞍座在成桥状态下的切点位置.由于解析法在成桥状态和空缆状态等特殊工况的应用尚可，用于一般施工工况则显得十分不方便.为了建立更为通用的计算方法，文献［4］提出采用一种两节点新单元——鞍座-索单元来模拟主缆与鞍座的约束关系，推导的切线刚度矩阵是精确的，计算效率较高，但其主要缺点在于切线刚度矩阵推导过于繁琐，对于具有平弯和竖弯的空间鞍座，由于与鞍座相接触主缆段的无应力长度的计算很复杂，整个单元的力状态量和位移状态量间的状态方程的显式表达式可能很难获得，使得应用受到限制.在此基础上，文献［5］中采用简便的方法建立了一种包含塔顶鞍座与两侧主缆在内的3节点7自由度组合单元，命名为“平面鞍座单元”，避免了冗长的切线刚度矩阵的推导，该方法通用性强，计算精度高，计算效率也较好，并可较方便地推广到空间鞍座.以上文献都是针对平面鞍座开展的研究，对空间鞍座与两侧主缆的接触非线性问题的研究则较少.

1 空间鞍座单元的分析

目前塔顶空间鞍座根据边、中跨主缆的线形一般设计为两种类型：

（1）当边、中跨主缆均为空间三维曲线，此时，鞍座鞍槽曲线在一个斜平面内，相当于平面鞍座绕桥轴线旋转一个角度［6］，如广州猎德桥的塔顶鞍座，称这种鞍座型式为斜平面圆曲线鞍座；

（2）当中跨主缆为空间曲线，边跨主缆仍为平面曲线，此时，主缆由平面曲线过渡到空间曲线，鞍槽曲线也相应为空间三维曲线，如杭州江东大桥的塔顶鞍座，称这种鞍座型式为空间曲线鞍座.

为了满足空间缆索悬索桥有限元计算的需要，本文在文献［5］的基础上，建立了“空间鞍座单元”，与平面鞍座单元相比，空间鞍座单元的建立要复杂得多，主要体现在：

（1）单元自由度由7自由度增加到12自由度，相应的切线刚度矩阵及单元内力计算更为复杂；

（2）空间鞍座鞍槽曲线为具有竖弯和平弯的空间三维曲线，主缆与鞍座脱离点的计算更为复杂；

（3）空间鞍座单元绕塔顶中心存在三向转动，鞍座的最终位置与转动次序相关.

下面以空间曲线鞍座为例来说明空间鞍座单元的建立方法.

1.1 空间鞍座单元的解析模型

如图1所示，空间鞍座单元包含了空间鞍座及两侧的空间主缆，是一个3节点12自由度的组合单元，单元的节点1（I点）取在左侧切点TPL到边跨侧靠近桥塔的第一个索夹之间的主缆上的任意一点，应确保该点不在鞍座上，该节点具有3个平动自由度.节点2取在塔顶中心，除了具有3个平动自由度外，还具有3个转动自由度θX、θY、θZ（分别对应绕总体坐标系X、Y、Z轴旋转），表征了塔顶截面的偏转引起的鞍座绕整体坐标系的转动.节点3（J点）的选取与节点1类同.鞍座的顶点D取成桥状态下过理论IP点的横桥向竖直面与主缆中心的交点，该点即是主缆与鞍座的永不脱离点，也是两侧主缆的分界点.连接节点2和顶点D用两根相互垂直的刚性杆AD、AK，鞍座的顶推可通过改变刚性拉杆AK的长度来实现.

需要指出的，两侧主缆与鞍座的切点TPL、TPR的位置相对鞍座是不断变化的，需要根据当前状态通过迭代计算得到，鞍座单元隐含了切点位置的确定，这是鞍座单元与普通单元的重要区别.

图1 空间鞍座单元模型Fig.1 Schematic diagram of space saddle element model

1.2 单元的已知条件

（1）3节点的坐标（Xi，Yi，Zi），i＝1，2，3，节点2绕整体坐标轴的转角θX，θY，θZ；

（2）节点1至鞍座顶点D间主缆无应力长度为S0L，节点3至鞍座顶点D间主缆无应力长度为S0R，S0L和S0R可根据文献［7］提出的成桥状态下的解析法计算得到；

（3）主缆的弹性模量E，左侧主缆面积AL、自重荷载集度qL，右侧主缆面积AR、自重荷载集度qR；

（4）在鞍座单元坐标系中，鞍座顶点坐标（xD，yD，zD）、平弯圆心坐标（xCH，zCH）及平弯半径RH、竖弯圆心坐标（xCV，yCV）及竖弯半径RV.刚性压杆AD的长度LC，刚性拉杆AK的长度LT.

1.3 单元的计算假定

（1）主缆为小应变理想柔性索，悬空段主缆（与鞍座分离的主缆段）的线形为悬链线［8-10］；

（2）鞍座为刚体，与鞍座相接触的主缆不发生径向变形.鞍座顶点D是两侧主缆的分界点，主缆在该点不会发生相对鞍座的滑动；

（3）在计算切点坐标时，假定主缆与鞍座在单元坐标系xoy投影面内相切，因鞍座构造所限而不保证两者在xoz面也相切；

（4）切点到鞍座顶点D间主缆段的索力分布简化如下：左侧主缆接触段（与鞍座相接触的主缆）索力的水平分力等于左切点TPL处的水平分力，右侧主缆接触段的水平分力等于右切点TPR处的水平分力.

2 空间鞍座单元的求解

空间鞍座单元计算的总体思路如下：

（1）根据主缆与鞍座的几何关系及两节点间主缆无应力长度保持不变，对单元进行当前状态的求解，确定主缆与鞍座的切点坐标和切点索力；

（2）由静力平衡条件及几何关系，计算单元各自由度对应的节点力；

（3）根据切线刚度矩阵的定义，逐一改变各自由度对应的节点位移，计算节点力增量，从而确定单元的切线刚度矩阵.

2.1 坐标转换矩阵

单元坐标系固定在鞍座上，鞍座各特征点在单元坐标系下的坐标保持不变，鞍座与主缆接触段长度、切点坐标及切点分力宜在单元坐标系下求得，然后转换到整体坐标系中.在施工过程中，塔顶截面发生绕整体坐标系X、Y、Z轴的转角位移θX、θY、θZ，相应的单元坐标系绕整体坐标系发生相同的转角位移.单元坐标系的最终位置和转动次序是相关的，严格意义上的大转动问题计算比较复杂［11-12］.事实上，在施工阶段及正常使用阶段悬索桥的塔顶截面的转角位移是很小的，可以假定塔顶截面的转角位移θX、θY、θZ是微小量，这样可以忽略转动次序的影响，同时总转角位移可用增量转角直接叠加计算.假定单元坐标系绕整体坐标系的旋转顺序是绕Y→Z→X依次旋转θY→θZ→θX角度，则可推导出整体坐标系和单元坐标系的转换关系，利用该矩阵可方便进行两坐标系间各物理量的转换.

2.2 单元的状态求解

空间鞍座单元的状态求解就是以两侧主缆的无应力长度不变为条件，根据主缆与空间鞍座的约束关系来确定当前状态的切点位置，当切点位置确定后，单元的各计算参数才完全确定.在单元坐标系中，鞍槽内主缆在坐标面xoy上的投影为一段竖弯圆弧线，在xoz面的投影为平弯圆弧段＋直线段，如图2所示.

图2 空间鞍座鞍槽内主缆的线形描述Fig.2 Alignment description of main cables in space saddle groove

其曲线方程可表述为

现以右侧主缆与鞍座的竖弯切点为例说明计算方法.

（1）根据节点2（即单元坐标系的原点）的转角位移θX，θY，θZ，可得到单元坐标系与整体坐标系间的转换矩阵及矢量的坐标变换关系；

（2）设切点处索力在整体坐标系的三向分力为FTXR、FTYR、FTZR，由矢量的坐标变换关系得到切点索力在单元坐标系中的三向分力fTxR、fTyR、fTzR，由此可计算出切点在单元坐标系中的纵向坐标xTR，代入式（1）可计算出切点的竖向及横向坐标yTR、zTR；

（3）根据接触段主缆的索力分布模式及主缆曲线方程，利用数值积分可计算出鞍座顶点D（xD，yD，zD）到右切点（xTR，yTR，zTR）间主缆的无应力长度S0cR，进而得到悬空段主缆的无应力长度S0hR；

（4）进行坐标转换得到切点在整体坐标系中的坐标（XTR，YTR，zTR），对悬空段主缆，将切点坐标、切点处索力及无应力长度代入索段状态方程中，可计算另一端点（即节点3）的坐标（X3，Y3，Z3）.

通过上述方法建立起了迭代变量FTXR、FTYR、FTZR与目标变量X3、Y3、Z3间的数学关系，可采用修正的影响矩阵法［7］进行求解，从而确定切点坐标及切点索力.

2.3 单元的切线刚度矩阵

在计算单元切线刚度矩阵的第i列元素时，令第i个自由度的位移发生微小增量，此时切点的位置也相应发生变化，需根据当前状态重新计算切点坐标，计算节点力增量，从而得到刚度矩阵的第i列元素.这样确定的单元切线刚度矩阵需通过迭代计算获得，是精确的，但难以得到显式的刚度矩阵.空间鞍座单元的切线刚度矩阵计算流程图如图3所示.

如果在计算刚度矩阵时不考虑切点位置的变化，可得到切线刚度矩阵的显式表达式.其详细推导方法可参考文献［5，13］.虽然按该方法得到的切线刚度矩阵的显式表达式由于没有考虑切点位置的变化存在一定的误差，但在非线性有限元计算中，只有不平衡力计算精度影响计算结果的精度，切线刚度矩阵只影响收敛速度和收敛性［14-15］，因此，本文创建的鞍座单元在切线刚度矩阵的推导上不必追求严格意义上的精确.

图3 空间鞍座单元切线刚度矩阵计算流程图Fig.3 Calculation flow chart of tangent stiffness matrix ofthe space saddle element

2.4 单元的节点力推导

当左、右侧切点位置确定后，根据主缆及鞍座几何关系及静力平衡条件，即可计算出各节点的节点力，节点力的计算图示如图4所示.

图4 空间鞍座单元节点力的计算Fig.4 Calculation schema of nodal forces of the space saddle element

由式（2）计算得到整体坐标系中的节点力为

由于切点坐标及切点索力计算采用解析法，故由此计算的节点力是精确值.另外，单元的节点力是在总体坐标系下推导，无需进行坐标转换.

3 空间鞍座单元的验证

基于上述理论，采用VC＋＋编制了空间鞍座单元计算模块.由于空间鞍座单元是本文首次开发的，难以用其他方法和有限元软件进行模拟和验证，为验证其正确性，采用作者自编的解析法程序SPCC［7］加以验证.以江东大桥为例，为简化计算取对称结构，如图5所示，计算其空缆状态和成桥状态，比较两种状态下SPCC程序和空间鞍座单元有限元程序的计算结果.表1为SPCC程序和空间鞍座单元有限元程序两种计算方法求得的切点坐标及切点索力.表2为边跨跨中、中跨1／4及中跨跨中点的坐标值.

由表1、2可以看出：两种计算方法无论是线形还是主缆索力几乎没有差异，说明数值解析解和有限元解是完全相同的.计算表明，本文开发的空间鞍座单元计算精度高，收敛速度较快，在每次状态求解时，迭代次数一般在12次以内.该单元适用于空间缆索悬索桥的计算，提高了悬索桥的计算精度，减少了计算工作量，具有较好的实用性.

图5 计算模型Fig.5 Calculation model

表1 主缆与鞍座的切点坐标及切点索力的比较Tab.1 Comparison of coordinate and cable force of tangent point between main cable and saddle

表2 主缆节点位移的比较Tab.2 Displacement comparison of main cable nodesm

4 结束语

空间鞍座单元是在平面鞍座单元的基础上发展而来的，用以解决空间鞍座与两侧主缆的接触非线性计算.针对斜平面圆曲线鞍座和空间曲线鞍座两种不同类型的鞍座，建立了空间鞍座单元的求解方法.在单元内部首先通过数值解析法迭代出主缆与鞍座的切点坐标，确定单元的当前状态，然后根据鞍座与主缆的几何约束关系及静力平衡条件，推导当前状态精确的单元节点力.采用工程中常用的以增量代替微分的方法计算切线刚度矩阵，根据切线刚度矩阵的定义计算单元刚度矩阵的元素.计算表明，该单元能准确计算主缆与鞍座的切点坐标和切点索力，计算精度高，收敛速度较快，能像常规单元用于空间缆索悬索桥的有限元计算中，通过修改一个参数即可方便实现鞍座顶推的计算.该单元的引入使计算模型更加接近实际结构，提高了空间缆索悬索桥的计算精度.

［1］ 潘永仁，杜国华，范立础.悬索桥恒载结构几何形状及内力的精细计算［J］.中国公路学报，2000，13（4）：33-36.PAN Yongren，DU Guohua，FAN Lichu.A fine calculation ofthegeometryandinternalforceof suspension bridge under dead load［J］.China Journal of Highway and Transport，2000，13（4）：33-36.

［2］ 唐茂林，强士中，宁晓骏.悬索桥索鞍理论位置的精确计算［J］.工程力学，1999（增刊）：339-344.

［3］ 李传习，王雷，刘光栋.悬索桥索鞍位置的分离计算法［J］.中国公路学报，2005，18（1）：63-68.LI Chuanxi，WANG Lei，LIU Guangdong.Separate calculationmethodonsuspensionbridgesaddle's position［J］.China Joural of Highway and Transport，2005，18（1）：63-68.

［4］ 罗喜恒，肖汝诚，项海帆.用于悬索桥非线性分析的鞍座—索单元［J］.土木工程学报，2005，38（6）：47-53.LUO Xiheng，XIAO Rucheng，XIANG Haifan.Saddlecable element for nonlinear analysis of suspension bridge［J］.China Civil Engineering Journal，2005，38（6）：47-53.

［5］ 齐东春，沈锐利，陈卫国.悬索桥结构分析中鞍座单元的研究及应用［J］.桥梁建设，2011（1）：35-38.QI Dongchun，SHEN Ruili，CHEN Weiguo.Research and application of saddle element in structural analysis ofsuspensionbridge［J］.BridgeConstruction，2011（1）：35-38.

［6］ 罗喜恒，肖汝诚，项海帆.空间缆索悬索桥的主缆线形分析［J］.同济大学学报：自然科学版，2004，32（10）：1394-1354.LUO Xiheng，XIAO Rucheng，XIANG Haifan.Cable shapeanalysisofsuspensionbridgewithspatial cables［J］.JournalofTongjiUniversity：Natrual Science，2004，32（10）：1394-1354.

［7］ 齐东春.大跨径悬索桥主缆精细化计算研究［D］.成都：西南交通大学，2012.

［8］ PEYROT A H，GOULOIS A M.Analysis of cable structures［J］.Computers＆Structures，1979，10（5）：805-813.

［9］ KIM H K，LEE M J，CHANG S P.Non-linear shapefindinganalysisofaself-anchoredsuspension bridge［J］.Engineer Structures，2002，24：1547-1559.

［10］ 沈锐利.悬索桥主缆系统设计及架设计算方法研究［J］.土木工程学报，1996，29（2）：3-9.SHEN Ruili.Calculation methods for design and erection of cable curve of suspension bridge［J］.China Civil Engineering Journal，1996，29（2）：3-9.

［11］ 潘永仁.悬索桥结构非线性分析理论与方法［M］.北京：人民交通出版社，2004：35-112.

［12］ ZUPAN D，SAJE M.Rotational invariants in finite elementformulationofthree-dimensionalbeam theories［J］.Computer and Structures，2004，82：2027-2040.

［13］ 齐东春，沈锐利，唐茂林.悬索桥的锚碇-锚跨单元研究及应用［J］.公路交通科技，2011，28（2）：82-92.QI Dongchun，SHEN Ruili，TANG Maolin.Study and applicationofanchorage-anchorspanelementfor suspensionbridges［J］.JournalofHighwayand TransportationResearchandDevelopment，2011，28（2）：82-92.

［14］ 唐茂林.大跨度悬索桥空间几何非线性分析与软件开发［D］.成都：西南交通大学，2003.

［15］ 丁顺泉，陈艾荣，项海帆.空间杆系结构实用几何非线性分析［J］.力学季刊，2001，22（3）：300-306.DING Shunquan，CHEN Airong，XIANG Haifan.Geometrically nonlinear practical analysis of spatial frames［J］.ChineseQuarterlyMechanics，2001，22（3）：300-306.

（中文编辑：秦 瑜 英文编辑：兰俊思）

3-Node Sptial Saddle Element for Finite Element Calculation of Suspension Bridge

QI Dongchun1，2， SHEN Ruili2， LIU Zhangjun1， TAN Yunzhi1
（1.College of Civil Engineering＆Architecture，China Three Gorges University，Yichang 443002，China；2.School of Civil Engineering，Southwest Jiaotong University，Chengdu 610031，China）

In order to solve the computational problem of nonlinear contact between main cable and saddle on tower top，a new 3-node spatial saddle element，including the saddle on tower top and both sides of the main cable，was produced.Based on the spatial catenary theory and the geometric relationship between main cable and saddle，the location of the tangent points and the cable force at the tangent points were obtained by solving the element state determination problem with known conditions.The accurate nodal force of the element was derived according to the static equilibrium.The elements of the tangent stiffness matrix were calculated based on its definition by replacing the differential with the increment.The new element could automatically satisfy the condition that the main cable is always tangent to the saddle，and thus the saddle jacking could be conveniently realized by modifying a parameter.Calculation shows that the new element has high calculation accuracy and convergence rate.The calculation results are the same with the numerical analytical solutions.The number of iterations is generally less than 12 in each element state solution.

suspension bridge；space saddle element；finite element method；contact nonlinearity；tangent stiffness matrix

U448.25

A

0258-2724（2014）06-0942-06

10.3969／j.issn.0258-2724.2014.06.002

2014-03-18

国家自然科学基金资助项目（51178396）；湖北省自然科学基金资助项目（2014CFB331）

齐东春（1978－），男，博士，研究方向为大跨度桥梁结构非线性理论及精细化分析，E-mail：qidongchun＠163.com

齐东春，沈锐利，刘章军，等.悬索桥有限元计算的三节点空间鞍座单元［J］.西南交通大学学报，2014，49（6）：942-947.