《中心差分法的基本理论与程序设计(共26页).docx》由会员分享，可在线阅读，更多相关《中心差分法的基本理论与程序设计(共26页).docx（26页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上中心差分法的基本理论与程序设计1 程序设计的目的与意义该程序通过用C语言（部分C+语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。2 程序功能及特点该程序采用C语言（部分C+语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。计算简便且在算法稳定的条件下，精度较高。3 中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：式中， 是系统结点位移向量，和分别是
2、系统的结点加速度向量和结点速度向量，和分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。通常的直接积分是基于两个概念，一是将在求解域内的任何时刻都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔的离散的时间点满足运动方程；二是在一定数目的区
3、域内，假设位移、速度、加速度的函数形式。 中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。在中心差分法中，加速度和速度可以用位移表示，即：时间的位移解答，可由时间的运动方程应得到满足，即由下式：而得到。为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：若已经求得和，则从上式可以进一步解出。所以上式是求解各个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。需要指
4、出的是，此算法有一个起步问题。因为当时，为了计算，除了知道初始条件已知的，还需要知道，所以必须用一专门的起步方法。根据以上加速度和速度的表达式可知：其中和可以从给定的初始条件中得到，而则可以利用时的运动方程得到，即：中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利用它求解具体问题时，时间步长必须小于由该问题求解方程性质所决定的某个临界值，否则算法将是不稳定的。中心差分法比较适用于由冲击、爆炸类型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。4 中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示：1. 初始计算(1) 形成刚度矩阵、质量矩阵
5、和阻尼矩阵；(2) 给定，和；(3) 选择时间步长，并计算积分常数， ；(4) 计算；(5) 形成有效质量矩阵；(6) 三角分解：。2. 对于每一时间步长（）(1) 计算时间的有效载荷(2) 求解时间的位移(3) 如果需要，计算时间的加速度和速度5 程序设计5.1 程序流程图1 程序流程图各子程序主要功能为：ArrayLU：LU三角分解；Inverse：求矩阵的转置矩阵；ArrayMVector：矩阵和向量的乘法；LUSolve：求解方程。5.2 输入数据及变量说明5.2.1 输入数据该程序的原始输入数据应包括三个部分：(1) 刚度矩阵，质量矩阵和阻尼矩阵；(2) 初始条件：时间时刻的位移，速
6、度，加速度；(3) 确定时间步长，其中为了保证该算法的稳定性，需要满足。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时间时刻的有效载荷；c0,c1,c2,c3积分常数；6 算例6.1 问题描述应用本程序计算一个三自由度系统，它的运动方程是：初始条件：当时，。已知此系统的固有频率为：，。相应的振动周期为：，。当时，利用公式，可
7、以计算得到；时间步长分别取和进行计算。6.2 理论计算6.2.1 中心差分法（理论解）(1) 时间步长当时，其积分常数为： 则起步条件为有效质量矩阵为：对于每一时间步长，需先计算有效载荷：再从下列方程计算时间的位移由上式得到的每一时间步长的位移结果如表1所示：表1 理论解时间23456789100.000.000.000.030.130.360.791.462.373.420.000.030.190.581.262.243.434.695.846.770.401.482.974.525.826.717.227.517.858.45(2) 时间步长按照相同的步骤，所得结果如下： 再计算下去，位移
8、将继续增大，这是不稳定的典型表现。其原因是在条件稳定的中心差分方法中采用了远大于的时间步长，所以不可能得到有意义的结果。6.2.2 振型叠加法（精确解）采用振型叠加法对上述问题进行计算，可以得到该问题的精确解。首先应求解的广义特征值问题为：按照一般的线性代数方法可以得到上式的解答为： 利用上式，可以将原问题转换为以为基向量的3个互不耦合的运动方程，即：原系统的初始条件是，经转换后为： 利用无阻尼情形的杜哈美积分公式，可以得到上述方程的解答为：最后利用振型叠加得到系统的位移为：根据上式计算得到的每一时间步长的位移值是系统响应的精确解，如表2 和表3所示。(1) 的精确位移值。表2 精确解时间23
9、456789100.000.000.010.040.140.370.791.452.323.340.000.040.210.591.252.213.384.635.806.760.391.452.924.485.816.747.297.607.928.46(2) 的精确位移值。表3 精确解时间23456789103.893.46-1.192.251.83-2.401.782.25-1.233.406.756.760.006.736.770.006.726.790.006.708.178.410.608.979.241.209.199.050.618.366.3 程序计算(1) 时间步长当时，计
10、算得到的起步条件为该程序的计算过程如图2所示，具体输出结果如表4所示。表4 数值解时间23456789100.000.000.010.030.130.360.791.462.373.420.000.040.190.581.262.243.434.695.846.770.401.482.974.525.826.717.227.517.858.45图2 时的程序计算过程及计算结果(2) 时间步长当时，计算得到的起步条件为该程序的计算过程如图3所示，具体输出结果如表5所示。表5 数值解时间23450.000.000.00图3 时的程序计算过程及计算结果6.4 结果讨论(1) 时间步长对于的情况，三者
11、的比较如表6和图4所示。表6 理论解，精确解和数值解比较结果理论解精确解数值解t(s)a1(m)a2(m)a3(m)a1(m)a2(m)a3(m)a1(m)a2(m)a3(m)0.726000.4000.39000.41.08900.031.4800.041.4500.041.481.45200.192.970.010.212.920.010.192.971.8150.030.584.520.040.594.480.030.584.522.1780.131.265.820.141.255.810.131.265.822.5410.362.246.710.372.216.740.362.246.
12、712.9040.793.437.220.793.387.290.793.437.223.2671.464.697.511.454.637.61.464.697.513.632.375.847.852.325.87.922.375.847.853.9933.426.778.453.346.768.463.426.778.45(a) 理论解，精确解和数值解比较结果()(b) 理论解，精确解和数值解比较结果()(c)理论解，精确解和数值解比较结果()图4 理论解，精确解和数值解比较结果由上述结果可知，当时，此时较小（），不论用理论计算还是程序计算，都与精确解吻合的非常好。说明了该程序用于计算实际问
13、题的准确性和有效性，能够满足精度要求。(2) 时间步长对于的情况，如果采用中心差分法，无论是理论计算还是程序计算，所得到的位移都会无限增大，说明此时的算法是不稳定的，与精确解并没有可比性。这也可以证明中心差分法是条件稳定算法，在利用其求解具体问题时，时间步长必须小于由该问题求解方程性质所决定的某个临界值，否则算法将是不稳定的。为了保证解的稳定性，中心差分法的步长必须满足下列条件：其中，是系统的最高阶固有振动频率，是系统的最小固有振动周期。原则上说，可以利用一般矩阵特征值问题的求解方法得到。实际上只需要求解系统中最小尺寸单元的最小固有振动周期，该数值可以通过一定的方法估计得到。7 总结本文中的程
14、序采用了C语言（部分C+语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。计算简便且在算法稳定的条件下，精度较高。中心差分法是对运动方程不进行方程形式的变换而直接进行逐步数值积分，属于直接积分法。中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。中心差分法的基本算法步骤为：输入刚度矩阵，质量矩阵和阻尼矩阵，并且给定初始时刻的位移，速度，加速度（可
15、通过计算得到）；选择时间步长，计算积分常数；计算起步条件；形成有效质量矩阵，并对其进行三角分解；对于每一个时间步长计算时间时刻的有效载荷并求解时间的位移，然后迭代求解。本文通过求解一个三自由度系统的运动方程验证了该程序的正确性。首先通过中心差分法，振型叠加法和该程序求解得到了问题的理论解，精确解和数值解。然后对这些结果进行了比较，结果表明，当时，理论解，精确解和数值解吻合的非常好，证明该算法的合理性与准确性；当时，理论解和数值解得到的位移结果都会无限增大，说明此时的算法是不稳定的，与精确解并没有可比性，这也间接证明中心差分法是条件稳定算法。8 心得体会我认为在编写程序的过程中，设计算法是核心部
16、分，需要对基础理论有很深刻的理解，清楚求解方法的每一个步骤，同时还要考虑程序语言的实现。在得到算法步骤以后，需要将其翻译成计算机程序设计语言，然后进行编辑，调试和运行，得到运行结果。得到运行结果并不意味着程序正确，还需要算例对其进行验证分析。附录/ 中心差分法程序#include #include #include #include using namespace std;/调用子程序void ArrayLU(double*, double*, double*, int); void Inverse(double*, double*, int); void ArrayMVector(doubl
17、e*, double*, double*, int); void LUSolve(double*, double*, double*, double*,int); void Centre(double*, double*, double*, double*, double*, double*, double*, int, int, double);int main()/设置矩阵维度，时间步，时间步长 int num=3; int timeStep=11; double dt=0.363; / double dt=18.14;/开辟矩阵空间 double *M = new double*num;
18、 for(int i=0; inum; i+) Mi = new doublenum; double *C = new double*num; for(int i=0; inum; i+) Ci = new doublenum; double *K = new double*num; for(int i=0; inum; i+) Ki = new doublenum; double *X = new double*timeStep+1; for(int i=0; itimeStep+1; i+) Xi = new doublenum; double *V = new double*timeSt
19、ep+1; for(int i=0; itimeStep+1; i+) Vi = new doublenum; double *A = new double*timeStep+1; for(int i=0; itimeStep+1; i+) Ai = new doublenum; double *P = new doublenum; /设置计算参数与初始条件 X10 = 0; X11 = 0; X12 = 0; V10 = 0; V11 = 0; V12 = 0; A10 = 0; A11 = 0; A12 = 6; M00 = 1; M01 = 0; M02 = 0;M10 = 0; M11
20、 = 3; M12 = 0; M20 = 0; M21 = 0; M22 = 1; C00 = 0; C01 = 0; C02 = 0; C10 = 0; C11 = 0; C12 = 0; C20 = 0; C21 = 0; C22 = 0; K00 = 2; K01 = -1; K02 = 0; K10 = -1; K11 = 4; K12 = -2; K20 = 0; K21 = -2; K22 = 2; P0=0; P1=0; P2=6;/开始计算，调用子程序Centre Centre(M,C,K,P,X,V,A,num,timeStep,dt); /输出计算结果 ofstream o
21、ut(out.txt); /输出结果文件out.txt if(!out.is_open() return; for(int i=0; i=timeStep; i+) couti*dt s:endl; /屏幕显示 out i*dt s:endl; /将结果写入结果文件 for(int j=0; jnum; j+) cout setw(20) Xij setw(20) Vij setw(20) Aij endl; outsetw(20) Xij setw(20) Vij setw(20) Aij endl; out.close(); return 0;/子程序Centrevoid Centre(d
22、ouble*M, double*C, double*K, double *P, double*X, double*V, double*A, int size, int timeStep, double dt) double c0,c1,c2,c3; double *PT=new doublesize; double *temp=new doublesize; double *tran=new doublesize; double *MM=new double*size; for(int i=0; isize; i+) MMi=new doublesize; double *L = new do
23、uble* size; for(int i=0; isize; i+) Li = new double size; double *U = new double* size; for(int i=0; isize; i+) Ui = new double size; /计算积分常数 c0=1/(dt*dt); c1=1/(2*dt); c2=2/(dt*dt); c3=(dt*dt)/2;/起步问题，计算 X00=X10-dt*V10+c3*A10; X01=X11-dt*V11+c3*A11; X02=X12-dt*V12+c3*A12; /形成有效质量矩阵for(int i=0; isiz
24、e; i+) for(int j=0; jsize; j+) MMij=c0*Mij+c1*Cij; /调用子程序ArrayLU对有效质量矩阵进行LU分解 ArrayLU(MM,L,U,size); /计算有效载荷for(int ti=1; ti=timeStep-1; ti+) for(int i=0; isize; i+) tempi = c2*Xtii-c0*Xti-1i; /调用ArrayMVector，得到矩阵和向量乘积 ArrayMVector(M,temp,tran,size); for(int i=0; isize; i+) PTi = Pi+trani; for(int i=
25、0; isize; i+) tempi = c1*Xti-1i; ArrayMVector(C,temp,tran,size); for(int i=0; isize; i+) PTi = PTi+trani; for(int i=0; isize; i+) tempi = Xtii; ArrayMVector(K,temp,tran,size); for(int i=0; isize; i+) PTi = PTi-trani;/计算时间的位移LUSolve(L,U,PT,temp,size); /输出结果：位移，速度，加速度for(int i=0; isize; i+) Xti+1i = t
26、empi; Atii = c0*(Xti-1i-2*Xtii+Xti+1i); Vtii = c1*(-Xti-1i+Xti+1i); /删除不用数组for(int i=0; isize; i+) delete MMi; delete MM; for(int i=0; isize; i+) delete Li; delete L; for(int i=0; isize; i+) delete Ui; delete U; delete PT; delete temp; delete tran; /子程序ArrayLU，对有效质量矩阵进行LU分解void ArrayLU(double*A, dou
27、ble*L, double*U, int size) double scale; double *B=new double*size; for(int i=0; isize; i+) Bi=new double size; for(int i=0; isize; i+) for(int j=0; jsize; j+) Bij = 0; Lij = 0; Uij = 0; Bii=1; for(int i=0; isize-1; i+) if(Aii=0) coutthe Array is singularendl; for(int j=i+1; jsize; j+) scale = Aji/A
28、ii; for(int k=i; ksize; k+) Ajk = Ajk-Aik*scale; for(int k=0; ksize; k+) Bjk = Bjk-Bik*scale; /调用子程序Inverse，对矩阵进行转置 Inverse(B, L, size); for(int i=0; isize; i+) for(int j=i; jsize; j+) Uij = Aij; for(int i=0; isize; i+) delete Bi; delete B; /子程序Inverse，对矩阵进行转置void Inverse(double*A, double*B, int siz
29、e) double scale; for(int i=0; isize; i+) Bii = 1; for(int i=0; isize-1; i+) if(Aii=0) cout the Array is singular endl; for(int j=i+1; jsize; j+) scale = Aji/Aii; for(int k=i; ksize; k+) Ajk = Ajk-Aik*scale; for(int k=0; k0; i-) for(int k=0; k=0; j-) scale = Aji/Aii; Aji = 0; for(int k=0; ksize; k+)
30、Bjk = Bjk-Bik*scale; for(int k=0; ksize; k+) B0k = B0k/A00; /子程序ArrayMVector，矩阵和向量的乘法void ArrayMVector(double *A, double *X, double *B, int size) for(int i=0; isize; i+) Bi = 0; for(int j=0; jsize; j+) Bi = Bi + Aij*Xj; /子程序LUSolve，求解方程void LUSolve(double*L, double*U, double*P, double*X, int size) double sum; double* Y = new doublesize; Y0=P0/L00; for(int i=1; isize; i+) sum = 0; for(int j=0; j=0; i-) sum=0; for(int j=i+1; jsize; j+) sum = sum + Uij*Xj; Xi=(Yi-sum)/Uii; delete Y; 专心-专注-专业