势能增量驻值原理与切线刚度矩阵的解构规则
陈常松1, 2, 陈政清1, 颜东煌2
(1. 中南大学 土木建筑学院, 湖南 长沙, 410075;
2. 长沙理工大学 桥梁与结构工程学院, 湖南 长沙, 410076)
摘要: 基于有限位移理论应变能密度函数的定义, 利用Kirchhoff应力张量和Green应变张量, 推出了非线性分析中增量形式的势能驻值公式, 并证明了由势能增量驻值原理得到的增量平衡方程形式与由虚位移原理所得的结果完全一致。 然后, 根据势能增量驻值原理并利用泛函驻值条件, 得到了有限位移分析中T.L格式和U.L格式的切线刚度矩阵的解构规则。 利用该解构规则形成切线刚度矩阵的方法简单、 直观, 通过平面梁元例题验证了可利用该方法求得任何单元或结构的切线刚度矩阵。
关键词: 有限位移理论; 势能增量驻值原理; Kirchhoff应力张量; Green应变张量; 切线刚度矩阵; 非线性分析; 解构规则
中图分类号:O242.21 文献标识码:A 文章编号: 1672-7207(2005)05-0892-07
Principle of stationary incremental potential energy and formation rules of tangent stiffness matrix
CHEN Chang-song1, 2, CHEN Zheng-qing1, YAN Dong-huang2
(1. School of Civil and Architectural Engineering, Central South University, Changsha 410075, China;
2. School of Bridge and Structure Engineering, Changsha University of Science and Technology, Changsha 410076, China)
Abstract: Based on definition of strain energy function, increment formula of stationary potential energy of finite displacement theory were derived in terms of Kirchhoff stress tensor and Green strain tensor. Incremental equilibrium equations of nonlinear analysis obtained by former increment formula were proved to be the same as the one obtained by principle of virtual work. Especially, in accordance with the stationary conditions of functional, the formation rules of tangent stiffness matrix may be deduced by principle of stationary incremental potential energy. These rules are simple and straightforward and confirmed by example of the tangent stiffness matrix of updated lagrangian formulation. It is a convenient method to derive tangent stiffness of any members and structures.
Key words: finite displacement theory; principle of stationary incremental potential energy; Kirchhoff stress tensor; Green strain tensor; tangent stiffness matrix; nonlinear analysis; formation rule
在有限位移理论中, 推导切线刚度矩阵时一般从虚功原理出发, 得到各增量状态下的增量平衡方程, 并由此得到单元或结构刚度矩阵[1-6]。 但是, 从势能增量驻值原理角度也可以得到相同的增量平衡方程, 并且从泛函变分的意义上可以给出对应于非线性分析时切线刚度矩阵的组成规则。 该规则物理意义更加明确, 方法简单明了, 适用于各类单元。 文献[7]仅对小变形线弹性桁段有限元进行了分析。 在此, 作者首先利用势能为标量的特点, 分析有限位移变形过程中Kirchhoff应力张量和Green应变张量的增量形式, 得到增量形式的势能驻值原理, 并且证明由势能增量驻值原理得到的T.L格式和U.L格式的增量平衡方程与由虚功原理得到的形式完全相同。 通过例题表明利用势能增量驻值原理及刚度矩阵的组成规则可以很方便地求得单元的切线刚度矩阵。
1 有限变形理论的总势能
在有限变形理论中, 一般假定变形是在等温或是绝热情况下产生的, 并且存在表示应力应变关系的函数[8], 且使得
Sij=Sij(E11, E12, …, E33), (i, j=1, 2, 3)。(1)
其中: Sij表示Kirchhoff应力张量分量, Eij表示Green应变张量分量。
当应变分量足够少时, 可将(1)式展开成Eij的幂级数, 并略去各高阶项, 得到线性化的应力-应变关系:
Sij=CijklEkl。(2)
显然, 由于应力张量和应变张量的对称性, (2)式中各系数之间存在下列关系:
Cijkl=Cjikl=Cijlk=Cklij。(3)
1.1 应变能U
设变形体中的任意微元为一个变形前为单位体积的长方体, 沿着一条加载路径, 其Green应变张量从变形前的零应变到(E11, E12, …, E33), 则应变能密度函数A为:
根据(3)式, 有:
则由上式可知, (4)式中dA为全微分, 且应变能密度函数A与加载途径无关, 仅依赖于最终的应变状态。
因此,
将(5)式代入(4)式可得:
所以, 变形体的应变能为:
其中, 0v表示初始时刻即0时刻物体的体积。
1.2 外力势能V
假设外加力系是保守的, 则外力可从位势函数导出。 设在变形过程中t时刻体积力的势函数为tΦ, 面积力的势函数为tΨ, 则体积力f和面积力T为:
式中: ei为Euler坐标系的坐标矢量; tui为t时刻i坐标方向(i=1, 2, 3)的位移。 在各变量符号中, 左上标表示变形平衡时刻, 右下标表示坐标轴方向。
外力势能V则可表示为:
式中: tv表示t时刻变形体的体积, tsσ表示t时刻变形体的应力条件表面积。
1.3 总势能
变形体的总势能等于应变能U与外力势能V的代数和, 根据(7)和(9)式可得:
1.4 势能增量驻值原理
在非线性分析中采用增量理论求解时, 通常是将变形体的变形加载途径分成若干段平衡状态, 每一平衡状态对应于变形过程中的一个平衡时刻, 设各平衡时刻的状态分别表示为:
Ω(0), Ω(1), …, Ω(N), Ω(N+1), …, Ω(f)。
其中: Ω(0)和Ω(f)分别表示变形体的初始状态和最终状态; Ω(N)和Ω(N+1)分别表示变形过程任一时刻t和t+Δt的中间平衡状态。根据势能驻值原理, 在时刻t的平衡状态下, 在所有的可能位移中只有真实位移使总势能
取驻值, 即
δt=0。(11)
式中,δt表示对t变分。
在t+Δt时刻的平衡状态, 设此时的总势能为t+Δt, 同理有:
δt+Δt=0。(12)
由于总势能为标量, 有:
t+Δt=t+Δ。(13)
式中, Δ为时刻t和t+Δt之间总势能的增量, 显然由(11)~(13)式可得:
δΔ=0。(14)
上式就是有限位移增量理论中的势能增量驻值原理。 这表明: 以t时刻的平衡位形为参考位形时, 在t+Δt时刻平衡状态下, 在所有的可能位移增量中, 只有真实的位移增量使势能增量Δ取驻值。
2 T.L格式和U.L格式下的势能增量驻值原理
在非线性分析中选择参考构形时采用较多的2种格式是T.L格式和U.L格式。 T.L格式以未变形状态的构形(也即初始构形, 用Ω(0)表示; 此时t=0)为参考构形, 所有的物理量都是相对这一参考构形的增量值。 U.L格式是以前一时刻已知的平衡状态的构形Ω(N)为参考构形, 所有的物理量都是相对这一时刻的增量。 每一增量步求解中采用的参考构形都必须重新调整, 以最近的平衡状态作为参考构形。
2.1 T.L格式
根据式(10), 时刻t和t+Δt的总势能分别为:
对于T.L格式, Kirchhoff应力、 Green应变、 质点位移、 体积力和面积力的增量关系分别如下:
由(15)~(17)可得势能增量为:
式(18)为T.L格式的势能增量表达式。在建立单元或整体结构的刚度方程时, 可令δΔ=0, 可得到由结点位移tui及其增量ui表示的驻值方程,按后面的解构规则可建立单元切线刚度矩阵或结构的整体切线刚度矩阵。
可以证明,利用上述势能增量驻值原理可得到与虚位移原理相同的T.L增量平衡方程。
由(18)式可得:
由(3)式可知, (19)式右边第一积分项中的被积函数为:
(19)式方程右边第二、 第三积分项中的被积函数分别为:
将式(20)~(22)代入式(19), 由势能增量驻值条件可得增量平衡方程:
式(23)左边第二、 第三项为外力t+Δtfi和t+ΔtTi在t+Δt时刻所做的虚功δt+ΔtW, 即:
将Green应变增量0Eij分解成线性部分0eij和非线性部分0ηij, 即:
其中:
上式中, 0uk,i表示uk/0xi, t0uk,i表示tuk/0xi, 其他符号规定与之相类似。
将(24)和(25)式代入(23)式可得:
上式与文献[6]和[8]中的T.L格式增量平衡方程完全相同。 此外,
则(26)式可写成:
式中, δΔW=δt+ΔtW-δtW, 为虚功的增量。 式(2)就是有限位移理论常用的虚功增量方程。
2.2 U.L格式
对于U.L格式, Kirchhoff应力和Green应变的增量关系分别如下:
其中: tσij为t时刻的Cauchy应力。 Green应变的增量tEij可分解成线性增量teij和非线性增量tηij:
其中:
设有t时刻变形体中一单位体积的长方体, 其在随后变形过程中的应变能表示为tA, 则根据式(4)和(28)有:
从时刻t到t+Δt, 总势能的增量Δ为:
将式(29)和(30)代入式(31), 利用
上式与文献[6]、 [8]和[9]中的U.L格式增量平衡方程完全相同。 方程右边第二项实际上就是时刻t的外力虚功, 即:
3 切线刚度矩阵的构成规则
下面根据势能增量驻值原理, 从泛函变分驻值条件来分析切线刚度矩阵中元素的构成规则。 设单元位移模式为:
{u}=[A]{U}, {U}=[Φ]{u}e。(34)
其中: {u}=[u1 u2 u3]T, 为单元内任一点的空间位移增量; {U}=[u1 u2 …um]T, 为单元中任一点的广义位移增量; {u}e=[u11 u21 … um1 u12 … umn]T, 为单元结点广义位移增量; m为单元的广义增量位移维数; n为单元的结点数。 则有
{u}=[A][Φ]{u}e=[N]{u}e。(35)
以U.L格式为例, 由(35)式可知, 单元内任一点位移是关于结点位移的线性插值函数, 所以, 利用式(2)和(29), 势能增量驻值方程(33)式各项可表示成结点位移增量与其变分项的线性表达式:
势能增量驻值方程可表示为式(36)的形式, 由δuij的任意性可知, 有下列m×n个驻值条件成立:
上式就是有限单元法中m×n个刚度方程, 因此, 式(36)中的系数k1ij, kl和k2ij, kl就是刚度矩阵中的元素。
根据势能增量驻值条件(37), 可以得到切线单元刚度矩阵的构形规则:
a. 各变分项结点位移增量的序号对应于单元刚度矩阵的行号, 而变分项系数表达式中各结点位移增量序号对应单元刚度矩阵的列号;
b. 设单元结点位移增量的在单元位移增量列向量中的序号为r=(i-1)m+j, s=(k-1)m+l, I1和I2积分式中δuij项的系数表达式中ukl前面的系数置于刚度矩阵的第r行第s列并累加, 将所有的系数累加完成即得切线刚度矩阵;
c. 分析式(36)的I1和I2积分式可知, I1中的系数单独累加将形成线弹性刚度矩阵, 将I2中的系数单独累加将得到几何刚度矩阵。
4 例 题
根据势能增量驻值原理可以很方便地得到各种单元的切线刚度矩阵, 下面推导常用的平面梁元的U.L格式的切线刚度矩阵。
在Euler坐标系的txtOty平面内弯曲的梁, 设梁元在t时刻横截面面积为tA, 抗弯惯性矩为tI, 单元两端结点编号分别为i和j, 结点位移增量为:
{u}e=[ui vi θi uj vj θj]T。
设单元形心轴上任一点tx(t时刻的位形)的位移增量为{u}=[u v]T, 横截面上任一点(tx, ty)的位移增量为{U}=[U V]T, 在Euler-Bernoulli的假设下, 考虑增量过程为小变形过程, 则有:
由于平面梁元为一维单元, 并设Green应变张量只计增量tE11, Kirchhoff应力张量只计tσ11和增量tS11, 则有:
采用多项式位移模式(其中u为一次式, v为三次式)[10, 11], 则有:
将式(40)代入式(39), 经整理得:
将式(41)~(43)代入式(33)左边第一积分式和第二积分式, 分别得:
根据切线刚度矩阵的构成规则, 式(44)和式(45)方程右边δui, δvi, δθi, δuj, δvj和δθj各变分项结点位移的序号对应于单元刚度矩阵的第1~6行, 而变分项系数表达式中ui, vi, θi, uj, vj和θj各结点位移序号对应单元刚度矩阵的第1~6列。 因此, 只要将各变分项系数表达式中的系数项按照其相应的行号和列号置于单元刚度矩阵相应的位置, 就可得到单元的刚度的矩阵。 由式(44)和(45)得到的单元切线刚度矩阵为:
[kT]e=[k0]e+[kσ]e。(46)
其中: [k0]e由(44)得出, 显然矩阵[k0]e为t时刻线弹性刚度矩阵, 其具体表达式见文献[10] 和[11]。 [kσ]e由式(45)可得:
上式表示的矩阵[kσ]e就是U.L.格式增量求解步中的初应力矩阵, 也称几何刚度矩阵。 当式(39)中第二式考虑(U/tx)2/2项时, [kσ]e中元素要更复杂一些, 但可类似推出, 具体讨论见参考文献[11]。
5 结 论
a. 利用应变能密度推导出应用于非线性分析的势能增量驻值原理, 应用于T.L格式和U.L格式, 得到了相应的势能增量驻值公式, 并且证明了势能增量驻值原理与虚功原理的一致性。
b. 采用势能增量驻值原理并利用泛函驻值条件, 可以得到大变形条件下切线刚度矩阵的构成规则, 根据此规则可以推导出单元或复杂结构的切线刚度矩阵。
参考文献:
[1]Mattiasson K, Bengtsson A, Samuelesson A. On the accuracy and efficiency of numerical algorithms for geometrically nonlinear structural analysis[A]. Proceeding of the Europe-US symposium[C]. Trondheim Norway: the Norwegain Institute of Technology, 1985.
[2]许红胜, 周绪红, 舒兴平. 空间框架几何非线性分析的一种新单元[J]. 工程力学, 2003(4): 39-44.
XU Hong-sheng , ZHOU Xu-hong , SHU Xing-ping. A new element for geometric nonlinear analysis of three-dimension steel frames[J]. Engineering Mechanics, 2003(4): 39-44.
[3]舒兴平, 沈蒲生. 空间钢框架结构的非线性全过程分析[J]. 工程力学, 1997(3): 36-44.
SHU Xing-ping, SHEN Pu-sheng. A full range nonlinear analysis for space frames[J]. Engineering Mechanics, 1997(3): 36-44.
[4]辛克贵, 刘钺强, 杨国平. 大跨度斜拉桥恒载非线性静力分析[J]. 清华大学学报(自然科学版), 2002, 42(6): 818-821.
XIN Ke-gui, LIU Yue-qiang, YANG Guo-ping. Nonlinear static analysis of long-span cable-stayed bridges under dead loads[J]. Journal of Tsinghua Univ(Sci & Tech), 2002, 42(6): 818-812.
[5]杨孟刚, 陈政清. 两节点曲线索单元精细分析的非线性有限元法[J]. 工程力学, 2003(1): 42-47.
YANG Meng-gang , CHEN Zheng-qing. Nonlinear analysis of cable structure using a two-node curved cable element of high precision[J]. Engineering Mechanics, 2003(1): 42-47.
[6]王勖成, 邵敏. 有限单元法基本原理和数值方法[M]. 北京: 清华大学出版社, 1997.
WAN Xu-cheng, SHAO Min. Fundamental Theory of Finite Element Method and Numerical Method[M]. Beijing: Tsinghua University Press, 1997.
[7]曾庆元, 杨平. 形成矩阵的“对号入座”法则与桁段有限元法[J]. 铁道学报, 1986(2) : 44-59.
ZENG Qing-yuan, YANG Ping. Rule of ‘number for site’ to form matrix and truss segment finite element method[J]. Journal of Railway, 1986(2) : 44-59.
[8]鹫津久一郎. 弹性和塑性力学中的变分法[M]. 北京: 科学出版社, 1984.
Washizu K. Variational Methods in Elasticity and Plasticity[M]. Beijing: Science Press, 1984.
[9]陈政清, 曾庆元, 颜全胜. 空间杆系结构大挠度问题内力分析的UL列式法[J]. 土木工程学报, 1992, 25(5): 34-44.
CHEN Zheng-qing, ZENG Qing-yuan, YAN Quan-sheng. A UL formulation for internal force analysis of special frame structures with large displacement[J]. China Civil Engineering Journal, 1992, 25(5): 34-44.
[10]普齐米尼斯其 J S. 矩阵结构分析理论[M]. 北京: 国防工业出版社, 1975.
Przemieniecki J S. Theory Matrix Structural Analysis[M]. Beijing: National Defense Industry Press, 1975.
[11]项海帆, 刘光栋. 拱结构的稳定与振动[M]. 北京: 人民交通出版社, 1991.
XIANG Hai-fan, LIU Guang-dong. Stability and Vibration of Arch Structure[M]. Beijing: Communication Press, 1991.
收稿日期: 2004-12-05
基金项目: 国家自然科学基金资助项目(59895410); 湖南省普通高校青年骨干教师培养计划项目([2004]252)
作者简介: 陈常松(1972-), 男, 湖南攸县人, 副教授, 博士研究生, 从事桥梁结构分析研究
论文联系人: 陈常松, 男, 副教授, 博士研究生; 电话: 0731-5219360(O); E-mail: changsongchen@vip.sina.com
