波动方程有限元解有关理论

1.1控制方程

在经典的线弹性理论当中，线弹性均匀介质的运动方程为

⋅

⋅=+∂∂i i j

ij u f x ρρσ (3,2,1=i ) (1-1)

其中，应力张量ij σ和应变张量kl ε满足Hooke 定律，即满足本构方程:

kl ijkl ij c εσ= (1-2)

应变张量kl ε可通过位移矢量i u 确定，从而给出几何方程:

)(2

1,,k l l k kl u u +=ε (k,l=1,2,3) (1-3)

各向同性线弹性材料的弹性特征可以用下列材料常数进行描述:弹性模量E,剪切 模量μ、泊松比v 、体积模量K 和Lame 常数λ等。这些常数中任意两个独立，而其它 几个常数能够用它们间接进行表示。在各向异性的情形下，则需要采用21个相互独立 的弹性常数ijkl C ，进行表示。 介质表面S 上的边界条件为

u i i S x t x u t x u ∈=),,(),( (1-4)

σσS x t x p n t x i j ij ∈=),,(),( (1-5)

由方程(2-1)-(2-3)这三个线弹性动力学的控制方程所描述的线弹性动力问题的初始 条件为

V x x u t x u i i ∈=),(),(00 (1-6)

V x x v t x u i i ∈=⋅

),(),(00 (1-7)

其中，)(x u i ，)(0x v i ，),(t x u i 和),(t x p i 是已知量。

由弹性动力学问题解的唯一性定理可知:若弹性体(体积为Y ，表面为S)的解能 满足方程(1-1)-式(1-7)，则其位移场、应力场和应变场的解答是唯一的。 将三个控制方程合并，可以得到用位移表示的运动方程:

⋅

⋅=+i i lj k ijkl u f u c ρρ, (3,2,1=i ) (1-8)

方程(1-8)所示的就是著名的Navie-Cauchy 运动方程。而在各向同性线弹性的情形 下，则有

)(jk il jl ik kl ij ijkl c δδδδμδλδ++= (1-9)

因此，各向同性材料的本构关系也为表示为

⎩⎨

⎧=≠=+=)

(,1)

(,0,2j i j i ij ij kk ij ij δμεελδσ (1-10)

1.2有限元方程

现在我们考虑一个任意形状的封闭区域，并将其用一定形式的网格划分为相应的有 限单元，各单元的单元矩阵可由如下所示的虚功方程求得。

0][][][][=Γ-Ω-Ω-Ω⎰⎰⎰⎰

Γ

Ω

Ω

Ω

d f u d f u d f u d B T

A T D T T δδδσδε (1-11)

其中，B f 为作用在边界上的外力矩阵; A f 和D f 为惯性力矩阵，其形式由式(1-12)和(1-13)

⋅

⋅-=u f A ρ (1-12)

⋅

-=u c f D (1-13)

给出。

在单元内任一点的位移e

u ，可以用单元形函数N 和单元节点位移u 近似表示为

Nu u e = (1-14)

由此，可相应的得到

⋅

⋅=u N u e

(1-15) ⋅

⋅⋅⋅=u N u e

(1-16)

u B u LN δδδε== (1-17) DBu =σ (1-18)

f N f T B = (1-19)

其中，D 为弹性矩阵，L 为微分算子。

将以上各个经过插值后的量(式(1-12)-式(1-19))，代入式(1-11)中，可得

0][][][][=Γ-Ω+Ω+Ω⎰⎰⎰⎰

Γ

⋅

⋅Ω

⋅

Ω

Ω

fd N u d u N N u d u cN N u DBud B u T T

T

T

T

T

T

T

δρδδδ(1-20)

由于虚位移u δ是任意的，因此整理可得

⎰⎰⎰⎰Γ

⋅

⋅Ω

⋅

⋅Ω

⋅

⋅

Ω

Γ=Ω+Ω+Ωfd N u d u N N u d u cN N u DBd B T

T

T

T

ρ (1-21) 方程(1-21)也可写为有限元方程的标准形式:

F u M u C Ku =++⋅

⋅⋅ (1-22)

其中，K 为刚度矩阵，C 为一致性阻尼矩阵(对于无阻尼问题，该矩阵为零阵)，M 为

一致性质量矩阵，F 为外力矩阵，分别可表示为

⎰Ω

Ω=DBd B K T (1-23)

Ω=⎰Ω

cNd N C T (1-24)

Ω=⎰Ω

Nd N M T ρ (1-25)

Γ=⎰Γ

fd N F T (1-26)

这就表明，对于该封闭区域，通过单元网格划分可以确定其总体刚度矩阵和总体质 量矩阵。因此，只要给定问题的边界条件和初始条件，就能对所列出的线性方程组(1-22) 进行求解，求解得到的结果即为该问题的解。由于材料的粘滞系数难以确定(一般只能 用实验方法测定)，在工程计算中通常采用Rayleigh 阻尼法确定阻尼矩阵，即假定阻尼 矩阵为刚度矩阵和质量矩阵的线性组合

K M C βα+= (1-27)

特别要注意的是，对于某些无阻尼问题，有限元方程的形式仅为

F u M Ku =+⋅

⋅ (1-28)

1.3求解算法

在建立了动力分析的有限元方程后，通过给定的边界条件就可以求解出体系的动力 反应。在对有限元方程进行求解时，常使用直接积分方法进行计算。

直接积分方法是指在积分运动方程之前不进行方程形式的变换，而直接进行数值积 分。这类方法的特点是求解域内仅在离散时间点(而非任意时刻)能够满足运动方程的要求。而在离散时间点间假定了节点加速度、节点速度和节点位移的函数形式，从而建立起在相邻离散时间点时各节点反应值间的相互关系，由此每一时刻的反应值便可逐步计算出来。直接积分方法可分为两类，分别是隐式方法和显式方法。 1有限元方程的隐式解法

隐式方法的研究历史较长，且大部分为无条件稳定，例如，线性加速度法、常平均加速度法、Newmark 法、Wilson-θ法等。隐式方法需要对藕联的线性方程组进行直接求解，通过求得的反应值建立下一迭代步的线性方程组(以下一时刻反应值作为未知数)，进而求解并得到该方程组中未知的反应值，以此类推。隐式方法的计算精度较高，算法稳定性好，但对于自由度数目较大或非线性程度较高的问题，隐式方法所需的计算量很大，计算成本较高。 (1) 线性加速度法

线性加速度法假设在相邻的离散时间点P 和P+1间，节点m 的加速度分量mi u ⋅

⋅以线 性规律变化，如下所示。

t t t t t u u u u P

mi

P mi P

mi mi ∆+≤≤-∆-+

=⋅⋅+⋅⋅⋅⋅⋅

⋅τττ),()(1

(1-37) 由式(1-37)积分并整理可得节点m 在P+1时刻的速度和位移的表达式

t u t u u d u u u P mi P mi P

mi t

t t

mi P

mi P mi ∆+∆+=+=+⋅⋅⋅⋅⋅∆+⋅

⋅⋅+⋅⎰1

12

121)(ττ (1-38)

1

2216

131)(+⋅⋅⋅⋅⋅∆+⋅

+∆+∆+∆+=+=⎰

P mi P mi P mi P

mi

mi t t t

P mi

P mi

u t u t u t u d u u u

ττ (1-39) 将式(1-38)和(1-39)代入P+1时刻的平衡方程，即可整理得到P+1时刻的节点加速