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

1.1控制方程在经典的线弹性理论当中，线弹性均匀介质的运动方程为iijijufx(3,2,1i)(1-1)其中，应力张量ij和应变张量kl满足Hooke定律，即满足本构方程:klijklijc(1-2)应变张量kl可通过位移矢量iu确定，从而给出几何方程:)(21,,kllkkluu(k,l=1,2,3)(1-3)各向同性线弹性材料的弹性特征可以用下列材料常数进行描述:弹性模量E,剪切模量、泊松比v、体积模量K和Lame常数等。这些常数中任意两个独立，而其它几个常数能够用它们间接进行表示。在各向异性的情形下，则需要采用21个相互独立的弹性常数ijklC，进行表示。介质表面S上的边界条件为uiiSxtxutxu),,(),((1-4)Sxtxpntxijij),,(),((1-5)由方程(2-1)-(2-3)这三个线弹性动力学的控制方程所描述的线弹性动力问题的初始条件为Vxxutxuii),(),(00(1-6)Vxxvtxuii),(),(00(1-7)其中，)(xui，)(0xvi，),(txui和),(txpi是已知量。由弹性动力学问题解的唯一性定理可知:若弹性体(体积为Y，表面为S)的解能满足方程(1-1)-式(1-7)，则其位移场、应力场和应变场的解答是唯一的。将三个控制方程合并，可以得到用位移表示的运动方程:iiljkijklufuc,(3,2,1i)(1-8)方程(1-8)所示的就是著名的Navie-Cauchy运动方程。而在各向同性线弹性的情形下，则有)(jkiljlikklijijklc(1-9)因此，各向同性材料的本构关系也为表示为)(,1)(,0,2jijiijijkkijij(1-10)1.2有限元方程现在我们考虑一个任意形状的封闭区域，并将其用一定形式的网格划分为相应的有限单元，各单元的单元矩阵可由如下所示的虚功方程求得。0][][][][dfudfudfudBTATDTT(1-11)其中，Bf为作用在边界上的外力矩阵;Af和Df为惯性力矩阵，其形式由式(1-12)和(1-13)ufA(1-12)ucfD(1-13)给出。在单元内任一点的位移eu，可以用单元形函数N和单元节点位移u近似表示为Nuue(1-14)由此，可相应的得到uNue(1-15)uNue(1-16)uBuLN(1-17)DBu(1-18)fNfTB(1-19)其中，D为弹性矩阵，L为微分算子。将以上各个经过插值后的量(式(1-12)-式(1-19))，代入式(1-11)中，可得0][][][][fdNuduNNuducNNuDBudBuTTTTTTTT(1-20)由于虚位移u是任意的，因此整理可得fdNuduNNuducNNuDBdBTTTT(1-21)方程(1-21)也可写为有限元方程的标准形式:FuMuCKu(1-22)其中，K为刚度矩阵，C为一致性阻尼矩阵(对于无阻尼问题，该矩阵为零阵)，M为一致性质量矩阵，F为外力矩阵，分别可表示为DBdBKT(1-23)cNdNCT(1-24)NdNMT(1-25)fdNFT(1-26)这就表明，对于该封闭区域，通过单元网格划分可以确定其总体刚度矩阵和总体质量矩阵。因此，只要给定问题的边界条件和初始条件，就能对所列出的线性方程组(1-22)进行求解，求解得到的结果即为该问题的解。由于材料的粘滞系数难以确定(一般只能用实验方法测定)，在工程计算中通常采用Rayleigh阻尼法确定阻尼矩阵，即假定阻尼矩阵为刚度矩阵和质量矩阵的线性组合KMC(1-27)特别要注意的是，对于某些无阻尼问题，有限元方程的形式仅为FuMKu(1-28)1.3求解算法在建立了动力分析的有限元方程后，通过给定的边界条件就可以求解出体系的动力反应。在对有限元方程进行求解时，常使用直接积分方法进行计算。直接积分方法是指在积分运动方程之前不进行方程形式的变换，而直接进行数值积分。这类方法的特点是求解域内仅在离散时间点(而非任意时刻)能够满足运动方程的要求。而在离散时间点间假定了节点加速度、节点速度和节点位移的函数形式，从而建立起在相邻离散时间点时各节点反应值间的相互关系，由此每一时刻的反应值便可逐步计算出来。直接积分方法可分为两类，分别是隐式方法和显式方法。1有限元方程的隐式解法隐式方法的研究历史较长，且大部分为无条件稳定，例如，线性加速度法、常平均加速度法、Newmark法、Wilson-法等。隐式方法需要对藕联的线性方程组进行直接求解，通过求得的反应值建立下一迭代步的线性方程组(以下一时刻反应值作为未知数)，进而求解并得到该方程组中未知的反应值，以此类推。隐式方法的计算精度较高，算法稳定性好，但对于自由度数目较大或非线性程度较高的问题，隐式方法所需的计算量很大，计算成本较高。(1)线性加速度法线性加速度法假设在相邻的离散时间点P和P+1间，节点m的加速度分量miu以线性规律变化，如下所示。tttttuuuuPmiPmiPmimi),()(1(1-37)由式(1-37)积分并整理可得节点m在P+1时刻的速度和位移的表达式tutuuduuuPmiPmiPmitttmiPmiPmi112121)((1-38)12216131)(PmiPmiPmiPmimitttPmiPmiutututuduuu(1-39)将式(1-38)和(1-39)代入P+1时刻的平衡方程，即可整理得到P+1时刻的节点加速度分量1Pmiu的表达式lijmlijmlijmPliPmjPmjPmilijmPmjPmjPmjPmjlijmPmjmtcktfututucutututuku636)366()366(21221..（1-40）通过式(1-40)即可通过P时刻的结果求得P+1时刻的节点加速度，将求得的加速度再分别代入式(1-38)和(1-39)，又可得到P+1时刻的节点速度分量和节点位移分量1njPu、1njPu通过依次迭代直至计算结束，便可实现有限元方程的求解。线性加速度方法是一种具有二阶计算精度的条件稳定的自起步逐步积分格式。(2)常平均加速度法常平均加速度法假设在相邻的离散时间点p和p+1间，节点m的加速度分量miu是一个常量，如式((2-41)所示。tttuuuPmiPmimi,2)(..1..(1-41)由式(2-41)积分并整理可得节点}r2在P+1时刻的速度和位移的表达式tutuuduuuPmiPmiPmitttmiPmiPmi1....1.2121)((1-42)212..14141)(tututuuduuuPmiPmiPmiPmitttmiPmiPmi(1-43)将式(1-42)和(1-43)代入P+1时刻的平衡方程，即可得到P+1时刻的节点加速度分量1mjPu的表认式lijmlijmlijmPliPmjPmilijmPmjPmjPmjlijmPmjmtcktfutucututuku424)2()44(2121..(1-44)通过式(1-44)即可通过P时刻的结果求得P+1时刻的节点加速度，将求得的加速度再分别代入式(1-42)和式(1-43)，又可得到P+1时刻的节点速度分量和节点位移分量1Pmju、1Pmju。依次迭代直至计算结束，便可实现有限元方程的求解。常平均加速度方法为条件稳定的自起步的逐步积分格式，计算精度为二阶。(3)Newmark法Newmark法是1959年Newmark提出的一种著名的数值求解方法。假定相邻离散时间点上节点m的速度和位移满足如下的差分公式tuuuuPmiPmiPmiPmi])1[(1.....1(1-45)21.....1])5.0[(tuuutuuPmiPmiPmiPmiPmi(1-46)其中，参数、用于控制积分的格式精度及稳定性。将式(1-45)和(1-46)代入P+1时刻的平衡方程，可得P+1时刻节点m的加速度满足如下差分公式lijmlijmlijmPliPmjPmjPmjlijmPmjPmjlijmPmjkttcmfututukutucu212..1..))5.0(())1(((1-47)通过式(1-47)即可通过P时刻的结果求得P+1时刻的节点加速度，将求得的加速度再分别代入式(1-45)和式(1-46)，又可得到P+1时刻的节点速度分量和节点位移分量1mjPu、1mjPu。依次迭代直至计算结束，由此便建立了一个求解动力学有限元方程的Newmark逐步积分格式群:当5.0时，Newmark法的计算精度为二阶，若不满足此条件则为一阶;当、“满足条件5.025.025.0）（时，Newmark法作为一种自起步的逐步积分格式，是无条件稳定的;当5.0，=1/6时，Newmark法在形式上将退化为线性加速度法;当5.0，=0.25时，Newmark法则将退化为常平均加速度法。2、有限元方程的显式解法显式解法的计算量和存储量较隐式方法要小得多，较隐式方法更适合于自由度数目较大或非线性程度较高的问题的求解。采用显式算法进行计算时，无需求解藕联的线性方程组，只需以一定格式进行递推迭代即可，因此大大简化了求解过程，提高了求解计算的速率。在对时间维进行离散时，使用较多的是中心差分格式，并由此衍生出一系列的显式求解算法。(1)中心差分法将速度和加速度用中心差分格式进行离散，有tuuuPPP211.(1-48)211..2tuuuuPPPP(1-49)取集中化的质量矩阵，同时将各节点的速度向量和加速度向量分别用式(1-48)和(1-49)进行离散并代入P+1时刻的平衡方程，可得中心差分法的公式:pnjjpnjjpmipmimpmipnjjpmimukuctuumtfuctumtmin1min121min1221)2(1211(1-50)当阻尼阵为非对角的非零矩阵时，式(1-50)仍是隐式的，用该式进行计算无法体现中心差分法的优越性。当阻尼矩阵为对角阵时，式(1-50)变为显式积分公式，具有二阶计算精度。)2(22211min221pmipmimpmimpnjjpmipmiuumutcuktftu(1-51)当阻尼矩阵为零阵时，式(1-50)转变为无阻尼问题的显式求解公式，具有二阶计算精度。1min121212pmipmipnjjmpmimpmiuuukmtfmtu(1-52)由式(1-51)和式(1-52)可知，该方法在进行求解时只需将前面的两个迭代步中计算得到的结果代入当前步中即可直接推演得到当前迭代步的位移值。这表明，相对于隐式积分方法来说显式方法的计算效率要高很多，这也是显式方法最大的优势所在。目前，中心差分法格式己被广泛的应用于工程计算领域。例如，著名的显式有限元分析软件LS-DYNA就采用了这种显式求解格式。(2)中心差分一单边差分法用中心差分法离散节点加速度(如式(1-49))，用单边差分法离散节点速度，如式(1-53)所示。)(11.PPPuutu(1-53)取集中化的质量矩阵，同时将各节点的加速度向量和速度