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

动力分析中平衡方程组的解法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2112211222(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21(25)其中tttttUCaMaUMaKRRˆ102因此，如果质量矩阵是对角形的，则解方程组（11.1）时就不需要进行矩阵的分解，即只需进行矩阵相乘便可求得右端项的有效荷载向量tRˆ，从而利用ii)i(t)i(ttmtRˆU2(26)可得出位移向量的各个分量，其中)i(ttU和)i(tRˆ分别表示向量ttU和tRˆ的第i个分量，而iim是质量矩阵的第i个对角线元素，并且假定0iim。如果对总刚度矩阵和质量矩阵都不需进行三角分解，也就不必形成总体的K和M。此时，求解式(23)可以在单元一级来解决，然后将每个单元的结果累加即可，即)UU(MtUKRRˆtittitiitt212(27)使用式(26)和(27)形式的中心差分法的优点是很明显的，因为它不需要计算总刚度矩阵和总质量矩阵，求解过程基本上是在单元一级上进行，所需要的内存比较少。如果所有相继的单元刚度矩阵和质量矩阵均相同，则该方法就显得更有效，因为这时只需计算或从后备存贮器上连续读出对应于系统中第一个单元的矩阵。至于中心差分法的缺点，必需承认，该过程的效果与对角形质量矩阵中采用和忽略通常依赖于速度的阻尼力有关，若只包合一个对角形阻尼矩阵，则仍然可保持在单元一级上进行求解的优点。从实用上看，只能用于对角形质量矩阵的这个缺点通常是不很严重的，因为可以采用足够精细的有限元离散化来使解有良好的精度。使用中心差分格式的另一个十分重要的考虑是，该积分方法要求时间步长t小于一个临界值crt，可由整个单元分割体的刚度和质量的性质来算出crt。更准确地说，要得到一个有效的解必须ncrTtt(28)其中nT是分析物体的最小周期，n是单元系统的阶。要求使用的时间步长t小于临界时间步长crt的差分格式，例如中心差分法，称为条件稳定的。如果使用一个大于crt的时间步长，则积分是不稳定的，这意味着由数值积分或在计算机上的舍入所导致的误差都会增大，并且在许多情形下会使响应的计算失去意义，积分稳定性的概念是非常重要的。上面我们讨论了中心差分方法的主要缺点：其格式是条件稳定的，许多其它积分方法也会是条件稳定的。由于在结构分析中条件稳定方法的有效使用会受到限制，因此在下面几节我们来考虑一些通常使用的无条件稳定的积分格式。无条件稳定的积分格式的有效性可从积分能取得较好精度这一事实显示出来，时间步长t可以不按式(28)的要求束选取，在许多情形下t可以大于式(28)所允许的t几个量级。所有下面讨论的积分方法都是隐式的，即在求解时要求对刚度矩阵K，或者说得更确切些是对有效刚度矩阵进行三角分解。（b）Houbolt法Houbolt积分格式与前面讨论的中心差分方法是有联系的，因为它也是根据标准的有限差分表示式用位移的各个分量来近似地表示加速度和速度的各个分量。Houbolt积分格式采用了下面的有限差分展开式ttttttttttttttttttUUUUtUUUUUtU222291811614521(29)为了得到在时刻tt的解，现在我们考虑在时刻tt的式(21)(而不是象中心差分法那样考虑在时刻t的情形)，它给出ttttUCtMtRUKCtMt35611222ttttUCtMtUCtMt222311234(30)利用上式求解ttU，需要先知道tU，ttU和ttU2，虽然知道000U,U,U对开始Houbolt积分格式是有用的，但用其它方法计算ttU,U2会更精确，即需要利用特殊的起始方法。用来积分式(21)以得出ttU,U2，是用一个不同的积分格式并尽可能是用以t的几分之一为时间步长，例如用中心差分格式。Houbolt法与中心差分法的一个基本差别在于，刚度矩阵K是作为所求位移ttU的因子出现的。ttKU项的出现是因为在式(21)中取时刻tt的平衡而不像中心差分法取时刻t的平衡，因此Houbolt方法是一个隐式积分格式，而中心差分方法则是一个显式过程。关于积分所用的时间步长t不存在临界的时间步长限制，一般t可以选得比式(28)所给的中心差分法的步长大一些。一个值得注意的特点是，如果忽略质量和阻尼的影响，基于Houbolt方法的逐步积分格式可直接简化为静力分析，但中心差分解法则不能这样做。换言之，若0C和0M，Houbolt求解方法就得出与时间有关的荷裁作用下的静力解。