有限元软件中通用梁单元的统一形成方法
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
最为复杂的部分。由于工程实际需求，必须要考虑多 种形式的载荷，例如：集中力、均布载荷、三角形分 布载荷、梯形分布载荷等，且载荷作用在纯弹性梁、 组合刚域-弹性梁、组合弹簧-弹性梁。按梁的形式 分工况形成载荷无疑是一种可能性，但我们在这里 再一次从简支梁岀发，给出一个简洁的方案。限于篇 幅，仅叙述集中力的作用下的等效载荷，其他分布载 荷情况可考虑为单位集中力的等效结点力的积分。
迹） ai
aj
\/L 1 -1/L 0 l/£ 0 -1/L 1
Wi
Or
=Su(e)
(1)
WjLeabharlann 其中，f（e） = ｛Qr Mi Qj Mj｝是梁单元两端的作 用力，而= ｛Mi M,｝是相应简支梁两端的作用 力矩，Qi = -Qj = + Mj\L反映了梁端作用力 的平衡关系。由此可知，通过这两个关系矩阵，任意 梁单元的变形和梁端作用力均可由一个相应的简支 梁表达岀来，当然我们也可以由简支梁得到普通梁 单元的刚度矩阵与等效载荷向量。
如（£）=处（0 + （1 - £）妙 + gw』 m =比（。+违兰
其中w（e）, 0（£）, ws（e）, 0s（e）分别以自然坐标e = x/l 表示梁单元和相应简支梁的挠度与转角。w（o, 0（e） 满足T-梁的控制方程，当且仅当地（£）, %（£）满足 该方程。记"⑷=｛妙,仇,吗,切｝为梁单元两端的挠 度与截面转角，“0 = ｛ai,aj｝为简支梁两端截面转 角，则有
第3期
陈璞等：有限元软件中通用梁单元的统一形成方法
449
较为方便地获得刚度矩阵的方法仍然是利用最 -(e)
小余能原理，修改不含转角弹簧的柔度矩阵Fs o例 如，用最小余能原理得到图4所示含右端转动弹簧 的柔度矩阵
耳2)=几+ 0
(11)
1/fc
其逆即为相应的刚度矩阵。
用柔度法可得两端的转角
〔警2 8 =L
第3期
陈璞等：有限元软件中通用梁单元的统一形成方法
447
在中外文献中关于梁单元的形成方法层出不穷,
这里不一一列举，但大部分文献或只考虑某一特定
的单元形式，或形成方案过于繁琐，不宜在工程软件
中获得应用。本文从工程实用的角度出发，着重讨论
前5类功能在长为L的二维平面梁单元中的统一实 现，读者不难将其推广至三维情形。
1.2等截面简支梁的柔度矩阵与刚度矩阵
记g = x/L为梁单元的自然坐标，（*），为（*）对 e的导数。由于简支梁是静定结构，我们可以通过考
虑剪切效应的最小余能原理求简支T-梁的柔度矩 阵
PT”
1 广 F'Tjv
~eF^ + t]0 GAS dC
其中的弯矩取线性插值
F(g) = (l—0
(3)
满足梁平衡关系。对于等截面梁，简支T_梁的柔 度矩阵可解析给岀
Abstract Based on the authors5 experience in the development of various finite element software, this paper presents the method of forming the stiffness matrix and the corresponding element load vectors for beam elements to meet the different engineering demands, focusing on the unified formation of Euler and Timoshenko beam elements and variable section beam elements based on the flexibility method, and the unified formation of the element stiffness matrix and the load vector under different beam end conditions. This paper provides a reference for developers in the engineering finite element software industry and to help the teaching of the finite element methods. Key words finite element analysis, Euler beam element, Timoshenko beam element, variable section beam, flexibility method
W
■ LN2 LN4 '
L0
LM2 LM4
(6)
其中
-1 J n2~
'h2'
1 +讥
£(1-£) 11 J/
m2
1 J 'H'2
i-ell
M4
1 + M H： +叶
1.4梁端刚域
在结构力学中梁-柱的计算是所谓的中心线计 算。在实际工程中，梁在其与柱的交汇区域内几乎没 有弯曲变形，这一部分称为刚域(图2),其长度一般 可取为柱的中线到柱外缘距离的75%o
-(e) 1 L — c —c
fs =
L -d L-d
= PT 於 e)
综合起来，中心线位置的简支梁柔度矩阵为
耳(Ze)X = pp-(e) pTT
(10)
1.5梁端部分刚接
梁端有时与柱等其他构件的连接不紧密，不能 完全传递力矩，此时一般称为部分刚接，它也可在弹 塑性计算中模拟理想塑性较。在有限元计算中，部分 刚接用梁端的附加扭转弹簧来模拟。
严)=ST厝)
另一方面，我们还可以利用最小势能原理获得 梁单元与相应简支梁单元的刚度矩阵变换关系以及 简支梁单元的刚度矩阵同。记h2（C） = £（£ - I）2, 日4（£） = C2（C - 1）为［0, 1］上标准的三次Hermite插
448
力学与实践
2021年第43卷
值函数中关于导数的两个，则满足无分布力简支T梁方程的位移插值形函数可以写为
一旦计入不变形的刚域，梁的弹性变形部分缩 短，刚度增大。如图3,仍记梁柱中心线处简支梁的 转角为潦)=｛%旳｝,又记刚域外缘弹性梁的转角 为五丁 = ｛乱,必｝,那么由线性化的几何关系得到
基于最小势能原理的单元刚度矩阵为
卅)=S'1"善 /
+
1 r1 1 L GAs(Ms - 7V^)t(Ms -
" S=
k
L ［12
刚接，与1.2节得到的刚度相质：
由以上讨论可以得到梁单元刚度矩阵形成方式: (e)
首先计算弹性简支梁的柔度矩阵Fs ,如果存在
刚域或部分刚接的情形，则需进行端头刚域与集
中弹簧修哼，获得黒叫最后通过变换K(Q =
ST ［耳S计算梁单元的刚度矩阵。
2等效结点力
2.1梁跨中载荷的等效结点力 梁载荷的等效结点力大约是工程有限元软件中
( fl M,tM' (13)
其中
b&
£ W a/L
Mp(E) = P
(14)
a(l -g), £ > a/L
是简支梁在集中力p作用下的弯矩。对于等截面简 支梁［5］可得到
对于等截面直梁，当2 T 0时，T竿
K
L
101.表示右端不能传递力矩，退化为完全较接;
当g T +OO时，K$)T 孕 2 1 ,退化为完全
关键词 有限元分析，Euler梁单元，Timoshenko梁单元，变截面梁，柔度法
中图分类号：0242.21 文献标识码：A doi: 10.6052/1000-0879-20-258
UNIFIED FORMULATION OF BEAM ELEMENTS IN FINITE ELEMENT SOFTWARE*12*)
本文于2020-06-16收到。 1) 国家重点研发计划(2018YFC0809700)和国家自然科学基金(11472014, 11521202, 11672362)资助项目。 2) 陈璞，教授，主要研究方向为计算力学。E-mail: chenpu@
引用格式：陈璞，杜晖，孙树立等.有限元软件中通用梁单元的统一形成方法.力学与实践，2021, 43(3): 446-452 Chen Pu, Du Hui, Sun Shuli, et al. Unified formulation of beam elements in fin计e element software. Mechanics in Engineering^ 2021, 43(3): 446-452
梁单元无疑是有限元分析中应用最为广泛的单 元之一。为了满足各种工程需求，梁单元常常可能要 求具备以下特性或功能：
(1) 剪切变形； (2) 端头刚域，梁-柱交汇处的区域； (3) 部分刚接，梁-柱交汇等处的刚度损失；
(4) 端头较接，梁-柱交汇等处的构造或刚度缺 失；
(5) 截面变化，例如桥梁等大跨结构； (6) 部分偏心，例如与墙柱相连的外立面梁； (7) 几何刚度，轴向力引起的附加刚度； (8) 轴线为曲线，例如：拱。
StK^S
(7)
由式(7)中的积分，可以获得与式(5)完全相同 的刚度矩阵。式⑷和式⑸互为逆，表明两端截面 转角、力矩和插值所得的跨内挠度、弯矩之间的准确 关系。
当耳=0时，剪切刚度无穷大，T-梁回归到E梁。作为£的函数，Ns就是Hermite插值H, Ms 是其微商H o此时，式(5)也退化到E-梁的刚度 矩阵。
1梁单元刚度矩阵的统一形成方式
以Timoshenko梁（T -梁）讨论刚度矩阵以及 载荷向量的形成方案，Euler梁（E-梁）可以作为 T-梁的一个特例。为了处理变截面梁、端头刚域以 及端头部分刚接，我们以柔度法为讨论重点。
1.1梁单元与相应的简支梁单元 考虑一个置于Oxy平面内长为L的梁单元,
其主轴与x轴重合，左端与原点O重合。与标 准的有限元教材［1-2］稍许不同，考虑到单元的通用 性，工程软件中梁单元的4x4刚度矩阵（这里 只考虑弯曲部分）往往不是按相应的公式直接形成 的，而是通过简支梁或悬臂梁的2x2刚度矩阵 扩张得到的。为此，对照梁单元与同跨度的简支梁 （图1）,它们的挠度和转角满足变换
L 4 +耳 一2 +耳
12EI 一2 +耳 4 +耳
⑷
其中无量纲参数T）= 12EI/（GASL2）,表示弯曲与剪 切刚度之比。其刚度矩阵可由式（4）求逆获得
（e）= El ［4 + r? 2呵
⑸
s £（1 + 7,） ［2-rj 4 + n_
其中，S是两组结点位移之间的几何关系矩阵，其转 置是相应的两端外力关系矩阵，即
CHEN Pu*® DU Hui* SUN Shuli* FU Xiangrong+
* (Department of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing 100871, China) t (College of Water Resources and Civil Engineering, China Agricultural University, Beijing 100083, China)
式(2) 一般可通过较高阶的Gauss积分或Lobatto积分获得。我们的经验是5点Gauss积分的精 度已足够工程应用。变截面梁的刚度矩阵如采用位 移法得到，则可能隐含较大的误差⑷。
现在，长为L-c-d的弹性简支梁两端的截面 转角於)与外力矩A)满足
其中弹性部分的柔度矩阵按式(2)计算，但长 度替换为L-c-do由刚域外缘与中心线的力矩平 衡得出
1.3变截面梁单元
当截面沿中性轴变化时，位移的三次Hermite 插值形函数不再满足梁的平衡方程，因此最小势能 原理所导岀的刚度矩阵不再是简支梁的精确刚度矩 阵。但作为静定结构，简支梁的弯矩分布仍然是线性 的，插值函数式(3)仍然给出结构力学意义下的精确 内力，式(2)是精确的表达式。由此，最小余能原理 所导出的柔度矩阵还是简支梁的精确柔度矩阵。在 这种意义下，变截面梁单元的刚度矩阵可以从柔度 矩阵求逆而准确得到。
第43卷第3期
力学与实践
2021年6月
有限元软件中通用梁单元的统一形成方法"
陈璞*，2)杜晖*孙树立*傅向荣t
*(北京大学工学院力学与工程科学系，北京100871) t(中国农业大学水利与土木工程学院，北京100083)
摘要本文根据作者参与多个有限元软件开发的经验，讨论并给出了满足多种工程需求的通用梁单元的 刚度矩阵及其单元载荷向量形成方法，特别是基于柔度法的Euler梁单元、Timoshenko梁单元和变截面梁单 元的统一形成方式，以及不同梁端条件下的单元刚度矩阵与载荷向量的统一形成方式，便于工程有限元软件开 发人员参考借鉴，并希望能对有限元方法的教学实践有所帮助。
迹） ai
aj
\/L 1 -1/L 0 l/£ 0 -1/L 1
Wi
Or
=Su(e)
(1)
WjLeabharlann 其中，f（e） = ｛Qr Mi Qj Mj｝是梁单元两端的作 用力，而= ｛Mi M,｝是相应简支梁两端的作用 力矩，Qi = -Qj = + Mj\L反映了梁端作用力 的平衡关系。由此可知，通过这两个关系矩阵，任意 梁单元的变形和梁端作用力均可由一个相应的简支 梁表达岀来，当然我们也可以由简支梁得到普通梁 单元的刚度矩阵与等效载荷向量。
如（£）=处（0 + （1 - £）妙 + gw』 m =比（。+违兰
其中w（e）, 0（£）, ws（e）, 0s（e）分别以自然坐标e = x/l 表示梁单元和相应简支梁的挠度与转角。w（o, 0（e） 满足T-梁的控制方程，当且仅当地（£）, %（£）满足 该方程。记"⑷=｛妙,仇,吗,切｝为梁单元两端的挠 度与截面转角，“0 = ｛ai,aj｝为简支梁两端截面转 角，则有
第3期
陈璞等：有限元软件中通用梁单元的统一形成方法
449
较为方便地获得刚度矩阵的方法仍然是利用最 -(e)
小余能原理，修改不含转角弹簧的柔度矩阵Fs o例 如，用最小余能原理得到图4所示含右端转动弹簧 的柔度矩阵
耳2)=几+ 0
(11)
1/fc
其逆即为相应的刚度矩阵。
用柔度法可得两端的转角
〔警2 8 =L
第3期
陈璞等：有限元软件中通用梁单元的统一形成方法
447
在中外文献中关于梁单元的形成方法层出不穷,
这里不一一列举，但大部分文献或只考虑某一特定
的单元形式，或形成方案过于繁琐，不宜在工程软件
中获得应用。本文从工程实用的角度出发，着重讨论
前5类功能在长为L的二维平面梁单元中的统一实 现，读者不难将其推广至三维情形。
1.2等截面简支梁的柔度矩阵与刚度矩阵
记g = x/L为梁单元的自然坐标，（*），为（*）对 e的导数。由于简支梁是静定结构，我们可以通过考
虑剪切效应的最小余能原理求简支T-梁的柔度矩 阵
PT”
1 广 F'Tjv
~eF^ + t]0 GAS dC
其中的弯矩取线性插值
F(g) = (l—0
(3)
满足梁平衡关系。对于等截面梁，简支T_梁的柔 度矩阵可解析给岀
Abstract Based on the authors5 experience in the development of various finite element software, this paper presents the method of forming the stiffness matrix and the corresponding element load vectors for beam elements to meet the different engineering demands, focusing on the unified formation of Euler and Timoshenko beam elements and variable section beam elements based on the flexibility method, and the unified formation of the element stiffness matrix and the load vector under different beam end conditions. This paper provides a reference for developers in the engineering finite element software industry and to help the teaching of the finite element methods. Key words finite element analysis, Euler beam element, Timoshenko beam element, variable section beam, flexibility method
W
■ LN2 LN4 '
L0
LM2 LM4
(6)
其中
-1 J n2~
'h2'
1 +讥
£(1-£) 11 J/
m2
1 J 'H'2
i-ell
M4
1 + M H： +叶
1.4梁端刚域
在结构力学中梁-柱的计算是所谓的中心线计 算。在实际工程中，梁在其与柱的交汇区域内几乎没 有弯曲变形，这一部分称为刚域(图2),其长度一般 可取为柱的中线到柱外缘距离的75%o
-(e) 1 L — c —c
fs =
L -d L-d
= PT 於 e)
综合起来，中心线位置的简支梁柔度矩阵为
耳(Ze)X = pp-(e) pTT
(10)
1.5梁端部分刚接
梁端有时与柱等其他构件的连接不紧密，不能 完全传递力矩，此时一般称为部分刚接，它也可在弹 塑性计算中模拟理想塑性较。在有限元计算中，部分 刚接用梁端的附加扭转弹簧来模拟。
严)=ST厝)
另一方面，我们还可以利用最小势能原理获得 梁单元与相应简支梁单元的刚度矩阵变换关系以及 简支梁单元的刚度矩阵同。记h2（C） = £（£ - I）2, 日4（£） = C2（C - 1）为［0, 1］上标准的三次Hermite插
448
力学与实践
2021年第43卷
值函数中关于导数的两个，则满足无分布力简支T梁方程的位移插值形函数可以写为
一旦计入不变形的刚域，梁的弹性变形部分缩 短，刚度增大。如图3,仍记梁柱中心线处简支梁的 转角为潦)=｛%旳｝,又记刚域外缘弹性梁的转角 为五丁 = ｛乱,必｝,那么由线性化的几何关系得到
基于最小势能原理的单元刚度矩阵为
卅)=S'1"善 /
+
1 r1 1 L GAs(Ms - 7V^)t(Ms -
" S=
k
L ［12
刚接，与1.2节得到的刚度相质：
由以上讨论可以得到梁单元刚度矩阵形成方式: (e)
首先计算弹性简支梁的柔度矩阵Fs ,如果存在
刚域或部分刚接的情形，则需进行端头刚域与集
中弹簧修哼，获得黒叫最后通过变换K(Q =
ST ［耳S计算梁单元的刚度矩阵。
2等效结点力
2.1梁跨中载荷的等效结点力 梁载荷的等效结点力大约是工程有限元软件中
( fl M,tM' (13)
其中
b&
£ W a/L
Mp(E) = P
(14)
a(l -g), £ > a/L
是简支梁在集中力p作用下的弯矩。对于等截面简 支梁［5］可得到
对于等截面直梁，当2 T 0时，T竿
K
L
101.表示右端不能传递力矩，退化为完全较接;
当g T +OO时，K$)T 孕 2 1 ,退化为完全
关键词 有限元分析，Euler梁单元，Timoshenko梁单元，变截面梁，柔度法
中图分类号：0242.21 文献标识码：A doi: 10.6052/1000-0879-20-258
UNIFIED FORMULATION OF BEAM ELEMENTS IN FINITE ELEMENT SOFTWARE*12*)
本文于2020-06-16收到。 1) 国家重点研发计划(2018YFC0809700)和国家自然科学基金(11472014, 11521202, 11672362)资助项目。 2) 陈璞，教授，主要研究方向为计算力学。E-mail: chenpu@
引用格式：陈璞，杜晖，孙树立等.有限元软件中通用梁单元的统一形成方法.力学与实践，2021, 43(3): 446-452 Chen Pu, Du Hui, Sun Shuli, et al. Unified formulation of beam elements in fin计e element software. Mechanics in Engineering^ 2021, 43(3): 446-452
梁单元无疑是有限元分析中应用最为广泛的单 元之一。为了满足各种工程需求，梁单元常常可能要 求具备以下特性或功能：
(1) 剪切变形； (2) 端头刚域，梁-柱交汇处的区域； (3) 部分刚接，梁-柱交汇等处的刚度损失；
(4) 端头较接，梁-柱交汇等处的构造或刚度缺 失；
(5) 截面变化，例如桥梁等大跨结构； (6) 部分偏心，例如与墙柱相连的外立面梁； (7) 几何刚度，轴向力引起的附加刚度； (8) 轴线为曲线，例如：拱。
StK^S
(7)
由式(7)中的积分，可以获得与式(5)完全相同 的刚度矩阵。式⑷和式⑸互为逆，表明两端截面 转角、力矩和插值所得的跨内挠度、弯矩之间的准确 关系。
当耳=0时，剪切刚度无穷大，T-梁回归到E梁。作为£的函数，Ns就是Hermite插值H, Ms 是其微商H o此时，式(5)也退化到E-梁的刚度 矩阵。
1梁单元刚度矩阵的统一形成方式
以Timoshenko梁（T -梁）讨论刚度矩阵以及 载荷向量的形成方案，Euler梁（E-梁）可以作为 T-梁的一个特例。为了处理变截面梁、端头刚域以 及端头部分刚接，我们以柔度法为讨论重点。
1.1梁单元与相应的简支梁单元 考虑一个置于Oxy平面内长为L的梁单元,
其主轴与x轴重合，左端与原点O重合。与标 准的有限元教材［1-2］稍许不同，考虑到单元的通用 性，工程软件中梁单元的4x4刚度矩阵（这里 只考虑弯曲部分）往往不是按相应的公式直接形成 的，而是通过简支梁或悬臂梁的2x2刚度矩阵 扩张得到的。为此，对照梁单元与同跨度的简支梁 （图1）,它们的挠度和转角满足变换
L 4 +耳 一2 +耳
12EI 一2 +耳 4 +耳
⑷
其中无量纲参数T）= 12EI/（GASL2）,表示弯曲与剪 切刚度之比。其刚度矩阵可由式（4）求逆获得
（e）= El ［4 + r? 2呵
⑸
s £（1 + 7,） ［2-rj 4 + n_
其中，S是两组结点位移之间的几何关系矩阵，其转 置是相应的两端外力关系矩阵，即
CHEN Pu*® DU Hui* SUN Shuli* FU Xiangrong+
* (Department of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing 100871, China) t (College of Water Resources and Civil Engineering, China Agricultural University, Beijing 100083, China)
式(2) 一般可通过较高阶的Gauss积分或Lobatto积分获得。我们的经验是5点Gauss积分的精 度已足够工程应用。变截面梁的刚度矩阵如采用位 移法得到，则可能隐含较大的误差⑷。
现在，长为L-c-d的弹性简支梁两端的截面 转角於)与外力矩A)满足
其中弹性部分的柔度矩阵按式(2)计算，但长 度替换为L-c-do由刚域外缘与中心线的力矩平 衡得出
1.3变截面梁单元
当截面沿中性轴变化时，位移的三次Hermite 插值形函数不再满足梁的平衡方程，因此最小势能 原理所导岀的刚度矩阵不再是简支梁的精确刚度矩 阵。但作为静定结构，简支梁的弯矩分布仍然是线性 的，插值函数式(3)仍然给出结构力学意义下的精确 内力，式(2)是精确的表达式。由此，最小余能原理 所导出的柔度矩阵还是简支梁的精确柔度矩阵。在 这种意义下，变截面梁单元的刚度矩阵可以从柔度 矩阵求逆而准确得到。
第43卷第3期
力学与实践
2021年6月
有限元软件中通用梁单元的统一形成方法"
陈璞*，2)杜晖*孙树立*傅向荣t
*(北京大学工学院力学与工程科学系，北京100871) t(中国农业大学水利与土木工程学院，北京100083)
摘要本文根据作者参与多个有限元软件开发的经验，讨论并给出了满足多种工程需求的通用梁单元的 刚度矩阵及其单元载荷向量形成方法，特别是基于柔度法的Euler梁单元、Timoshenko梁单元和变截面梁单 元的统一形成方式，以及不同梁端条件下的单元刚度矩阵与载荷向量的统一形成方式，便于工程有限元软件开 发人员参考借鉴，并希望能对有限元方法的教学实践有所帮助。