中心差分法的基本理論與程序設(shè)計(jì)

1、中心差分法的基本理論與程序設(shè)計(jì)1 程序設(shè)計(jì)的目的與意義該程序通過用C語言（部分C+語言）編寫了有限元中用于求解動(dòng)力學(xué)問題的中心差分法，鞏固和掌握了中心差分法的基本概念，提高了實(shí)際動(dòng)手能力，并通過實(shí)際編程實(shí)現(xiàn)了中心差分法在求解某些動(dòng)力學(xué)問題中的運(yùn)用，加深了對(duì)該方法的理解和掌握。2 程序功能及特點(diǎn)該程序采用C語言（部分C+語言）實(shí)現(xiàn)了用于求解動(dòng)力學(xué)問題的中心差分法，可以求解得到運(yùn)動(dòng)方程的解答，包括位移，速度和加速度。計(jì)算簡(jiǎn)便且在算法穩(wěn)定的條件下，精度較高。3 中心差分法的基本理論在動(dòng)力學(xué)問題中，系統(tǒng)的有限元求解方程（運(yùn)動(dòng)方程）如下所示：式中， 是系統(tǒng)結(jié)點(diǎn)位移向量，和分別是系統(tǒng)的結(jié)點(diǎn)加速度向量和結(jié)點(diǎn)

2、速度向量，和分別是系統(tǒng)的質(zhì)量矩陣、阻尼矩陣、剛度矩陣和結(jié)點(diǎn)載荷向量，并分別由各自的單元矩陣和向量集成。與靜力學(xué)分析相比，在動(dòng)力分析中，由于慣性力和阻尼力出現(xiàn)在平衡方程中，因此引入了質(zhì)量矩陣和阻尼矩陣，最后得到的求解方程不是代數(shù)方程組，而是常微分方程組。常微分方程的求解方法可以分為兩類，即直接積分法和振型疊加法。中心差分法屬于直接積分法，其對(duì)運(yùn)動(dòng)方程不進(jìn)行方程形式的變換而直接進(jìn)行逐步數(shù)值積分。通常的直接積分是基于兩個(gè)概念，一是將在求解域內(nèi)的任何時(shí)刻都應(yīng)滿足運(yùn)動(dòng)方程的要求，代之僅在一定條件下近似地滿足運(yùn)動(dòng)方程，例如可以僅在相隔的離散的時(shí)間點(diǎn)滿足運(yùn)動(dòng)方程；二是在一定數(shù)目的區(qū)域內(nèi)，假設(shè)位移、速度、加速

3、度的函數(shù)形式。 中心差分法的基本思路是用有限差分代替位移對(duì)時(shí)間的求導(dǎo)，將運(yùn)動(dòng)方程中的速度和加速度用位移的某種組合表示，然后將常微分方程組的求解問題轉(zhuǎn)換為代數(shù)方程組的求解問題，并假設(shè)在每個(gè)小的時(shí)間間隔內(nèi)滿足運(yùn)動(dòng)方程，則可以求得每個(gè)時(shí)間間隔的遞推公式，進(jìn)而求得整個(gè)時(shí)程的反應(yīng)。在中心差分法中，加速度和速度可以用位移表示，即：時(shí)間的位移解答，可由時(shí)間的運(yùn)動(dòng)方程應(yīng)得到滿足，即由下式：而得到。為此將加速度和速度的表達(dá)式代入上式中，即可得到中心差分法的遞推公式：若已經(jīng)求得和，則從上式可以進(jìn)一步解出。所以上式是求解各個(gè)離散時(shí)間點(diǎn)的解的遞推公式，這種數(shù)值積分方法又稱為逐步積分法。需要指出的是，此算法有一個(gè)起步問

4、題。因?yàn)楫?dāng)時(shí)，為了計(jì)算，除了知道初始條件已知的，還需要知道，所以必須用一專門的起步方法。根據(jù)以上加速度和速度的表達(dá)式可知：其中和可以從給定的初始條件中得到，而則可以利用時(shí)的運(yùn)動(dòng)方程得到，即：中心差分法避免了矩陣求逆的運(yùn)算，是顯式算法，且其為條件穩(wěn)定算法，利用它求解具體問題時(shí)，時(shí)間步長(zhǎng)必須小于由該問題求解方程性質(zhì)所決定的某個(gè)臨界值，否則算法將是不穩(wěn)定的。中心差分法比較適用于由沖擊、爆炸類型載荷引起的波傳播問題的求解，而對(duì)于結(jié)構(gòu)動(dòng)力學(xué)問題則不太合適。4 中心差分法的有限元計(jì)算格式利用中心差分法逐步求解運(yùn)動(dòng)方程的算法步驟如下所示：1. 初始計(jì)算(1) 形成剛度矩陣、質(zhì)量矩陣和阻尼矩陣；(2) 給定，

5、和；(3) 選擇時(shí)間步長(zhǎng)，并計(jì)算積分常數(shù)， ；(4) 計(jì)算；(5) 形成有效質(zhì)量矩陣；(6) 三角分解：。2. 對(duì)于每一時(shí)間步長(zhǎng)（）(1) 計(jì)算時(shí)間的有效載荷(2) 求解時(shí)間的位移(3) 如果需要，計(jì)算時(shí)間的加速度和速度5 程序設(shè)計(jì)5.1 程序流程圖1 程序流程圖各子程序主要功能為：ArrayLU：LU三角分解；Inverse：求矩陣的轉(zhuǎn)置矩陣；ArrayMVector：矩陣和向量的乘法；LUSolve：求解方程。5.2 輸入數(shù)據(jù)及變量說明5.2.1 輸入數(shù)據(jù)該程序的原始輸入數(shù)據(jù)應(yīng)包括三個(gè)部分：(1) 剛度矩陣，質(zhì)量矩陣和阻尼矩陣；(2) 初始條件：時(shí)間時(shí)刻的位移，速度，加速度；(3) 確定時(shí)

6、間步長(zhǎng)，其中為了保證該算法的穩(wěn)定性，需要滿足。5.2.2 變量說明該程序的各個(gè)變量含義如下：(1) num，timeStep，dtnum矩陣維度；timeStep時(shí)間步數(shù)；dt時(shí)間步長(zhǎng)；(2) M，C，K，X，V，A，P，MM，PT，c0，c1，c2，c3M質(zhì)量矩陣；C阻尼矩陣；K剛度矩陣；X位移矩陣；V速度矩陣；A加速度矩陣；P載荷向量；MM有效質(zhì)量矩陣；PT時(shí)間時(shí)刻的有效載荷；c0,c1,c2,c3積分常數(shù)；6 算例6.1 問題描述應(yīng)用本程序計(jì)算一個(gè)三自由度系統(tǒng)，它的運(yùn)動(dòng)方程是：初始條件：當(dāng)時(shí)，。已知此系統(tǒng)的固有頻率為：，。相應(yīng)的振動(dòng)周期為：，。當(dāng)時(shí)，利用公式，可以計(jì)算得到；時(shí)間步長(zhǎng)分別取

7、和進(jìn)行計(jì)算。6.2 理論計(jì)算6.2.1 中心差分法（理論解）(1) 時(shí)間步長(zhǎng)當(dāng)時(shí)，其積分常數(shù)為： 則起步條件為有效質(zhì)量矩陣為：對(duì)于每一時(shí)間步長(zhǎng)，需先計(jì)算有效載荷：再?gòu)南铝蟹匠逃?jì)算時(shí)間的位移由上式得到的每一時(shí)間步長(zhǎng)的位移結(jié)果如表1所示：表1 理論解時(shí)間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) 時(shí)間步長(zhǎng)按照相同的步驟，所得結(jié)果如下： 再計(jì)算下去，位移將繼續(xù)增大，這是不穩(wěn)定的典

8、型表現(xiàn)。其原因是在條件穩(wěn)定的中心差分方法中采用了遠(yuǎn)大于的時(shí)間步長(zhǎng)，所以不可能得到有意義的結(jié)果。6.2.2 振型疊加法（精確解）采用振型疊加法對(duì)上述問題進(jìn)行計(jì)算，可以得到該問題的精確解。首先應(yīng)求解的廣義特征值問題為：按照一般的線性代數(shù)方法可以得到上式的解答為： 利用上式，可以將原問題轉(zhuǎn)換為以為基向量的3個(gè)互不耦合的運(yùn)動(dòng)方程，即：原系統(tǒng)的初始條件是，經(jīng)轉(zhuǎn)換后為： 利用無阻尼情形的杜哈美積分公式，可以得到上述方程的解答為：最后利用振型疊加得到系統(tǒng)的位移為：根據(jù)上式計(jì)算得到的每一時(shí)間步長(zhǎng)的位移值是系統(tǒng)響應(yīng)的精確解，如表2 和表3所示。(1) 的精確位移值。表2 精確解時(shí)間23456789100.000

9、.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 精確解時(shí)間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 程序計(jì)算(1) 時(shí)間步長(zhǎng)當(dāng)時(shí)，計(jì)算得到的起步條件為該程序的

10、計(jì)算過程如圖2所示，具體輸出結(jié)果如表4所示。表4 數(shù)值解時(shí)間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 時(shí)的程序計(jì)算過程及計(jì)算結(jié)果(2) 時(shí)間步長(zhǎng)當(dāng)時(shí)，計(jì)算得到的起步條件為該程序的計(jì)算過程如圖3所示，具體輸出結(jié)果如表5所示。表5 數(shù)值解時(shí)間23450.000.000.00圖3 時(shí)的程序計(jì)算過程及計(jì)算結(jié)果6.4 結(jié)果討論(1) 時(shí)間步長(zhǎng)對(duì)于的情況，三者的比較如表6和圖4所示。表

11、6 理論解，精確解和數(shù)值解比較結(jié)果理論解精確解數(shù)值解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.712.9040.793.

12、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) 理論解，精確解和數(shù)值解比較結(jié)果()(b) 理論解，精確解和數(shù)值解比較結(jié)果()(c)理論解，精確解和數(shù)值解比較結(jié)果()圖4 理論解，精確解和數(shù)值解比較結(jié)果由上述結(jié)果可知，當(dāng)時(shí)，此時(shí)較?。ǎ?，不論用理論計(jì)算還是程序計(jì)算，都與精確解吻合的非常好。說明了該程序用于計(jì)算實(shí)際問題的準(zhǔn)確性和有效性，能夠滿

13、足精度要求。(2) 時(shí)間步長(zhǎng)對(duì)于的情況，如果采用中心差分法，無論是理論計(jì)算還是程序計(jì)算，所得到的位移都會(huì)無限增大，說明此時(shí)的算法是不穩(wěn)定的，與精確解并沒有可比性。這也可以證明中心差分法是條件穩(wěn)定算法，在利用其求解具體問題時(shí)，時(shí)間步長(zhǎng)必須小于由該問題求解方程性質(zhì)所決定的某個(gè)臨界值，否則算法將是不穩(wěn)定的。為了保證解的穩(wěn)定性，中心差分法的步長(zhǎng)必須滿足下列條件：其中，是系統(tǒng)的最高階固有振動(dòng)頻率，是系統(tǒng)的最小固有振動(dòng)周期。原則上說，可以利用一般矩陣特征值問題的求解方法得到。實(shí)際上只需要求解系統(tǒng)中最小尺寸單元的最小固有振動(dòng)周期，該數(shù)值可以通過一定的方法估計(jì)得到。7 總結(jié)本文中的程序采用了C語言（部分C+語

14、言）實(shí)現(xiàn)了用于求解動(dòng)力學(xué)問題的中心差分法，可以求解得到運(yùn)動(dòng)方程的解答，包括位移，速度和加速度。計(jì)算簡(jiǎn)便且在算法穩(wěn)定的條件下，精度較高。中心差分法是對(duì)運(yùn)動(dòng)方程不進(jìn)行方程形式的變換而直接進(jìn)行逐步數(shù)值積分，屬于直接積分法。中心差分法的基本思路是用有限差分代替位移對(duì)時(shí)間的求導(dǎo)，將運(yùn)動(dòng)方程中的速度和加速度用位移的某種組合表示，然后將常微分方程組的求解問題轉(zhuǎn)換為代數(shù)方程組的求解問題，并假設(shè)在每個(gè)小的時(shí)間間隔內(nèi)滿足運(yùn)動(dòng)方程，則可以求得每個(gè)時(shí)間間隔的遞推公式，進(jìn)而求得整個(gè)時(shí)程的反應(yīng)。中心差分法的基本算法步驟為：輸入剛度矩陣，質(zhì)量矩陣和阻尼矩陣，并且給定初始時(shí)刻的位移，速度，加速度（可通過計(jì)算得到）；選擇時(shí)間步

15、長(zhǎng)，計(jì)算積分常數(shù)；計(jì)算起步條件；形成有效質(zhì)量矩陣，并對(duì)其進(jìn)行三角分解；對(duì)于每一個(gè)時(shí)間步長(zhǎng)計(jì)算時(shí)間時(shí)刻的有效載荷并求解時(shí)間的位移，然后迭代求解。本文通過求解一個(gè)三自由度系統(tǒng)的運(yùn)動(dòng)方程驗(yàn)證了該程序的正確性。首先通過中心差分法，振型疊加法和該程序求解得到了問題的理論解，精確解和數(shù)值解。然后對(duì)這些結(jié)果進(jìn)行了比較，結(jié)果表明，當(dāng)時(shí)，理論解，精確解和數(shù)值解吻合的非常好，證明該算法的合理性與準(zhǔn)確性；當(dāng)時(shí)，理論解和數(shù)值解得到的位移結(jié)果都會(huì)無限增大，說明此時(shí)的算法是不穩(wěn)定的，與精確解并沒有可比性，這也間接證明中心差分法是條件穩(wěn)定算法。8 心得體會(huì)我認(rèn)為在編寫程序的過程中，設(shè)計(jì)算法是核心部分，需要對(duì)基礎(chǔ)理論有很深刻

16、的理解，清楚求解方法的每一個(gè)步驟，同時(shí)還要考慮程序語言的實(shí)現(xiàn)。在得到算法步驟以后，需要將其翻譯成計(jì)算機(jī)程序設(shè)計(jì)語言，然后進(jìn)行編輯，調(diào)試和運(yùn)行，得到運(yùn)行結(jié)果。得到運(yùn)行結(jié)果并不意味著程序正確，還需要算例對(duì)其進(jìn)行驗(yàn)證分析。附錄/ 中心差分法程序#include <stdio.h>#include <iostream.h>#include <iomanip.h>#include <fstream.h> using namespace std;/調(diào)用子程序void ArrayLU(double*, double*, double*, int); void

17、Inverse(double*, double*, int); void ArrayMVector(double*, double*, double*, int); void LUSolve(double*, double*, double*, double*,int); void Centre(double*, double*, double*, double*, double*, double*, double*, int, int, double);int main()/設(shè)置矩陣維度，時(shí)間步，時(shí)間步長(zhǎng) int num=3; int timeStep=11; double dt=0.363

18、; / double dt=18.14;/開辟矩陣空間 double *M = new double*num; for(int i=0; i<num; i+) Mi = new doublenum; double *C = new double*num; for(int i=0; i<num; i+) Ci = new doublenum; double *K = new double*num; for(int i=0; i<num; i+) Ki = new doublenum; double *X = new double*timeStep+1; for(int i=0;

19、 i<timeStep+1; i+) Xi = new doublenum; double *V = new double*timeStep+1; for(int i=0; i<timeStep+1; i+) Vi = new doublenum; double *A = new double*timeStep+1; for(int i=0; i<timeStep+1; i+) Ai = new doublenum; double *P = new doublenum; /設(shè)置計(jì)算參數(shù)與初始條件 X10 = 0; X11 = 0; X12 = 0; V10 = 0; V11

20、= 0; V12 = 0; A10 = 0; A11 = 0; A12 = 6; M00 = 1; M01 = 0; M02 = 0;M10 = 0; M11 = 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; P

21、2=6;/開始計(jì)算，調(diào)用子程序Centre Centre(M,C,K,P,X,V,A,num,timeStep,dt); /輸出計(jì)算結(jié)果 ofstream out("out.txt"); /輸出結(jié)果文件out.txt if(!out.is_open() return; for(int i=0; i<=timeStep; i+) cout<<i*dt<<" s:"<<endl; /屏幕顯示 out<< i*dt<<" s:"<<endl; /將結(jié)果寫入結(jié)果文件

22、 for(int j=0; j<num; j+) cout << setw(20) << Xij << setw(20) << Vij << setw(20) << Aij << endl; out<<setw(20) << Xij << setw(20) << Vij << setw(20) << Aij << endl; out.close(); return 0;/子程序Centrevoid Centre(doubl

23、e*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; i<size; i+) MMi=new doublesize; double *L = new do

24、uble* size; for(int i=0; i<size; i+) Li = new double size; double *U = new double* size; for(int i=0; i<size; i+) Ui = new double size; /計(jì)算積分常數(shù) c0=1/(dt*dt); c1=1/(2*dt); c2=2/(dt*dt); c3=(dt*dt)/2;/起步問題，計(jì)算 X00=X10-dt*V10+c3*A10; X01=X11-dt*V11+c3*A11; X02=X12-dt*V12+c3*A12; /形成有效質(zhì)量矩陣for(int i

25、=0; i<size; i+) for(int j=0; j<size; j+) MMij=c0*Mij+c1*Cij; /調(diào)用子程序ArrayLU對(duì)有效質(zhì)量矩陣進(jìn)行LU分解 ArrayLU(MM,L,U,size); /計(jì)算有效載荷for(int ti=1; ti<=timeStep-1; ti+) for(int i=0; i<size; i+) tempi = c2*Xtii-c0*Xti-1i; /調(diào)用ArrayMVector，得到矩陣和向量乘積 ArrayMVector(M,temp,tran,size); for(int i=0; i<size; i+

26、) PTi = Pi+trani; for(int i=0; i<size; i+) tempi = c1*Xti-1i; ArrayMVector(C,temp,tran,size); for(int i=0; i<size; i+) PTi = PTi+trani; for(int i=0; i<size; i+) tempi = Xtii; ArrayMVector(K,temp,tran,size); for(int i=0; i<size; i+) PTi = PTi-trani;/計(jì)算時(shí)間的位移LUSolve(L,U,PT,temp,size); /輸出結(jié)果

27、：位移，速度，加速度for(int i=0; i<size; i+) Xti+1i = tempi; Atii = c0*(Xti-1i-2*Xtii+Xti+1i); Vtii = c1*(-Xti-1i+Xti+1i); /刪除不用數(shù)組for(int i=0; i<size; i+) delete MMi; delete MM; for(int i=0; i<size; i+) delete Li; delete L; for(int i=0; i<size; i+) delete Ui; delete U; delete PT; delete temp; dele

28、te tran; /子程序ArrayLU，對(duì)有效質(zhì)量矩陣進(jìn)行LU分解void ArrayLU(double*A, double*L, double*U, int size) double scale; double *B=new double*size; for(int i=0; i<size; i+) Bi=new double size; for(int i=0; i<size; i+) for(int j=0; j<size; j+) Bij = 0; Lij = 0; Uij = 0; Bii=1; for(int i=0; i<size-1; i+) if(A

29、ii=0) cout<<"the Array is singular"<<endl; for(int j=i+1; j<size; j+) scale = Aji/Aii; for(int k=i; k<size; k+) Ajk = Ajk-Aik*scale; for(int k=0; k<size; k+) Bjk = Bjk-Bik*scale; /調(diào)用子程序Inverse，對(duì)矩陣進(jìn)行轉(zhuǎn)置 Inverse(B, L, size); for(int i=0; i<size; i+) for(int j=i; j<size; j+) Uij = Aij; for(int i=0; i<size; i+) delete Bi; delete B; /子程序Inverse，對(duì)矩陣進(jìn)行轉(zhuǎn)置void Inverse(double*A, double*B, int size) double scale; for(int i=0; i<size; i+) Bii = 1; for(int i=0; i<size-1;