中科大計(jì)算流體力學(xué)CFD

1、CFD 實(shí)驗(yàn)報(bào)告姓名：學(xué)號(hào)、題目： 利用中心差分格式近似導(dǎo)數(shù)d2y/dx2，數(shù)值求解常微分方程(0 x 1)如=sin2 (0 x 1)dx2x=x=1sin 2T步長(zhǎng)分別取 Ax =0.05, 0.01, 0.001, 0.0001。二、報(bào)告要求：列出全部計(jì)算公式和步驟；表列出程序中各主要符號(hào)和數(shù)組意義；繪出數(shù)值計(jì)算結(jié)果的函數(shù)曲線,并與精確解比較；比較不同差分格式和不同網(wǎng)格步長(zhǎng)計(jì)算結(jié)果的精度和代價(jià)附源程序。三、相關(guān)差分格式二階導(dǎo)數(shù)d2y/dx2的三點(diǎn)差分格式有向前差分、向后差分和中心差分，表達(dá)式分別如下：一階向前差分： 竺上= t2j+i *j + O(Ax) TOC o 1-5 h z d

2、x2Ax2一階向后差分：列=.j-1 + j + O (Ax)dx2Ax2二階中心差分：紐=+廠乞+1 + O (Ax 2)dx2Ax2代入微分方程可以得到差分方程,表達(dá)式分別如下：一階向前差分一階向后差分二階中心差分u 2u + u一階向前差分一階向后差分二階中心差分/1j = sin 2 x HYPERLINK l bookmark12 o Current Document Ax2ju 2u + ui2jj = sin 2 x HYPERLINK l bookmark16 o Current Document Ax2ju 2u + ujji = sin 2 x HYPERLINK l bo

3、okmark20 o Current Document Ax2j對(duì)于三種差分格式，差分格式可以改寫(xiě)成AY = b的形式，其中A是相同的,非齊次項(xiàng)b不同，如下所示:b=1b=3-211-2 11 -211b=1b=3-211-2 11 -211-2系數(shù)矩陣ysin（2 xAx 21sin（2 x） Ax 2sin（2x） Ax2 - yk - 21sin（2x ）Ax2 - ysin（2x ）Ax23sin（2 x） Ax 2sin （2x ）-Ax2 - y I k1sin（2x ） Ax2 - ysin（2x ）Ax22sin（2 x） Ax 2kv 2sin（2x ）Ax2 - yk-11

4、一階向前差分一階向后差分二階中心差分yk - 2人丿。求解AY = byk - 2人丿。12四、計(jì)算公式和步驟；關(guān)于精確解的推導(dǎo)：已知士 = si 2 x，對(duì)X進(jìn)行兩次積分得到y(tǒng) = -4sin2x + C1x + C2，再結(jié)合x(chóng)=1邊界條件yx=0 = 0和y“ =1 -晉得到相對(duì)應(yīng)的-和Cx=1關(guān)于數(shù)值求解方法：對(duì)于方程組AY = b可直接求解，也可以使用追趕法求解，下面介紹簡(jiǎn)單追 趕法求解三對(duì)角方程組的過(guò)程。對(duì)于三對(duì)角方程組：dcyb1111adcb22 2 22adcybn1n1n1n1n1adybn n n n系數(shù)矩陣A中僅三對(duì)角線上的數(shù)值不全為0,其余位置上的數(shù)值全為0,是典 型的

5、對(duì)角占優(yōu)的三對(duì)角矩陣，利用高斯消去法，經(jīng)過(guò)n-1次消元可以化為同解方 程組：1 u11u2 y1 y2q1巳21 uyqn1n1n11ynqn如上所示，求這些值的過(guò)程（即消元過(guò)程）稱(chēng)為追：u1uiqiu1uiqi=c / d , q = b / d ,1 / 1 1 =c / Id u a 丿=(b q a )/ (dii 1 ii(i = 2,., n 1)u a )i 1 i(i = 2,., n )再利用回代過(guò)程求出方程組的各變量：y = q ,(i (i = n 1,n 2,.,1)y = q u yi i i i+1這一逆序求變量的過(guò)程（回代過(guò)程）稱(chēng)為趕。五、計(jì)算結(jié)果與分析根據(jù)題意需

6、要分別取步長(zhǎng)Ax =0.05, 0.01, 0.001, 0.0001進(jìn)行計(jì)算，因此采用MATLAB 進(jìn)行編程運(yùn)算，然后將數(shù)值解與精確解進(jìn)行比較，如下圖所示：步長(zhǎng) Ax =0.05步 長(zhǎng) Ax =0.01步長(zhǎng) Ax =0.001步長(zhǎng) Ax =0.0001從上面四個(gè)圖可以看出，這幾個(gè)格式的精確解和數(shù)值解之間符合度很好，其 中中心差分格式精確度更高，并且隨著步長(zhǎng)的減小與精確解的符合度越高。下面 將計(jì)算結(jié)果用表格的形式列出：追趕法求解Ax =0.05Ax =0.01Ax =0.001Ax =0.0001計(jì)算耗時(shí)、/.、八前差0.0068530.0148030.0489311.111667后差0.00

7、66450.0141160.0486621.103086中心差0.0055260.0158000.0596941.072931表1不同格式不同步長(zhǎng)計(jì)算耗時(shí)（單位s）追趕法求解Ax =0.05Ax =0.01Ax =0.001Ax =0.0001計(jì)算精度、/.、八前差4.6976e-048.5578e-058.3727e-068.3542e-07后差3.6746e-048.1482e-058.3317e-068.3501e-07中心差6.8119e-085.4443e-105.4441e-135.3257e-16表2 不同格式不同步長(zhǎng)計(jì)算精度由不同格式不同步長(zhǎng)計(jì)算耗時(shí)和計(jì)算精度的結(jié)果對(duì)比可知：1

8、、從時(shí)間代價(jià)來(lái)看，三種格式采用追趕法求解所需要的時(shí)間基本一樣，從計(jì)算 精度來(lái)看，中心差分格式的計(jì)算精度明顯優(yōu)于向前差分和向后差分格式。2、隨著步長(zhǎng)的減小，雖然求解的計(jì)算精度增高，但所需要的時(shí)間卻大大增加。3、對(duì)于這三種格式而言，中心差分格式明顯優(yōu)于另兩種格式。4、從上可知，提高計(jì)算精度可以通過(guò)減少步長(zhǎng)和采用高階格式的方法，由于減 少步長(zhǎng)需要用時(shí)間作為代價(jià)，所以采用高階格式相對(duì)而言比較好。六、源程序及其說(shuō)明符號(hào)說(shuō)明：long變量取值范圍dx步長(zhǎng)n+1節(jié)點(diǎn)數(shù)y0始端邊界值yl末端邊界值A(chǔ)系數(shù)矩陣b非齊次項(xiàng)向量Y數(shù)值解向量Y e精確解向量T E節(jié)點(diǎn)截?cái)嗾`差Q精度源程序如下：clear; %清空l(shuí)on

9、g=l;%變量取值范圍dx=0.05; %設(shè)置步長(zhǎng)n=long/dx; % 此處n=節(jié)點(diǎn)數(shù)一1 y0=0;%設(shè)定始端邊界條件 yl=1-sin(2)/4; %設(shè)定末端邊界條件A=zeros(n-1,n-1);%設(shè)置系數(shù)矩陣 AA(1,1:2)=-2 1;%始端邊界for i=2:n-2A(i,i-1:i+1)=1,-2,1;endA(n-1,n-2:n-1)=1,-2;% 末端邊界b=zeros(n-1,1);%設(shè)置矩陣 bfor i=1:n-1%b(i)=dxA2*sin(i-1)*2*dx); %階向前差分%b(i)=dxA2*sin(i+1)*2*dx); %階向后差分b(i)=dxA2

10、*sin(i*2*dx); %中心差分endb(1)=b-yO;%設(shè)置邊界條件b(n-1)=b(n-1)-yl;x=dx*(0:n);%求解精確解Y_e=x-sin(2*x)/4;);disp(SC16OO5O41 李文建);disp( );tic; % 追趕法求解并開(kāi)始計(jì)時(shí)a=1; c=-2;u(1)=a/c; q(1)=b(1)/c; for i=2:n-2u(i)=a/(c-u(i-1)*a);endfor i=2:n-1q(i)=(b(i)-q(i-1)/(c-u(i-1)*a);end y(n-1)=q(n-1); for i=n-2:-1:1y(i)=q(i)-u(i)*y(i+1);endY(2:n)=y(l:n-1);