有限元程序对拉索的模拟计算方法探讨

有限元程序对索单元的模拟计算方法探讨邓科，李翠，郝海龙，王灿东，郑德江（中交公路规划设计院有限公司,北京100088）摘要：本文讨论了目前有限元软件对索单元模拟计算的考虑垂度的Ernst公式刚度折减和考虑初应力效应的几何刚度方法，并给出了精确描述索单元力学行为的解析方法。根据目前桥梁工程中主缆和斜拉索常用的力学和尺寸参数进行计算，分析两种模拟计算方法与解析解之间的误差，最后总结了悬索桥和斜拉桥计算时对主缆和斜拉索应采用的合理的模拟计算方法。关键词：索单元，垂度，几何刚度DiscussiononFEMsimulatingmethodofcableDENGKe,LICui,HAOHai-long,WANGCan-dong,ZHENGDe-jiang(CCCCHighwayConsultantsCO.,Ltd.)Abstract:Thispaperdiscussedthesimulationofcableelementsandanalyzedwithfiniteelementsoftware.ThesimulationofcableelementsconcludedErnststiffnessreductionmethodconsideringsagofcable,geometricstiffnessmethodconsideringtheeffectofinitialstressandanaccuratemethodofdescriptionforthecablemechanicalbehaviorwasproposed.Accordingtothecalculationwithmechanicsanddimensionalparameterscommonlyusedinmaincableandstayed-cable,errorsbetweenthesimulationmethodandtheanalyticsolutionwereanalyzed.Finallythereasonablesimulationmethodforthemaincableandstayed-cablewhileanalyzingsuspensionbridgeandcable-stayedbridgewasproposed.Keywords:cable,sag,geometricstiffness1概述索是指仅有轴向拉压刚度而没有抗弯刚度的受力构件，在桥梁工程中大量应用这种构件，如悬索桥中的主缆、吊索，斜拉桥中的斜拉索。在桥梁设计时，目前的有限元程序都是将主缆、吊杆和斜拉索采用桁架单元模拟。一方面，索单元由于其不存在抗弯刚度，在重力作用下存在垂度；另一方面，结构建模都是基于成桥状态模拟结构构形，此时索单元内部存在初应力，这使得索单元的非线性行为非常复杂。在斜拉桥计算时，常采用Ernst公式对斜拉索刚度进行折减；而悬索桥计算时又采用有限位移理论，考虑主缆和吊杆的初应力来提高缆索系统的刚度。同样是索单元，何时考虑刚度折减何时考虑刚度提高，这是一个让桥梁设计工程师倍感困惑的问题。本文通过理论分析和数值验证，给出这一问题的解答。2索单元的模拟计算方法2.1考虑垂度的刚度折减索单元在自重影响下会产生一定的垂度，这一垂度的大小与索力有关，而且垂度与索力呈非线性关系，推导如下。αfTTmqcosαqsinα图1考虑垂度的索单元假设索的重力集度为q，弹性模量为E，索原长S0，变形后端点之间直线长度为S，索与水平面的夹角为α，索内的轴力为T，索中部m的垂度为f，取索的一半为隔离体，对索端部取矩，不难得到：αqSTfcos81=2(1)当垂度较小时，可认为索的形状为抛物线，那么索的形状长度为：2232224cos+=38+=TαlqSSfSP(2)又根据索的轴力为T，那么其形状长度为：01SEATP(3)由(2)和(3)，可以得到：22320024cos=Δ=TαSqSEATSSS(4)根据弹性模量的定义：dTSdASASSddTεdσdEeffΔ=Δ==00ETSαASqE3023212cos+1=(5)近似认为S≈S0，其引起的误差可忽略不计，将(5)整理后，得到ETαASqEEeff322212cos+1=(6)这就是考虑垂度效应后对拉索轴向刚度进行修正的的Ernst公式，从这个公式不难发现，无论索两端轴力为何值，其等效刚度Eeff总小于弹性刚度E，而且随着拉力的减小，其刚度减小更多。2.2考虑初应力的刚度提高众所周知，构件受到轴向压力后，其侧向刚度减小，当轴向压力达到构件的压曲荷载时，其侧向刚度为零，此时构件发生失稳。而构件受到轴向拉力时，其侧向刚度提高，抵抗侧向荷载的能力加强。这个特性是由于结构的几何刚度改变所引起。索的几何刚度矩阵可以简单推导如下。TTTTViVjijΔiΔj图2考虑初拉力的索单元如图2所示，拉索两端受到拉力为T，若拉索两端产生横向位移Δi和Δj，在小位移情形，假设拉索轴力保持T不变，则对i点取矩可得：TlVijj(7)其中为l拉索单元长度。又根据竖向受力平衡，不难知道Vi=Vj，由此可以得到拉索抵抗侧向荷载的平衡方程为：jijilTVV1111(8)其中1111lT便是拉索的几何刚度矩阵。不难看出，拉索初拉力的存在使得其产生横向位移需要施加侧向力，而且随着拉索拉力T的增加，使拉索两端产生横向位移所需的侧向力相应加大，意味着拉索横向刚度的提高。2.3拉索力学行为的解析解事实上，索的受力行为具有高度的非线性，对其力学行为的真实描述不能仅考虑垂度和初拉力，应结合索的基本属性，通过受力分析加以确定。对索单元，引入以下假定：(1)忽略索的抗弯刚度；(2)材料满足胡克定律；(3)索的抗拉刚度变形前后不发生变化。设索弹性模量为E，横截面积为A，延米自重为q，水平长度为l，竖向高度为h，无应力长度为S0，变形后形状长度为P0，左端的水平力和竖向力分别为H、Vi，右端的水平力和竖向力分别为H、Vj，其中由于自重作用方向竖直向下，因此索段处处水平力分量相同，均为H。HViHVjlhpsx,yMO图3索单元解析解设索段上一点M，它与左端点O之间的索段的无应力长度为s，变形后形状长度为p，在笛卡尔坐标系中的坐标为x、y，该处的轴力，水平分力和竖向分力分别为T、H、V。那么根据索段长度定义，可知122dpdydpdx(9)根据受力平衡：TqsVTVdpdyTHdpdxi和(10)根据胡克定律，索段的伸长1EATdsdp(11)将(9)和(10)联立，不难得到22qsVHTi(12)(10)、(11)和(12)联立，可知22qsVHHEAHdsdpdpdxdsdxi(13)分离变量并积分，可得EAHSl02020lnqSVHqSVqHii22lniVHVi(14)类似的22qsVHpsVEApsVdsdpdpdydsdyiii(15)分离变量并积分，可得22202020122iVHqSVHqEASVqShii(16)在l和h确定的前提下，H、Vi和S0三个变量只有一个是独立的，当其中一个确定后，剩下两个变量可以通过(14)和(16)联立求得。dl，dh，dH和dVi的相互关系推导如下。iiiiiiiiiiiidVBdHBdVVHqSVHHdHVHVHVqSVHqSVHqSVHHldVVldHHldl2122202222220220202111111(17)iiiiiiiiiidVBdHBdVVHVqSVHqSVqEASdHVHqSVHHdVVhdHHhdh43222020022202111(18)2.4各种模拟方法的对比研究为了对索的三种模拟计算方法进行对比研究，我们取用工程实践中常采用的平行钢丝拉索或索股的力学参数进行计算研究。索的弹性模量取200GPa；索单元水平投影长度l取10~500m；索单元与水平面的夹角α取10°~80°，拉索上端固定，下端水平约束，竖向自由，计算其竖向发生单位位移时所需力，即为其下端的竖向刚度，如图4所示。HViHVjlhΔhαΔViT,σT,σ图4索单元受力模式若采用桁架单元不考虑垂度效应，也不考虑初拉力提高刚度，那么拉索下端的竖向刚度为：lEAK20sincos(19)对于考虑垂度影响的Ernst公式计算拉索下端的竖向刚度：lAEKeff21sincos(20)对于考虑初应力刚度提高的拉索模拟方法，下端的竖向刚度：lTlEAK322cossincos(21)对于索单元采用解析方法，根据（17）（18），可知：324113BBBBBdhdVKi(22)图5~图8给出了当索内初应力为600MPa、400MPa、200MPa，不同水平投影长度及倾角的K0，K1，K2与K3的比值。01002003004005006000.850.900.951.001.051.101.15Ki/K3l(m)K0/K3α=10°K1/K3α=10°K2/K3α=10°K0/K3α=20°K1/K3α=20°K2/K3α=20°图5竖向刚度比(σ=600MPa,α=10、20°)01002003004005006000.850.900.951.001.051.101.15Ki/K3l(m)K0/K3α=30°K1/K3α=30°K2/K3α=30°K0/K3α=40°K1/K3α=40°K2/K3α=40°图6竖向刚度比(σ=600MPa,α=30、40°)01002003004005006000.850.900.951.001.051.101.15Ki/K3l(m)K0/K3α=50°K1/K3α=50°K2/K3α=50°K0/K3α=60°K1/K3α=60°K2/K3α=60°图7竖向刚度比(σ=600MPa,α=50、60°)01002003004005006000.850.900.951.001.051.101.15Ki/K3l(m)K0/K3α=70°K1/K3α=70°K2/K3α=70°K0/K3α=80°K1/K3α=80°K2/K3α=80°图8竖向刚度比(σ=600MPa,α=70、80°)01002003004005006000.91.01.11.21.31.4Ki/K3l(m)K0/K3α=10°K1/K3α=10°K2/K3α=10°K0/K3α=20°K1/K3α=20°K2/K3α=20°图9竖向刚度比(σ=400MPa,α=10、20°)01002003004005006000.91.01.11.21.31.4Ki/K3l(m)K0/K3α=30°K1/K3α=30°K2/K3α=30°K0/K3α=40°K1/K3α=40°K2/K3α=40°图10竖向刚度比(σ=400MPa,α=30、40°)01002003004005006000.91.01.11.21.31.4Ki/K3l(m)K0/K3α=50°K1/K3α=50°K2/K3α=50°K0/K3α=60°K1/K3α=60°K2/K3α=60°图11竖向刚度比(σ=400MPa,α=50、60°)01002003004005006000.91.01.11.21.31.4Ki/K3l(m)K0/K3α=70°K1/K3α=70°K2/K3α=70°K0/K3α=80°K1/K3α=80°K2/K3α=80°图12竖向刚度比(σ=400MPa,α=70、80°)01002003004005006001.01.52.02.5