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

中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。

2程序功能及特点该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。

计算简便且在算法稳定的条件下，精度较高。

3中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：()()()()Ma t Ca t Ka t Q t++=式中，()a t分别是系统的结点加速度向a t是系统结点位移向量，()a t和()量和结点速度向量，,,M C K和()Q t分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。

与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。

常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。

中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。

通常的直接积分是基于两个概念，一是将在求解域0t T内的任何时刻t都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔t∆的离散的时间点满足运动方程；二是在一定数目的t∆区域内，假设位移a、速度a、加速度a的函数形式。

中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。

在中心差分法中，加速度和速度可以用位移表示，即：()212t t t t t a a a a t -∆+∆=-+∆ ()12t t t t a a a t-∆+∆=-+∆时间t t ∆+的位移解答t t a +∆，可由时间t 的运动方程应得到满足，即由下式：t t t t Ma Ca Ka Q ++=而得到。

为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：2221121122t t t tt t M C a Q K M a M C a t t t t t +∆-∆⎛⎫⎛⎫⎛⎫+=---- ⎪ ⎪ ⎪∆∆∆∆∆⎝⎭⎝⎭⎝⎭若已经求得t a 和t t a -∆，则从上式可以进一步解出t t a +∆。

所以上式是求解各个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。

需要指出的是，此算法有一个起步问题。

因为当0t =时，为了计算t a ∆，除了知道初始条件已知的0a ，还需要知道t a -∆，所以必须用一专门的起步方法。

根据以上加速度和速度的表达式可知：20002tt a a ta a -∆∆=-∆+其中0a 和0a 可以从给定的初始条件中得到，而0a 则可以利用0t =时的运动方程0000Ma Ca Ka Q ++=得到，即：()10000a M Q Ca Ka -=--中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利用它求解具体问题时，时间步长t ∆必须小于由该问题求解方程性质所决定的某个临界值cr t ∆，否则算法将是不稳定的。

中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。

4 中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示： 1. 初始计算(1) 形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C ； (2) 给定0a ，0a 和0a ； (3) 选择时间步长t ∆，2ncr nT tt ωπ∆∆==，并计算积分常数021c t =∆，112c t=∆，202c c =， 321c c =；(4) 计算0030t a a ta c a -∆=-∆+；(5) 形成有效质量矩阵01ˆM c M c C =+； (6) 三角分解ˆM：ˆM LU =。

2. 对于每一时间步长（0,,2t t t =∆∆⋅⋅⋅） (1) 计算时间t 的有效载荷()()201ˆt t t t tQ Q K c M a c M c C a -∆=---- (2) 求解时间t t +∆的位移ˆt t tLUa Q +∆= (3) 如果需要，计算时间t 的加速度和速度()02t t t t t t a c a a a -∆+∆=-+()1t t t t t a c a a -∆+∆=-+5程序设计5.1程序流程图1 程序流程图各子程序主要功能为：=；ArrayLU：LU三角分解A LUAInverse：求矩阵的转置矩阵；ArrayMVector：矩阵和向量的乘法；=。

LUSolve：求解方程LUX P5.2 输入数据及变量说明5.2.1 输入数据该程序的原始输入数据应包括三个部分： (1) 刚度矩阵K ，质量矩阵M 和阻尼矩阵C ；(2) 初始条件：时间0t =时刻的位移0a ，速度0a ，加速度0a ；(3) 确定时间步长t ∆，其中为了保证该算法的稳定性，需要满足2ncr nT tt ωπ∆∆==。

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 问题描述应用本程序计算一个三自由度系统，它的运动方程是：100210003014200010226a a -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥+--=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦初始条件：当0t =时，00a =，00a =。

已知此系统的固有频率为：1ω=2ω=3ω=110.89T =，2 4.444T =，3 3.628T =。

当0t =时，利用公式()10000a M Q Ca Ka -=--，可以计算得到[]0006Ta =；时间步长分别取3100.363t T ∆==和3518.14t T ∆==进行计算。

6.2 理论计算6.2.1 中心差分法（理论解）(1) 时间步长3100.363t T ∆==当3100.363t T ∆==时，其积分常数为：0217.589c t ==∆ 111.3772c t==∆ 20215.178c c == 321 6.5882c c e ==- 则起步条件为0030000000.36300.0659000060.3953ta a ta c a -∆⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-∆+=-+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦有效质量矩阵ˆM为：011000007.5900ˆ7.59030 1.38000022.770001000007.59M c M c C ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+=+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦对于每一时间步长，需先计算有效载荷：()()201ˆ013.18107.59000141.532022.77060213.18007.59t t t t tt t tQ Q K c M a c M c C a a a -∆-∆=----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦再从下列方程计算t t +∆时间的位移t t a +∆7.5900ˆ022.770007.59t t a Q +∆⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦由上式得到的每一时间步长的位移结果如表1所示：表1 3100.363t T ∆==理论解时间t ∆ 2t ∆ 3t ∆ 4t ∆ 5t ∆ 6t ∆ 7t ∆ 8t ∆ 9t ∆ 10t ∆ 1a 0.00 0.00 0.00 0.03 0.13 0.36 0.79 1.46 2.37 3.42 2a 0.00 0.03 0.19 0.58 1.26 2.24 3.43 4.69 5.84 6.77 3a0.401.482.974.525.826.717.227.517.858.45(2) 时间步长3518.14t T ∆==按照相同的步骤，所得结果如下：2009.8710t a ∆⎡⎤⎢⎥=⎢⎥⎢⎥⨯⎣⎦ 52502.07106.4610t a ∆⎡⎤⎢⎥=⨯⎢⎥⎢⎥⨯⎣⎦ 78387.13102.36105.6610t a ∆⎡⎤⨯⎢⎥=⨯⎢⎥⎢⎥⨯⎣⎦再计算下去，位移将继续增大，这是不稳定的典型表现。

其原因是在条件稳定的中心差分方法中采用了远大于()cr t T π∆=的时间步长()3518.14t T ∆==，所以不可能得到有意义的结果。

6.2.2 振型叠加法（精确解）采用振型叠加法对上述问题进行计算，可以得到该问题的精确解。

首先应求解的广义特征值问题为：2210100142030022001φωφ-⎡⎤⎡⎤⎢⎥⎢⎥--=⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦按照一般的线性代数方法可以得到上式的解答为：2113ω= 15123Tφ⎡⎤=⎢⎥⎣⎦222ω= []1201Tφ=-233ω= []1112Tφ=--利用上式，可以将原问题转换为以123,,φφφ为基向量的3个互不耦合的运动方程，即：()()1119103x t x t +=()()22265x t x t += ()()33332x t x t +=-原系统的初始条件是00a =，00a =，经转换后为：00i t x == 00i t x == ()1,2,3i =利用无阻尼情形的杜哈美积分公式，可以得到上述方程的解答为：())1101cos 27x t ⎡⎤=-⎣⎦ ())231cos 5x t ⎡⎤=-⎣⎦())311cos 2x t ⎡⎤=--⎣⎦最后利用振型叠加得到系统的位移为：())))101cos 27121353011cos 521211cos 2a t ⎡⎤⎡⎤-⎢⎥⎣⎦--⎡⎤⎢⎥⎢⎥⎢⎥⎡⎤=-⎢⎥⎣⎦⎢⎥⎢⎥-⎢⎥⎣⎦⎡⎤⎢⎥--⎣⎦⎢⎥⎣⎦根据上式计算得到的每一时间步长的位移值是系统响应的精确解，如表2 和表3所示。

(1) 3100.363t T ∆==的精确位移值。

表2 3100.363t T ∆==精确解时间t ∆ 2t ∆ 3t ∆ 4t ∆ 5t ∆ 6t ∆ 7t ∆ 8t ∆ 9t ∆ 10t ∆ 1a 0.00 0.00 0.01 0.04 0.14 0.37 0.79 1.45 2.32 3.34 2a 0.00 0.04 0.21 0.59 1.25 2.21 3.38 4.63 5.80 6.76 3a0.391.452.924.485.816.747.297.607.928.46(2) 3518.14t T ∆==的精确位移值。