振动问题的有限元分析法
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
离为 x 点的位移为
u(x, t) = N1 (x)u1(t) + N 2 (x)u2 (t)
式中 N1 (x), N 2 (x) 为形函数也称为插值函数，它们应使点的位移满足单元的边界条件，即
u(0,t) = u1(t) ， u(l,t) = u2 (t)
177
PDF created with pdfFactory Pro trial version
可得到
∂u
εx
ε
y
∂x ∂v
∂y ∂w
ε
=
εγ
z x
=
∂u
∂z +
∂v
γ γ
y z
∂y ∂x
∂v ∂∂wz ∂x
+ +
∂w
∂y ∂u
∂z
ε = Bq e = [B1 B2 B3 L]qe
D 为应力应变关系矩阵，又称为弹性矩阵，所以单元体的应变能为
（7-5） （7-6）
(7-7)
U
=
∫∫∫ V
1 2
ε
Tσ
dV
=
∫∫∫ 1 (q e )T BT DBq e
V2
dV
（7-8）
设单元体振动时，受有正比于速度的阻尼力，阻尼系数为 c，则单元体积上所受的阻尼力为
− cq 。单元体上阻尼力所消耗的能量为
常应变状态和刚体位移。单元的形状越复杂，形函数的阶次就越高，单元适应能力就越强。 将式（7-23）写成矩阵形式，有
所以形函数矩阵为
u1
v1
q
=
u v w
=
N1
0
0
0 N1 0
0 0 N1
N2 0 0
0 N2 0
0 0 N2
L L L
有限元法中采用的单元类型有二维、三维等，其形状、大小可以变化，各单元互相之间也容易 连接，因此它能适应复杂的结构，如板、壳、杆等组合的结构物，也适用于各种不同的边界条件。 有限元法中的分析顺序是比较固定的，因此便于计算机计算，并安排标准程序和通用程序。
7.1 单元体的运动方程式
取一单元体，设单元体的动能为 T，应变能为 U，阻尼消耗的能量为Wd ，外力的势能为We 。
∫∫∫ ∫∫∫ Wd = − V
1 cqTq dV = −
2
V
1 c(q e )T N T Nqe dV 2
（7-9）
单元体上所受的外力分为两部分，即体积力 FV = (Fx , Fy , Fz ) 和表面力 FS = (Fx , Fy , Fz ) 。
它们的势能分别为We1 ，We2 。
∫∫∫ ∫∫∫ We1 = q T FV dV = (qe )T N T FV dV
N
=
1 −
x l
x l
形函数是用局部坐标在单元中定义的。如将坐标原点取在单元中点如图（7-2）所示，以 ξ = x l
2
做局部自然坐标，则可得形函数为
N1
(ξ)
=
1
− 2
ξ
，
N
2
(ξ)
=
1
+ 2
ξ
对于二维单元（ ξ, η ）如图（7-3）所示。正方形单元有 4 个节点，则其形函数为
N1
其中子矩阵为
∂N i
∂x
0
0 ∂N i
0
0
∂y
0
Bi
=
∂N
i
∂y
0
∂N i ∂x
∂N i ∂z
0
0
∂N i
∂N
i
∂z ∂y
∂N ∂z
i
0
∂N i ∂x
（7-26） （7-27）
如果形函数是用 ξ, η, ς 局部坐标给出的，则应根据偏微分法则先对 ξ, η, ς 求偏导数，然后再对
∫∫ + (q e )T N T FS d S S
（7-12）
由哈密尔顿原理，将其在时间区间 (t1 , t2 ) 上对 L 积分，并使其变分等于零，考虑到 D 的对称性后，
有
∫ ∫∫∫ ∫∫∫ δL = t2 [δ(qe )T ( BT DB dV )qe − δ(q& e )T ( ρN T N dV )q& e −
t1
V
V
V
− (∫∫∫ N T FV dV ) − (∫∫ N T FS d S)]d t = 0
V
S
（7-16）
令 K eq 、 M eq 、 Ceq 、 Feq 分别表示单元体的刚度矩阵、质量矩阵、阻尼矩阵和载荷矩阵，则
有 由此得到
K eq = ∫∫∫ BT DB dV V
M eq = ∫∫∫ ρN T N dV V
176
PDF created with pdfFactory Pro trial version
7.2.1 形函数矩阵 N
在有限元中，形函数的作用十分重要，因为单元形状和相应的形函数确定以后，其它运算可依 照标准步骤和普遍公式进行。单元上任一点的位移用节点的位移表示为
式中 N i 为形函数，ui ， vi ， wi 是第 i 个节点的位移。形函数 N i 是单元内部坐标的连续函数，它所
满足的条件是：在节点 i 处， N i =1;在其它节点处 N i =0。用它定义的未知量 u、v、w 要保证相邻单
元间的连续性。为了保证收敛于精确解，形函数应包含任意线性项和满足 ΣN i = 1，使单元包含有
于是可得到
N1(0) = 1， N1(l) = 0 ， N2 (0) = 0 ， N2 (l) = 1
图 7-1 一维单元
图 7-2 一维单元
设形函数 Ni (x) = ai + bi x ，包含常数项和线性项。将边界条件代入可得
N1 ( x)
=1−
x l
， N2 (x) =
x l
满足 ΣN i = 0 。所以对于一维单元，形函数矩阵为
t1
V
V
∫∫∫ ∫∫∫ ∫∫ − δ(q& e )T ( cN T N dV )q e − δ(q e )T ( N T FV dV ) − δ(q e )T ( N T FS d S )]d t = 0
V
V
S
（7-13）
175
PDF created with pdfFactory Pro trial version
（i = 1，2，3，4）
图 7-3 二维单元
图 7-4 三维单元
对于三维单元（ ξ, η, ς ）如图（7-4）所示。正六面体 −1 ≤ ξ ≤ + 1 ，−1 ≤ η ≤ + 1 ，−1 ≤ ς ≤ + 1。
将坐标原点取在单元形心上，单元边界是六个平面，ξ = ±1,η = ±1, ς = ±1单元的节点是八个角点，
=
(1
−
ξ)(1 4
−
η)
，
N
2
=
(1 + ξ)(1 − η) 4
， N3
=
(1
−
ξ)(1 4
+
η)
，
N
4
=
(1 + ξ)(1 + η) 4
引入新的变量 ξ0 = ξi ξ ， η0 = ηiη 。其中 ξi , ηi 为节点 i 的坐标，于是上面的四个形函数可合
并表示为
Ni
=
(1 + ξ0 )(1 + η0 ) 4
则其形函数为
178
PDF created with pdfFactory Pro trial version
Ni
=
(1 + ξ0 )(1 + η0 )(1 + ς0 ) 8
（i = 1，2，3，4，5，6，7，8）
7.2.2 应变位移关系矩阵 B
将位移函数矩阵表达式（7-24），代入由弹性力学的应变与位移的关系
∫ ∫∫∫ ∫ ∫∫∫ t2 [δ(q& e )T ( cN T N dV )q e ]d t − t2 [δ(qe )T ( cN T N dV )q& e ]d t
t1
V
t1
V
（7-15）
于是式（7.13）成为
∫ ∫∫∫ ∫∫∫ ∫∫∫ δL = t2 δ(qe )T [( B T DB dV )qe + ( ρN T N dV )q&&e + ( cN T N dV )q& e
第 7 章 振动问题的有限元法
有限元法是力学模型系统上近似的数值计算方法。它先将拟分析的工程结构模型，假想地分割 成有限个单元，组成离散化模型。各个单元之间在单元的外节点处互相连接起来，单元节点可为铰 接或其它形式。然后导出各单元体的运动方程式，最后将这些单元体的运动方程式叠加而得到离散 了的工程结构的有限元运动方程式。因此，有限元法中分析的结构，是一个由有限个单元组成的与 原结构非常接近的离散系统，有限元法是离散化方法。计算所得结果的精确程度取决于单元体的划 分。
应用分布积分公式，上式的第二项有
∫ ∫∫∫ t2 [δ(q& e )T ( t1
ρN T N dV )q& e ]d t
V
∫∫∫ ∫ ∫∫∫ = [δ(qe )T ( V
ρN T N
d
V
)q&
e
]t2 t1
−
t2 t1
[δ(qe )T (
V
ρN T N dV )q&&e ]d t
（7-14）
式中 δqe (t1 ) = 0 ， δq e (t2 ) = 0 ，则只剩下第二项。用同样的方法可得到
q = Nqe
用 u、v、w 表示一点在空间沿 x、y、z 方向的位移，则
u = ΣN iui = N1u1 + N 2u2 + N3u3 + L
v = ΣN ivi = N1v1 + N 2v2 + N3v3 + L
（7-23）
w = ΣN i wi = N1w1 + N 2 w2 + N3 w3 + L
建立拉格朗日函数为
L = T − U − Wd − We
（7-1）
设 q 为单元体中任一点的位移矢量， qe 为单元体上各节点的位移矢量，它是时间 t 的函数。令
单元体中任一点的位移矢量 q 用单元体上各节点的位移矢量 qe 表示为
q = Nqe
式中 N 为形函数矩阵，它是坐标 x、y、z 的函数。 单元体中任一点的位移矢量 q 又可表示为
Ceq = ∫∫∫ cN T N dV V
∫∫∫ ∫∫ Feq = ( N T FV dV ) + ( N T FS d S)
V
S
∫ δL =
t2 t1
δ(qe )T (K eq qe + M eqq&&e + C eqq& e − Feq ) d t = 0
（7-17） （7-18） （7-19） （7-20）
uw21
v2
w2
M
N1 0 0 N 2 0 0 L
N
=
0
N1
0
0 N2 0 L
0 0 N1 0 0 N 2 L
（7-24） （7-25）
在一维单元中，有二个节点，节点的位移为 u1 (t),u2 (t) ，如图（7-1）所示。设单元上距 1 点距
T
=
∫∫∫ V
1 2
ρq& Tq& dV
=
∫∫∫ V
1 2
ρ (q& e )T
N T Nq& e
dV
式中 ρ 为单元体积的质量。
按照弹性力学中的公式，应变与节点位移之间的关系为
ε = Bq e
式中 B 为应变位移关系矩阵，它是几何矩阵，与 t 无关。
单元体上的应力σ 为
σ = Dε = DBq e
有限元法的基本过程为，（1）将结构离散化，即把结构划分成离散的单元。（2）考虑单元的性 质，建立单元的质量矩阵、刚度矩阵、阻尼矩阵、载荷矩阵，推导出单元体的运动方程式。（3）组 合与叠加各单元的质量矩阵、刚度矩阵、阻尼矩阵，得到整个离散系统的运动方程式。（4）解特征 方程，求出频率与振型，或求解动力响应问题和动应力问题。
q = [u(t) v(t) w(t)]T
（7-2） （7-3）
式中 u(t) 、 v(t) 、 w(t) 分别表示该点沿 x、y、z 方向的位移，它们都是时间 t 的函数，如用“·”表 示“ ∂ ”，则
∂t
于是可求得单元体的动能为
q& = [u& v& w& ]T = Nq& e
（7-4）
174
PDF created with pdfFactory Pro trial version
（7-21）
由于单元位移的变分 δ(q e )T 是任取的，所以可由式（7-21）得到单元的运动方程为
M eq q&&e + C eqq& e + K eq qe = Feq
（7-22）
7.2 单元体的特性分析
在有限元法中，将单元节点的位移作为基本量，单元体中各点的位移、应变、应力等量都要表 示为节点位移的函数，单元特性分析就是研究如何得到这些关系。
其整体坐标 x，y，z 求偏导数。 对于平面二维单元，可求得 B 为
V
V
e )T N T FS d S
S
S
于是拉格朗日泛函为
（7-10） （7-11）
∫∫∫ L = 1
2V
[ρ (q e )T N T Nq& e − (q e )T B T DBq e − c(q& e )T N Tq e + 2(qe )T N T FV ]dV
u(x, t) = N1 (x)u1(t) + N 2 (x)u2 (t)
式中 N1 (x), N 2 (x) 为形函数也称为插值函数，它们应使点的位移满足单元的边界条件，即
u(0,t) = u1(t) ， u(l,t) = u2 (t)
177
PDF created with pdfFactory Pro trial version
可得到
∂u
εx
ε
y
∂x ∂v
∂y ∂w
ε
=
εγ
z x
=
∂u
∂z +
∂v
γ γ
y z
∂y ∂x
∂v ∂∂wz ∂x
+ +
∂w
∂y ∂u
∂z
ε = Bq e = [B1 B2 B3 L]qe
D 为应力应变关系矩阵，又称为弹性矩阵，所以单元体的应变能为
（7-5） （7-6）
(7-7)
U
=
∫∫∫ V
1 2
ε
Tσ
dV
=
∫∫∫ 1 (q e )T BT DBq e
V2
dV
（7-8）
设单元体振动时，受有正比于速度的阻尼力，阻尼系数为 c，则单元体积上所受的阻尼力为
− cq 。单元体上阻尼力所消耗的能量为
常应变状态和刚体位移。单元的形状越复杂，形函数的阶次就越高，单元适应能力就越强。 将式（7-23）写成矩阵形式，有
所以形函数矩阵为
u1
v1
q
=
u v w
=
N1
0
0
0 N1 0
0 0 N1
N2 0 0
0 N2 0
0 0 N2
L L L
有限元法中采用的单元类型有二维、三维等，其形状、大小可以变化，各单元互相之间也容易 连接，因此它能适应复杂的结构，如板、壳、杆等组合的结构物，也适用于各种不同的边界条件。 有限元法中的分析顺序是比较固定的，因此便于计算机计算，并安排标准程序和通用程序。
7.1 单元体的运动方程式
取一单元体，设单元体的动能为 T，应变能为 U，阻尼消耗的能量为Wd ，外力的势能为We 。
∫∫∫ ∫∫∫ Wd = − V
1 cqTq dV = −
2
V
1 c(q e )T N T Nqe dV 2
（7-9）
单元体上所受的外力分为两部分，即体积力 FV = (Fx , Fy , Fz ) 和表面力 FS = (Fx , Fy , Fz ) 。
它们的势能分别为We1 ，We2 。
∫∫∫ ∫∫∫ We1 = q T FV dV = (qe )T N T FV dV
N
=
1 −
x l
x l
形函数是用局部坐标在单元中定义的。如将坐标原点取在单元中点如图（7-2）所示，以 ξ = x l
2
做局部自然坐标，则可得形函数为
N1
(ξ)
=
1
− 2
ξ
，
N
2
(ξ)
=
1
+ 2
ξ
对于二维单元（ ξ, η ）如图（7-3）所示。正方形单元有 4 个节点，则其形函数为
N1
其中子矩阵为
∂N i
∂x
0
0 ∂N i
0
0
∂y
0
Bi
=
∂N
i
∂y
0
∂N i ∂x
∂N i ∂z
0
0
∂N i
∂N
i
∂z ∂y
∂N ∂z
i
0
∂N i ∂x
（7-26） （7-27）
如果形函数是用 ξ, η, ς 局部坐标给出的，则应根据偏微分法则先对 ξ, η, ς 求偏导数，然后再对
∫∫ + (q e )T N T FS d S S
（7-12）
由哈密尔顿原理，将其在时间区间 (t1 , t2 ) 上对 L 积分，并使其变分等于零，考虑到 D 的对称性后，
有
∫ ∫∫∫ ∫∫∫ δL = t2 [δ(qe )T ( BT DB dV )qe − δ(q& e )T ( ρN T N dV )q& e −
t1
V
V
V
− (∫∫∫ N T FV dV ) − (∫∫ N T FS d S)]d t = 0
V
S
（7-16）
令 K eq 、 M eq 、 Ceq 、 Feq 分别表示单元体的刚度矩阵、质量矩阵、阻尼矩阵和载荷矩阵，则
有 由此得到
K eq = ∫∫∫ BT DB dV V
M eq = ∫∫∫ ρN T N dV V
176
PDF created with pdfFactory Pro trial version
7.2.1 形函数矩阵 N
在有限元中，形函数的作用十分重要，因为单元形状和相应的形函数确定以后，其它运算可依 照标准步骤和普遍公式进行。单元上任一点的位移用节点的位移表示为
式中 N i 为形函数，ui ， vi ， wi 是第 i 个节点的位移。形函数 N i 是单元内部坐标的连续函数，它所
满足的条件是：在节点 i 处， N i =1;在其它节点处 N i =0。用它定义的未知量 u、v、w 要保证相邻单
元间的连续性。为了保证收敛于精确解，形函数应包含任意线性项和满足 ΣN i = 1，使单元包含有
于是可得到
N1(0) = 1， N1(l) = 0 ， N2 (0) = 0 ， N2 (l) = 1
图 7-1 一维单元
图 7-2 一维单元
设形函数 Ni (x) = ai + bi x ，包含常数项和线性项。将边界条件代入可得
N1 ( x)
=1−
x l
， N2 (x) =
x l
满足 ΣN i = 0 。所以对于一维单元，形函数矩阵为
t1
V
V
∫∫∫ ∫∫∫ ∫∫ − δ(q& e )T ( cN T N dV )q e − δ(q e )T ( N T FV dV ) − δ(q e )T ( N T FS d S )]d t = 0
V
V
S
（7-13）
175
PDF created with pdfFactory Pro trial version
（i = 1，2，3，4）
图 7-3 二维单元
图 7-4 三维单元
对于三维单元（ ξ, η, ς ）如图（7-4）所示。正六面体 −1 ≤ ξ ≤ + 1 ，−1 ≤ η ≤ + 1 ，−1 ≤ ς ≤ + 1。
将坐标原点取在单元形心上，单元边界是六个平面，ξ = ±1,η = ±1, ς = ±1单元的节点是八个角点，
=
(1
−
ξ)(1 4
−
η)
，
N
2
=
(1 + ξ)(1 − η) 4
， N3
=
(1
−
ξ)(1 4
+
η)
，
N
4
=
(1 + ξ)(1 + η) 4
引入新的变量 ξ0 = ξi ξ ， η0 = ηiη 。其中 ξi , ηi 为节点 i 的坐标，于是上面的四个形函数可合
并表示为
Ni
=
(1 + ξ0 )(1 + η0 ) 4
则其形函数为
178
PDF created with pdfFactory Pro trial version
Ni
=
(1 + ξ0 )(1 + η0 )(1 + ς0 ) 8
（i = 1，2，3，4，5，6，7，8）
7.2.2 应变位移关系矩阵 B
将位移函数矩阵表达式（7-24），代入由弹性力学的应变与位移的关系
∫ ∫∫∫ ∫ ∫∫∫ t2 [δ(q& e )T ( cN T N dV )q e ]d t − t2 [δ(qe )T ( cN T N dV )q& e ]d t
t1
V
t1
V
（7-15）
于是式（7.13）成为
∫ ∫∫∫ ∫∫∫ ∫∫∫ δL = t2 δ(qe )T [( B T DB dV )qe + ( ρN T N dV )q&&e + ( cN T N dV )q& e
第 7 章 振动问题的有限元法
有限元法是力学模型系统上近似的数值计算方法。它先将拟分析的工程结构模型，假想地分割 成有限个单元，组成离散化模型。各个单元之间在单元的外节点处互相连接起来，单元节点可为铰 接或其它形式。然后导出各单元体的运动方程式，最后将这些单元体的运动方程式叠加而得到离散 了的工程结构的有限元运动方程式。因此，有限元法中分析的结构，是一个由有限个单元组成的与 原结构非常接近的离散系统，有限元法是离散化方法。计算所得结果的精确程度取决于单元体的划 分。
应用分布积分公式，上式的第二项有
∫ ∫∫∫ t2 [δ(q& e )T ( t1
ρN T N dV )q& e ]d t
V
∫∫∫ ∫ ∫∫∫ = [δ(qe )T ( V
ρN T N
d
V
)q&
e
]t2 t1
−
t2 t1
[δ(qe )T (
V
ρN T N dV )q&&e ]d t
（7-14）
式中 δqe (t1 ) = 0 ， δq e (t2 ) = 0 ，则只剩下第二项。用同样的方法可得到
q = Nqe
用 u、v、w 表示一点在空间沿 x、y、z 方向的位移，则
u = ΣN iui = N1u1 + N 2u2 + N3u3 + L
v = ΣN ivi = N1v1 + N 2v2 + N3v3 + L
（7-23）
w = ΣN i wi = N1w1 + N 2 w2 + N3 w3 + L
建立拉格朗日函数为
L = T − U − Wd − We
（7-1）
设 q 为单元体中任一点的位移矢量， qe 为单元体上各节点的位移矢量，它是时间 t 的函数。令
单元体中任一点的位移矢量 q 用单元体上各节点的位移矢量 qe 表示为
q = Nqe
式中 N 为形函数矩阵，它是坐标 x、y、z 的函数。 单元体中任一点的位移矢量 q 又可表示为
Ceq = ∫∫∫ cN T N dV V
∫∫∫ ∫∫ Feq = ( N T FV dV ) + ( N T FS d S)
V
S
∫ δL =
t2 t1
δ(qe )T (K eq qe + M eqq&&e + C eqq& e − Feq ) d t = 0
（7-17） （7-18） （7-19） （7-20）
uw21
v2
w2
M
N1 0 0 N 2 0 0 L
N
=
0
N1
0
0 N2 0 L
0 0 N1 0 0 N 2 L
（7-24） （7-25）
在一维单元中，有二个节点，节点的位移为 u1 (t),u2 (t) ，如图（7-1）所示。设单元上距 1 点距
T
=
∫∫∫ V
1 2
ρq& Tq& dV
=
∫∫∫ V
1 2
ρ (q& e )T
N T Nq& e
dV
式中 ρ 为单元体积的质量。
按照弹性力学中的公式，应变与节点位移之间的关系为
ε = Bq e
式中 B 为应变位移关系矩阵，它是几何矩阵，与 t 无关。
单元体上的应力σ 为
σ = Dε = DBq e
有限元法的基本过程为，（1）将结构离散化，即把结构划分成离散的单元。（2）考虑单元的性 质，建立单元的质量矩阵、刚度矩阵、阻尼矩阵、载荷矩阵，推导出单元体的运动方程式。（3）组 合与叠加各单元的质量矩阵、刚度矩阵、阻尼矩阵，得到整个离散系统的运动方程式。（4）解特征 方程，求出频率与振型，或求解动力响应问题和动应力问题。
q = [u(t) v(t) w(t)]T
（7-2） （7-3）
式中 u(t) 、 v(t) 、 w(t) 分别表示该点沿 x、y、z 方向的位移，它们都是时间 t 的函数，如用“·”表 示“ ∂ ”，则
∂t
于是可求得单元体的动能为
q& = [u& v& w& ]T = Nq& e
（7-4）
174
PDF created with pdfFactory Pro trial version
（7-21）
由于单元位移的变分 δ(q e )T 是任取的，所以可由式（7-21）得到单元的运动方程为
M eq q&&e + C eqq& e + K eq qe = Feq
（7-22）
7.2 单元体的特性分析
在有限元法中，将单元节点的位移作为基本量，单元体中各点的位移、应变、应力等量都要表 示为节点位移的函数，单元特性分析就是研究如何得到这些关系。
其整体坐标 x，y，z 求偏导数。 对于平面二维单元，可求得 B 为
V
V
e )T N T FS d S
S
S
于是拉格朗日泛函为
（7-10） （7-11）
∫∫∫ L = 1
2V
[ρ (q e )T N T Nq& e − (q e )T B T DBq e − c(q& e )T N Tq e + 2(qe )T N T FV ]dV