有限元动力学分析方程及解法

动力分析中平衡方程组的解法1前言 描述结构动力学特征的基本力学变量和方程与静力问题类似，但所有的变量都是时间的函数。

基本变量三大类变量(,)i u t ξ、(,)ij t εξ和(,)ij t σξ是坐标位置(,,)x y z ξ和时间t 的函数，一般将其记为()()()i ij ij u t t t εσ。

基本方程(1) 平衡方程利用达朗贝尔原理将惯性力和阻尼力等效到静力平衡方程中，有,()()()()0ij j i i i t b t u t u t σρν+--= (1)其中ρ为密度，ν为阻尼系数。

(2) 几何方程）,,1()(()())2ij i j j i t u t u t ε=+ (2)(3) 物理方程 ()()ij ijkl kl t D t σε= (3)其中ijkl D 为弹性系数矩阵。

(4) 边界条件位移边界条件()BC u 为，()()i i u t u t = 在u S 上 (4)力的边界条件()BC p 为，()()ij j i t n p t σ= 在p S 上 (5)初始条件0(,0)()i i u t u ξξ== (6){ 0(,0)()i i u t u ξξ== (7)虚功原理基于上述基本方程，可以写出平衡方程及力边界条件下的等效积分形式，,()()0pij j i i i ij j i S u u b u d n p dA δσρνδσΩ∏=---+Ω+-=⎰⎰ (8)对该方程右端第一项进行分部积分，并应用高斯-格林公式，整理得，()()0pijkl ij kl i i i i i i i i S D u u u u d b u d p u dA εδερδνδδδΩΩ-++Ω-Ω+=⎰⎰⎰ (9) 有限元分析列式单元的节点位移列阵为，111222()[(),(),(),(),(),()(),(),()]e t k k k U t u t v t w t u t v t w t u t v t w t = (10)单元内的插值函数为，(,)()()e t u t N U t ξξ= (11)%其中()N ξ为单元的形状函数矩阵，与相应的静力问题单元的形状函数矩阵完全相同，ξ为单元中的几何位置坐标。

基于上面的几何方程和物理方程及(11)式，将相关的物理量表达为节点位移的关系，有，(,)[](,)[]()()()()e e t t t u t N U t B U t εξξξξ=∂=∂= (12)(,)()()()()e e t t t D DB U t S U t σξεξξ=== (13)(,)()()e t u t N U t ξξ= (14)(,)()()e t u t N U t ξξ= (15)将(12)-(15)供稿到虚功方程(9)中，有，[()()()()]()0e e e e e e eT e t t t t t M U t C U t K U t R t U t δδ∏=++-= (16)由于()e t U t δ具有任意性，消去该项并简写有，e e e e et t t t U C U KU R ++= (17)其中，】e e T M NNd ρΩ=Ω⎰ (18) ee T C N Nd νΩ=Ω⎰ (19) ee T K B DBd Ω=Ω⎰ (20)e M 为单元质量矩阵，e C 为单元阻尼矩阵，e K 为单元刚度矩阵。

同样，将单元的各个矩阵进行组装，可形成系统的整体有限元方程，即，MU CU KU R ++= (21)其中M ，C 和K 分别是系统的质量、阻尼和刚度矩阵，R 是外荷载向量，U ,U和U 分别是有限元分割体的加速度、速度和位移向量。

方程(21)是通过考虑在时刻t 的静力平衡而推导出来的。

对静力或动力分析的选择(即在分析中是考虑或忽略与速度及加速度有关的力)，一般取决于工程上的判断，其目的在于减少所需要的分析工作量。

但是，应该认识到，一个静力分析的假定，应该有理由说明它是正确的，否则，分析的结果就是无意义的。

确实，在非线性分析中，采用忽略惯性力和阻尼力的假定，可能严重到难以求得甚至无法求得解答。

在数学上，方程(21)是一个二阶线性微分方程组，原则上可用求解常系数微分方程组的标准过程来求得方程组的解。

但是，如果矩阵的阶数很高，则采用求解一般微分方程组的过程可能要付出很高的费用，除非特别利用系数矩阵K ，C 和M 的特殊性质。

因此，在实用的有限元分析中，主要对几种有效的方法感兴趣，下面将集中介绍这几种方法。

我们所考虑的基本过程，可分为两种求解方法：直接积分法和振型叠加法。

初看起来，这两种方法似乎完全不同，但事实上它们有着密切的关系，至于选择这种或那种方法，只取决于它们的数值效果。

2直接积分法在直接积分中对方程(21)是逐步地进行数值积分的，“直接”的意思是，进行数值积分前没有进行把方程变为另一种形式的变换。

实质上，直接积分是基于下面的两个想法，第一个想法是只在相隔t ∆的一些离散的时间区间上而不是试图在任一时刻t 上满足方程(21)即包含有惯性力和阻尼力作用的(静力)平衡是在求解区间上的一些离散时刻点上获得的。

因此，似乎在静力分析中使用过的所有求解方法，在直接积分法中或许也能有效地使用；第二个想法是假定位移、速度和加速度在每一时间区间t ∆内变化。

下面假设分别用000U ,U ,U 来表示初始时刻)t (0=的位移、速度和加速度向量为已知，要求出方程(21)从0=t 到T t =的解。

在求解时，把时间全程T 划分为几个相等的时间区间t ∆(即n /T t =∆)，所用的积分格式是在时刻t ,∆0，T ,,t t ,t ,,t ∆+∆2上确定方程的近似解。

由于计算下一个时刻的解的算法要考虑到前面各个时刻的解，因此假定在时刻t ,,t , ∆0的解为已知，来推导出求时刻t t ∆+的解的算法。

计算时刻t t ∆+的解对于计算自此以后t ∆的时刻上的解是有代表意义的，这样就可建立用来计算在所有离散时间点上解的一般算法。

（a ） %（b ） 中心差分法若把式(21)的平衡关系看作是一个常系数常微分方程组，便可以用任一有限差分表达式通过位移来近似表示加速度和速度。

因此，在理论上，许多不同的有限差分表达式均可使用。

但是，我们要求求解格式必须是有效的，这样便只需考虑少数几种计算格式。

对某些问题求解是非常有效的一个过程是中心差分法，这个方法假定{}{}t t t t t t t t t t t U U t U U U U t U ∆+∆-∆+∆-+-∆=+-∆=21212 (22)将式(22)代入t 时刻的式(21)，可得t t t t t t U C t M tU M t K R U C t M t ∆-∆+⎪⎭⎫ ⎝⎛∆-∆-⎪⎭⎫ ⎝⎛∆--=⎪⎭⎫ ⎝⎛∆+∆2112211222 (23)从式(23)我们可以求出t t U ∆+。

应该注意，t t U ∆+的解是基于利用在时刻t 的平衡条件。

因此，该积分过程称为显式积分方法，且这样的积分格式在逐步解法中不需要对(有效)刚度矩阵进行分解。

另一方面，以后所考虑的Houbolt ，Wilson θ及Newmark 方法，要利用在t t ∆+上的平衡条件，因而称为隐式积分方法。

另外还应注意到，应用中心差分法时，t t U ∆+的计算包含有t U 和t t U ∆-，因此，计算在时刻t ∆的解，必需用一个具体的起始过程。

由于000U ,U ,U 都是已知的，由关系式(22)可求02002U t U t U U t ∆+∆-=∆- (24) 具体计算步骤为A ． 初始计算1．·2． 形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。

3．计算初始值000U ,U ,U 。

4．选取时间步长t ∆，要求cr t t ∆≤∆（临界值）。

5．计算系数201ta ∆=，t a ∆=211，022a a =，231a a =。

6．计算0300U a U t U U t +∆-=∆-。

7．形成有效质量矩阵C a M a M ˆ10+=。

8．对M ˆ作三角分解：T LDL M ˆ=B ． 每一时间步长内的计算1．计算在时刻t 的有效荷载：()()tt t t t U C a M a U M a K R R ˆ∆-----=102。

2．计算时刻t t ∆+的位移：tt t T R ˆU LDL =∆+。

3．.4． 必要时，按照式（）计算时刻t 速度和加速度。

假设所考虑的系统没有物理阻尼，即C 是零矩阵，在这种情形下式(23)可简化为t t t R ˆMU t=∆∆+21 (25) 其中()()tt t t t U C a M a U M a K R R ˆ∆-----=102 因此，如果质量矩阵是对角形的，则解方程组（）时就不需要进行矩阵的分解，即只需进行矩阵相乘便可求得右端项的有效荷载向量tR ˆ，从而利用 ⎪⎪⎭⎫ ⎝⎛∆=∆+ii )i (t )i (t t m t R ˆU2 (26)可得出位移向量的各个分量，其中)i (t t U ∆+和)i (t R ˆ分别表示向量t t U ∆+和tR ˆ的第i 个分量，而ii m 是质量矩阵的第i 个对角线元素，并且假定0>ii m 。

如果对总刚度矩阵和质量矩阵都不需进行三角分解，也就不必形成总体的K 和M 。

此时，求解式(23)可以在单元一级来解决，然后将每个单元的结果累加即可，即)U U (M t U K R R ˆt it t i t i i t t 212-∆--=∑∑∆- (27) 使用式(26)和(27)形式的中心差分法的优点是很明显的，因为它不需要计算总刚度矩阵和总质量矩阵，求解过程基本上是在单元一级上进行，所需要的内存比较少。

如果所有相继的单元刚度矩阵和质量矩阵均相同，则该方法就显得更有效，因为这时只需计算或从后备存贮器上连续读出对应于系统中第一个单元的矩阵。

|至于中心差分法的缺点，必需承认，该过程的效果与对角形质量矩阵中采用和忽略通常依赖于速度的阻尼力有关，若只包合一个对角形阻尼矩阵，则仍然可保持在单元一级上进行求解的优点。

从实用上看，只能用于对角形质量矩阵的这个缺点通常是不很严重的，因为可以采用足够精细的有限元离散化来使解有良好的精度。

使用中心差分格式的另一个十分重要的考虑是，该积分方法要求时间步长t ∆小于一个临界值cr t ∆，可由整个单元分割体的刚度和质量的性质来算出cr t ∆。

更准确地说，要得到一个有效的解必须πncr T t t =∆≤∆ (28)其中n T 是分析物体的最小周期，n 是单元系统的阶。

要求使用的时间步长t ∆小于临界时间步长cr t ∆的差分格式，例如中心差分法，称为条件稳定的。