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

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

三大类变量(,)iut、(,)ijt和(,)ijt是坐标位置(,,)xyz和时间t的函数，

一般将其记为()()()iijijuttt。 基本方程 (1) 平衡方程 利用达朗贝尔原理将惯性力和阻尼力等效到静力平衡方程中，有

,()()()()0ijjiiitbtutut (1)

其中为密度，为阻尼系数。 (2) 几何方程

,,1()(()())2ijijjitutut (2) (3) 物理方程 ()()ijijklkltDt (3)

其中ijklD为弹性系数矩阵。 (4) 边界条件 位移边界条件()BCu为，

()()iiutut 在uS上 (4)

力的边界条件()BCp为， ()()ijjitnpt 在pS上 (5)

初始条件 0(,0)()iiutu (6)

0(,0)()iiutu (7) 虚功原理 基于上述基本方程，可以写出平衡方程及力边界条件下的等效积分形式，

,()()0pijjiiiijjiSuubudnpdA

 (8)

对该方程右端第一项进行分部积分，并应用高斯-格林公式，整理得， ()()0pijklijkliiiiiiiiSDuuuudbudpudA (9)

有限元分析列式 单元的节点位移列阵为，

111222()[(),(),(),(),(),()(),(),()]etkkkUtutvtwtutvtwtutvtwt (10)

单元内的插值函数为， (,)()()etutNUt (11)

其中()N为单元的形状函数矩阵，与相应的静力问题单元的形状函数矩阵完全相同，为单元中的几何位置坐标。 基于上面的几何方程和物理方程及(11)式，将相关的物理量表达为节点位移的关系，有，

(,)[](,)[]()()()()eetttutNUtBUt (12)

(,)()()()()eetttDDBUtSUt (13)

(,)()()etutNUt (14)

(,)()()etutNUt (15)

将(12)-(15)供稿到虚功方程(9)中，有， [()()()()]()0eeeeeeeTetttttMUtCUtKUtRtUt (16)

由于()etUt具有任意性，消去该项并简写有， eeeeettttUCUKUR (17)

其中，

eeTMNNd

(18)

eeTCNNd

(19) eeTKBDBd

(20)

eM为单元质量矩阵，eC为单元阻尼矩阵，eK为单元刚度矩阵。同样，将单元

的各个矩阵进行组装，可形成系统的整体有限元方程，即， MUCUKUR (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到Tt的解。在求解时，把时间全程T划分为几个相等的时间区间t(即n/Tt)，所用的积分格式是在时刻t,0，T,,tt,t,,t2上确定方程的近似解。由于计算下一个时刻的解的算法要考虑到前面各个时刻的解，因此假定在时刻t,,t,0的解为已知，来推导出求时刻tt的解的算法。计算时刻tt的解对于计算自此以后t的时刻上的解是有代表意义的，这样就可建立用来计算在所有离散时间点上解的一般算法。 （a） 中心差分法 若把式(21)的平衡关系看作是一个常系数常微分方程组，便可以用任一有限差分表达式通过位移来近似表示加速度和速度。因此，在理论上，许多不同的有限差分表达式均可使用。但是，我们要求求解格式必须是有效的，这样便只需考虑少数几种计算格式。对某些问题求解是非常有效的一个过程是中心差分法，这个方法假定



tttttttttttUUtUUUUtU





21212





(22)

将式(22)代入t时刻的式(21)，可得 ttttttUCtMtUMtKRUCtMt



2112211

222 (23)

从式(23)我们可以求出ttU。应该注意，ttU的解是基于利用在时刻t的平衡条件。因此，该积分过程称为显式积分方法，且这样的积分格式在逐步解法中不需要对(有效)刚度矩阵进行分解。另一方面，以后所考虑的Houbolt，Wilson及Newmark方法，要利用在tt上的平衡条件，因而称为隐式积分方法。

另外还应注意到，应用中心差分法时，ttU的计算包含有tU和ttU，因此，

计算在时刻t的解，必需用一个具体的起始过程。由于000U,U,U都是已知的，由关系式(22)可求

02002UtUtUUt (24) 具体计算步骤为 A． 初始计算 1． 形成刚度矩阵K、质量矩阵M和阻尼矩阵C。

2． 计算初始值000U,U,U。 3． 选取时间步长t，要求crtt（临界值）。 4． 计算系数201ta，ta211，022aa，231aa。 5． 计算0300UaUtUUt。 6． 形成有效质量矩阵CaMaMˆ10。 7． 对Mˆ作三角分解：TLDLMˆ B． 每一时间步长内的计算 1． 计算在时刻t的有效荷载： tttttUCaMaUMaKRRˆ102。

2． 计算时刻tt的位移：tttTRˆULDL。 3． 必要时，按照式（11.3）计算时刻t速度和加速度。 假设所考虑的系统没有物理阻尼，即C是零矩阵，在这种情形下式(23)可简化为

tttRˆMUt2

1 (25)

其中 tttttUCaMaUMaKRRˆ102

因此，如果质量矩阵是对角形的，则解方程组（11.1）时就不需要进行矩阵的分

解，即只需进行矩阵相乘便可求得右端项的有效荷载向量tRˆ，从而利用



ii)i(t)i(ttmt

RˆU2

(26)

可得出位移向量的各个分量，其中)i(ttU和)i(tRˆ分别表示向量ttU和tRˆ的第i个分量，而iim是质量矩阵的第i个对角线元素，并且假定0iim。 如果对总刚度矩阵和质量矩阵都不需进行三角分解，也就不必形成总体的K和M。此时，求解式(23)可以在单元一级来解决，然后将每个单元的结果累加即可，即

)UU(MtUKRRˆtittitiitt212 (27)