中心差分法的基本理论与程序设计

中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。2程序功能及特点该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。计算简便且在算法稳定的条件下，精度较高。3中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：MatCatKatQt式中，at是系统结点位移向量，at和at分别是系统的结点加速度向量和结点速度向量，,,MCK和Qt分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。通常的直接积分是基于两个概念，一是将在求解域0tT内的任何时刻t都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔t的离散的时间点满足运动方程；二是在一定数目的t区域内，假设位移a、速度a、加速度a的函数形式。中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。在中心差分法中，加速度和速度可以用位移表示，即：212tttttaaaat12ttttaaat时间tt的位移解答tta，可由时间t的运动方程应得到满足，即由下式：ttttMaCaKaQ而得到。为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：2221121122ttttttMCaQKMaMCattttt若已经求得ta和tta，则从上式可以进一步解出tta。所以上式是求解各个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。需要指出的是，此算法有一个起步问题。因为当0t时，为了计算ta，除了知道初始条件已知的0a，还需要知道ta，所以必须用一专门的起步方法。根据以上加速度和速度的表达式可知：20002ttaataa其中0a和0a可以从给定的初始条件中得到，而0a则可以利用0t时的运动方程0000MaCaKaQ得到，即：10000aMQCaKa中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利用它求解具体问题时，时间步长t必须小于由该问题求解方程性质所决定的某个临界值crt，否则算法将是不稳定的。中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。4中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示：1.初始计算(1)形成刚度矩阵K、质量矩阵M和阻尼矩阵C；(2)给定0a，0a和0a；(3)选择时间步长t，2ncrnTtt，并计算积分常数021ct，112ct，202cc，321cc；(4)计算0030taataca；(5)形成有效质量矩阵01ˆMcMcC；(6)三角分解ˆM：ˆMLU。2.对于每一时间步长（0,,2ttt）(1)计算时间t的有效载荷201ˆtttttQQKcMacMcCa(2)求解时间tt的位移ˆtttLUaQ(3)如果需要，计算时间t的加速度和速度02ttttttacaaa1tttttacaa5程序设计5.1程序流程给定初始时刻的位移，速度和及速度选择时间步长，计算积分常数计算起步条件LU分解M=LU(ArrayLU子程序)计算时间t的有效载荷（ArrayMVector子程序）求解时间t+Δt的位移LUa(t+Δt)=Q(LUSolve子程序)输出时间t+Δt的位移计算时间t的速度和加速度启动输入原始数据（M,C,K）输出结果结束图1程序流程图各子程序主要功能为：ArrayLU：LU三角分解ALUA；Inverse：求矩阵的转置矩阵；ArrayMVector：矩阵和向量的乘法；LUSolve：求解方程LUXP。5.2输入数据及变量说明5.2.1输入数据该程序的原始输入数据应包括三个部分：(1)刚度矩阵K，质量矩阵M和阻尼矩阵C；(2)初始条件：时间0t时刻的位移0a，速度0a，加速度0a；(3)确定时间步长t，其中为了保证该算法的稳定性，需要满足2ncrnTtt。5.2.2变量说明该程序的各个变量含义如下：(1)num，timeStep，dtnum——矩阵维度；timeStep——时间步数；dt——时间步长；(2)M，C，K，X，V，A，P，MM，PT，c0，c1，c2，c3M——质量矩阵；C——阻尼矩阵；K——刚度矩阵；X——位移矩阵；V——速度矩阵；A——加速度矩阵；P——载荷向量；MM——有效质量矩阵；PT——时间t时刻的有效载荷；c0,c1,c2,c3——积分常数；6算例6.1问题描述应用本程序计算一个三自由度系统，它的运动方程是：100210003014200010226aa初始条件：当0t时，00a，00a。已知此系统的固有频率为：113，22，33。相应的振动周期为：110.89T，24.444T，33.628T。当0t时，利用公式10000aMQCaKa，可以计算得到0006Ta；时间步长分别取3100.363tT和3518.14tT进行计算。6.2理论计算6.2.1中心差分法（理论解）(1)时间步长3100.363tT当3100.363tT时，其积分常数为：0217.589ct111.3772ct20215.178cc3216.5882cce则起步条件为0030000000.36300.0659000060.3953taataca有效质量矩阵ˆM为：011000007.5900ˆ7.590301.38000022.770001000007.59McMcC对于每一时间步长，需先计算有效载荷：201ˆ013.18107.59000141.532022.77060213.18007.59ttttttttQQKcMacMcCaaa再从下列方程计算tt时间的位移tta7.5900ˆ022.770007.59ttaQ由上式得到的每一时间步长的位移结果如表1所示：表13100.363tT理论解时间t2t3t4t5t6t7t8t9t10t1a0.000.000.000.030.130.360.791.462.373.422a0.000.030.190.581.262.243.434.695.846.773a0.401.482.974.525.826.717.227.517.858.45(2)时间步长3518.14tT按照相同的步骤，所得结果如下：2009.8710ta52502.07106.4610ta78387.13102.36105.6610ta再计算下去，位移将继续增大，这是不稳定的典型表现。其原因是在条件稳定的中心差分方法中采用了远大于crtT的时间步长3518.14tT，所以不可能得到有意义的结果。6.2.2振型叠加法（精确解）采用振型叠加法对上述问题进行计算，可以得到该问题的精确解。首先应求解的广义特征值问题为：2210100142030022001按照一般的线性代数方法可以得到上式的解答为：211315123T2221201T2331112T利用上式，可以将原问题转换为以123,,为基向量的3个互不耦合的运动方程，即：1119103xtxt22265xtxt33332xtxt原系统的初始条件是00a，00a，经转换后为：00itx00itx1,2,3i利用无阻尼情形的杜哈美积分公式，可以得到上述方程的解答为：1101cos1327xtt231cos25xtt311cos32xtt最后利用振型叠加得到系统的位移为：101cos1327121353011cos2521211cos32tattt根据上式计算得到的每一时间步长的位移值是系统响应的精确解，如表2和表3所示。(1)3100.363tT的精确位移值。表23100.363tT精确解时间t2t3t4t5t6t7t8t9t10t1a0.000.000.010.040.140.370.791.452.323.342a0.000.040.210.591.252.213.384.635.806.763a0.391.452.924.485.816.747.297.607.928.46(2)3518.14tT的精确位移值。表33518.14tT精确解时间t2t3t4t5t6t7t8t9t10t1a3.893.46-1.192.251.83-2.401.782.25-1.233.402a6.756.760.006.736.770.006.726.790.006.703a8.178.410.608.979.241.209.199.050.618.366.3程序计算(1)时间步长3100.363tT当3100.363tT时，计算得到的起步条件为000.395307Tta该程序的计算过程如图2所示，具体输出结果如表4所示。表43100.363tT数值解时间t2t3t4t5t6t7t8t9t10t1a0.000.000.010.030.130.360.791.462.373.422a0.000.040.190.581.262.243.434.695.846.773a0.401.482.974.525.826.717.227.517.858.45图23100.363tT时的程序计算过程及计算结果(2)时间步长3518.14tT当3518.14tT时，计算得到的起步条件为00987.179Tta该程序的计算过程如图3所示，具体输出结果如表5所示。表53518.14tT数值解时间t2t3t4t5t1a0.000.0077.1310111.2410141.59102a0.0052.171082.3610112