第二章 桥梁结构有限元法及可视化软件的开发1

第二章桥梁结构分析有限元法及可视化软件开发
2.1 前言
有限元法是伴随着计算机技术进步而发展起来的一种数值分析方法，广泛应用于分析各类复杂科学与工程问题。

由于它采用矩阵形式表达基本公式，能灵活地处理和求解各种复杂问题，便于应用于计算机编程运算等优点，有限元法显示了其强大的生命力。

目前，各种有限元软件广泛应用于科学研究和工程应用领域。

但许多大型通用软件涉及领域广泛，在具体问题的操作中难于满足工程设计人员各方面的要求；一些有限元应用软件程序庞大，不利于普通技术人员的掌握和使用。

为此，结合工程实际编制易学、易用的可视化专用有限元软件能极大的提高工作效率。

本文开发的基于Visual Basic语言的大跨桥梁结构施工分析系统(BSCAS)，采用桥梁结构施工控制中普遍采用的正装法(前进分析法)和倒装法(倒拆分析法)实现了程序的结构化分析和施工模拟。

VB友好的用户界面，使设计过程生动直观，具有有限元结构分析初步知识的工程技术人员在界面的引导下，可以畅通无阻地进行应用设计。

2.2 桥梁结构有限元法概述
有限元法是将连续体分成有限个单元，单元间相互有结点连接的理想结点系统。

分析时，先进行单元分析，用结点位移表示单元内力，然后将单元再合成结构，进行整体分析，建立整体平衡关系，由此求出结点位移。

有限元求解程序的内部过程如图2.2.1所示：
2.2.1 桥梁结构有限元分析一般步骤
采用有限元法进行桥梁施工控制中的结构分析计算，首先要建立数据文件。

数据文件准备一般分为四步：
(1) 桥梁结构的模型化。

桥梁结构的模型化就是将实际结构理想化为有限个单元的集合。

根据结构的受力特性与工作行为选择恰当的单元形式来模拟实际结构以及选择正确的约束模拟形式，选择正确的约束模拟形式尤为重要。

就结构分析模型来看，与一般的已成桥梁分析不同的是施工控制中的结构分析模型一般是随着施工的不断推进而不断变化的，这是由于桥梁在形成过程中的结构体系是在不断变化的。

实际工作中，通常可考虑建立一个统一的模型，而对某个施工状态的结构模拟则可通过某些单元的是否激活来实现。

计算模型中单元的选择应以能准确描述施工过程中结构受力与变形状态为准。

对于桥梁结构而言，最主要的是结构纵向的受力分析计算，而将纵向分析模
型近似地处理成杆件系统，往往是一种可取的简化模型。

一个模型可以是由一种单元组成，也可是由几种单元组成，对一些特殊的施工工艺需要采用特殊的单元
来描述。

如连续刚构桥的墩梁固结处，则假定墩顶零号块是一个具有无限大刚度的刚性杆件，应用程序的刚臂单元来处理。

(2) 桥梁结构的离散化。

桥梁结构的离散化就是在模型化处理后，将结构离散为带有有限个自由度的结构，这些单元在结点上相互联结，构成一个与真实结构等价的计算模型，即将真实的复杂空间结构模拟为平面的杆件系统。

单元大小与结点位置确定应充分考虑结构受力情况与施工单元的划分。

结构离散化必须遵循以下原则[5]：1) 计算模型应尽量符合实际结构的构造特点和受力特点，以保证计算分析结果的真实性；2) 保证体系的几何不变性，特别是在复杂的体系转换过程中更应注意，同时要避免与实际结构受力不符的多余联结；3) 在合理模拟的前提下，应尽量减少结点数目，减少未知量数目，以缩小计算规模，节省时间和计算机空间；4) 本构关系的确定应尽量反映材料的性质。

杆件单元的划分，应根据结构的构造特点、实际问题的需要以及计算精度
的要求来决定，用来划分单元的结点应包括构件的转折点、交接点、截面变化点。

由于杆件结构有限元程序的计算结果是以结点位移及内力形式输出的，故在进行结点划分时，应根据设计需要验算的截面以及求影响线的作用点的要求，来确定所需增加的中间结点。

(3) 按所用软件的输入要求形成数据文件
(4) 检查校正数据文件。

计算模型最终体现为数据文件，数据文件正确方能保证计算模型的正确，乃至才能保证计算结果的正确性。

其次，运行分析软件。

对于桥梁施工控制中的结构分析，由于计算模型随着施工过程的改变，同时要求分析跟踪进行，应采用具有施工控制跟踪、仿真分析功能的软件。

最后，对分析结果进行分析和处理。

2.2.2 桥梁结构几何非线性有限元分析的单元刚度矩阵
在结构发生大位移时，单元的外力和位移之间的线性关系不再适用。

可以将荷载分成几步逐渐施加，从而把一非线性问题用一系列线性步骤来处理，求出每一荷载增量引起的位移增量，即增量法。

在每一步中刚度矩阵可以分为两部分，一部分是在每一步按当时的几何构造计算的标准的弹性刚度矩阵[]E K ，另一部分是几何刚度矩阵[]G K 。

几何刚度矩阵不但与几何形状有关，而且与正在开始计算的这步开始时单元的内力有关，有时也称为初应力刚度矩阵。

梁单元同时考虑轴向变形和弯曲，其切线刚度矩阵[]T K 由两部分组成，即：
[][][]G E T K K K +=。

关于弹性刚度矩阵[]E K 、几何刚度矩阵[]G K 的推导，很多文
献中都有非常详尽的过程，这里不再重复。

弹性刚度矩阵为：
[]
⎥
⎥⎥⎥⎥⎥⎥
⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡-
-
-
---
--=
I l
I I
l
I l I l I l
I l I A A I l
I I l I l I l I l I l I A A
l E K e
460
2606120612000002604606120612000002
22
2
几何刚度矩阵为：
[]⎥⎥⎥⎥⎥
⎥⎥⎦
⎤⎢⎢⎢⎢
⎢⎢⎢⎣⎡--------=2222
430
30
33603360000000
30430
3360336000000030l l
l
l
l l l l l l l l l P K G
式中：P 为梁单元的轴向力(拉为正，压为负)，A 为梁单元的截面面积，I 为梁单元的截面惯性矩，E 为梁单元材料的弹性模量。

2.3 Visual Basic 软件简介及BSCAS 系统实现
2.3.1 Visual Basic 软件简介
Visual Basic 是Microsoft 公司在Basic 语言的基础上开发的Windows 应用程序开发工具，自1991年正式发布1.0版以来，已发展到7.0版[33]。

Visual Basic 是一种面向对象的编程语言，具有简单易用、开发效率高的特点，被认为是当前最为流行的“快速应用工具”。

Visual Basic 程序的开发思想是面向对象进行编程，对象是程序中的主角，既是过程的载体，也是被处理的实体。

一个Windows 应用程序中的窗体和窗体上的控件都是对象，这些对象不需要编写大量代码来描述其界面元素的外观和位置，而只要在设计程序时把VB 预先建立好的对象拖放到窗体中即可，并且这些对象具有属性、方法、事件，开发者所要做的是为对象编写事件驱动代码。

用VB 进行程序设计时，人机交互界面的设计非常方便，它包括：(1)窗口本身及里面所有自制项的设计；(2)设置每个对象的属性；(3)为某些对象编写相关的程序码；(4)对于试运行通过的软件在Windows 环境下运行。

在有限元结构分析中，是以结点位移为基本未知量的，在杆件结构中又称为矩阵位移法。

用这种位移有限元法分析结构时，所需输入的原始数据如控制信息、结构几何性质、材料弹性常数、支撑约束条件等，涉及到的数据类型多、数据量大、算法复杂，影响求解结果准确性的因素较多。

应用VB 开发可以使数据输入简便，操作方便，用户界面友好，并能及时反馈和纠正数据输入中出现的错误，能很好地满足设计计算要求。

2.3.2 程序开发思路
有限元软件，一般包括前处理、分析计算与后处理三大模块。

基于VB 开发的用户界面友好的有限元软件的基本思想是[60]：VB 用于人机交互界面设计和数值计算、AutoCAD 软件包用于前后处理的图形显示。

开发前处理模块，主要包括：用VB 设计数据输入界面，并采用模块化程序设计方法，开发输入数据的有效性和准确性检查模块、有限元网格自动生成模块；用AutoCAD 软件包实现生
成有限元网格的图形显示。

在后处理中，用VB设计有限元计算结果显示界面；用AutoCAD软件包实现有限元部分计算结果的图形显示；将重要有限元计算结果存入数据库，实现后处理模块中各子模块的集成。

并用VB对上述各模块进行集成。

BSCAS系统在编制过程中，采用现今流行的面向对象的程序设计思想，引入对象、属性、方法和事件等概念，将其分为几个部分：一为系统控制部分，为MDI主控窗口；二为系统的输入、输出模块；三为结构计算分析专用程序，接收到数据后，在后台运行计算程序；四为程序前后处理模块。

遵循以上程序设计思想，BSCAS系统的编制中，主要实现了桥梁结构的分析及施工模拟分析部分，前处理中实现了人机交互界面化设计、数据模块检查、结构图形生成和后处理中的结果数据文件显示及简单的结果图形显示。

限于时间，后处理的结果图形显示功能仍需进一步的强化，这也是今后需要完善的地方。

2.3.3 程序可视化实现
为了使设计计算过程充分可视化，软件界面中应用标签(Label)和文本框(TextBox)控件来输人参数与显示计算过程中的数据和计算结果，应用组合框(ComboBox)和单选框(OptionButton)控件进行参数与类型选择，并且采用AutoCAD软件包实现有限元部分计算结果的图形显示。

所有控件均按不同设计功能用框架(Frame)分类布置，使操作界面简洁明了。

Windows应用程序的基本载体是窗口，在窗口中做到了数量上的控制，主窗口能完成参数插入、缺省值选择、计算并显示结果等大部分设计工作，减小了用户的操作难度。

而次级窗口则可以进行一些主窗口无法进行或无需在主窗口中进行的设计操作与数据处理，如图2.3.2示，使应用程序窗口主次分明，功能各异。

同时，在程序可视化的实现中，特别注意到了程序界面的友好性，方便快捷的工具栏按钮上，设计了工具提示文本，实现了快速访问命令按钮。

窗体下部的状态栏可以在程序运行的任何时间，告知用户当前的运行状态。

主窗体可以随意
新建文本文档，且不同的文档中可以打开不同的文件，用程序窗口工具条中的纵铺、平铺功能，可同时打开输入数据模块和结果数据文件，进行对照检查(如图2.3.4)。

程序运行时每日一帖对话框和程序计算完毕时的提示框，如图2.3.3,
生动形象，充分展现了程序的人性化友好界面。

同时，程序的界面化过程力求取得与人们所熟悉的Windows环境在风格上的一致，便于工程设计人员使用。

2.4 BSCAS系统有限元编程思想简述
2.4.1 BSCAS计算程序和单元分析
BSCAS 系统适用于可作为平面杆系处理的桥梁结构体系。

能够逐段形成结构体系，自动完成施工阶段静力体系的转换。

对连续刚构桥梁施工中使用的悬臂施工工艺，能够从体系形成到使用阶段连续进行分析。

下面介绍程序的部分设计原理和设计思想。

(1) 结构总刚度矩阵的一维变带宽存储
总刚度矩阵是耗费内存最多的一个数组，一个结构的规模如超出计算机的容量往往都是由于该数组过大所致。

由于一般工程结构的总刚度矩阵都是对称带状矩阵，其非零元素集中在主对角线两侧的斜带形区域内。

如果采用一维变带宽存储对结构刚度矩阵的下三角部分采用变带宽按行存储，就能有效地节省计算机的存储量。

当矩阵采用一维变带宽存储时，只要知道每一行主对角元素在一维数组中的位置，就可以找出全部元素在一维数组中的地址。

为了确定对角元素的地址，必须先求出结构刚度矩阵各行的半带宽。

在矩阵的下三角部分，从某一行的第一个非零元素到主对角元素的个数称为半带宽。

由于稀疏矩阵每行的带宽是变化的，故用一维数组存放时，还必须有保存每行带宽的信息。

如果用数组()n V 来记录对角元素在一维数组中的序号，且约定()0V =0的话，则第i 行的半带宽为：()()11---i V i V ，(1=i ，2,n )，于是按下三角存放的系数矩阵的任一带内元素ij K 在一维数组R 中的地址转换公式为：
()j i i V p +-=
式中：p 为元素ij K 在一维数组中的地址；i 、j 分别为元素ij K 在结构刚度矩阵中的行号和列号。

总刚度矩阵K 与一维数组R 元素间的对应关系为：ij K =()p R 。

(2) 特殊单元的分析和处理
1)刚臂单元的处理。

众所周知，在计算钢筋混凝土刚架桥等一类结构物时，理论上简化的结构计算图式往往会与实际的结构物不相符合，其中有两个问题较为突出：一是结点上所有的杆件轴线未必会相交与一点，常常有几根杆的轴线偏离结点，如图2.4.1所示，这时理论的计算图式如果取一个结点计算势必失真，
如果取几个结点计算则增加了方程的阶数；二是杆件进入结点附近时，常常是和刚性很大的结点块联在一起，这时杆件的刚度只能是指结点刚性块以外的部分，如果把杆件的长度计算为结点之间的距离，势必把杆件的刚度算得太小了。

如连续刚构桥墩顶零号块往往假定是具有无限大刚度的刚性杆体，应充分考虑其刚域效应。

这种情况下，采用图2.4.1所示的计算图式，杆件与结点有刚臂相连接显然比较符合实际结构的受力，这就是带刚臂的杆单元问题。

在本程序的实现上，带刚臂单元的刚度矩阵是由无刚臂的单元刚度矩阵推导得出的。

在处理带刚臂单元杆件的子程序中，考虑到了程序的通用性，对于带或不带刚臂的杆件都可实现调用和分析计算。

根据图2.4.2所示的几何关系[29]，可以建立杆端位移和刚端位移的关系式：
⎪⎪⎪⎪⎭
⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣
⎡--=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧rk
rk
rk rj rj rj
k j
lk
lk lk lj lj lj d d θ
υ
μθυμθυμθυμ10
100
001
1
00
10001
(2.4.1)
或简写成： {}[]{}e r e l H δδ= (2.4.1) 式中： {}e
l δ——杆端位移列阵；
[]H ——变换矩阵；
{}e r δ——刚端位移列阵。

同理，由图2.4.3可得杆件的刚端内力{}e r F 与杆端内力{}e
l F 的关系：
{}e r F []H ={}e
l F (2.4.2)
此外，根据杆端、刚端结点力和结点位移关系：
[]{}e
l e
l K δ={}e
l F (2.4.3)
[]{}{}e r e r e r F K =δ (2.4.4)
由式(2.4.1)～(2.4.4)，经变换得到杆端刚度矩阵与刚端刚度矩阵的关系式：
[][][][]H K H K e
l T
e r = (2.4.5)
利用[]e r K 形成总刚度，通过解方程求得结点位移和每根杆的刚端位移{}e r δ，然后分别由式(2.3.1)和式(2.3.3)求杆端位移和杆端力。

以上带刚臂的单元刚度矩阵的实现，只需调用矩阵乘法子程序即可。

带刚臂单元实质上是在单元两端增设了两个刚性杆元(∞==EI EA )，这种处理方法可以避免在刚度矩阵中过多地引入刚性很大的单元，使刚度矩阵病态化，造成解的精度降低。

带刚臂单元主要应用于结构中刚性结点尺寸与杆件单元长度相比不容忽略及杆件单元轴线不交会于一点时，或考虑交结点刚域效应的情况。

2) 中间铰的处理。

桥梁结构的构件大都是刚性联结，但也有个别用铰接联结的。

对带铰杆件的处理，程序用特殊杆元，即两端铰接刚性杆(∞=EA ，0=EI )模拟,该特殊单元的长度可为大于零的任意值，容重ρ=0的单元，两端可带刚臂。

如图2.4.4所示。

图中，a )为仅传递竖向力的中间铰；b)为普通中间铰，用两个特殊单元使i 、
j 两结点的平移位移相同，而转角相互不约束模拟；c)中墩(或台)上支座为活动
支座或固定支座时处理方法不同，固定时处理与b)类似，活动时与a)类似。

程序的实现上，采用铰结端多编一个结点的方法。

如果杆件之间联结是铰接，如图2.4.5示，则在同一结点上多编几个结点号，这些结点号必须连续，并且坐标相同。

这种方法增加了未知
量数目，但结点位移分量编号规律性未被打乱，程序处理上简单，便于编写，且在计算分析中可实时判断杆件是否为铰接。

在铰接处连续条件的处理上，当i 结点与1+i 结点为铰接时，这两点的水平位移和竖向位移相等，而转角不相等，即变形连续条件为：
1+=i i u u ， 1+=i i v v
处理该连续条件的简单方法是，在两个结点间加一根水平连杆和一根垂直连杆，连杆的刚度取结构中最大刚度l EI 的300倍，在总刚度矩阵中把两根连杆的单元刚度矩阵也加进去。

另外程序还使用了两端固结刚性杆(∞=EA ，∞=EI )特殊单元，灵活应用这两个特殊杆元，可以方便解决实际结构中非刚性中间结点的约束模拟问题，使计算图式归结为全部刚结的形式，给程序使用中建立结构计算分析模型提供了方便。

(3) 逐段形成的结构体系的处理
对于大跨径预应力混凝土桥梁，其最终体系一般总是通过一系列的施工阶段而逐步形成的。

在各个阶段，结构可能具有不同的静力体系，因此在施工中存在
结构体系形成与转换的过程。

这一过程可能是多种多样的，不同的施工方法决定了不同的体系形成过程，反映在结点位移自由度、总刚度矩阵等都随施工阶段而发生变化。

显然，不同的施工方法，其施工阶段的受力状态各不相同，而施工阶段的受力状态又往往有异于桥梁建成后的使用阶段。

为了保证施工安全和长期正常使用，在桥梁结构分析时，就必须计算每一个受力阶段结构各部分的应力变形状态。

使用上述带刚臂的杆件单元和两种特殊杆元组成一个包括全部永久和临时杆件元组的结构体系作为综合计算图式(体系形成诸阶段计算图式的总和或包络，结构总刚度矩阵的长度由此来确定)，而后通过转换杆件元截面刚度特性EA、EI的办法来实现逐阶段计算图式的转换。

为了生成某一施工阶段的计算图式，需判别杆件单元在该阶段是否存在。

当该阶段某些杆件元不存在时，则将其截面刚度EA、EI充零。

如此根据该阶段存在的诸杆件元组集成的总刚度矩阵，显然代表了该阶段实际结构体系的刚度特性。

当某一杆件元被拆除后，在以后的各个阶段中，其截面刚度特性将永远充零，即认定该杆件元不复存在。

若杆件元不是拆除而是撤换，则用撤换的杆件元，需作为另一单元编号计入。

当撤除临时杆件元时，则将到撤除前一阶段为止的临时杆件元内力作为撤除阶段的外力施加之。

程序用以上方法，方便处理各种错综复杂的实际桥梁施工过程。

施工过程简例：如图2.4.6、图2.4.7所示为两跨连续梁桥悬臂施工过程中结构体系转换过程中，实际施工阶段和对应程序内计算图式：
图2.4.6 实际施工阶段图2.4.7 程序内计算图式
表2.4.1 程序施工过程处理
(4) 施工过程计算程序设计
在桥梁悬臂施工过程中，往往需要计算各施工状态的结构变形和内力。

由于各施工状态的结构一般都不相同，通常随着施工的进行，结构在不断增大，如果要将各施工状态逐一进行受力计算，并进行变形和内力的累计，则需解决各状态变形和内力的累加问题。

如果对桥梁的最终结构进行单元和结点编号时，就严格按施工的先后顺序来编制。

只要将最终结构的基本数据一次输入后，各施工状态的计算就可按一个一个工况来进行，每一状态的基本数据只要从一次输入的总的基本数据中取各数组的前面部分即可。

结点位移、单元内力和约束内力的累加可这样来进行：开始将他们全部清零，然后在各状态计算后，将这些已存在的结构相应的增量进行叠加即可。

施工过程计算主程序框图如图2.4.8所示。

(5) BSCAS计算分析通用程序流程图2.4.9所示：
(6) 程序功能介绍
以上对BSCAS 大跨连续刚构桥结构静力计算程序系统的原理、构成以及部分功能的力学本质进行了说明，现将程序的主要功能介绍如下：
1) 所计算的桥梁可以是预应力混凝土、钢筋混凝土或砖石结构，也可以是这些结构的组合，结构的不同构件可以采用不同的材料类型。

2) 结构体系可以分阶段形成，各阶段可以具有不同的静力图式；对连续刚构通常使用的悬臂施工工艺，程序能对结构体系进行施工阶段和使用阶段综合分析。

3) 构件可分阶段形成组合截面。

4) 程序采用了两端带刚臂的偏心梁单元，因而考虑了多根杆件交汇于一个结点的结点刚域效应。

5) 在施工阶段，程序考虑了如下各种荷载： ⑴体系调整(包括脱离工作的杆件恢复工作)；
⑵张拉预应力钢束。

预应力荷载的计算有多种方法，如等效荷载法[1]，分段求和法[5]。

本程序考虑到连续刚构节段施工中，单元划分较细，主梁单元长度最长不超过施工节段长，预应力按单元长分段，每一段上作用一对沿力束方向的集中力y N (如图4.4.10)，y N 的大小可取单元范围内的平均有效预加力。

然后将单元两端的y N 移至各自对应结点上，直接作为结点荷载，求出整体坐标系下的结点荷载的三个分量(x P 、y P ，M )，如
i 端的结点力分量为：
⎪
⎩⎪
⎨⎧⋅⋅=⋅=⋅=θ
θ
θ
cos sin cos e N M N P N P y y y y x ⑶拆除临时钢束； ⑷结构自重； ⑸永久集中力荷载； ⑹永久分布力荷载； ⑺永久强迫位移； ⑻临时集中力荷载； ⑼临时分布力荷载； ⑽临时强迫位移； ⑾拆除杆件元荷载；
⑿混凝土收缩徐变影响(正在开发)；
⒀混凝土收缩徐变引起的预应力损失(正在开发)； ⒁斜拉索张拉荷载。

6) 程序的计算内容：
⑴关心截面位移和累计位移； ⑵关心截面内力和累计内力；
2.5 BSCAS 系统结果可靠性验证
2.5.1 刚臂单元可靠性验算
算例1[15]： 如图2.5.1所示的半跨刚架拱中有两个结点块(3，5结点处)，相关杆件轴线在该处不交于一点，且结点处刚度较大，故采用刚臂连结。

实腹段(5，6，7号单元)为变截面，各单元按本段内平均等截面处理。

截面几何特性见表：
容重ρ=253m KN ，弹性模量KPa E 7103⨯=，线膨胀系数50.1-=E α，试对下面两种工况进行结构分析：
(1) 跨中作用100KN 的垂直力，3号单元中点作用50KN 的垂直力(对称)计结构自重；
(2) 各单元温度均匀下降C 15。

表2.5.1 截面几何特性表 序号 单元 2
m
A 4
m
I
序号 单元 2
m
A 4
m
I
1 1，3 0.219 0.0084 4 5 0.429 0.0576
2 2 0.140 0.0019 5 6 0.314 0.0234
3 4
0.245
0.0100
6
7
0.210
0.0080
计算结果比较如表2.5.2所示(限于篇幅，只列出结点1、2、8处支座反力以作对
照)：
表2.5.2
2.5.2 梁单元可靠性验证
算例2： 如图2.5.2所示桥梁结构，已知中跨跨径为40m ，两边跨均为20m ，墩高15m ，主梁25.3m A =，40.2m I =，墩25.3m A =，45.0m I =，主梁和墩均为混凝土材料，容重325m KN =ρ，弹性模量KPa E 7103⨯=。

计算仅当结点1。

限于篇幅，只列出关心截面5(中墩底)和9(中跨跨中)截面的内力做比较，计算结果与SAP2000解见表2.5.3：
表2.5.3 中墩底截面和跨中截面内力计算结果
算例3： 如图2.5.3所示为一带有外伸端，第一跨抗弯刚度分段不等，且既承受结点荷载，又承受非结点荷载的四跨连续梁。

试对该结构进行内力分析，并求各结点的位移。

已知GPa E 210=，44100.4cm I ⨯=。