中心差分法编程C语言,中心差分法的基本理论与程序设计资料.docx

中心差分法的基本理论与程序设计程序设计的目的与意义该程序通过用C语言(部分C++语言)编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。程序功能及特点该程序采用C语言(部分C++语言)实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。计算简便且在算法稳定的条件下，精度较高。中心差分法的基本理论在动力学问题中，系统的有限元求解方程(运动方程)如下所示：

式中， 是系统结点位移向量，和分别是系统的结点加速度向量和结点速度向量，和分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成???与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。通常的直接积分是基于两个概念，一是将在求解域内的任何时刻都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔的离散的时间点满足运动方程；二是在一定数目的区域内，假设位移、速度、加速度的函数形式。 中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。在中心差分法中，加速度和速度可以用位移表示，即：

时间的位移解答，可由时间的运动方程应得到满足，即由下式：

而得到。为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：

若已经求得和，则从上式可以进一步解出。所以上式是求解各个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。需要指出的是，此算法有一个起步问题。因为当时，为了计算，除了知道初始条件已知的，还需要知道，所以必须用一专门的起步方法。根据以上加速度和速度的表达式可知：

其中和可以从给定的初始条件中得到，而则可以利用时的运动方程得到，即：

中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利用它求解具体问题时，时间步长必须小于由该问题求解方程性质所决定的某个临界值，否则算法将是不稳定的。中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示：初始计算形成刚度矩阵、质量矩阵和阻尼矩阵；给定，和；选择时间步长，，并计算积分常数，，， ；计算；形成有效质量矩阵；三角分解：。对于每一时间步长()计算时间的有效载荷

求解时间的位移

如果需要，计算时间的加速度和速度

程序设计程序流程

图1 程序流程图各子程序主要功能为：ArrayLU：LU三角分解；Inverse：求矩阵的转置矩阵；ArrayMVector：矩阵和向量的乘法；LUSolve：求解方程。输入数据及变量说明输入数据该程序的原始输入数据应包括三个部分：刚度矩阵，质量矩阵和阻尼矩阵；初始条件：时间时刻的位移，速度，加速度；确定时间步长，其中为了保证该算法的稳定性，需要满足。变量说明该程序的各个变量含义如下：num，timeStep，dtnum——矩阵维度；timeStep——时间步数；dt——时间步长；M，C，K，X，V，A，P，MM，PT，c0，c1，c2，c3M——质量矩阵；C——阻尼矩阵；K——刚度矩阵；X——位移矩阵；V——速度矩阵；A——加速度矩阵；P——载荷向量；MM——有效质量矩阵；PT——时间时刻的有效载荷；c0,c1,c2,c3——积分常数；算例问题描述应用本程序计算一个三自由度系统，它的运动方程是：

初始条件：当时，，。已知此系统的固有频率为：，，。相应的振动周期为：，，。当时，利用公式，可以计算得到；时间步长分别取和进行计算。理论计算中心差分法(理论解)时间步长当时，其积分常数为： 则起步条件为

有效质量矩阵为：

对于每一时间步长，需先计算有效载荷：

再从下列方程计算时间的位移

由上式得到的每一时间步长的位移结果如表1所示：表1 理论解时间23456789100.000.000.000.030.130.360.791.462.37