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

1、有限差分法的一個(gè)算例計(jì)算流體力學(xué)大作業(yè)作者：郝柏函2010011545指導(dǎo)：李嵩1. 題目編程計(jì)算熱傳導(dǎo)方程, 邊界條件：初始條件：（1） 用FTCS格式分別在滿足和不滿足穩(wěn)定性條件兩種情況下計(jì)算，給出結(jié)果比較和分析。（2）自選一種其他格式編程計(jì)算，并給出結(jié)果和分析注：原題中給出的初始條件與邊界條件是矛盾的，所以將其改為2. FTCS格式2.1. 計(jì)算方法2.1.1. 差分格式及其相容性對(duì)于方程，采用FTCS差分格式，即 其中，。以下討論這一格式的相容性。 所以因此，該格式與原微分方程是相容的，而且對(duì)于時(shí)間精度是一階的，對(duì)于空間，精度是二階的。2.1.2. 穩(wěn)定性與收斂性對(duì)于適定的線性微分方程

2、，格式如果差分格式，那么穩(wěn)定和收斂是等價(jià)的。所以只需要討論穩(wěn)定性就可以了。設(shè)則式可寫為 即 放大因子 所以為保證,應(yīng)有 只要滿足式，差分格式就是穩(wěn)定的。2.1.3. 初始、邊界條件處理以及全部計(jì)算過(guò)程初始條件，差分格式為 邊界條件，為了保證空間的二階精度，采用二次多項(xiàng)式來(lái)構(gòu)造差分格式，結(jié)果為，其中是方向上位置的格點(diǎn)數(shù)。給出了初始、邊界條件，以及之前的差分格式，就可以給出完整地算法：（1）首先用式計(jì)算出第一個(gè)時(shí)層的溫度；（2）然后使用式就算出下一時(shí)層的溫度值，但是，此時(shí)還沒(méi)有就算出，然后利用邊界條件求出（3）不斷使用第（2）步，直至計(jì)算出所要求時(shí)層所對(duì)應(yīng)的溫度值注：為了保證計(jì)算效率，不應(yīng)過(guò)小。如

3、果要求計(jì)算結(jié)果是穩(wěn)定的，應(yīng)滿足式，如果要求不穩(wěn)定，應(yīng)不滿足式。2.2. 計(jì)算結(jié)果與分析本文采用matlab編程，程序見(jiàn)于第4小節(jié)。 在不穩(wěn)定的差分格式下，計(jì)算結(jié)果是不可采信的，如圖 1所示。圖 1 不穩(wěn)定格式計(jì)算得到溫度分布，s=1.04,其中，時(shí)間采用2400步，空間采用100步采用不穩(wěn)定格式雖然也能得到比較光滑的溫度分布圖，但是，根據(jù)本算例的物理意義，左端為恒定溫度0，右端為絕熱壁面，所以計(jì)算結(jié)果應(yīng)該是，溫度始終大于0，別且距左端越近，溫度越低。可見(jiàn)，非穩(wěn)定格式的計(jì)算結(jié)果是定性錯(cuò)誤的。而穩(wěn)定格式的計(jì)算結(jié)果是可以采信的。如圖 2，圖 3，圖 4所示。圖 2 穩(wěn)定格式計(jì)算得到溫度分布，s=0.

4、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é)果相對(duì)于之前的不穩(wěn)定計(jì)算結(jié)果，只是增大改變了時(shí)間步數(shù)。使得式得以滿足。但僅僅是這一條件的改變，使得之前所描述的定性結(jié)果是正確的。但是，這三個(gè)計(jì)算結(jié)果也是有微弱的差別的。僅僅看時(shí)，處，三個(gè)解算結(jié)果溫度值是不同的，分別為0.5226，0.5380，0.5395。盡管，這三個(gè)數(shù)值之間有微弱的差別，但是整體上來(lái)說(shuō)趨近于0.54這個(gè)數(shù)值。而且

5、溫度分布的整體趨勢(shì)、數(shù)值之間的差距也幾乎為0.這說(shuō)明計(jì)算確實(shí)是穩(wěn)定的、收斂的。而且，時(shí)間步數(shù)越多，結(jié)果會(huì)越精確。3. FTCS隱格式3.1. 計(jì)算方法3.1.1. 差分格式及其相容性采用FTCS隱格式，即 其中，。以下討論這一格式的相容性。 所以其中，均在處取值因此，該格式與原微分方程是相容的，而且對(duì)于時(shí)間精度是一階的，對(duì)于空間，精度是二階的。3.1.2. 穩(wěn)定性與收斂性對(duì)于適定的線性微分方程，格式如果差分格式，那么穩(wěn)定和收斂是等價(jià)的。所以只需要討論穩(wěn)定性就可以了。設(shè)則式可寫為 即 放大因子 所以該差分格式是無(wú)條件穩(wěn)定的。3.1.3. 初始、邊界條件處理以及全部計(jì)算過(guò)程初始條件，差分格式為 邊

6、界條件，為了保證空間的二階精度，采用二次多項(xiàng)式來(lái)構(gòu)造差分格式，結(jié)果為，其中是方向上位置的格點(diǎn)數(shù)。給出了初始、邊界條件，以及之前的差分格式，就可以給出完整地算法：（1）首先用式計(jì)算出第一個(gè)時(shí)層的溫度；（2）然后使用式以及，組成一個(gè)m元一次方程組，對(duì)于這個(gè)方程組，以第n時(shí)層的溫度值為初值，采用迭代法計(jì)算，直至誤差足夠小。即其中，k為迭代次數(shù)。（3）不斷使用第（2）步，直至計(jì)算出所要求時(shí)層所對(duì)應(yīng)的溫度值。3.2. 計(jì)算結(jié)果與分析為了與之前的顯格式進(jìn)行比較，這里取空間步數(shù)為100，時(shí)間步數(shù)則分別取50,000， 500,000， 5,000,000，得到如所示的結(jié)果。每一個(gè)時(shí)間步長(zhǎng)的計(jì)算精度都是。這三

7、個(gè)結(jié)果都滿足：正溫度；越靠近原點(diǎn)處，溫度越低；原點(diǎn)處溫度為0；x=2處絕熱等定性的條件。而且計(jì)算結(jié)果并沒(méi)有很大的區(qū)別。三種計(jì)算結(jié)果下，在0.5s時(shí)，x=2處的溫度值分別為，0.5396，0.5396，0.5396，幾乎沒(méi)有什么區(qū)別。在其他點(diǎn)處，計(jì)算結(jié)果也是幾乎完全相同。其實(shí)，對(duì)于隱格式，時(shí)間步長(zhǎng)也不需要到50,000這么多。時(shí)間步數(shù)為5,000，在0.5s時(shí)，x=2處的溫度值已經(jīng)是0.5397，步數(shù)為500，該溫度值已經(jīng)是0.5399。圖 5 隱格式下t=0.5s時(shí)刻的溫度分布，時(shí)間步長(zhǎng)50,000圖 6 隱格式下t=0.5s時(shí)刻的溫度分布，時(shí)間步長(zhǎng)500,000圖 7 隱格式下t=0.5s時(shí)

8、刻的溫度分布，時(shí)間步長(zhǎng)5,000,000相較于顯格式，隱格式很快收斂到了0.5396，然而，顯格式即使步數(shù)為5,000,000，也只能收斂到0.5397。由此可見(jiàn)，隱格式具有如下優(yōu)越性 （1）無(wú)條件穩(wěn)定，不需要考慮時(shí)間步長(zhǎng)和空間步長(zhǎng)的關(guān)系；（2）收斂步數(shù)短。另外，我們從以上數(shù)據(jù)還可以看出，對(duì)本算例而言，有如下特點(diǎn)：顯格式的計(jì)算結(jié)果相對(duì)于精確值偏小，而隱格式的計(jì)算結(jié)果偏大。4. 本文程序在壓縮包中，應(yīng)該有.m文件的程序，fluid_conpution.m為顯格式的計(jì)算程序，fluid_conpution2.m為隱格式的計(jì)算程序。4.1. FTCS格式程序% FTCS格式計(jì)算熱傳導(dǎo)方程%clccl

9、ear all% 1.定義參數(shù)alpha=2; % 熱傳導(dǎo)微分方程的系數(shù)T=0.5; % 計(jì)算時(shí)長(zhǎng)，可修改D=2; % 計(jì)算區(qū)域?qū)挾萀1=500000; % 時(shí)間網(wǎng)格數(shù)L2=100; % 空間網(wǎng)格數(shù)s=1; % 差分格式的系數(shù)while s>1/2 % 減小時(shí)間步長(zhǎng)，以保證s<=0.5，穩(wěn)定情況下適用 L1=L1*10; dt=T/L1; % 時(shí)間步長(zhǎng) dx=D/L2; % 空間步長(zhǎng) s=alpha*dt/(dx)2end % dt=T/L1;% dx=D/L2;% s=alpha*dt/(dx)2 % 2.給出初始條件x=0:dx:D;u=sin(pi*x/4);plot(x,u

10、,'k-'),hold on % 3.迭代計(jì)算下一時(shí)層的溫度f(wàn)or ii=1:L1 for jj=2:L2 u(jj)=s*(u(jj+1)-u(jj)-(u(jj)-u(jj-1)+u(jj); % 下一時(shí)層幾乎所有格點(diǎn)處的溫度值 end u(1)=0; % 下一時(shí)層x=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'),yla

11、bel('u/oC')legend('t=0時(shí)刻的溫度分布','t=0.5s時(shí)刻的溫度分布')4.2. FTCS隱格式程序% FTCS隱格式計(jì)算熱傳導(dǎo)方程%clcclear all% 1.定義參數(shù)alpha=2; % 熱傳導(dǎo)微分方程的系數(shù)T=0.5; % 計(jì)算時(shí)長(zhǎng)，可修改D=2; % 計(jì)算區(qū)域?qū)挾萀1=5e2; % 時(shí)間網(wǎng)格數(shù)L2=100; % 空間網(wǎng)格數(shù)eps=1e-16; % 每一時(shí)層的精度要求dt=T/L1; % 時(shí)間步長(zhǎng)dx=D/L2; % 空間步長(zhǎng)s=alpha*dt/(dx)2% 2.給出初始條件x=0:dx:D;u=sin(pi*x/4);plot(x,u,'k-'),hold on % 3.迭代法計(jì)算下一時(shí)層的溫度 for ii=1:L1 epsilon=1; % 每一時(shí)層的精度 u(1)=0; % 下一時(shí)層x=0處的溫度值 w=u; while epsilon>eps v=w; for jj=2:L2 w(jj)=(s*(w(jj+1)+w(jj-1)+u(jj)/(1+2*s); % 下一時(shí)層幾乎所有格點(diǎn)處的溫度值 end w(L2+1)=(4*w(L2)-w(L2-1)/3; % 下一時(shí)層x=D處的溫度值 epsilon=max(abs(w-v);