第十一章桥梁结构几何非线性

第十一章第十一章 桥梁结构几何非线性桥梁结构几何非线性 计算理论计算理论 n11.1 11.1 概概 述述n11.211.2大跨度桥梁几何非线性分析的有限元方法大跨度桥梁几何非线性分析的有限元方法 n11.311.3桥梁结构分析常用单元的切线刚度矩阵桥梁结构分析常用单元的切线刚度矩阵 n11.411.4桥梁结构几何非线性分析特殊问题的讨论桥梁结构几何非线性分析特殊问题的讨论 n11.511.5非线性方程组的求解非线性方程组的求解 n11.611.6小结小结 11.1 11.1 概概 述述 Oden说过“我们生活在一个非线性世界里”。早在十九世纪未，科学家就发现，固体力学的经典线性理论在许多情况下并不适用，于是开始了对非线性力学问题的研究。二十世纪中，科学家奠定了非线性力学的理论基础。但由于计算繁复，许多非线性微分方程的边值问题无法求解，用解析法解决非线性工程问题仍显得无能为力。直到二十世纪六十年代末，有限元法与计算机相结合，才使工程中的非线性问题逐步得以解决。固体力学中有三组基本方程，即本构方程、几何运动方程和平衡方程。经典线性理论基于三个基本假定，即材料的应力、应变关系满足广义虎克定律；位移是微小的；约束是理想约束，这些假定使得三组基本方程成为线性。只要研究对象不能满足线性问题基本假定中任何一个时，就转化为各种非线性问题。表11.1给出了非线性问题的分类及基本特点。非线性问题 定 义 特点 桥梁工程中 的典型问题 材料非线性 由材料的非线性应力、应变关系引起基本控制方程的非线性问题。材料不满足虎克定律。砼徐变、收缩和弹塑性问题。几何非线性 放弃小位移假设，从几何上严格分析单元体的尺寸、形状变化，得到非线性的几何运动方程，由此造成基本控制方程的非线性问题。几何运动方程为非线性。平衡方程建立在结构变形后的位置上，结构刚度除了与材料及初始构形有关外，与受载后的应力、位移整体也有关。柔性结构的恒载状态确定问题，柔性结构的恒、活载计算问题；桥梁结构的稳定分析问题。接触问题 不满足理想约束假定而引起的边界约束方程的非线性问题。受力后的边界条件在求解前未知。悬索桥主缆与鞍座的接触状态；支架上预应力梁张拉后的部分落架现象。非线性问题的分类及基本特点 表11.1 由表11.1可知，几何非线性理论将平衡方程建立在结构变形后位置上。事实上，任何结构的平衡只有在其变形后的位置上满足，才是真实意义上平衡的。线性理论之所以能得以广泛应用，只是因为一般结构的受力状态不因变形而发生明显改变。而有些问题则不然，以图11.1所示结构为例，按线性理论求解就无法找到平衡位置，按几何非线性分析方法处理，在P力作用下，B点产生竖向位移，当位移达到一定值时，AB、BC两杆件中轴力的竖向分力与P平衡，即为B点位移的解。可见，受力状态因变形而发生明显改变时，就必须用几何非线性方法进行分析。图11.1受集中力的二力杆 几何非线性分析理论在桥梁工程中的发展，起因于桥跨的长大化和柔性结构的应用。早在1888年，Melan就在悬索桥结构分析中提出了几何非线性的挠度理论，在考虑主缆拉力二阶影响的基础上将悬索桥的平衡方程建立在变形后的位置上，但忽略了吊杆伸长、结构水平位移及加劲梁剪切变形的影响。挠度理论从1908年开始应用于纽约的Manhattan大桥设计，大大节省了工程造价，充分显示了它的优越性。此后的数十年中，挠度理论为悬索桥和大跨径拱桥的发展作出了巨大贡献。但是，挠度理论平衡微分方程的求解仍是十分复杂的。Timoshenko于1928年提出了三角级数解，Godard通过忽略后期荷载对结构刚度的影响提出了线性挠度理论，我国李国豪教授于1941年提出了用于悬索桥分析的等代梁法，将挠度理论中的非线性项等代于偏心受拉梁的弯矩减小系数，揭示了悬索桥受力的本质。现代桥梁工程的发展和跨径的增大，使得结构越来越柔，越来越复杂，结构分析中梁柱效应、索的伸长、结构水平位移及后期荷载的二阶影响变得不可忽略，对各种复杂结构，建立挠度理论的平衡微分方程及其求解也越来越困难。为此，工程界渴望出现更精确、方便的理论和方法。六十年代初，M.J.Turner、Brotton等开始发表求解结构大位移、初应力问题的研究成果，Poskitti、Saffan等也在此领域里作出了贡献。这些理论方法都可归入几何非线性力学的有限位移理论。在建立以杆系结构有限位移理论为基础的大跨径桥梁结构几何非线性分析平衡方程时，一般考虑了三方面因素的几何非线性效应：1)单元初内力对单元刚度矩阵的影响。一般情况下是指单元轴力对弯曲刚度的影响，有时也考虑弯矩对轴向刚度的影响。常通过引入稳定函数或单元几何刚度矩阵的方法来考虑。在大跨径桥梁结构分析中遇到的初应力（或初应变）问题，就是指结构现有内力引起的结构刚度变化对本期荷载响应的影响问题。2)大位移对建立结构平衡方程的影响。在这个问题上，目前流行的T.L列式法和U.L列式法各有不同的处理方法。前者将参考座标选在未变形的结构上，通过引入大位移单元刚度矩阵来考虑大位移问题；后者将参考座标选在变形后的位置上，让节点座标跟随结构一起变化，从而使平衡方程直接建立在变形后的位置上。3)用杆单元近似模拟索类构件，由索垂度引起的单元刚度变化。简单的处理方法是引入Ernst公式，通过等效模量法来近似修正垂度效应。也可以通过导出索元切线刚度矩阵，用索单元直接描述索类构件。今天，有限位移理论一般用有限元方法通过计算机程序来求解。因此，程序的编制也应看成是非线性计算理论和方法不可分割的一部分。七十年代未，国外相继推出了ADINA,ANSYS,MARC,NASTRAN,ASKA,NON-SAP 等结构分析综合程序。它们可用于桥梁结构的部分非线性计算和局部应力分析。但由于缺少许多必备的功能，这些程序无法完整地完成桥梁设计计算。国内学者根据规范要求和实际情况，开发了桥梁通用程序，如同济大学桥梁系开发的BAP系统、交通部公规院开发的QJS系统，有的已具备非线性计算功能。随着计算机技术的发展，桥梁结构分析软件也得到了迅速发展，经历了从单一化结构分析到将数据管理、用户接口、图形加工与管理、面向对象的软件设计和可视化技术融为一体的发展过程。本章结合程序计算流程，讨论桥梁结构有限位移分析的理论与方法。11.2 11.2 大跨度桥梁几何非线性大跨度桥梁几何非线性 分析的有限元方法分析的有限元方法本节以杆系结构为对象，讨论拉格朗日列式的大跨度桥梁几何非线性有限元方法。11.2.1 11.2.1 变形体的运动描述变形体的运动描述 任何变形体在空间都占据一定的区域，构成一定的形状，这种几何形状简称为构形，物体在问题求解开始时的构形称为初始构形，在任一瞬时的构形称为现时构形，物体位移的改变叫运动。在下面讨论中，字母的左上标表示构形所处时刻。图11.2中有一物体，在t=时，物体有初始构形 。物体中一点 的坐标为 ;在t=时，物体运动有构形 ，点 运动到 ，在 时，物体运动有构形 ，点 运动到 。变形体及其上质点的运动状态，随不同坐标选取有以下几种描述方法：0t0A0P000123(,)xxxntnA0PnP1nntttt1nA0P1nP图11.2变形体的运动 (1)物质描述：独立变量为 和 ，即给出任意时刻物体中各质点的位置。这种描述在连续介质力学与有限元中很少使用。(2)参照描述：独立变量为任意选择的参照构形中质点的当前坐标与时刻t。这种描述法称为拉格朗日法(Lagrangian Formulation)。当选择t=0时的构形为参照构形时，称总体拉格朗日描述(T.L.Formulation)。(3)相关描述：以 为独立变量。参照构形与时间有关，取 为非线性增量求解时增量步的开始时刻，则称为更新的拉格朗日描述(U.L Formulation)。(4)空间描述：独立变量是质点当前位置 与时间 。这种描述称为欧拉描述。在欧拉描述中，有限元网络在空间中是固定的，材料流过这些网络。这种描述适用于流体及定常状态。0(1,2,3)ix i 0tntnt1nx1nt11.2.2 11.2.2 总体拉格朗日列式法总体拉格朗日列式法(Total Lagrangian(Total Lagrangian Formulation)Formulation)在整个分析过程中，以t=0时的构形作为参考，且参考位形保持不变，这种列式称为总体拉格朗日列式。对于任意应力应变关系与几何运动方程，杆系单元的平衡方程可由虚功原理推导得到：BdVfTV0BdVfTV0(11-1)B 式中：单元的应力向量；f 单元杆端力向量；V 单元体积分域，对T.L列式V是变形前的单元体积域；B 应变矩阵，是单元应变与节点位移的关系矩阵。即：杆端位移向量。(11-2)在有限位移情况下B是位移 的函数。后面将看到，B矩阵可分解为与杆端位移无关的部分 和与杆端位移有关的部分 两部分，即：0BLB0 LBBB11-3 直接按式(111)建立单元刚度方程并建立结构有限元列式，称为全量列式法。在几何非线性分析中，按全量列式法得到的单元刚度阵和结构刚度阵往往是非对称的，对求解不利。因此多采用增量列式法。将式(111)写成微分形式：dBdVdfTV()0 (114)dBdVBddVdfVTVT 或(115)根据式(113)和(115)等式左边第一项可写成：dBdVdBdVkdVTVLT 0(116)另一方面，单元的应力、应变增量关系可表示成：dDd(117)式中：D 弹性矩阵。当材料满足线弹性时：00()D（118）式中：0 单元的初应变向量；0单元的初应力向量。由(112)，(113)代入(117)得：d=D(B0+BL)d （119）于是，式(115)左边第二项可表示为：()BddVBDBdVBDBdVBDBdVBDBdV dTVTVTLVLTVLTLV00000000 kBDBdVTV000 kBDBdVBDBdVLTLVLTV BDBdVLTLV (1110)记 (1111)(1112)则式(115)最后可表达为：)(00000fddkdkkkTL (1113)式(1113)就是增量形式T.L列式的单元平衡方程。式中 是三个刚度阵之和，称为单元切线刚度矩阵，它表示荷载增量与位移增量之间的关系，也可理解为单元在特定应力、变形下的瞬时刚度。单元节点位移无关，是单元弹性刚度矩阵，称为单元初位移刚度矩阵或单元大位移刚度矩阵，是由大位移引起的结构刚度变化，是d 的函数。称为初应力刚度矩阵，它表示初应力对结构刚度的影响，当应力为压应力时，单元切线刚度减小，反之单元切线刚度增加。0TK0OK0LK0K将各单元切线刚度方程按节点力平衡条件组集成结构增量刚度方程，即有：0PddKT (1114)式中：为结构切线刚度矩阵，可以由单元切线刚度矩阵按常规方法进行组集形成；为荷载增量，由于荷载增量一般取为有限值而不可能取成微分形式，结构在求得的位移状态下，抗力与总外荷载之间有一差量，即失衡力，结构必须产生相对位移以改变结构的抗力来消除这个失衡力。在计算中，一般通过迭代法来求解。0TK d P11.2.3 11.2.3 更新的拉格朗日列式法更新的拉格朗日列式法(U.L(U.L列式列式)在建立t+t时刻物体平衡方程时，如果我们选择的参照构形不是未变形状态t=0时的构形，而是最后一个已知平衡状态，即以本增量步起始时的t时刻构形为参照构形，这种列式法称为更新的拉格朗日列式法(U.L列式)。由于采用了U.L列式，平衡方程式(115)中的积分须在t时刻单元体积内进行，且 的积分式是 的一阶或二阶小量，此点在下一节将作进一步说明。因此，代表 的积分式可以略去。这是U.L列式与T.L列式的一个重要区别。最后增量形式的U.L列式平衡方程可写成:tlk tok Lk)(0Pddkktt )(0Pddkktt (1115)11.2.4 T.L11.2.4 T.L列式与列式与U.LU.L列式的异同及适用范围列式的异同及适用范围 T.L列式与U.L列式是不同学派用不同的简化方程及理论导出的不同方法，但是，它们在相同的荷载增量步内其线性化的切线刚度矩阵应该相同，这一点已得到多个实际例题的证明。T.L列式与U.L列式的不同点由表112给出。T.L列式与U.L列式的不同点 比较内容T.L列式U.L列式注意点计算单刚的积分域在初始构形的体积域内进行 在变形后的t时刻体积域内进行U.L列式必须保留各节点座标值精度保留了刚度阵中所有线性与非线性项忽略了高阶非线性项U.L列式的荷载增量不能过大 单刚组集成总刚用初始时刻各单元结构总体座标系中的方向余弦形成转换阵，计算过程中不变用变形后t时刻单元在结构总体座标中的方向余弦形成转换阵，计算过程中不断改变U.L列式中组集载向量也必须注意方向余弦的改变本构关系的处理在大应变时，非线性本构关系不易引入比较容易引入大应变非线性本构关系U.L方法更适用于砼徐变分析表11-2 从理论上讲，这两种方法都可以用于各种几何非线性分析，但通过表112的对比可以发现T.L列式适用于大位移、中等转角和小应变的几何非线性问题，而U.L列式除了适应于上述问题外，还适用于非线性大应变分析、弹塑性、徐变分析。可以追踪变形过程的应力变化。目前，国内使用的桥梁非线性分析程序，一般都采用U.L列式方法。11.3 11.3 桥梁结构分析常用单元的切桥梁结构分析常用单元的切线刚度矩阵线刚度矩阵 由前面讨论可知，T.L列式下单元切线刚度阵可分为三个部分，即弹性刚度阵 、初位移刚度阵 和几何刚度阵 ，而U.L列式下单元切线刚度阵只有 和 两部分，本节进一步讨论桥梁结构分析中常用单元切线刚度阵的具体表达形式。00 k0 Lk0 k0 tk tk11.3.1 平面桁架单元的切线刚度矩阵平面桁架单元的切线刚度矩阵图11.3所示的桁架单元ij，杆长为l，截面积为A，在外荷载作用下，i、j端发生了位移。将桁架单元轴向应变由线性形式改写成大位移形式：xdudxdvdx122()(1116)dudx122()dvdx式中：x的线性表达式；-大位移情况下x由竖向位移引出的附加项。取桁架单元的形函数矩阵为：Nxlxlxlxl100010 (1117)则单元上任意点的位移可写成：uvN (1118)式中：单元的节点位移向量。式(1118)代入式(1116)整理得：xllll 10101201012()(1119)dlldlllldx 101001010101 和()BBdL0 (1120)其中：Bll01010 (1121)考虑到：01011llluvji()(1122)BllL0101 (1123)在T.L列式下，单元的局部座标建立在变形前的初始状态，将(1121)代入式(1111)得：将式(1121)和(1123)代入式(1112)计算。得：kEAl01010000010100000 (1124)kEAlL00002222 (1125)由式(1123)得：jjiiijTLvuvulvvlllB10100000101000001)(110102 (1126)式(1126)代入式(116)，并注意到：d VN lV (1127)得：dlNdVBdVTL1010000010100000 (1128)即：最后，T.L列式下桁架单元的切线刚度矩阵为：1010000010100000lNk (1129)kkkkTL0 (1130)在U.L列式下，单元的局部坐标系建立在t时刻，即变形后的位置上，其单元切线刚度矩阵为：tTttkkk0 (1131)下面我们来考察TL列式与UL列式两种切线刚度阵的等价性。图11.4了给出桁架单元发生变形前、后的位置，其中 坐标系为变形前的单元局部坐标系，为变位后的单元局部坐标系，两者夹角为，xy为整体坐标系。用U.L列式法建立单元切线平衡方程时，局部坐标系是建立在 上的。00 x yttx yttx y现将弹性刚度矩阵 从 座标系变换到变形前的位置 上，考虑到角较小，则sin，cos1，其变换矩阵为：ttx y0 tk00 x y100100001001cossin00sincos0000cossin00sincosT (1132)于是：0TkTktTT222221111lEA222200000000010100000101lEAlEALkk00 (1133)由此说明，T.L列式和U.L列式的单元切线刚度矩阵具有等价性。T.L列式下单元初位移矩阵的实质是让单元在变形后的位置上发挥其作用，以满足平衡方程必须建立在结构变形后位置上这一重要条件，而U.L列式则通过节点坐标的不断迁移来实现这一目标。11.3.2 平面柔索单元的切线刚度矩阵平面柔索单元的切线刚度矩阵 桥梁结构分析中，常遇到柔索构件。斜拉桥的斜拉索、悬索桥的主缆、施工中用的缆风和扣索等都可抽象成柔索。柔索的特点是抗弯刚度小，索的自重对结构平衡影响不可忽略，用拉压杆模拟柔索会引起误差。因此，有必要建立柔索单元的刚度方程。为讨论方便，且不影响计算精度，作如下假定：(1)柔索仅能承受张力而不承受弯曲内力（抗弯刚度为0）；(2)柔索仅受索端集中力和沿索长均匀分布的荷载作用，荷载合力效应为q；(3)柔索材料符合虎克定律；(4)局部座标系取在柔索荷载合力平面内。考察图11.5中所示柔索，无应力索长为 ，索的荷载集度q向下为正。0S图11.5 柔索单元图11.5中参变量之间有如下关系：0220220,TFFTTFFTqSFFHFFyyhxxljjiiijjiyxjyxiyyxxijij (1134)易导得各力素与几何变量之间的关系如下：qTTTTEAqqTTqSFEAShFTFTqEASFlijijijyyjyixijii)(21)2()ln1(22000 (1135)对式(1135)取全微分有:iiiiiiiiyyxxyyxxFFhFFhhFFlFFll (1136)于是，i端力和位移的增量关系可写成：hCAlCBFhCAlCBFiiyx1122 (1137)端位移和力的增量关系可写成：hCAlCBFhCAlCBFjjyx1122 (1138)式中:luuhvvjiji (1139)122102121)(1)11()(11BABACTFTFqEASBTTqFBATFTFqFAiyjyijxiyjyxijiiji (1140)将式(1137)(1139)合并整理后写成矩阵形式：FkeTe (1141)式中：()FFFFFexyxjyjTii (1142)44343324232214131211kSymmkkkkkkkkkkeT (1143)其中：TjjiivuvuCAkkkCAkkkkCBkkk),(14424222342314122334411 (1144)式(1141)即为柔索单元切线刚度方程，在索端平衡力已知的情况下，可直接计算柔索切线刚度矩阵。在索端平衡力未知的情况下，首先按单根柔索计算索端力，求解时先初估一个 和 ，若(1135)式自然满足，初估值即为真实值。否则，设估算值使式(1135)产生的误差为：xiFyiFhhelleyx00 (1145)式中：0l由估算的 和 计算出的柔索水平投影长 xiFyiF0h由估算的 和 计算出的柔索垂直投影长 xiFyiFlh下一次计算希望通过和的修正，使误差趋于零，即：00yxehel (1146)将式(1136)代入(1146)易得：CAeBeFCBeAeFyxyxyxii1122 (1147)yixiFF,用 修正 xiF和 yiF，再按照图11.6所示的流程迭代，就可以在已知l h q E A S,0的情况下，求出所有索端力。图11-6求索端力的计算流程 最后，计算切线刚度阵流程见图11.7。用直杆代替柔索计算是常用的近似方法，柔索的垂度效应可用 Ernst公式对弹性模量进行修正，这种方法在小位移、高应力水平下，具有较高精度。但如果索工作在大位移状态或应力水平不高的情况下，就会出现很大的误差。因此，采用近似方法计算时应慎重。11.3.3 平面梁单元的切线刚度矩阵平面梁单元的切线刚度矩阵 在桥梁结构分析中还经常会用到平面梁单元、空间梁单元等。与113.1中平面桁架单元一样，通过给出合适的几何运动方程。由式(114)出发作相应的推导，就能得到各自的单元切线刚度矩阵，限于篇幅，这里不作推导地给出平面刚架单元U.L列式下的单元切线刚度矩阵。平面受弯曲的梁单元，横截面积为A，长度为l，惯矩为I，轴向位移用一次多项式插值，梁竖向位移和转角用三次多项式插值，可导出平面梁单元的弹性刚度阵为：2222223461200260461206120000IlSymmIlIAlIlIlIlIlIIlIAlAllEket (1148)平面梁单元的初应力刚度阵为:lSymmlllllllNkt433600030433603360000003022 (1149)N是t时刻梁单内的轴向力，受拉为正。11.3.4 11.3.4 算例算例下面给出两个算例的计算结果，供读者参考。例11.1 如图11.8所示,杆AB连线长度L=100m,与水平线夹角=45,具有初始力T0=71.9208kN,杆的截面刚度系数EA=1000kN,分成10个等长的杆单元,在中点C受一竖直力P=52.8374kN的作用。数值与理论计算结果基本一致，见表113。图 11.8 有初轴力的直杆受竖直力作用方 法位 移内 力u uv vu uv v本文解6.559977.8388574.9276112.9214理论解6.559977.8388574.9276112.9214有初轴力直杆分析结果有初轴力直杆分析结果 (单位单位:M-kN:M-kN)表表11113 3 例11.2 用柔索单元计算图11.9所示的柔索索端力与索端位置的关系。图11.9 自重作用下索端力与索端位置的关系 图11.9给出了柔索索端力与索端位置的关系，图中I,J表示索的两上端点。11.4 11.4 桥梁结构几何非线性分析特殊桥梁结构几何非线性分析特殊问题的讨论问题的讨论11.4.111.4.1稳定函数与几何刚度矩阵稳定函数与几何刚度矩阵如图11.10所示压杆的M、Q和位移为正，其挠曲平衡微分方程为：图 11.10 压杆变形、内力)()1(xyNlxMlxMyEIji (1150)方程的解为：xlxNMlxNMlxClxCyji)1(cossin21 (1151)其中：EINl (1152)引入边界条件：lyylxx|,0|0 得：NMCNMMCiij21,sincos (1153)yMNxlxli(sin()sin)11MNxlxlxi(sinsin)(1154)于是：lxjjxiidxdydxdy|)(|)(0 (1155)()(JijJiicMsMEIlsMcMEIl (1156)其中：)1sin(1)1(122stgc (1157)如果轴力为拉力，则：)sinh1(1)1(122stghc (1158)c、s为轴力影响下，杆端单位力矩引起的杆端角变形，c为力矩作用端的角变形，s为另一端的角变形。最后，可导出有初轴力的杆单元刚度方程：ijjiiIJjjiiQQSCSClEIQSCSClEIMSCSClEIM)22()()()(22 (1159)或：jjiiJjiiuuCSymmlSClSCSlSCClSClSClSClSClEIMQMQ)(22)22(22222222 (1160)其中,2222,scsSsccC 是轴力N的函数，称为稳定函数，其值随的变化而变化。以稳定函数表达的刚度系数包含了轴力对弯曲刚度的影响，相当于前面切线刚度阵中弹性刚度系数与几何刚度系数之和。图11.11给出了稳定函数与单元切线刚度系数随变化的情况，说明1时，二者基本一样，当 3时，二者区别也不大。图11.11 稳定函数与单元切线刚度系数随的变化曲线 事实上，将C展开成级数形式有：.504019301(442C (1161)而)301(4)301(4)(2233033EINlkkEIl (1162)对比表明，几何刚度阵系数就是稳定函数忽略高阶项的轴力影响系数。从以上讨论可以看到：当 3时，随着的增大，几何 刚度矩阵的误差也增大。但由于与l 成正比，有限元分析中，只要减小单元长度，就可避免使用几何刚度阵产生的这种误差。11.4.2 11.4.2 弯矩对轴向刚度的影响弯矩对轴向刚度的影响 杆单元的弯曲将引起杆件计算长度(杆件两端节点的距离)的改变，从而影响杆件的轴向刚度。在杆件微段上，弯曲引起的杆件轴线计算长度的改变量 d为：dxdxdydxdsd22 (1163)22)(211)(1/yydxd (1164)dxyl02)(21 (1165)在外力作用下，杆件总的缩短量为：dxyEANll02)(21)/(1166)因此 lEAN (1167)式中：弯矩引起的轴向刚度修正系数。由于:)sin/sin()1sin)/1(sin(lxlxNMlxlxNMyji (1168)NlMMlxlNMlxlNMMyjiijisincossincos (1169)式中：由式(1152)定义，式(1169)代入式(1165)积分，由式(1166)、(1167)可得)(4/(1123MlNEA (1170)其中：当 P0(压杆)时)cot1(csc2)(2)csc)(cot()(222jijijiMMMMMMM (1171)当P0(拉杆)时)coth1(csc2)(2)csc)(coth()(2222jijijiMMhMMhMMM (1172)11.4.311.4.3活载的几何非线性分析活载的几何非线性分析 为了讨论方便,将非线性状态下荷载最不利加载区域称为影响区。活载几何非线性分析，会遇到如下问题：线性叠加原理失效，无法再用传统的影响线加载法进行活载分析。确定影响区本身是一个非线性问题。仅用恒载初始状态计算活载，会带来影响区范围改变和不正确载位引起的误差。单位强迫变位产生的等效力很大，用机动法求解影响区将破坏指定状态结构影响区的真实形状。考虑到结构在确定初始状态下，影响线函数的大小，仍代表单位荷载作用对关心截面计算参数的影响。因此，可以按如下方法计算非线性活载最不利响应：以结构恒载状态为初态，计算影响区函数，用相应的最不利活载作为试探，求出第一次近似，将前一次试探活载与恒载共同作用时的状态代替初态，重新计算影响区和最不利荷载。当本次活载效应与上次活载效应的误差落在某一允许范围时，计算收敛。这样,每叠代一次，都加强了活载对计算状态的影响。这一计算可归结为如下步骤：1)将结构恒载受力状态作为求解影响区的初始状态，计算出初始影响函数。2)用动态规划加载法，找出最不利加载位置，并作好记录。3)以恒载受力状态为计算初态，将活载按最不利载位一次性作用于结构，分析恒、活载共同作用下的结构受力状态和关心截面力学量。4)将恒、活载共同作用下的结构状态作为求解下一步影响区函数新的初态，重复1)3)的计算。经过数次迭代计算就可以得到活载作用下关心力学量的最值。求解活载影响区可用机动法,但单位强迫变位必须取用一个很小的数(如 ），这样，可以保证确定的影响区不失真，使动态规划法找到的载位即为相应内力状态下的最不利载位。51011.4.4 11.4.4 几何非线性调值计算几何非线性调值计算 桥梁结构在设计和施工计算中，常常遇到如下问题：设计阶段：需要通过调整部分构件的内力或支点位移，来优化整个结构体系的受力状态。比如，调整斜拉索力来优化斜张桥结构恒载内力状态；调整预应力索张拉力来改善预应力体系关心截面应力状态；调整连续梁支座来改善全桥正负弯矩分配等。施工阶段：由于结构模型噪声、施工精度及其它随机因素的影响，使结构状态偏离控制目标或出现超应力状态。需要通过调整少量构件的内力或改变部分约束节点的位移来减少已建结构与控制目标的偏差，并使后期结构的施工状态最大限度地逼近目标。竣工后的结构也会遇到与前相仿的问题。如通过拱脚位移顶推来调整拱圈内力;通过斜拉索再张拉，来调整使用多年的结构由于徐变、收缩引起的线形和内力改变等。上述问题可由桥梁结构关心截面内力、位移、应力调值计算来解决，一般情况下，要使关心截面中n个独立参量调整为指定值，就必须改变施调截面中的n个独立参量。随着新技术和施工新工艺的应用与发展，调值计算问题在工程中越来越多，调值计算成了大型结构计算和结构优化设计必不可少的计算方法 当几何非线性表现突出时，基于线性叠加原理的调值计算方法无法直接用于非线性结构的计算，下面讨论计入几何非线性影响的调值计算求解策略。1)影响矩阵及其形成方法为了方便讨论，首先作如下定义：受调向量：结构物中关心截面上m个独立元素所组成的列向量。这些元素一般由构件中的截面内力、应力或位移组成。它们在调值过程中接受调整，以期达到某种期望状态。受调向量记为：12(,)TmDd dd (1173)施调向量：结构物中指定可实施调整以改变受调向量的个独立元素（）所组成的列向量，记为：12(,)TlXx xx (1174)施调元素多为杆件内力或支座变位。影响向量：施调向量中第j个元素 发生单位变化，引起受调向量D的变化向量，记为：jx12(,)TjjjmjAaaa (1175)影响矩阵：个施调向量分别发生单位变化，引起的个影响向量依次排列形成的矩阵，记为：.AAAAaaaaaaaaalmmml1211121121222112 (1176)在影响矩阵中，元素可能是内力、应力、位移等力学量中的一个，影响矩阵是这些力学量混合组成的。从理论上讲，只要将单位施调向量逐一加到结构上，分别求出相应的影响向量，便能形成结构的影响矩阵。但当受调向量为内力时，由于内力无法直接加在结构上，一般是通过先将相应构件从结构中“断开”，并在断开处施以一对大小相等方向相反的单位力来进行计算的。显然，这样做破坏了原有的结构形式，用有限元方法计算，则每计算一个影响向量，就要形成和分解一次结构刚度阵，是很不经济的。为了减少形成影响矩阵的计算量，可先将内力元素的影响向量用相应位置和方向上杆件的单位强迫变形影响向量来代替，这样就不必将构件断开，而可在同一力学模型上进行影响向量的计算。用有限元法计算影响矩阵，可归结为如下步骤来进行：a）形成调值计算阶段结构总刚，并作 分解；b）对施调元循环；c）令第j号施调元调值量为1(如果是内力，则用单位强迫变形代替)形成相应的结构荷载列阵；d）回代求相应的节点位移。e）对受调元i循环，计算相应的受调元素 ，如果 是位移，就从求解得到的位移向量中取出该调值点对应方向的位移作为j列影响向量的第i个元素 ，如果 是内力，就由位移向量求出相应调值截面对应方向的内力值作为 。如果 是应力，则由内力求出相应点的应力值作为 。如果调值点又是受调点本身，在 中还需扣除单位强迫变位引起的单元的内力变化。重复bd各步，就可形成所有影响向量，从而形成影响矩阵A。TLDLidididijaidijaijaija2)等变量的调值计算 设结构中n个关心截面上期望的内力、应力、位移组成的向量为E，关心截面中现有相应向量为 ，活载引起的相应向量为 ，调值计算就是通过改变n个施调元的力学量，使结构状态在关心截面处达到E。此时，结构受调向量为：dF LF DEFFdL (1177)受调向量与施调向量的元素相等，这种调值计算称为等变量的调值计算。当结构满足线性叠加时，有：AXD (1178)式(1178)可唯一求得施调向量X，向量x表示:要使关心截面力学量达到E，必须使施调变量产生X的相应变化。3)几何非线性调值计算几何非线性对结构受力状态影响较大时，用式(1178)求得的施调向量X作用于结构物，并不能使被调向量达到期望值。这就是几何非线性的调值计算问题。几何非线性调值计算可以通过迭代法完成，具体步骤如下:a）以调值前结构状态为初态，计算广义影响矩阵A。b）求解方程(1178)，得到初始施调向量X0。c）将施调向量X0同时作用于结构，进行几何非线性分析，求结构实际状态。d）提取相应的初始被调向量D0，并求出受调差值向量：DDD0 e）以第三步中得到的结构状态替代调值前初态，重复a）d）的计算，直到D的范数|D|小于指定误差为止。(1179)以上计算收敛较快，一般迭代cd次就能满足工程精度要求。11.5 11.5 非线性方程组的求解非线性方程组的求解 11.5.1 11.5.1 求解方法概述求解方法概述 用有限元法进行结构非线性分析，其控制方程最终是一组非线性代数方程。非线性代数方程组的求解方法很多，其选择往往与物理问题的性质、特点、非线性程度、对计算结果的要求以及计算机的容量、计算速度等因素有关。这就要求研究者对非线性问题的求解过程以及程序设计有较全面的了解。以下介绍几种常用的求解方法。1)直接求解法 直接求解法是基于全量列式的求解过程，应用最多的是直接迭代法，由虚功原理建立的非线性有限元平衡方程为：)(PK (1180)0当设定位移向量的初值 后，改进的近似解可由下式得到：100KP (1181)整个迭代过程可用下式表示：PKnn11 (1182)当迭代结果满足预定的收敛准则时，就得到了所要求的节点位移向量。0图11.12(a)为取0 =0时单自由度问题的迭代过程取得收敛的示意图。直接迭代法应用简单，运算速度一般也较快，可应用于具有轻微非线性的问题。这一求解过程的成功与否很大程度上取决于对初值位移0 的正确估计。图11.12(b)表示的是直接迭代法迭代过程发散时的情形。为改善收敛性和收敛速度，可以采用将荷载分成若干级的做法。图 11.12 直接迭代法收敛和发散过程 2)增量法 增量形式的有限元列式方法具有一个共同的特点：将整个荷载变形过程划分为一连串增量段，每一增量段中结构的荷载反应被近似地线性化。简单增量法将每一级增量荷载下直接求得的状态变量视作结构平衡状态，计算相应的切线刚度阵，进而作下一级荷载计算，并不 断累加其位移增量。图11.13描述了简单增量法的求解过程。图 11.13 简单增量法的求解过程 几何非线性问题的有限元分析最初多采用简单增量法进行，虽然这种求解方法对每一级荷载作用时的计算速度较快，但由于每一级荷载作用前结构并未精确地到达平衡位置，所求得的解答会随着增量过程的继续而越来越偏离真实的荷载变形过程。为了保证计算精度，常常将增量区间划分得相当小。此外，为了评价解的精度，一般要对同一问题在进一步细分增量区间后再次求解，通过两次解的比较判定是否收敛。这样就需要消耗大量的计算时间。作为对这一方法的改进，可将不平衡力作为一种修正荷载并入下一级荷载增量。这就是有一阶自校正的增量法。一阶自校正增量法求解过程的示意图如图11.14所示。一阶自校正增量法具有较高的求解速度，同时也比简单增量法的计算精度高。这一方法在求解非线性问 题特别是求解塑性问题时得到广泛的应用。图 11.14 一阶自校正增量法的收敛过程11.5.2 Newton-Raphson11.5.2 Newton-Raphson法法对于式(1183)给出的单自由度非线性平衡方程()x 0 (1183)()xxn采用牛顿下山法，将在 点展开成 Taylor级数，取线性近似公式：()()()()xxddxxxnnn (1184)求非线性方程式(1183)的根x 可按如下公式进行迭代计算：xxddxxxxnnnnnn 111()/()(1185)单自由度非线性刚度方程一般为：kR()(1186)或 ()()()kRFR0 (1187)式(1187)代入式(1184)，有：()()()ddRFnnnn 1 (1188)对式(1187)求导，并注意到 dRd 0得：dddFdkT()()(1189)式(1189)即为体系在 处的切线刚度表达式，求 相应的迭代公式为：kRFRTnnnnnn111()(1190)RnN R为失衡力，式(1190)即为Newton-Raphson 法)求解结构非线性问题的最简单形式，其收敛过程如图11.15(a)所示。法(图11.15(a)对于多自由度体系，同样可以导出相应的迭代公式：)()(111nnnnnnTnRFRK (1191)这就是求解结构非线性平衡方程组的N R方法。由式(91)可见，N R法在每次迭代后都要重新形成 KT，对于大跨度桥梁结构进行这 一过程很费机时。为了减少形成总刚及其三角化分解的次数，有时用()KT0代替()KnT ，这样，仅进行一次切线刚度阵和三角化分解计算，后面的迭代只是线性方程组 的回代，这种 方法称为修正的 N R法(M N R法)。图11.15(b)给出了该方法的迭代过 程。M N R法在每次迭代中均用同一斜率，收敛较 N R差。图11.16给出了 N R法 和 M N R 法求解非线性方程组的流程，编程时可将这两种方法结合使用。图 11.15(b)N R法的收敛过程11.5.311.5.3收敛准则收敛准则 在迭代计算中，为了中止迭代过程，必须确定一个收敛标准。在实际应用中，可以从结构的不平衡力向量和位移增量向量两方面来判断迭代计算的敛散性。数的大小可以用其绝对值来衡量，而对于一个结构，无论其节点力还是节点位移都是向量，其大小一般用该向量的范数来表示。设列向量v=，该向量的范数可以定义为：123(,)TnV V VV(1)各元数绝对值之和：VViin11 (1192)(2)各元素平方和的根 VViin221 21()/(1193)(3)元素中绝对值最大者 VVni max (1194)这3个范数记为),.,2,1(PVP，应用中可任选其中的一种。有了列向量的范数，无论是结点力向量还是结点位移向量，其“大小”均可按其范数的大小来判断。所谓足够小就是指其范数已小于预先指定的某个小数。取位移增量为衡量收敛标准的准则称为位移准则，若满足下列条件就认为迭代收敛：11idiiuuu (1195)式中:d位移收敛容差1iu位移增量向量的某种范数。实践证明，对有些问题，前后两次迭代所得到的位移向量范数之比值会出现剧烈跳动，以导致收敛不可靠。取不平衡结点力为衡量收敛标准的准则称为平衡力准则，若满足下列条件就认为迭代收敛：iPPP (1196)式中:P为外荷载向量;Pi为不平衡力向量;p为不平衡力收敛容差。上面公式中取哪一种范数，按道理可以任选。有学者认为，在用平衡力准则时，取 P2比较好；在用位移准则时，取 u更为方便。在非线性比较严重的 问题中，用位移准则更合适。有的学者还用能量 PuT作为收敛标准，综合了力与位移两个方面，但要增加更多的计算量。例：图11.17结构在图示荷载作用下杆件材料仍将处在线弹性工作阶段，弹性模量E=2.0 kN/，杆件的横截面积A=4 。试采用带有一阶自校正的增量法和Newton-Raphson 方法对其进行几何非线性分析。4102cm2cm图 1117 结构示意 解：根据对称性，该结构仅有一个线位移自由度，即C点的竖向位移，可记为。将CB杆视作一个单元，其长度 cm，将该单元局部坐标系 02.1002100220l的原点设在C点，单元的方向角 02100 02 arcsin(.)1)带有一阶自校正的增量法。这一求解方法在每级增量荷载之后对结构进行一次平衡修正，把不平衡力并入下一级荷载增量中，由(1124)、(1125)和(1129)可得:20002002,2,sin2lEAKlNKlEAKL01.0cos00l由式(1121)和(1123)可知：1110lBBBL将上式代入式(114)，并注意到N=A，C点局部坐标系中的杆端力为：XYN111由此可计算C节点竖直方向的节点合力：)02.0(2)cos(sin2)cossin(2000101NNYXF采用P=P/10的荷载增量步，按一阶自校正法求解的过程如表11-3所示。表114 一阶自校正法求解的过程 由表114可以看出：(1)本例采用一阶自校正的增量法所得节点位移和杆件轴向力分别比精确值大5.32%和8.28%，比采用简单增量法(相应误差分别为24.32%和40.31%)有很大的改善。(2)利用一阶自校正的增量法可以得到比较精确的荷载变位曲线，而且计算工作量不大，因此一阶自校正增量法很合适于求解塑性问题等路径有关的非线性问题。(3)一阶自校正增量法要求在荷载变位反应非线性程度较大的区段采用较小的增量步。2)Newton-Raphson 方法 采用Newton-Raphson 方法求解时，结构切线刚度矩阵和不平衡力与一阶自校正增量法计算时相同，但现在要求通过迭代使结构达到平衡。由于几何非线性解与加载路径无关，本例将全部荷载一次作用于原结构。取位移收敛条件 /103则Newton-Raphson 方法的求解过程如表115所示。，由表115可见：(1)采用Newton-Raphson 方法经6次迭代运算已满足预定的位移收敛条件，此时节点位移和杆件轴向力的误差分别为0.08%和0.12%。(2)迭代收敛时结构初位移矩阵KL和初应力矩阵K的值已远大于线性刚度矩阵K0的值，其中以KL的值为最大。(3)若需要得到精确的荷载位移曲线，可以将荷载分成若干级进行迭代，使结构在每 一级增量荷载下都达到平衡。表115 Newton-Raphson 方法的求解过程11116 6 小结小结 本章介绍了桥梁结构非线性问题的分类及几何非线性分析理论的发展过程。通过变形体的运动描述，介绍了总体拉格朗日列式法和更新的拉格朗日列式法求解非线性问题的思想及其异同和适用范围；通过平面桁架单元的切线刚度矩阵的推导，说明了TL列式与UL列式两者的统一性。本章还讨论了稳定函数与几何刚度矩阵关系、弯矩对轴向刚度的影响、活载几何非线性分析方法、几何非线性调值计算方法等桥梁结构几何非线性分析中的特殊问题，最后，给出了非线性方程的求解和收敛准则的确定方法，为大跨径桥梁几何非线性分析奠定了理论基础。习习 题题 111 图11.16所示结构，假定杆件始终处于弹性工作阶段，材料的弹性模量 E=2.0 kN/，杆件的横截面积A=4 。4102cm2cm 试采用NewtonRaphson 方法计算结构在C点的竖向位移以及杆件的轴向力。112 试分别采用一阶自校正增量法和NewtonRaphson方法计算图11.17所示结构中点的水平位移。已知杆件为绝对刚性；弹簧铰的转动刚度为k=20kNm/弧度。113 试用NewtonRaphson方法编写平面桁架大位移分析的计算程序。