西南交大李乔教授：桥梁结构有限元建模中的若干问题

在结构分析领域，有限元法是目前应用最为广泛的数值方法。 按照未知数的特性划分，有限元法分为位移法、力法和混合法，其中位移法因其高度的灵活性和适应性、强大的模拟能力和求解能力，成为主要的结构分析方法。

桥梁结构分析属于结构分析家族中的一员，当然离不开有限元法(图1)。目前从桥梁的设计计算到施工模拟计算、从理论研究到应用研究，都高度依赖有限元法。可以说，有限元法于桥梁结构分析而言，近乎于“万能”的。有限元法对桥梁工程技术人员如此重要，如同一个忠实的仆人，每天为你服务。但你对她真的很熟悉吗？你知道她也有强烈的个性吗？你如果不尊重她的个性，她就会悄悄地给你设置陷阱，让你吃过期的食品或者给你用假冒伪劣商品。

在对桥梁结构进行有限元建模时，首先要了解桥梁结构分析的理论，了解有限元法原理和特性，只有这样才能建立一个能够反映结构实际特征、符合力学原理的有限元数值模型 。忽略力学原理，只把结构计算当作一堆没有物理意义的数字对待，这是目前不少学生容易犯的毛病。本文拟对常规的桥梁结构有限元建模进行讨论，主要讨论其中一些容易被忽略的问题以及建模中的技巧，供大家参考。为了简单，以下均以平面单元和梁单元为例，对于其他单元，其原理是一样的。

（c）

图1  桥梁结构有限元分析（图片a、b引自互联网）

有些桥梁设计师或年轻桥梁工程师在使用有限元软件的初期，容易有一个潜意识里的错误认识，即认为有限元是一个精确的分析方法。

有限元是一个强有力的、有较高精度的数值方法，但并不是一个精确的方法。简单回顾一下有限元法理论，就能容易地知道这一点。下面以位移法有限元为例来说明。

有限元法是把连续的（具有无穷多个质点的）求解域划分为含有有限个子域（单元）的离散网格形式（见图2）。这些单元之间通过有限数量的网格节点相互连接，用这些离散的节点处的位移表征实际结构中连续的位移场。在这些节点上，位移是连续的，但位移的各阶导数并非都是连续的，取决于节点位移未知数中是否包含这些导数。这是有限元法的第一个近似之处。在单元内部，由于不知道真实的位移场函数，所以有限元法采用近似的数学表达式（如多项式）通过节点位移进行插值来描述，该表达式称为插值函数。这是第二个近似之处。在节点之间的连续边界上，位移及其各阶导数是否连续则取决于单元位移模式和插值函数的形式，一般最多使一阶导数连续。这是第三个近似之处。

图2  有限元网格

在位移法有限元中，采用能量原理建立以节点位移为基本未知数的平衡方程组，所以求解有限元问题归结为求解一个大型线性方程组的问题（即使对于非线性问题，在进行每一迭代步或者增量步求解时，也是一个线性的方程组）。该线性方程组一般为式（1）所示的形式。其中[K]为方程组的系数矩阵，称为总刚度矩阵，{δ}为方程组的基本未知数，即结构的节点位移列向量，{P}为荷载列向量。方程组可能含有成千上万乃至几百万个未知数，那么由于计算机中数字的字长（通俗地说就是有效数字位数）是有限的，因此在求解过程的大量计算中，会由于舍入误差而产生计算误差，这是第四个近似之处。这些只是有限元法本身的近似之处，还不包括力学理论、模型与实际结构差异、计算参数与实际结构差异等引起的误差。

上面的描述虽然不太严谨，但已经能够清楚地看出，有限元法是一种近似的数值方法，是有误差的方法。  

如上所述，既然有限元是一种离散的、近似的数值方法，那么她的计算精度显然跟网格的划分密切相关，并非随便怎么划分都能得到准确的结果。

网格密度： 这是个经验性的东西，当没有经验可借鉴时，则需要通过不同网格化分的数值计算比较（也有称之为数值试验的）来确定。虽然从纯理论上说网格越密，计算结果越逼近精确解，但过密的网格会导致计算量加大，也会增加舍入误差。所以实际计算并不是网格越密越好。而网格过稀则同样会导致计算精度的降低。

单元的形状： 如上所述，单元内部的位移场是用节点位移通过插值函数来近似表达的。以一个平面四边形单元为例，如图3a、b所示，如果要在平面内的两个方向上通过插值得到对位移场较好的近似，显然 单元 应该尽量接近正方形，而 不能过于歪斜 。如果采用图3c所示形状的单元，则插值出来的位移场跟实际会有较大的偏差，从而导致较大的计算误差。

单元的尺度： 首先讨论相邻单元之间的相对尺度。式（1）中的总刚度矩阵[K]是由各单元的刚度矩阵叠加得到的，而各单元刚度矩阵则与单元的几何构形、材料特性以及插值函数等相关。如果相邻单元的尺度相差过大，就会使相邻单元刚度矩阵的数值相差悬殊，因而导致系数矩阵[K]中相邻的元素数值相差过大。当这种差别的量级接近计算的有效数字位数时，就会使方程组成为数学上的所谓病态方程，从而使计算结果产生误差，严重的病态方程会得出完全错误的结果。因此，应注意， 相邻单元尺度不能相差过大 ，要逐步过渡。有限元软件一般采用8字节长度的浮点数字，按十进制大约有15位有效数字。因此，一般不过于悬殊的尺寸差别（比如梁单元长度相差10倍以内）不会导致太大的误差。

相邻单元材料的弹性模量： 这是最容易导致病态方程的问题。有时为了刚性连接两个节点，就在这两个节点之间设置一个弹性模量比常规单元材料高十几个甚至几十个数量级的单元。这会导致严重的病态方程和错误结果。所以请注意， 不要使相邻单元刚度相差过大 。

在很多有限元软件中，都有所谓的主从节点和带刚臂单元，利用它们可以模拟很多结构的连接和特殊构造，这在下面各节的讨论中会看到它们的用途。那么它们究竟是什么工作原理呢？

主从节点： 如图4所示，设有两个节点 j A 和 j B ，二者可以坐标相同（即位置重合），也可以相距一段距离(Lx, Ly)。如果节点 j B 的部分位移{dB}不是独立的，而是跟 j A 的对应位移{dA}直接相关（即如果知道了{dA}，就可以求得{dB}）。这时就称这两个节点为主从节点，其中 j A 称为主节点， j B 称为从节点，{dA}与{dB}则称为有主从关系的自由度。根据线性方程组的特性，必须在式（1）中去掉不独立的未知数{dB}及其对应的方程，否则方程组就没有解。假如在有限元软件里面不考虑主从节点之间的距离效应，处理方法就是直接令{dB}={dA}。如果考虑这种距离效应，就必须考虑主从节点间转动与平动的相互关系。

图4   主从节点

（图中两个节点之间的水平位移有主从关系）

带刚臂单元： 图5所示为一个平面梁单元，在其两端 i’ 和 j’ 处各带有一个刚度无穷大的刚臂 i’-i 和 j’-j 。单元是通过刚臂与节点 i 和 j 相连接的，而不是梁端点直接连接节点。这样的单元就叫作带刚臂单元。此时 i 和 j 是整个有限元模型中的独立节点，它们的位移是独立的未知数。而 i’ 和 j’ 不是独立节点，它们只是带刚臂单元的内部端点。显然，它跟普通梁单元不同之处在于多了刚臂，因此在每个刚臂的两端（ i’ 与 i ， j’ 与 j ）之间的转动与平动相互关联。例如，不但 i 的平动会引起 i’ 的平动， i 的转动也会引起 i’ 的平动，反之亦然。

图5   带刚臂单元

如图6a所示，节点 j B 为连续梁支点处的节点，位于梁轴线上，节点 j P 为对应的桥墩顶部节点。显然 j B 与 j P 间距离 h BP 为梁轴线到其下缘的距离加上支座高度。要正确模拟这种连接，有三种方法可供使用：

1）在两个节点之间设置一个刚度较大的单元（图6b中的 e P1 ）。为避免重复计算重量，该单元自重应该设为零。同时该单元刚度应该比对应的桥墩单元大（1~3个数量级），但不能过大，否则会出现前面说的病态问题。如果墩梁之间是铰支座，那么在点 j P 处应该设置两个坐标重合且具有主从关系的节点 j P 1 和 j P 2 （图6b），其中 j P 1 为单元 e P1 的下端， j P 2 为单元 e P2 的上端。如果是连续刚构，则只用一个节点 j P 作为上下单元 e P1 和 e P2 的公共节点即可；

2）采用考虑了距离效应的主从节点模拟。如可把 j P 作为主节点， j B 作为从节点，二者之间的所有自由度都有主从关系（连续刚构），或者2个平动自由度有主从关系（固定铰支座），或者只有竖向自由度有主从关系（活动铰支座）。这时要特别注意不同的软件对距离效应的处理方式是不同的；

3）采用带刚臂单元模拟（图6c），这种情况只适合连续刚构。 e P2 为上端 j P 带有一个刚臂的单元，连接到主梁节点 j B 。注意此时的 j P 不再是一个独立的节点（相当于图4中的 i’ 或者 j’ ）。

在计算梁拱组合体系桥梁时，拱肋在拱脚处跟主梁的连接部位与上面的梁墩连接类似，也可以用前述的方法进行模拟。

多数斜拉桥的索-梁及索-塔锚固点坐标都跟梁轴线或者塔轴线不重合，即都具有偏心距，如图7所示。此时斜拉索的拉力除了对梁或塔施加一个索轴线方向的拉力外，还会有一个附加力矩（拉力乘以偏心距）。这种附加力矩对结构的影响不可忽视，所以不能简单地把锚固点直接搬到梁轴线或者塔轴线上。正确的模拟应该是把斜拉索也作为带刚臂的索单元，通过刚臂连接到梁及塔的轴线上（见图7）。

图7   斜拉索偏心距模拟

如图8a所示，钢-混凝土组合梁是由混凝土桥面板和钢梁组成，在大跨度桥梁中使用这种结构时，比如斜拉桥中，一般是先架设钢梁，然后进行斜拉索初张拉，再架设预制桥面板，然后通过剪力连接件和湿接缝把混凝土桥面板和钢梁连接起来，形成组合截面。显然，钢梁自重、斜拉索初张力、预制混凝土桥面板重量以及湿接缝重量都是作用在钢梁上的，只有在湿接缝达到强度后施加的索力和荷载才由组合截面承担。

在进行这种结构的有限元建模时，如果简单地直接把两种材料构成的组合截面梁作为一个整体来模拟，即每个单元都包含混凝土桥面板和钢梁（图8a），采用换算截面几何特性（把两种材料换算成一种材料），单元轴线为组合截面轴线（图8a中的红线），就等于忽略了钢梁和混凝土桥面板分别架设的过程，而这些过程可能是钢梁受力的最不利阶段；这种模拟也无法正确计算混凝土桥面板的收缩徐变效应，因为即使把截面换算成混凝土，所计算的收缩徐变也是全截面的，而不仅仅是桥面板的，同时由于忽略了分阶段安装过程，构件受力产生差别，也导致收缩徐变计算不正确。

图8 组合梁的模拟

那么怎样才能正确模拟组合梁桥呢？方法之一是采用上下两个单元分别模拟混凝土桥面板和钢梁，从而可以分别在不同时间安装它们，并可正确计算混凝土桥面板的收缩徐变效应。该方法的关键不在这里，而在于如何考虑桥面板轴线与钢梁轴线之间的距离h（图8a）。切不可把上下单元端部节点直接共享，这将导致严重的错误，因为这忽略了距离h的作用，使桥面板和钢梁轴力差构成的那部分弯矩被忽略，而这部分弯矩是组合梁弯矩中的主要组成部分。有效的办法是钢梁采用普通梁单元模拟，而混凝土桥面板采用前述的带刚臂单元模拟（见说桥17），它通过刚臂连接到钢梁单元节点上，如图8b所示。当然这种模拟方法也有缺点。因混凝土桥面板和钢梁只在单元两端连接，中间部位没有连接，这与实际的组合截面梁有所差异。当单元中间部位没有荷载作用时尚能正确模拟，但当中间有荷载时，就会产生误差，这可以从力学和有限元理论得到证明。所以，采用这种方法时，要加密单元划分以缩小单元长度，减小误差。实际计算表明，只要单元不过长，这种方法能够正确模拟组合截面桥梁。

模拟组合梁的方法之二是把上述的刚臂换成短的大刚度梁单元（图8c中的连接单元），即通过连接单元把混凝土桥面板单元和钢梁单元连接起来。其优点是当所用的有限元软件没有带刚臂单元时也能计算。缺点则是增加了单元和节点数量，同时为了使桥面板单元和钢梁单元满足组合截面的平截面特性，必须使其间的连接单元刚度比桥面板单元和钢梁单元大几个数量级，且由于这样的单元数量众多，会在整个结构刚度矩阵中形成大量的刚度差异区域，如前所述，会影响计算精度。

模拟组合梁的方法之三是在形成组合截面之前与方法之一和之二相同，而在形成组合截面之后，把原来的钢梁单元拆掉，换成组合截面梁单元。也可以通过改变梁单元截面几何与材料特性，演变成组合截面单元。这种方法的问题在于如何考虑截面由钢梁变为组合截面后，轴线位置发生变化（轴线上移），且变化前后的截面应力如何继承与累加。如果软件不具备自动处理这种关系的功能，那么可能得到的计算结果会有较大的误差。  

钢管混凝土结构的模拟较为简单，类似于上述组合梁，仍然把钢管和管内混凝土分别建立单元（见图9），并在不同阶段安装。此时钢管单元与混凝土单元轴线重合，端节点重合，不需要设置刚臂。在混凝土浇筑后且强度没有达到要求之前，混凝土重量及其他荷载由钢管承受。而在形成组合截面后施加的荷载，由钢管和混凝土共同承担。这里重点谈谈如此模拟到底忽略了什么，以及应该注意什么。

图9 钢管混凝土结构模拟

首先，与上述组合梁的模拟类似，钢管单元和混凝土单元只在节点处连接，因此单元长度不能太大。第二，没有考虑钢管对管内混凝土的套箍作用，因此不能考虑由于套箍作用而使混凝土强度提高，也不能考虑由于套箍作用引起的混凝土单元变形与单向受力单元变形的差异。好在这种近似是偏于保守的。第三，忽略了由于上述套箍作用而使混凝土的徐变效应发生的变化，这种变化导致钢管和混凝土之间承担荷载比例发生变化，也导致整个结构徐变次内力发生变化。这到底是偏于保守还是偏于不安全不能一概而论，对不同的构件产生的影响是不同的。

很多连续刚构桥采用双薄壁墩结构，如图10a所示。其墩梁连接的模拟已在“说桥17”中讨论，此处讨论桥墩与承台及桩基础组成的体系如何模拟（图10a中圆圈内的部位）。

图10 双薄壁墩-承台-桩基础的模拟

对于这种连接的模拟可以采用多种方法。例如，可以按照图10b所示方法进行模拟。因为承台的主要受力特征为受弯和受剪，因此图中的承台单元为水平放置的梁单元。墩单元和桩单元可以采用带刚臂单元与承台节点连接，也可以不用带刚臂单元，而采用主从节点与承台节点相连。对于基础的弹性支承作用，则可近似采用竖向及水平弹簧支承于桩单元节点上（图中未画出）。

如果想简化计算，也可以如图10 c所示那样，把所有桩基础及所受约束简化为一个弹簧支承，该弹簧具有水平、竖向及转动约束。这样做还有一个好处，即当需要计算基础不均匀沉降引起的次内力时，能够通过一个支承使同一基础上的双壁墩的两个壁（即两个墩柱）同时沉降，比较符合实际情况。而若采用图10b的模型，则各根桩都有独立的约束，因此很难进行最不利的基础不均匀沉降组合。如果把承台取消，直接在两个墩柱底各设置一个约束，当软件不能自动把这种情况识别为一个基础时，也会使两个墩柱分别沉降，从而引起巨大的次内力，这是不符合实际情况的。

这个问题乍看上去有点儿莫名其妙，是不是太业余了？恒载位移就是恒载引起的位移，还能是什么？先别急着下结论，看了下面的例子再说。

图11所示为一等截面悬臂梁，由同一种材料构成。从结构力学可知，该梁在其自重作用下，自由端的挠度为 f =(qL ⁴ )/(8EI) 。其中q为自重荷载线集度，L为梁长，EI为截面抗弯刚度。但如果采用有限元软件计算，并且分阶段安装各单元，比如划分为4个单元（图12a），并分4个阶段依次安装1~4号单元，则多数有限元软件输出的自由端挠度 f 都会比上述的 (qL ⁴ )/(8EI) 小很多。而如果你把所有单元在一个阶段内一次安装上去，输出的自由端挠度就是 (qL ⁴ )/(8EI) 了。可是从结构力学原理可以知道，对于这样一个结构，若荷载不变，一次安装与分次安装的挠度好像不应该有什么区别。那么问题来了，软件输出的恒载挠度到底是什么涵义呢？如此看来，本节提出的问题也并非毫无意义了。

图11 自重作用下的悬臂梁

在有限元软件中，当一个单元安装时，例如图12b所示的3号单元安装，该单元与先前已经存在的节点3相连接的左边端点，其初始位移只能等于节点3已经发生的位移V 3 。但本次新安装的节点4，程序对其初始位移有两种处理方法，一是新安装单元3沿着前面已存在的节点3处的切线方向放置（图12b中的红线），从而得到新安装节点的初始位置4’，这种方法称为“切线安装法”。另一种方法是将新安装节点4的初始位移设定为零（图12 b中蓝线及点4”），这称为“零位置安装法”。

图12 两种分阶段安装方法

显然在分阶段安装过程中，软件采用上述不同的处理方法所得到的恒载位移数值是不一样的。采用切线安装法时，分阶段安装和一次安装的最终结果是一样的，对于上述例子，自由端挠度与结构力学公式结果相同。但当采用零位置安装法时，因为每个节点的初始位移为零，所以最终位移数值是从其安装开始到最终状态之间的位移增量。对于上面的例子，自由端挠度值只是4号单元安装后由于其重量引起的挠度。与切线位置安装法相比，少了安装时由于前面已安装单元发生的位移和转角而引起的刚体位移。图12b中的V 40 就是4号节点的这种差异。

综上所述，在使用有限元软件计算恒载位移时，要搞清楚软件采用的是什么安装法，不能想当然。目前的大多数软件在默认情况下都是采用零位置安装法的，因为这种方法有很多方便之处，如编程方便、计算安装线形方便、模拟施工过程各节点位移方便（从吊装或者立模标高起算）等。如果要想计算切线安装法的位移，可以先一次性把体系转换前的所有单元都安装上去，但各个单元的重量为零，然后分阶段施加单元重量和其他荷载。这样得到的体系转换前的位移就是切线安装法的位移，因为此时没有施加重量和其他荷载的单元随着前面的单元一起发生刚体位移。而在体系转换之后的位移，仍然可以计算切线安装法位移，但涉及到的问题较为复杂，且已经超出了本文的研讨范围，限于篇幅，这里不进行讨论。

以上讨论了采用有限元法进行桥梁结构建模时常遇到的一些问题，由于建摸问题的复杂性和多样性，本文只是涉及到了其中很少的一部分，期望能够为初期使用有限元建摸的人员提供一些参考。另外，文中只是以平面有限元为例进行说明，至于三维有限元建摸，其道理是类似的，只是更麻烦一些而已。