梁杆结构几何非线性有限元的数值实现方法

转换矩阵 ， 并将它 代 入 连续体的 三维 虚功 增量 方 程， 引入相应的变形假定(如梁截面变形的直线假定) 可以使虚功方程的三维积分降维，最终导出以广义 位移为未知量的虚功方程。经有限元离散后，可编 程实现数值计算。1979 年 Bathe 等[1]提出了应用三 维 虚功 增量 方程建立 三维 梁结构 的 大 位移 问题的 UL 型增量有限元分析方法。 Bathe 建立的有限元平 衡方程仍然是三维积分形式，需要通过三维数值积 分计算梁单元刚度矩阵，输出结果是积分点应力， 而不是梁截面内力， 因此只能适用矩形截面和圆(环) 形截面组成的梁结构。在计算效率和适用范围上都 无法满足工程实用要求。1989 年，作者发现 Bathe 方法没有注意到空间梁结构的广义位移本质上是一
工
程
力
学
43
维的。依据梁的这一特点，作者从三维虚功增量方 程出发，推出了适用梁结构大位移分析的一维虚功 增量方程，进而由它建立起了一套完整的三维梁杆 结构几何非线性分析算法，仍然满足严格的几何非 线性理论要求。 与 Bathe 方法相比， 新方法可完全消 除数值积分，计算效率是不言而喻的，并且在计算 格式上与已有的三维梁结构线性有限元程序完全一 致，直接计算截面内力，从而可应用于任何复杂截 面的实际工程结构。由于新方法直接计算出结构设 计人员需要的梁截面内力弯矩，因此这一方法当时 命名为空间杆系结构大挠度问题内力分析方法，并 于 1992 年在《土木工程学报》发表 ，经过改进后 1993 年在《Computers and Structures》上发表[3]。作 者据此方法编制的 NACS 程序从 1990 年开始至今， 已 直 接 应 用 于 自 我 国第 一 座 悬索桥以来 多 个 悬索 桥、斜拉桥的工程计算。后来，杨孟刚与陈政清又 进一步将上述方法扩展到索结构，提出了两种曲线 索单元的非线性分析方法并开发了计算程序，论文 于 2003 年分别发表在《工程力学》 与《土木工程 学报》[5]上。作者后来发现，Wen R K 等[6]、Yang Y B 等[7]、Spillers W R[8]等也在 1980 年代建立了类似 的方法。2002 年以来，Yang Y B 等[9
同一个固定的空间坐标系， 变形体在 Ω (0) 、Ω ( N ) 和 Ω ( N + 1) 状态下所对应的空间与边界分别为 V0 与 S0 、 Vt 与 St 、 Vt + ∆t 与 St +∆t 。 以 Ω ( Ν ) 状态为参考时， 变形体在 t+Δt 时刻的 克西霍夫应力为：
t + ∆t ij tT
第 31 卷第 6 期 2014 年 6 月
Vol.31 No.6 June 2014
工
程
力
学 42
ENGINEERING MECHANICS
文章编号：1000-4750(2014)06-0042-11
梁杆结构几何非线性有限元的数值实现方法
陈政清
(湖南大学土木工程学院，长沙 410082)
摘 要：梁杆结构几何非线性有限元方法主要包括两个部分， 建立虚功方程和实现数值求解。 该文运用对比方法， 分析了采用 UL 型增量理论的梁杆结构几何非线性有限元法求解过程与连续体求解过程的主要不同点，特别是论 述了确定加载步末的内力状态的重要性和方法。 关键词：梁杆结构；几何非线性；有限元法；增量法；非线性分析； 中图分类号：TU375.1 文献标志码：A doi: 10.6052/j.issn.1000-4750.2013.05.ST08
Vt
t + ∆t
Vt
bkφk β dvt + ∫∫
t +∆t
St
pk dst
(17) (18)
A1′ = ∫ Cijkl ∆∗ε kl δ∆ ∗ε ij dvt
Vt
ψ β = ∫∫∫ tσ ijφi β , j dvt
′ = ∫ tσ ij δ∆∗ηij dvt A2
Vt
式 中 特别用重 积分 符号 强 调了计算 每 个元 素 的 积 分区域，用左上标表明了他们所处的状态( 时刻) ，
Abstract:
The geometrically nonlinear finite element method for beam structures consists of two parts: the
development of virtual work equations and the corresponding numerical implementation. In this study, the solution procedure of the geometrically nonlinear finite element analysis using the Updated Lagrangian (UL) incremental approach for beam structures and three-dimensional continuum is presented, and major differences in the solution process are identified. Particularly, emphasis has been placed on the importance and the methodology of determining internal forces at the end of each load step. Key words: beam structures; geometric nonlinear; finite element method; incremental method; nonlinear analysis 大 跨度的 梁杆 索 结构 的施 工 过 程 以及使用寿 命期内在特殊荷载( 如台风、地震)作用下的变形与 内力都具有明显的几何非线性特性，需要精确的几 何非线性分析方法。长期以来，梁杆索结构几何非 线性问题的研究都是分门别类，逐个问题地加以研 究解决。例如，研究压弯杆件导致梁柱理论、研究 悬索桥面内荷载导致挠度理论等。研究的方法也往 往随具体问题而变，无一定法则可循。20 世纪 80 年代末，计算力学总结出了完善的连续体三维虚功 增量方程。利用它可以导出连续体任何一种非线性 单元模式，并建立相应的增量求解方法。对于梁、 板、壳等这一类以中线或者中面广义位移为未知量 的结构单元，只要建立结构质点位移与广义位移的
= tσ ij + ∆∗T ij = ∆∗ Eij
(1) (2)
而它在 t+Δt 时刻柯西应变就等于其增量：
t + ∆t t Eij
式中， ∆ Eij 为：
∗
∆∗ Eij = ∆∗ε ij + ∆∗ηij 1 ∆∗ε ij = (∆ui ,j + ∆u j ,i ) 2 1 ∆∗ηij = ∆uk ,i ∆uk ,j 2
式中各项分别为： (6)
t E kαβ = ∫∫∫ Cijklφiα ,k φ j β ,l dvt Vt Vt G kαβ = ∫∫∫ tσ ijφkα ,iφk β , j dvt
将以上关系代入虚功原理可以导出 UL 法的三维虚 (7) (8) (9) (10) (11)
t
t
t + ∆t
Pβ = ∫∫∫
E 以上诸式中， kαβ 所有变量都以 Ω ( N ) 状态为参考。 G 称为弹性刚度矩阵， kαβ 称为几何刚度矩阵或应力
′= A3
t +∆t tP Vt
′ = ∫ tσ ij δ∆∗ε ij dvt A4
′导 在有限元离散中， A1′ 导出线性刚度矩阵 K 0 ， A2 ′ 是 Ω ( N + 1) 出初应力阵或称几何刚度矩阵 Kσ ，A3 ′ − A4 ′ 仍是荷载增量与不平衡荷载 状态外力虚功，A3
[2]
1
1.1
连续体几何非线性有限元及其 数值实现方法
连续体大变形的有限元增量平衡方程 连续体大变形问题的求解一般采用增量方法。
连续体 大 变形的有限元 增量平衡 方程 的推 导需 要 张量的字母下标记法和连续介质力学的知识，本文 不加证明的介绍基本推导过程，需要了解的读者可 以参考文献[12]。 如图 1 所示， 一个连续体在已知体力荷载 b j 和 面力荷载 p j 作用下经历如下的大变形过程：
(3) (4) (5)
44
工
程
力
学
E G [ t kαβ ]{∆qα } = {t+ ∆t Pβ − tψ β } + t kαβ
仍然假定变形体的应变增量是小应变，应 力应变增量关系可以记为：
(14) (15) (16)
′ ∆∗ε kl ∆∗T ij = Cijki
功增量方程如下： ′ = A3 ′ − A4 ′ A1′ + A2 式中：
NUMERICAL IMPLEMENTATION OF GEOMETRICALLY NONLINEAR FINITE ELEMENT METHOD FOR BEAM STRUCTURES
CHEN Zheng-qing
(College of Civil Engineering, Hunan University, Changsha 410082, China)
[4]
又提出了直
接 由 梁元 的 大 旋 转 特性 导 出 非线性平衡 方程 的 方 法。近年来，有学者提出了通过强迫单元服从梁壳 结构变形的假定直接由多个连续体单元的结合建立 梁壳结构非线性单元的方法，称为 CB 梁单元和 CB 壳单元，可见庄茁翻译的文献[11]。 梁杆 索 结构几何非线性有限元方法主要包括 两个部分： 1) 建立虚功方程。 将以质点位移为基本未知量 的 虚功方程 转 化为 以 广义位移为 基 本未知 量的 并 按广义坐标积分的虚功方程，以梁元为例，要转化 为 以 梁 截 面的 位移 与 转 角 为未知 量的 沿梁 的 中性 轴积分的虚功方程； 2) 实现数值求解。 其中的关键步骤是确定加载 步末的几何与内力状态。由于广义内力至今没有严 格的变形前后转换关系式，因此确定加载步末的精 确内力状态只能通过平衡迭代来实现，需要设计比 连续体大变形分析更为复杂的增量法求解过程。此 外，在大旋转条件下跟踪单元局部坐标系在三维空 间的变化也比线性小变形分析复杂，要在每个加载 步求解过程中应用欧勒角理论处理。因此，梁杆结 构几何非线性有限元 的 实现要 难 于三维弹 性 体 几
———————————————
收稿日期：2013-05-01；修改日期：2014-03-06 基金项目：国家自然科学基金项目(91215302) 作者简介： 陈政清(1947―)， 男， 湖南湘潭人， 教授， 博士， 湖南大学风工程研究中心主任， 主要从事结构振动与控制研究(E-mail: zqchen@hnu.edu.cn).
的虚功之和。在 UL 法中，每一步的参考状态都要 改变，包括了每个变量的当前坐标改变和积分域改 变，因为式(7)中的积分与微分都是对 Ω ( N ) 状态的 空间坐标 X i ( 它也是 Ω ( N + 1) 状态的物质坐标) 进 行的。 任何线性有限元 使用的 单 元 模式 都可以用来 构建非线性有限元方程，例如图 2 所示的最简单的 12 个自由度的 4 面体单元，它有 4 个节点，每个节 点有三个位移自由度。一般地，设块体单元内的位 移 增量 ∆ui 由单 元 插 值 函 数 φiα 和 节 点 位移 增量
Ω (0) ，Ω (1) ，L ，Ω ( N ) ，Ω ( N + 1) ，L ，ຫໍສະໝຸດ Baidu ( f ) 设 0、 t、 t+Δt 分别表示 Ω (0) 、 Ω (N ) 和 Ω (N + 1)
状态所对应的时刻，并用左下标表示参考状态，左 上标表示当前状态。例如， 0t ui 表示以 Ω (0) 为参考
t ij 时， Ω (N ) 状态的位移； t + ∆0 T 表示以 Ω (0) 状态为 参考的在 Ω (N + 1) 状态的 Kirchhoff 应力。相对于
―10]
何非线性有限元。在有限元发展史上，梁杆结构小 变形的线性分析是最先实现的，但大变形的非线性 分析反而在三维连续体之后，也许原因正在这里。 现 已 发 表 的 关 于梁杆结构几何非线性分析的 学 术 论文，很少或没有论述数值实现方法，难以指导开 发相应的计算软件。 梁杆索结构都适合运用 UL 型增量理论，本文 运用对比方法，分析了采用 UL 型增量理论的梁杆 结构几何非线性有限元法求解 过 程与 连续体 求解 过程的主要不同点，特别是论述了确定加载步末的 内力状态的重要性和方法。