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

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

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

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

4、05,其中，時間采用50,000步，空間采用100步圖 3 穩(wěn)定格式計算得到溫度分布，s=0.005,其中，時間采用500,000步，空間采用100步圖 4 穩(wěn)定格式計算得到溫度分布，s=0.0005,其中，時間采用5,000,000步，空間采用100步這三個計算結(jié)果相對于之前的不穩(wěn)定計算結(jié)果，只是增大改變了時間步數(shù)。使得式得以滿足。但僅僅是這一條件的改變，使得之前所描述的定性結(jié)果是正確的。但是，這三個計算結(jié)果也是有微弱的差別的。僅僅看時，處，三個解算結(jié)果溫度值是不同的，分別為0.5226，0.5380，0.5395。盡管，這三個數(shù)值之間有微弱的差別，但是整體上來說趨近于0.54這個數(shù)值。而且

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

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

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

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

9、ear all% 1.定義參數(shù)alpha=2; % 熱傳導(dǎo)微分方程的系數(shù)T=0.5; % 計算時長，可修改D=2; % 計算區(qū)域?qū)挾萀1=500000; % 時間網(wǎng)格數(shù)L2=100; % 空間網(wǎng)格數(shù)s=1; % 差分格式的系數(shù)while s>1/2 % 減小時間步長，以保證s<=0.5，穩(wěn)定情況下適用 L1=L1*10; dt=T/L1; % 時間步長 dx=D/L2; % 空間步長 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.迭代計算下一時層的溫度for ii=1:L1 for jj=2:L2 u(jj)=s*(u(jj+1)-u(jj)-(u(jj)-u(jj-1)+u(jj); % 下一時層幾乎所有格點處的溫度值 end u(1)=0; % 下一時層x=0處的溫度值 u(L2+1)=(4*u(L2)-u(L2-1)/3; % 下一時層x=D處的溫度值end % 4.繪制計算結(jié)果 plot(x,u,'k-')title('穩(wěn)定格式下t=0.5s時刻的溫度分布，s=0.0005')xlabel('x/m'),yla

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