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

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

中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言(部分C语言)编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。2程序功能及特点该程序采用C语言(部分C语言)实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。计算简便且在算法稳定的条件下，精度较高。3中心差分法的基本理论在动力学问题中，系统的有限元求解方程(运动方程)如下所示MATCTKATQ式中，是系统结点位移向量，和分别是系统的结点加速度AT向量和结点速度向量，和分别是系统的质量矩阵、阻尼矩阵、刚,T度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。通常的直接积分是基于两个概念，一是将在求解域内的任何时刻都应满足运动方程的要求，代之仅在一定条件下近似地0TTT满足运动方程，例如可以仅在相隔的离散的时间点满足运动方程；二是在一T定数目的区域内，假设位移、速度、加速度的函数形式。TAA中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。在中心差分法中，加速度和速度可以用位移表示，即21TTTAATT时间的位移解答，可由时间的运动方程应得到满足，即由下式TTATTTTMCAKQ而得到。为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式22211TTTTAAMCATT若已经求得和，则从上式可以进一步解出。所以上式是求解各TTT个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。需要指出的是，此算法有一个起步问题。因为当时，为了计算，除了知道初始0TTA条件已知的，还需要知道，所以必须用一专门的起步方法。根据以上加0ATA速度和速度的表达式可知200TTA其中和可以从给定的初始条件中得到，而则可以利用时的运动方程0A00T得到，即0MCKQ1000AMQCAK中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利用它求解具体问题时，时间步长必须小于由该问题求解方程性质所决定的T某个临界值，否则算法将是不稳定的。中心差分法比较适用于由冲击、爆CRT炸类型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。4中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示1初始计算1形成刚度矩阵、质量矩阵和阻尼矩阵；KMC2给定，和；0A03选择时间步长，，并计算积分常数，T2NCRTT021CT，，；12CT203214计算；TATCA5形成有效质量矩阵；01ˆMC6三角分解。ˆLU2对于每一时间步长()0,2TT1计算时间的有效载荷T201ˆTTTTQKCMACCA2求解时间的位移TˆTTLUQ3如果需要，计算时间的加速度和速度T02TTTTACA1TTT5程序设计51程序流程给定初始时刻的位移，速度和及速度选择时间步长，计算积分常数计算起步条件LU分解MLUARRAYLU子程序计算时间T的有效载荷(ARRAYMVECTOR子程序)求解时间TΔT的位移LUATΔTQLUSOLVE子程序输出时间TΔT的位移计算时间T的速度和加速度启动输入原始数据(M,C,K)输出结果结束图1程序流程图各子程序主要功能为ARRAYLULU三角分解；ALUINVERSE求矩阵的转置矩阵；ARRAYMVECTOR矩阵和向量的乘法；LUSOLVE求解方程。XP52输入数据及变量说明521输入数据该程序的原始输入数据应包括三个部分1刚度矩阵，质量矩阵和阻尼矩阵；KMC2初始条件时间时刻的位移，速度，加速度；0T0A00A3确定时间步长，其中为了保证该算法的稳定性，需要满足。2NCRTT522变量说明该程序的各个变量含义如下1NUM，TIMESTEP，DTNUM矩阵维度；TIMESTEP时间步数；DT时间步长；2M，C，K，X，V，A，P，MM，PT，C0，C1，C2，C3M质量矩阵；C阻尼矩阵；K刚度矩阵；X位移矩阵；V速度矩阵；A加速度矩阵；P载荷向量；MM有效质量矩阵；PT时间时刻的有效载荷；TC0,C1,C2,C3积分常数；6算例61问题描述应用本程序计算一个三自由度系统，它的运动方程是10210346AA初始条件当时，，。0T00已知此系统的固有频率为，，。相应的振动周期为1323，，。1089T24368T当时，利用公式，可以计算得到T1000AMQCAK；时间步长分别取和进行计算。06TA36TT35184TT62理论计算621中心差分法(理论解)1时间步长3106TT当时，其积分常数为0217589CT1372CT201578C321658CE则起步条件为03006359635TATCA有效质量矩阵为ˆM01007590ˆ7593182759MCC对于每一时间步长，需先计算有效载荷201ˆ01387590452638759TTTTTTQKCMACCAA再从下列方程计算时间的位移TTA7590ˆ2759TQ由上式得到的每一时间步长的位移结果如表1所示表1理论解306TT时间T2T345T6T7T8T9T10T1A00000000000301303607914623734220000030190581262243434695846773A0401482974525826717227517858452时间步长35184TT按照相同的步骤，所得结果如下209871TA507164TA7831026TA再计算下去，位移将继续增大，这是不稳定的典型表现。其原因是在条件稳定的中心差分方法中采用了远大于的时间步长，所以CRTT35184TT不可能得到有意义的结果。622振型叠加法(精确解)采用振型叠加法对上述问题进行计算，可以得到该问题的精确解。首先应求解的广义特征值问题为22101043按照一般的线性代