計(jì)算流體力學(xué)算例

1、有限差分法的一個(gè)算例計(jì)算流體力學(xué)大作業(yè)作者：郝柏函2010011545指導(dǎo)：李嵩1. 題目編程計(jì)算熱傳導(dǎo)方程注：原題中給出的初始條件心o=sin才與邊界條件"L0是矛盾的，所以將其改uutx2，t 0,x0,2,-2小u0邊界條件：ux 00,x|x 2x初始條件：ut=0-si n4(1) 用FTCS格式n 1Ujnsuj !12s u；nsuj !un 1ujn2 u j 1x2u；nu 1t分別在滿足和不滿足穩(wěn)定性條件兩種情況下計(jì)算，給出結(jié)果比較和分析。(2)自選一種其他格式編程計(jì)算，并給出結(jié)果和分析為 Ut=0 =sin2. FTCS格式2.1 計(jì)算方法 2.1.1.差分格式

2、及其相容性對于方呈汁2胃，采用FTCS差分格式xn 1 nUj u jtx2nUj 12u；nUj 1 ，其中，sn 1nujsuj 112snUjnsUj 1(1.1)以下討論這一格式的相容性。n 1UjnUj12x所以nUj 12u：nUj2-2 oxnUj1 nUj tx2nnUj 1 2UjnUj 12U2xx2因此，間，精度是二階的。該格式與原微分方程是相容的,而且對于時(shí)間精度是一階的,對于空2.1.2 穩(wěn)定性與收斂性對于適定的線性微分方程，格式如果差分格式，那么穩(wěn)定和收斂是等價(jià)的。 所以只需要討論穩(wěn)定性就可以了設(shè) u； rnei j則(1.1)式可寫為i jrn 1ei jsrne

3、2s rne ji j 1srnerni se1 2si se放大因子rn 1se1 2si sern2s2scos2s 1 cos24ssin 所以1 4s G 1為保證1，應(yīng)有(1-2)只要滿足（1.2）式，差分格式就是穩(wěn)定的。2.1.3.初始、邊界條件處理以及全部計(jì)算過程初始條件u t=0 =sin差分格式為(1.3)邊界條件u % o 0,X x 20，為了保證空間的二階精度,采用二次多項(xiàng)式來構(gòu)造差分1 .5 =sinnn格式，結(jié)果為u：0,um應(yīng)，其中m是x方向上位置的格點(diǎn)數(shù)。3給出了初始、邊界條件，以及之前的差分格式，就可以給出完整地算法:（1）首先用（1.3）式計(jì)算出第一個(gè)時(shí)層的

4、溫度；（2）然后使用（1.1）式就算出下一時(shí)層的溫度值，但是，此時(shí)口畀山焉1還沒有就算出,n 1 n 1然后利用邊界條件求出u；11 o,um1 如1 吐3（3）不斷使用第（2 ）步，直至計(jì)算出所要求時(shí)層所對應(yīng)的溫度值注：為了保證計(jì)算效率，s 一2不應(yīng)過小。如果要求計(jì)算結(jié)果是穩(wěn)定的，應(yīng)滿足 （1.2）x式，如果要求不穩(wěn)定，應(yīng)不滿足（1.2）式。2.2 計(jì)算結(jié)果與分析本文米用matlab編程，程序見于第4小節(jié)。在不穩(wěn)定的差分格式下，計(jì)算結(jié)果是不可采信的，如圖1所示圖1不穩(wěn)定格式計(jì)算得到溫度分布，s=1.04,其中，時(shí)間采用 2400步，空間采用100步采用不穩(wěn)定格式雖然也能得到比較光滑的溫度分布

5、圖，但是，根據(jù)本算例的物理意義，左端為恒定溫度0，右端為絕熱壁面，所以計(jì)算結(jié)果應(yīng)該是，溫度始 終大于0,別且距左端越近，溫度越低?？梢姡欠€(wěn)定格式的計(jì)算結(jié)果是定性錯(cuò) 誤的。而穩(wěn)定格式的計(jì)算結(jié)果是可以采信的。如 圖2，圖3，圖4所示。圖2 穩(wěn)定格式計(jì)算得到溫度分布，S=0.05,其中，時(shí)間采用 50,000步，空間采用100步圖3穩(wěn)定格式計(jì)算得到溫度分布，S=0.005,其中，時(shí)間采用 500,000步，空間采用100步圖4 穩(wěn)定格式計(jì)算得到溫度分布，s=0.0005,其中，時(shí)間采用 5,000,000步，空間采用100步這三個(gè)計(jì)算結(jié)果相對于之前的不穩(wěn)定計(jì)算結(jié)果，只是增大改變了時(shí)間步數(shù)。使得(1

6、.2)式得以滿足。但僅僅是這一條件的改變，使得之前所描述的定性結(jié)果是正確的。但是，這三個(gè)計(jì)算結(jié)果也是有微弱的差別的。僅僅看t 0.5s時(shí)，x 2m處，三個(gè)解算結(jié)果溫度值是不同的，分別為0.5226，0.5380，0.5395。盡管，這三個(gè)數(shù)值之間有微弱的差別，但是整體上來說趨近于0.54這個(gè)數(shù)值。而且溫度分布的整體趨勢、數(shù)值之間的差距也幾乎為0.這說明計(jì)算確實(shí)是穩(wěn) 定的、收斂的。而且，時(shí)間步數(shù)越多，結(jié)果會越精確。3. FTCS隱格式3.1.計(jì)算方法3.1.1.差分格式及其相容性采用FTCS隱格式n 1 nUj u jtn 1nUj 1 2Uj xn 1Uj 1，即n 1Ujn 11 2s SU

7、j 1nn 1sUj 1(1.4)其中，s 2。以下討論這一格式的相容性。xn 1 nUjt12x所以n 1Uj 12u；n 1Uj 1其中,n 1Ujx2n 1Uj 1n 1Uj 1U2U因此，間，精度是二階的。,znU；1處取值t x該格式與原微分方程是相容的,均在而且對于時(shí)間精度是一階的，對于空3.1.2 穩(wěn)定性與收斂性對于適定的線性微分方程，格式如果差分格式，那么穩(wěn)定和收斂是等價(jià)的 所以只需要討論穩(wěn)定性就可以了。設(shè) u： rnei j則(1.1)式可寫為rn 1e '1sr：1 2s1eij1i ji j 1rnesrn 1e即1 2s2 s cosr： 1r：放大因子r： 1

8、1G 4r：1 2s2 s cos11 2s1 cos所以1 2s該差分格式是無條件穩(wěn)定的。3.1.3.初始、邊界條件處理以及全部計(jì)算過程初始條件U=sin-，差分格式為t 02(1.5)Uj sin 亠邊界條件u0,0,為了保證空間的二階精度，采用二次多項(xiàng)式來構(gòu)造差分X x 24 nn格式，結(jié)果為u ； 0,u m 如1 也，其中m是x方向上位置的格點(diǎn)數(shù)。3給出了初始、邊界條件，以及之前的差分格式，就可以給出完整地算法：(1) 首先用(1.3)式計(jì)算出第一個(gè)時(shí)層的溫度；n 1 n 1(2) 然后使用(1.4)式以及u； 1 0,um1m1 Um2，組成一個(gè)m元一次方程組,3對于這個(gè)方程組，以

9、第 n時(shí)層的溫度值為初值，采用迭代法計(jì)算，直至誤差足夠小。即n 1,0Uju；,j 1,2,L ,mn 1,kq0n 1,k 1Uj1n 1,knn 1,k 11 2s sUj 1 ujsuj 1, j 1,2,L ,m 1n 1,k 1Umn 1n 14um 1 um 23其中，k為迭代次數(shù)。(3)不斷使用第(2)步，直至計(jì)算出所要求時(shí)層所對應(yīng)的溫度值。3.2 計(jì)算結(jié)果與分析為了與之前的顯格式進(jìn)行比較，這里取空間步數(shù)為100，時(shí)間步數(shù)則分別取50,000，500,000，5,000,000，得到如所示的結(jié)果。每一個(gè)時(shí)間步長的計(jì)算精度都是1 10 16。這三個(gè)結(jié)果都滿足：正溫度；越靠近原點(diǎn)處，

10、溫度越低；原點(diǎn) 處溫度為0; x=2處絕熱等定性的條件。而且計(jì)算結(jié)果并沒有很大的區(qū)別。三種計(jì)算結(jié)果下，在0.5s時(shí)，x=2處的溫度值分別為，0.5396，0.5396，0.5396，幾乎沒有什么區(qū)別。在其他點(diǎn)處，計(jì)算結(jié)果也是幾乎完全相同。其實(shí)，對于隱格式，時(shí)間步長也不需要到50,000這么多。時(shí)間步數(shù)為5,000，在0.5s時(shí)，x=2處的溫度值已經(jīng)是0.5397 ,步數(shù)為500，該溫度值已經(jīng)是0.5399。0.90.806/030.20.60.812162x/m0.9D.80.7/06703D.20.60.812162x/m% D5°D°D< D5隱格式下冃晁時(shí)剖的溫

11、度分布、時(shí)間步長印,00C!匸0時(shí)刻的溫度分布匸m空時(shí)刻的溫屢好布圖6隱格式下t=0.5s時(shí)刻的溫度分布，時(shí)間步長500,000圖5隱格式下t=0.5s時(shí)刻的溫度分布，時(shí)間步長50,000隱格式下t=CL£g時(shí)刻的溫度弁布*時(shí)間步長500.0001珂）時(shí)刻的溫度分布 匸U.別時(shí)刻的溫度分布圖7 隱格式下t=0.5s時(shí)刻的溫度分布，時(shí)間步長5,000,000相較于顯格式，隱格式很快收斂到了0.5396，然而，顯格式即使步數(shù)為5,000,000 ,也只能收斂到0.5397。由此可見，隱格式具有如下優(yōu)越性（1）無條件穩(wěn)定，不需要考慮時(shí)間步長和空間步長的關(guān)系；（2）收斂步數(shù)短。另外，我們從以

12、上數(shù)據(jù)還可以看出，對本算例而言，有如下特點(diǎn)：顯格式的計(jì)算結(jié)果相對于精確值偏小，而隱格式的計(jì)算結(jié)果偏大。4. 本文程序在壓縮包中，應(yīng)該有 .m文件的程序， fluid_conpution.m 為顯格式的計(jì)算程序，fluid_conpution2.m為隱格式的計(jì)算程序。4.1. FTCS格式程序% % FTCS格式計(jì)算熱傳導(dǎo)方程% % clcclear all% 1.定義參數(shù) alpha=2;%熱傳導(dǎo)微分方程的系數(shù)T=0.5;%計(jì)算時(shí)長，可修改D=2;L1=500000;L2=100;s=1;while s>1/2適用L1=L1*10;dt=T/L1;dx=D/L2;s=alpha*dt/(

13、dxF2end% 計(jì)算區(qū)域?qū)挾? 時(shí)間網(wǎng)格數(shù)% 空間網(wǎng)格數(shù)% 差分格式的系數(shù)% 時(shí)間步長% 空間步長% 減小時(shí)間步長，以保證 s<=0.5 ，穩(wěn)定情況下% dt=T/L1;% dx=D/L2;% s=alpha*dt/(dx)A2 % 2. 給出初始條件 x=0:dx:D;u=sin(pi*x/4);plot(x,u,'k-'),hold on % 3. 迭代計(jì)算下一時(shí)層的溫度for ii=1:L1for jj=2:L2u(jj)=s*(u(jj+1)-u(jj)-(u(jj)-u(jj-1)+u(jj);% 下一時(shí)層幾乎所有格點(diǎn)處的溫度值endu(1)=0;%下一時(shí)層x

14、=0處的溫度值u(L2+1)=(4*u(L2)-u(L2-1)/3;% 下一時(shí)層 x=D處的溫度值end % 4. 繪制計(jì)算結(jié)果plot(x,u, 'k-' )title( ' 穩(wěn)定格式下 t=0.5s 時(shí)刻的溫度分布， s=0.0005' )xlabel( 'x/m' ),ylabel( 'u/AoC')legend( 't=0 時(shí)刻的溫度分布 ','t=0.5s 時(shí)刻的溫度分布 ' )4.2. FTCS 隱格式程序% FTCS急格式計(jì)算熱傳導(dǎo)方程%clcclear allon% 熱傳導(dǎo)微分方程的

15、系數(shù)% 計(jì)算時(shí)長，可修改% 計(jì)算區(qū)域?qū)挾? 時(shí)間網(wǎng)格數(shù)% 空間網(wǎng)格數(shù)% 每一時(shí)層的精度要求% 時(shí)間步長% 空間步長% 1. 定義參數(shù)alpha=2;T=0.5;D=2;L1=5e2;L2=100;eps=1e-16;dt=T/L1;dx=D/L2;s=alpha*dt/(dxF2% 2. 給出初始條件 x=0:dx:D;u=sin(pi*x/4);plot(x,u,'k-' ),hold % 3. 迭代法計(jì)算下一時(shí)層的溫度for ii=1:L1epsilon=1;% 每一時(shí)層的精度u(1)=0;%下一時(shí)層x=0處的溫度值w=u;while epsilon>epsv=w;for jj=2:L2w(jj)=(s*(w(jj+1)+w(jj-1)+u(jj)/(1+2*s);% 下一時(shí)層幾乎所有格點(diǎn)處的溫度值endw(L2+1)=(4*w(L2)-w(L2-1)/3;