空氣動(dòng)力學(xué)數(shù)值方法：有限體積法(FVM)：二維氣體動(dòng)力學(xué)方程的FVM求解1空氣動(dòng)力學(xué)數(shù)值方法：有限體積法(FVM)：二維氣體動(dòng)力學(xué)方程的FVM求解1.1緒論1.1.1有限體積法的起源與應(yīng)用有限體積法(FiniteVolumeMethod,FVM)是一種廣泛應(yīng)用于流體力學(xué)、熱力學(xué)和空氣動(dòng)力學(xué)領(lǐng)域的數(shù)值方法。它的起源可以追溯到20世紀(jì)50年代，最初是為了求解流體動(dòng)力學(xué)中的偏微分方程而發(fā)展起來的。FVM的核心思想是基于控制體的概念，將計(jì)算域劃分為一系列非重疊的體積，然后在每個(gè)體積上應(yīng)用守恒定律，從而將連續(xù)的偏微分方程轉(zhuǎn)化為離散的代數(shù)方程組。這種方法不僅能夠保證守恒性，還具有較好的穩(wěn)定性和準(zhǔn)確性，因此在處理復(fù)雜的流體動(dòng)力學(xué)問題時(shí)表現(xiàn)出色。在空氣動(dòng)力學(xué)中，F(xiàn)VM被用于求解二維和三維的氣體動(dòng)力學(xué)方程，包括歐拉方程和納維-斯托克斯方程。這些方程描述了氣體在不同條件下的運(yùn)動(dòng)，如速度、壓力、溫度和密度的變化。通過FVM，工程師和科學(xué)家能夠模擬飛機(jī)、火箭等飛行器在不同飛行條件下的氣動(dòng)特性，為設(shè)計(jì)和優(yōu)化提供關(guān)鍵數(shù)據(jù)。1.1.2維氣體動(dòng)力學(xué)方程簡(jiǎn)介二維氣體動(dòng)力學(xué)方程通常包括連續(xù)性方程、動(dòng)量方程和能量方程，它們是描述氣體在二維空間中運(yùn)動(dòng)的基本方程。這些方程可以寫成守恒形式，即：?其中，U是狀態(tài)向量，包含了密度ρ、動(dòng)量ρu、ρv和總能量E；F和G分別是x和在有限體積法中，這些方程被應(yīng)用于每個(gè)控制體，通過積分形式轉(zhuǎn)化為離散方程。例如，對(duì)于一個(gè)控制體，連續(xù)性方程可以轉(zhuǎn)化為：d其中，Ω是控制體的體積，Γ是控制體的邊界，nx和n1.2有限體積法求解二維氣體動(dòng)力學(xué)方程1.2.1網(wǎng)格劃分在應(yīng)用FVM求解二維氣體動(dòng)力學(xué)方程之前，首先需要對(duì)計(jì)算域進(jìn)行網(wǎng)格劃分。網(wǎng)格可以是結(jié)構(gòu)化的（如矩形網(wǎng)格）或非結(jié)構(gòu)化的（如三角形網(wǎng)格）。下面是一個(gè)使用Python和matplotlib庫進(jìn)行簡(jiǎn)單矩形網(wǎng)格劃分的例子：importnumpyasnp

importmatplotlib.pyplotasplt

#定義網(wǎng)格參數(shù)

nx=50#x方向網(wǎng)格數(shù)

ny=50#y方向網(wǎng)格數(shù)

Lx=1.0#x方向計(jì)算域長(zhǎng)度

Ly=1.0#y方向計(jì)算域長(zhǎng)度

#創(chuàng)建網(wǎng)格

x=np.linspace(0,Lx,nx+1)

y=np.linspace(0,Ly,ny+1)

X,Y=np.meshgrid(x,y)

#繪制網(wǎng)格

plt.figure(figsize=(8,6))

plt.plot(X,Y,'k-',lw=0.5)

plt.plot(X.T,Y.T,'k-',lw=0.5)

plt.xlabel('x')

plt.ylabel('y')

plt.title('2DStructuredMesh')

plt.show()1.2.2離散化過程離散化過程是將連續(xù)的氣體動(dòng)力學(xué)方程轉(zhuǎn)化為離散形式的關(guān)鍵步驟。在FVM中，這通常通過在每個(gè)控制體上應(yīng)用積分守恒定律來實(shí)現(xiàn)。下面是一個(gè)離散化連續(xù)性方程的例子：假設(shè)我們有一個(gè)矩形控制體，其邊界為xi?Δx/2到d其中，ρi,j是控制體i,j內(nèi)的平均密度，ρui1.2.3時(shí)間推進(jìn)在離散化方程后，需要使用時(shí)間推進(jìn)方法來求解方程。常見的方法包括顯式歐拉法、隱式歐拉法和Runge-Kutta法。下面是一個(gè)使用顯式歐拉法進(jìn)行時(shí)間推進(jìn)的例子：假設(shè)我們已經(jīng)離散化了連續(xù)性方程，動(dòng)量方程和能量方程，并將它們寫為：U其中，Ui,jn是狀態(tài)向量在時(shí)間n和位置i,j的值，Δt是時(shí)間步長(zhǎng)，F(xiàn)i+1.2.4邊界條件處理邊界條件是任何數(shù)值模擬中不可或缺的一部分，它們定義了計(jì)算域邊緣的物理行為。在空氣動(dòng)力學(xué)中，常見的邊界條件包括固壁邊界、進(jìn)氣邊界、排氣邊界和對(duì)稱邊界。下面是一個(gè)處理固壁邊界條件的例子：假設(shè)我們有一個(gè)固壁邊界，其速度為0。在離散化方程時(shí)，我們需要將固壁邊界上的速度通量設(shè)為0。例如，對(duì)于x方向上的固壁邊界，我們有：ρ對(duì)于y方向上的固壁邊界，我們有：ρ1.2.5求解算法求解算法是將離散方程和邊界條件結(jié)合起來，通過迭代求解狀態(tài)向量的過程。常見的算法包括迭代法、直接法和多網(wǎng)格法。下面是一個(gè)使用迭代法求解離散方程的例子：假設(shè)我們已經(jīng)將所有方程和邊界條件離散化，并將它們寫為一個(gè)代數(shù)方程組。我們可以通過迭代法求解這個(gè)方程組，直到滿足收斂條件。例如，我們可以使用以下偽代碼：#初始化狀態(tài)向量

U=initialize(U)

#設(shè)置迭代參數(shù)

max_iterations=1000

tolerance=1e-6

residual=1.0

#迭代求解

foriterationinrange(max_iterations):

#計(jì)算界面通量

F,G=calculate_fluxes(U)

#應(yīng)用時(shí)間推進(jìn)

U=time_advance(U,F,G,dt)

#計(jì)算殘差

residual=calculate_residual(U)

#檢查收斂

ifresidual<tolerance:

break

#輸出結(jié)果

output(U)在這個(gè)例子中，initialize函數(shù)用于初始化狀態(tài)向量，calculate_fluxes函數(shù)用于計(jì)算界面通量，time_advance函數(shù)用于應(yīng)用時(shí)間推進(jìn)，calculate_residual函數(shù)用于計(jì)算殘差，output函數(shù)用于輸出結(jié)果。1.3結(jié)論有限體積法在空氣動(dòng)力學(xué)數(shù)值模擬中扮演著重要角色，它能夠準(zhǔn)確、穩(wěn)定地求解復(fù)雜的氣體動(dòng)力學(xué)方程。通過網(wǎng)格劃分、離散化過程、時(shí)間推進(jìn)、邊界條件處理和求解算法，工程師和科學(xué)家能夠模擬和分析飛行器在不同飛行條件下的氣動(dòng)特性，為飛行器設(shè)計(jì)和優(yōu)化提供關(guān)鍵數(shù)據(jù)。雖然這里只提供了一個(gè)簡(jiǎn)化的例子，實(shí)際的FVM求解過程可能涉及更復(fù)雜的方程和算法，但基本原理和步驟是相同的。2有限體積法基礎(chǔ)2.1控制體積的概念在空氣動(dòng)力學(xué)中，有限體積法（FVM）是一種廣泛使用的數(shù)值方法，用于求解流體動(dòng)力學(xué)方程。FVM的核心思想是將連續(xù)的流體域離散化為一系列控制體積，每個(gè)控制體積內(nèi)流體的物理量（如速度、壓力、密度）被視為常數(shù)。這種方法基于守恒定律，確保了質(zhì)量、動(dòng)量和能量在每個(gè)控制體積內(nèi)的守恒。2.1.1舉例說明考慮一個(gè)二維流體域，我們可以將其劃分為許多小的矩形區(qū)域，每個(gè)矩形區(qū)域即為一個(gè)控制體積。在每個(gè)控制體積內(nèi)，我們計(jì)算流體的平均物理量。例如，對(duì)于一個(gè)控制體積，其位置由xi,yj表示，我們可以定義其內(nèi)的平均密度為ρi2.2守恒定律與積分形式方程流體動(dòng)力學(xué)中的守恒定律包括質(zhì)量守恒、動(dòng)量守恒和能量守恒。在有限體積法中，這些守恒定律被表示為積分形式的方程，即在每個(gè)控制體積內(nèi)，物理量的總變化等于通過控制體積邊界流入或流出的物理量的凈變化。2.2.1質(zhì)量守恒方程質(zhì)量守恒方程描述了流體質(zhì)量在控制體積內(nèi)的變化。在二維情況下，質(zhì)量守恒方程可以表示為：?在有限體積法中，這個(gè)方程被轉(zhuǎn)化為積分形式：d其中，V是控制體積，S是控制體積的邊界，n是邊界上的外法向量。2.2.2動(dòng)量守恒方程動(dòng)量守恒方程描述了流體動(dòng)量在控制體積內(nèi)的變化。在二維情況下，沿x和y方向的動(dòng)量守恒方程可以表示為：??其中，p是壓力，fx和f在有限體積法中，這兩個(gè)方程也被轉(zhuǎn)化為積分形式：dd2.2.3能量守恒方程能量守恒方程描述了流體能量在控制體積內(nèi)的變化。在二維情況下，能量守恒方程可以表示為：?其中，E是總能量，q是熱傳導(dǎo)率。在有限體積法中，這個(gè)方程也被轉(zhuǎn)化為積分形式：d2.2.4數(shù)值離散化將上述積分方程離散化，我們得到每個(gè)控制體積的離散方程。例如，對(duì)于質(zhì)量守恒方程，離散化后的方程可以表示為：d其中，Δx和Δy是控制體積在x和y方向上的尺寸，ρu2.2.5Python代碼示例下面是一個(gè)使用Python和NumPy庫實(shí)現(xiàn)的簡(jiǎn)單二維有限體積法求解質(zhì)量守恒方程的示例。在這個(gè)例子中，我們假設(shè)流體在一個(gè)二維矩形域內(nèi)流動(dòng)，域被劃分為10×importnumpyasnp

#定義流體域的尺寸和控制體積的數(shù)目

Lx=1.0#流體域在x方向上的長(zhǎng)度

Ly=1.0#流體域在y方向上的長(zhǎng)度

nx=10#控制體積在x方向上的數(shù)目

ny=10#控制體積在y方向上的數(shù)目

#計(jì)算控制體積的尺寸

dx=Lx/nx

dy=Ly/ny

#初始化密度和速度

rho=np.zeros((nx,ny))

u=np.zeros((nx+1,ny))

v=np.zeros((nx,ny+1))

#設(shè)置初始條件

rho[:,:]=1.0#初始密度為1.0

u[0,:]=1.0#左邊界速度為1.0

v[:,0]=1.0#下邊界速度為1.0

#定義時(shí)間步長(zhǎng)和迭代次數(shù)

dt=0.01

nt=100

#迭代求解質(zhì)量守恒方程

forninrange(nt):

#更新密度

foriinrange(nx):

forjinrange(ny):

rho[i,j]+=dt*((u[i+1,j]-u[i,j])/dx+(v[i,j+1]-v[i,j])/dy)

#輸出最終的密度分布

print(rho)在這個(gè)例子中，我們使用了顯式歐拉方法來更新密度。在實(shí)際應(yīng)用中，可能需要使用更復(fù)雜的數(shù)值方法來確保解的穩(wěn)定性和準(zhǔn)確性。2.2.6結(jié)論有限體積法是一種強(qiáng)大的數(shù)值方法，用于求解流體動(dòng)力學(xué)方程。通過將流體域劃分為控制體積，并在每個(gè)控制體積內(nèi)應(yīng)用守恒定律，我們可以得到一系列離散方程，這些方程可以使用數(shù)值方法求解。在本教程中，我們介紹了控制體積的概念，以及如何將質(zhì)量、動(dòng)量和能量守恒方程轉(zhuǎn)化為積分形式，并進(jìn)行了數(shù)值離散化。我們還提供了一個(gè)使用Python和NumPy庫實(shí)現(xiàn)的簡(jiǎn)單示例，用于求解二維流體域內(nèi)的質(zhì)量守恒方程。請(qǐng)注意，上述代碼示例僅用于說明目的，實(shí)際應(yīng)用中可能需要更復(fù)雜的邊界條件處理和數(shù)值穩(wěn)定性控制。此外，對(duì)于動(dòng)量和能量守恒方程的求解，可能需要引入額外的物理模型，如狀態(tài)方程和湍流模型。3離散化技術(shù)3.1網(wǎng)格生成與類型在有限體積法(FVM)中，網(wǎng)格的生成是求解二維氣體動(dòng)力學(xué)方程的基礎(chǔ)。網(wǎng)格可以分為結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格。3.1.1結(jié)構(gòu)化網(wǎng)格結(jié)構(gòu)化網(wǎng)格通常由規(guī)則的四邊形或六面體組成，每個(gè)網(wǎng)格單元的形狀和大小在計(jì)算域內(nèi)是已知的。這種網(wǎng)格便于編程和處理，但在復(fù)雜幾何形狀的邊界附近可能需要大量的網(wǎng)格單元來準(zhǔn)確表示。3.1.2非結(jié)構(gòu)化網(wǎng)格非結(jié)構(gòu)化網(wǎng)格由不規(guī)則的四邊形、三角形、五面體或四面體組成，適用于復(fù)雜幾何形狀的邊界。這種網(wǎng)格的靈活性高，但在處理時(shí)需要更多的數(shù)據(jù)結(jié)構(gòu)和算法支持。3.2通量的離散化方法在有限體積法中，通量的離散化是關(guān)鍵步驟，它涉及到如何在網(wǎng)格單元之間傳遞信息。主要的離散化方法包括中心差分法和上風(fēng)差分法。3.2.1中心差分法中心差分法基于網(wǎng)格單元中心點(diǎn)的值來計(jì)算通量。這種方法簡(jiǎn)單，但在非光滑流場(chǎng)中可能會(huì)產(chǎn)生數(shù)值振蕩。3.2.1.1示例代碼#假設(shè)我們有網(wǎng)格單元i和i+1的保守變量U

defcentral_difference_flux(U_i,U_i_plus_1,dx):

"""

計(jì)算網(wǎng)格單元i和i+1之間的中心差分通量。

:paramU_i:網(wǎng)格單元i的保守變量

:paramU_i_plus_1:網(wǎng)格單元i+1的保守變量

:paramdx:網(wǎng)格單元的寬度

:return:通量

"""

#計(jì)算通量的中心差分

flux=(U_i_plus_1-U_i)/(2*dx)

returnflux3.2.2上風(fēng)差分法上風(fēng)差分法考慮了流體的流動(dòng)方向，使用上風(fēng)側(cè)的網(wǎng)格單元值來計(jì)算通量。這種方法可以減少數(shù)值振蕩，但在逆風(fēng)方向可能會(huì)產(chǎn)生數(shù)值擴(kuò)散。3.2.2.1示例代碼#假設(shè)我們有網(wǎng)格單元i和i+1的保守變量U，以及流體速度v

defupwind_difference_flux(U_i,U_i_plus_1,v,dx):

"""

計(jì)算網(wǎng)格單元i和i+1之間的上風(fēng)差分通量。

:paramU_i:網(wǎng)格單元i的保守變量

:paramU_i_plus_1:網(wǎng)格單元i+1的保守變量

:paramv:流體速度

:paramdx:網(wǎng)格單元的寬度

:return:通量

"""

ifv>0:

#正向流動(dòng)，使用網(wǎng)格單元i的值

flux=(U_i_plus_1-U_i)/dx

else:

#反向流動(dòng)，使用網(wǎng)格單元i+1的值

flux=(U_i-U_i_plus_1)/dx

returnflux3.2.3數(shù)據(jù)樣例假設(shè)我們有以下網(wǎng)格單元的保守變量U和流體速度v：網(wǎng)格單元Uvi-11.00.5i1.50.3i+12.00.23.2.3.1中心差分法計(jì)算通量#網(wǎng)格單元寬度

dx=1.0

#計(jì)算網(wǎng)格單元i和i+1之間的通量

flux_i_to_i_plus_1=central_difference_flux(1.5,2.0,dx)

print("網(wǎng)格單元i和i+1之間的中心差分通量:",flux_i_to_i_plus_1)3.2.3.2上風(fēng)差分法計(jì)算通量#計(jì)算網(wǎng)格單元i和i+1之間的通量

flux_i_to_i_plus_1=upwind_difference_flux(1.5,2.0,0.2,dx)

print("網(wǎng)格單元i和i+1之間的上風(fēng)差分通量:",flux_i_to_i_plus_1)通過上述示例，我們可以看到中心差分法和上風(fēng)差分法在計(jì)算通量時(shí)的不同處理方式，以及如何在Python中實(shí)現(xiàn)這些方法。在實(shí)際應(yīng)用中，選擇合適的離散化方法對(duì)于獲得準(zhǔn)確的數(shù)值解至關(guān)重要。4數(shù)值通量計(jì)算4.1Godunov方法Godunov方法是一種保守的數(shù)值方法，用于求解雙曲型守恒律方程。在空氣動(dòng)力學(xué)中，它被廣泛應(yīng)用于求解氣體動(dòng)力學(xué)方程。Godunov方法的核心思想是使用精確或近似的Riemann問題解來計(jì)算界面通量，從而確保數(shù)值解的保守性和穩(wěn)定性。4.1.1原理Godunov方法基于以下步驟：網(wǎng)格離散化：將計(jì)算域劃分為一系列控制體積，每個(gè)控制體積由網(wǎng)格節(jié)點(diǎn)定義。狀態(tài)重構(gòu)：在每個(gè)時(shí)間步內(nèi)，使用控制體積內(nèi)的平均狀態(tài)來近似網(wǎng)格節(jié)點(diǎn)處的狀態(tài)。Riemann問題求解：在每個(gè)網(wǎng)格界面，基于左右兩側(cè)的狀態(tài)，求解Riemann問題，得到界面通量。更新狀態(tài)：使用界面通量和控制體積內(nèi)的狀態(tài)，通過時(shí)間離散化更新下一個(gè)時(shí)間步的狀態(tài)。4.1.2示例假設(shè)我們有以下一維氣體動(dòng)力學(xué)方程：?其中，U=ρ,ρu,ET是保守變量向量，ρ是密度，u是速度，E使用Godunov方法，我們首先需要在每個(gè)網(wǎng)格界面求解Riemann問題。以下是一個(gè)使用Python實(shí)現(xiàn)的Godunov方法的簡(jiǎn)化示例：importnumpyasnp

defgodunov_flux(U_left,U_right,gamma):

"""

計(jì)算Godunov界面通量

:paramU_left:左側(cè)狀態(tài)向量

:paramU_right:右側(cè)狀態(tài)向量

:paramgamma:比熱比

:return:界面通量

"""

rho_left,u_left,E_left=U_left

rho_right,u_right,E_right=U_right

p_left=(gamma-1)*(E_left-0.5*rho_left*u_left**2)

p_right=(gamma-1)*(E_right-0.5*rho_right*u_right**2)

#求解Riemann問題，這里簡(jiǎn)化為使用左右狀態(tài)的最小值和最大值

flux=np.zeros(3)

flux[0]=min(u_left,0)*rho_left+max(u_right,0)*rho_right

flux[1]=min(u_left,0)*(rho_left*u_left+p_left)+max(u_right,0)*(rho_right*u_right+p_right)

flux[2]=min(u_left,0)*(E_left+p_left)+max(u_right,0)*(E_right+p_right)

returnflux

#示例數(shù)據(jù)

U_left=np.array([1.0,0.5,2.5])#左側(cè)狀態(tài)向量

U_right=np.array([0.5,0.7,2.0])#右側(cè)狀態(tài)向量

gamma=1.4#比熱比

#計(jì)算界面通量

flux=godunov_flux(U_left,U_right,gamma)

print("Godunov界面通量:",flux)4.1.3描述上述代碼示例中，我們定義了一個(gè)godunov_flux函數(shù)，它接受左側(cè)和右側(cè)的狀態(tài)向量以及比熱比作為輸入，返回界面通量。這里我們簡(jiǎn)化了Riemann問題的求解，僅使用了左右狀態(tài)的最小值和最大值來近似界面通量，這在實(shí)際應(yīng)用中可能需要更復(fù)雜的Riemann求解器。4.2Roe平均與Roe通量Roe平均是一種用于計(jì)算界面通量的平均狀態(tài)方法，它基于Riemann問題的線性化，通過計(jì)算左右狀態(tài)的Roe平均狀態(tài)來簡(jiǎn)化通量計(jì)算。Roe通量是基于Roe平均狀態(tài)計(jì)算的界面通量。4.2.1原理Roe平均狀態(tài)的計(jì)算基于以下步驟：Roe平均速度和壓力：計(jì)算左右狀態(tài)的平均速度和壓力。Roe平均矩陣：基于平均狀態(tài)，構(gòu)建Roe平均矩陣，該矩陣用于線性化方程。特征分解：對(duì)Roe平均矩陣進(jìn)行特征分解，得到特征值和特征向量。計(jì)算界面通量：使用特征值和特征向量，以及左右狀態(tài)的差值，計(jì)算界面通量。4.2.2示例以下是一個(gè)使用Python實(shí)現(xiàn)的Roe平均與Roe通量計(jì)算的簡(jiǎn)化示例：defroe_average(U_left,U_right):

"""

計(jì)算Roe平均狀態(tài)

:paramU_left:左側(cè)狀態(tài)向量

:paramU_right:右側(cè)狀態(tài)向量

:return:Roe平均狀態(tài)向量

"""

rho_left,u_left,E_left=U_left

rho_right,u_right,E_right=U_right

p_left=(gamma-1)*(E_left-0.5*rho_left*u_left**2)

p_right=(gamma-1)*(E_right-0.5*rho_right*u_right**2)

#計(jì)算Roe平均速度和壓力

u_roe=(u_left+u_right)/2

p_roe=np.sqrt(p_left*p_right)

rho_roe=np.sqrt(rho_left*rho_right)

#構(gòu)建Roe平均矩陣

A_roe=np.array([[0,rho_roe,0],

[gamma*p_roe/rho_roe,0,-gamma*p_roe],

[0,rho_roe*u_roe,0]])

returnA_roe,u_roe,p_roe,rho_roe

defroe_flux(U_left,U_right,A_roe,u_roe,p_roe,rho_roe):

"""

計(jì)算Roe界面通量

:paramU_left:左側(cè)狀態(tài)向量

:paramU_right:右側(cè)狀態(tài)向量

:paramA_roe:Roe平均矩陣

:paramu_roe:Roe平均速度

:paramp_roe:Roe平均壓力

:paramrho_roe:Roe平均密度

:return:界面通量

"""

#計(jì)算左右狀態(tài)的差值

delta_U=U_right-U_left

#計(jì)算界面通量

flux=0.5*(np.dot(A_roe,delta_U)+np.abs(u_roe)*delta_U)

returnflux

#示例數(shù)據(jù)

U_left=np.array([1.0,0.5,2.5])#左側(cè)狀態(tài)向量

U_right=np.array([0.5,0.7,2.0])#右側(cè)狀態(tài)向量

gamma=1.4#比熱比

#計(jì)算Roe平均狀態(tài)

A_roe,u_roe,p_roe,rho_roe=roe_average(U_left,U_right)

#計(jì)算Roe界面通量

flux=roe_flux(U_left,U_right,A_roe,u_roe,p_roe,rho_roe)

print("Roe界面通量:",flux)4.2.3描述在Roe平均與Roe通量的計(jì)算中，我們首先通過roe_average函數(shù)計(jì)算Roe平均狀態(tài)，包括平均速度、壓力和密度，以及Roe平均矩陣。然后，通過roe_flux函數(shù)，使用Roe平均矩陣和左右狀態(tài)的差值，計(jì)算界面通量。這個(gè)示例簡(jiǎn)化了實(shí)際的Roe方法，實(shí)際應(yīng)用中需要更精確的特征分解和通量計(jì)算。以上示例展示了Godunov方法和Roe平均與Roe通量計(jì)算的基本原理和實(shí)現(xiàn)，它們是有限體積法求解二維氣體動(dòng)力學(xué)方程的重要組成部分。在實(shí)際應(yīng)用中，這些方法需要與時(shí)間離散化和空間重構(gòu)技術(shù)結(jié)合，以獲得高精度和穩(wěn)定的數(shù)值解。5邊界條件處理5.1物理邊界條件的實(shí)施在空氣動(dòng)力學(xué)數(shù)值模擬中，物理邊界條件的實(shí)施是確保模擬結(jié)果準(zhǔn)確性和物理意義的關(guān)鍵步驟。邊界條件反映了流體與邊界之間的相互作用，包括但不限于壁面、進(jìn)氣口、排氣口等。在有限體積法(FVM)中，邊界條件的處理通常涉及對(duì)邊界面上的通量進(jìn)行計(jì)算，以反映流體在邊界上的行為。5.1.1壁面邊界條件壁面邊界條件通常假設(shè)流體在壁面上的速度為零（無滑移條件），并且壁面是絕熱的（無熱流）。在FVM中，這意味著在壁面邊界上，通量的法向分量為零。5.1.2進(jìn)氣口邊界條件進(jìn)氣口邊界條件通常設(shè)定為給定的流體速度和溫度。在FVM中，這可以通過在進(jìn)氣口邊界上直接指定速度和溫度的值來實(shí)現(xiàn)。5.1.3排氣口邊界條件排氣口邊界條件通常設(shè)定為自由流出，即流體可以自由地離開計(jì)算域，而不受任何額外的力的影響。在FVM中，這通常意味著在排氣口邊界上，壓力被設(shè)定為大氣壓力，而速度則由內(nèi)部流場(chǎng)決定。5.2數(shù)值邊界條件的處理數(shù)值邊界條件的處理涉及到如何在計(jì)算網(wǎng)格的邊界上應(yīng)用物理邊界條件。在FVM中，邊界條件的數(shù)值處理通常包括對(duì)邊界面上的通量進(jìn)行計(jì)算，以及如何處理邊界上的離散方程。5.2.1通量計(jì)算在FVM中，邊界面上的通量計(jì)算是基于流體動(dòng)力學(xué)方程的。例如，對(duì)于二維Euler方程，邊界面上的通量可以表示為：#假設(shè)我們有速度u,v和壓力p

#以及邊界法向量n_x和n_y

u,v,p=...#這些是速度和壓力的值

n_x,n_y=...#這些是邊界法向量的分量

#計(jì)算邊界面上的通量

flux_x=u*p*n_x+0.5*(u**2+v**2)*n_x

flux_y=v*p*n_y+0.5*(u**2+v**2)*n_y5.2.2離散方程處理在邊界上，離散方程的處理需要特別注意，以確保數(shù)值穩(wěn)定性。例如，對(duì)于壁面邊界，由于速度為零，我們可能需要修改離散方程，以避免除以零的錯(cuò)誤。#假設(shè)我們有邊界上的離散方程

#以及邊界條件（例如，壁面邊界條件）

discrete_eq=...#這是離散方程

boundary_condition=...#這是邊界條件

#處理邊界上的離散方程

ifboundary_condition=='wall':

#修改離散方程以反映無滑移條件

discrete_eq=discrete_eq.subs(u,0).subs(v,0)5.2.3數(shù)據(jù)樣例考慮一個(gè)二維計(jì)算網(wǎng)格，其中包含一個(gè)壁面邊界。我們可以通過以下方式處理壁面邊界上的離散方程：#假設(shè)我們有網(wǎng)格數(shù)據(jù)和邊界信息

grid_data=...#這是網(wǎng)格數(shù)據(jù)

boundary_info=...#這是邊界信息

#遍歷邊界信息，處理壁面邊界

forboundaryinboundary_info:

ifboundary['type']=='wall':

#獲取邊界上的網(wǎng)格點(diǎn)

boundary_points=boundary['points']

#遍歷邊界上的網(wǎng)格點(diǎn)，修改離散方程

forpointinboundary_points:

#獲取網(wǎng)格點(diǎn)上的速度和壓力

u,v,p=grid_data[point]['u'],grid_data[point]['v'],grid_data[point]['p']

#修改離散方程以反映壁面邊界條件

grid_data[point]['discrete_eq']=grid_data[point]['discrete_eq'].subs(u,0).subs(v,0)通過上述方法，我們可以確保在壁面邊界上，流體的速度為零，從而滿足無滑移條件。同樣，對(duì)于進(jìn)氣口和排氣口邊界，我們可以通過直接指定速度和壓力的值，或者設(shè)定為自由流出條件，來處理邊界上的離散方程。5.2.4結(jié)論邊界條件的處理在空氣動(dòng)力學(xué)數(shù)值模擬中至關(guān)重要。通過正確實(shí)施物理邊界條件，并在數(shù)值上處理邊界上的離散方程，我們可以確保模擬結(jié)果的準(zhǔn)確性和物理意義。在FVM中，邊界面上的通量計(jì)算和邊界上的離散方程處理是實(shí)現(xiàn)這一目標(biāo)的關(guān)鍵步驟。6時(shí)間推進(jìn)方法在空氣動(dòng)力學(xué)數(shù)值模擬中，時(shí)間推進(jìn)方法是解決瞬態(tài)問題的關(guān)鍵。有限體積法(FVM)在處理二維氣體動(dòng)力學(xué)方程時(shí)，通常需要選擇合適的時(shí)間推進(jìn)策略來確保解的穩(wěn)定性和準(zhǔn)確性。本教程將詳細(xì)介紹兩種主要的時(shí)間推進(jìn)方法：顯式時(shí)間推進(jìn)和隱式時(shí)間推進(jìn)。6.1顯式時(shí)間推進(jìn)6.1.1原理顯式時(shí)間推進(jìn)方法基于當(dāng)前時(shí)間步的解直接計(jì)算下一個(gè)時(shí)間步的解，無需求解線性方程組。這種方法簡(jiǎn)單直觀，但其穩(wěn)定性受到時(shí)間步長(zhǎng)的限制，通常需要滿足CFL條件（Courant-Friedrichs-Lewy條件）。6.1.2內(nèi)容在二維氣體動(dòng)力學(xué)方程的FVM求解中，顯式時(shí)間推進(jìn)可以表示為：u其中，un是當(dāng)前時(shí)間步的解，un+1是下一個(gè)時(shí)間步的解，6.1.3示例假設(shè)我們正在解決二維Euler方程，使用顯式時(shí)間推進(jìn)方法。以下是一個(gè)Python代碼示例，展示如何使用顯式時(shí)間推進(jìn)更新網(wǎng)格上的解：importnumpyasnp

defexplicit_time_advance(u,flux,dt,dx,dy):

"""

顯式時(shí)間推進(jìn)更新解

:paramu:當(dāng)前時(shí)間步的解，形狀為(Nx,Ny,4)，其中4表示密度、動(dòng)量、能量等

:paramflux:通量函數(shù)，計(jì)算基于當(dāng)前解的通量

:paramdt:時(shí)間步長(zhǎng)

:paramdx:x方向的網(wǎng)格間距

:paramdy:y方向的網(wǎng)格間距

:return:更新后的時(shí)間步解

"""

#計(jì)算通量

flux_x=flux(u,axis=0)

flux_y=flux(u,axis=1)

#更新解

u_new=u-dt*(np.diff(flux_x,axis=0)/dx+np.diff(flux_y,axis=1)/dy)

#邊界條件處理

u_new[0,:]=u[0,:]#左邊界

u_new[-1,:]=u[-1,:]#右邊界

u_new[:,0]=u[:,0]#下邊界

u_new[:,-1]=u[:,-1]#上邊界

returnu_new

#假設(shè)的通量函數(shù)

defflux(u,axis):

"""

計(jì)算通量

:paramu:解

:paramaxis:計(jì)算通量的方向

:return:通量

"""

#這里使用一個(gè)簡(jiǎn)單的示例通量函數(shù)，實(shí)際應(yīng)用中需要根據(jù)Euler方程計(jì)算通量

returnu*0.5#簡(jiǎn)化示例

#初始化解

Nx,Ny=100,100

u=np.zeros((Nx,Ny,4))

#設(shè)置時(shí)間步長(zhǎng)和網(wǎng)格間距

dt=0.01

dx=0.1

dy=0.1

#更新解

u=explicit_time_advance(u,flux,dt,dx,dy)6.2隱式時(shí)間推進(jìn)6.2.1原理隱式時(shí)間推進(jìn)方法在計(jì)算下一個(gè)時(shí)間步的解時(shí)，考慮了未來狀態(tài)的影響。這種方法通常需要求解線性方程組，但可以使用較大的時(shí)間步長(zhǎng)，從而提高計(jì)算效率。6.2.2內(nèi)容隱式時(shí)間推進(jìn)可以表示為：I其中，I是單位矩陣，A是包含空間導(dǎo)數(shù)的矩陣，Δt6.2.3示例使用隱式時(shí)間推進(jìn)更新解通常涉及求解線性方程組。以下是一個(gè)使用Python和SciPy庫的示例，展示如何使用隱式時(shí)間推進(jìn)更新網(wǎng)格上的解：importnumpyasnp

fromscipy.sparseimportdiags

fromscipy.sparse.linalgimportspsolve

defimplicit_time_advance(u,flux,dt,dx,dy):

"""

隱式時(shí)間推進(jìn)更新解

:paramu:當(dāng)前時(shí)間步的解，形狀為(Nx,Ny,4)

:paramflux:通量函數(shù)

:paramdt:時(shí)間步長(zhǎng)

:paramdx:x方向的網(wǎng)格間距

:paramdy:y方向的網(wǎng)格間距

:return:更新后的時(shí)間步解

"""

#計(jì)算通量

flux_x=flux(u,axis=0)

flux_y=flux(u,axis=1)

#構(gòu)建矩陣A

A_x=diags([-1,1],[-1,1],shape=(u.shape[0],u.shape[0]))/dx

A_y=diags([-1,1],[-1,1],shape=(u.shape[1],u.shape[1]))/dy

A=np.kron(np.eye(u.shape[1]),A_x)+np.kron(A_y,np.eye(u.shape[0]))

#構(gòu)建右側(cè)向量

b=u.flatten()+dt*(np.diff(flux_x,axis=0)/dx+np.diff(flux_y,axis=1)/dy).flatten()

#求解線性方程組

u_new_flat=spsolve(I-dt*A,b)

#重塑解

u_new=u_new_flat.reshape(u.shape)

#邊界條件處理

u_new[0,:]=u[0,:]#左邊界

u_new[-1,:]=u[-1,:]#右邊界

u_new[:,0]=u[:,0]#下邊界

u_new[:,-1]=u[:,-1]#上邊界

returnu_new

#假設(shè)的通量函數(shù)

defflux(u,axis):

"""

計(jì)算通量

:paramu:解

:paramaxis:計(jì)算通量的方向

:return:通量

"""

#這里使用一個(gè)簡(jiǎn)單的示例通量函數(shù)，實(shí)際應(yīng)用中需要根據(jù)Euler方程計(jì)算通量

returnu*0.5#簡(jiǎn)化示例

#初始化解

Nx,Ny=100,100

u=np.zeros((Nx,Ny,4))

#設(shè)置時(shí)間步長(zhǎng)和網(wǎng)格間距

dt=0.1#注意這里的時(shí)間步長(zhǎng)比顯式方法大

dx=0.1

dy=0.1

#更新解

u=implicit_time_advance(u,flux,dt,dx,dy)在上述示例中，我們使用了SciPy庫中的spsolve函數(shù)來求解線性方程組，其中A是包含空間導(dǎo)數(shù)的矩陣，b是基于當(dāng)前解和通量計(jì)算的右側(cè)向量。通過隱式時(shí)間推進(jìn)，我們可以使用較大的時(shí)間步長(zhǎng)，從而減少計(jì)算迭代次數(shù)，提高效率。7高精度格式與限幅器7.1高精度重構(gòu)技術(shù)在空氣動(dòng)力學(xué)數(shù)值模擬中，有限體積法(FVM)是一種廣泛使用的方法，它基于守恒定律，通過在網(wǎng)格上離散化連續(xù)方程來求解流體動(dòng)力學(xué)問題。為了提高FVM的精度，特別是在處理復(fù)雜的流場(chǎng)結(jié)構(gòu)時(shí)，高精度重構(gòu)技術(shù)變得至關(guān)重要。這些技術(shù)通過在每個(gè)網(wǎng)格單元內(nèi)進(jìn)行更高階的插值，來更準(zhǔn)確地估計(jì)界面通量，從而減少數(shù)值擴(kuò)散和振蕩。7.1.1原理高精度重構(gòu)技術(shù)的核心在于如何在網(wǎng)格單元內(nèi)部進(jìn)行插值，以獲得更精確的界面狀態(tài)。常見的高精度重構(gòu)方法包括：線性重構(gòu)：假設(shè)網(wǎng)格單元內(nèi)的狀態(tài)是線性變化的，使用網(wǎng)格單元中心的值和其梯度來估計(jì)界面狀態(tài)。高階重構(gòu)：使用多項(xiàng)式擬合或更復(fù)雜的函數(shù)來描述網(wǎng)格單元內(nèi)的狀態(tài)變化，以提高精度。7.1.2內(nèi)容在二維氣體動(dòng)力學(xué)方程的FVM求解中，高精度重構(gòu)技術(shù)通常應(yīng)用于狀態(tài)變量（如密度、速度、壓力等）的界面值估計(jì)。例如，使用二階重構(gòu)，我們可以基于網(wǎng)格單元中心的值和其梯度，來估計(jì)界面處的狀態(tài)變量值。7.1.2.1示例假設(shè)我們有一個(gè)二維網(wǎng)格，每個(gè)網(wǎng)格單元內(nèi)部的狀態(tài)變量（如密度）可以表示為：ρ其中，ρc是網(wǎng)格單元中心的密度值，?ρ?x和在界面x=xiρ7.1.3代碼示例以下是一個(gè)使用Python實(shí)現(xiàn)的簡(jiǎn)單二階重構(gòu)示例，用于估計(jì)二維網(wǎng)格單元界面處的密度值：importnumpyasnp

deflinear_reconstruction(rho_c,grad_rho_x,grad_rho_y,x_i,y_i,x_c,y_c):

"""

使用線性重構(gòu)技術(shù)估計(jì)界面處的密度值。

參數(shù):

rho_c:網(wǎng)格單元中心的密度值。

grad_rho_x:密度在x方向上的梯度。

grad_rho_y:密度在y方向上的梯度。

x_i:界面在x方向的位置。

y_i:界面在y方向的位置。

x_c:網(wǎng)格單元中心在x方向的位置。

y_c:網(wǎng)格單元中心在y方向的位置。

返回:

rho_i:界面處的密度值。

"""

rho_i=rho_c+grad_rho_x*(x_i-x_c)+grad_rho_y*(y_i-y_c)

returnrho_i

#示例數(shù)據(jù)

rho_c=1.225#中心密度值

grad_rho_x=-0.01#x方向的密度梯度

grad_rho_y=0.02#y方向的密度梯度

x_i=0.5#界面x位置

y_i=0.5#界面y位置

x_c=0.0#中心x位置

y_c=0.0#中心y位置

#計(jì)算界面處的密度值

rho_i=linear_reconstruction(rho_c,grad_rho_x,grad_rho_y,x_i,y_i,x_c,y_c)

print(f"界面處的密度值:{rho_i}")7.2通量限幅器的使用在使用高精度重構(gòu)技術(shù)時(shí)，可能會(huì)遇到數(shù)值振蕩的問題，特別是在激波或其它強(qiáng)不連續(xù)性附近。為了解決這一問題，通量限幅器被引入，它是一種非線性限制器，用于控制高階重構(gòu)方案的振蕩，確保數(shù)值解的穩(wěn)定性。7.2.1原理通量限幅器通過調(diào)整界面通量的計(jì)算，來限制數(shù)值解的振蕩。它基于局部梯度信息，動(dòng)態(tài)調(diào)整重構(gòu)方案的斜率，以確保數(shù)值解不會(huì)產(chǎn)生不物理的振蕩。常見的限幅器包括：VanLeer限幅器Superbee限幅器Minmod限幅器7.2.2內(nèi)容在二維氣體動(dòng)力學(xué)方程的FVM求解中，通量限幅器通常應(yīng)用于計(jì)算界面通量的階段。例如，使用VanLeer限幅器，我們可以調(diào)整界面處的密度梯度，以減少振蕩。7.2.2.1示例假設(shè)我們使用VanLeer限幅器來調(diào)整界面處的密度梯度。VanLeer限幅器的公式為：?其中，r是上下游網(wǎng)格單元中心密度梯度的比值。7.2.3代碼示例以下是一個(gè)使用Python實(shí)現(xiàn)的VanLeer限幅器示例，用于調(diào)整界面處的密度梯度：defvan_leer_limiter(r):

"""

使用VanLeer限幅器調(diào)整密度梯度。

參數(shù):

r:上下游網(wǎng)格單元中心密度梯度的比值。

返回:

phi:限幅后的梯度調(diào)整因子。

"""

phi=(r+abs(r))/(1+abs(r))

returnphi

#示例數(shù)據(jù)

grad_rho_x_left=-0.02#左側(cè)網(wǎng)格單元中心的x方向密度梯度

grad_rho_x_right=-0.01#右側(cè)網(wǎng)格單元中心的x方向密度梯度

r=grad_rho_x_right/grad_rho_x_left#比值

#使用VanLeer限幅器調(diào)整梯度

phi=van_leer_limiter(r)

print(f"限幅后的梯度調(diào)整因子:{phi}")通過上述代碼示例，我們可以看到如何在二維氣體動(dòng)力學(xué)方程的FVM求解中，應(yīng)用高精度重構(gòu)技術(shù)和通量限幅器來提高數(shù)值解的精度和穩(wěn)定性。這些技術(shù)是空氣動(dòng)力學(xué)數(shù)值模擬中不可或缺的組成部分，能夠幫助我們更準(zhǔn)確地理解和預(yù)測(cè)流體動(dòng)力學(xué)現(xiàn)象。8穩(wěn)定性與收斂性8.1穩(wěn)定性分析穩(wěn)定性是數(shù)值方法中一個(gè)關(guān)鍵的概念，它確保了計(jì)算過程不會(huì)因?yàn)槲⑿〉臄_動(dòng)而產(chǎn)生巨大的誤差。在有限體積法(FVM)中，穩(wěn)定性分析通常涉及到對(duì)時(shí)間步長(zhǎng)和網(wǎng)格尺寸的限制，以確保數(shù)值解的正確性和可靠性。8.1.1CFL條件CFL條件（Courant-Friedrichs-Lewy條件）是判斷穩(wěn)定性的一個(gè)重要準(zhǔn)則。它基于波在時(shí)間步長(zhǎng)內(nèi)不能超過一個(gè)網(wǎng)格單元的距離這一物理直覺。CFL數(shù)定義為：C其中，u是流體的速度，Δt是時(shí)間步長(zhǎng)，Δ8.1.2穩(wěn)定性判據(jù)除了CFL條件，還有其他穩(wěn)定性判據(jù)，如矩陣穩(wěn)定性分析和頻譜分析。這些方法通常用于更復(fù)雜的情況，如非線性方程或不規(guī)則網(wǎng)格。8.2收斂性判斷與加速收斂性是指隨著網(wǎng)格細(xì)化和時(shí)間步長(zhǎng)減小，數(shù)值解逐漸接近真實(shí)解的性質(zhì)。判斷和加速收斂性是有限體積法中優(yōu)化計(jì)算效率的重要步驟。8.2.1收斂性判斷收斂性可以通過監(jiān)測(cè)殘差（即迭代過程中方程的不平衡量）來判斷。當(dāng)殘差低于一個(gè)預(yù)設(shè)的閾值時(shí)，可以認(rèn)為解已經(jīng)收斂。8.2.2收斂加速技術(shù)為了加速收斂，可以采用多種技術(shù)，如多網(wǎng)格方法、預(yù)條件技術(shù)、以及時(shí)間步長(zhǎng)的自適應(yīng)調(diào)整。其中，多網(wǎng)格方法通過在不同尺度的網(wǎng)格上迭代求解，可以有效減少高頻率誤差，從而加速收斂。8.2.3示例：CFL條件檢查假設(shè)我們有一個(gè)簡(jiǎn)單的二維氣體動(dòng)力學(xué)問題，其中流體的速度為u=1，網(wǎng)格間距Δx=Δ#定義流體速度和網(wǎng)格參數(shù)

u=1.0

dx=0.1

dy=0.1

CFL=0.5

#計(jì)算時(shí)間步長(zhǎng)

dt=CFL*min(dx,dy)/u

#輸出時(shí)間步長(zhǎng)

print(f"時(shí)間步長(zhǎng):{dt}")8.2.4示例：殘差監(jiān)測(cè)在迭代求解過程中，監(jiān)測(cè)殘差是判斷收斂性的重要手段。以下是一個(gè)使用Python和NumPy庫監(jiān)測(cè)二維氣體動(dòng)力學(xué)方程殘差的示例：importnumpyasnp

#定義殘差閾值

residual_threshold=1e-6

#初始化殘差

residual=1.0

#迭代次數(shù)

iteration=0

#迭代求解過程

whileresidual>residual_threshold:

#假設(shè)這里是求解方程的代碼

#...

#計(jì)算殘差

residual=np.linalg.norm(solution-previous_solution)

#更新迭代次數(shù)

iteration+=1

#輸出當(dāng)前迭代次數(shù)和殘差

print(f"迭代次數(shù):{iteration},殘差:{residual}")

#保存當(dāng)前解作為下一次迭代的前一次解

previous_solution=solution.copy()通過這些示例和理論介紹，我們能夠更好地理解和應(yīng)用有限體積法中的穩(wěn)定性與收斂性概念。9案例分析9.1維激波管問題9.1.1問題描述二維激波管問題是一個(gè)經(jīng)典的測(cè)試案例，用于驗(yàn)證有限體積法(FVM)在解決二維氣體動(dòng)力學(xué)方程時(shí)的準(zhǔn)確性和穩(wěn)定性。該問題模擬了在二維空間中，由一個(gè)初始條件突變引起的激波和膨脹波的傳播。激波管通常被分為兩個(gè)區(qū)域，每個(gè)區(qū)域具有不同的初始狀態(tài)，如壓力、密度和速度。9.1.2數(shù)學(xué)模型二維氣體動(dòng)力學(xué)方程包括連續(xù)性方程、動(dòng)量方程和能量方程，可以表示為：?其中，U是保守變量向量，F(xiàn)和G是沿x和y方向的通量向量。9.1.3離散化在有限體積法中，我們將計(jì)算域劃分為一系列控制體積，然后在每個(gè)控制體積上應(yīng)用積分形式的氣體動(dòng)力學(xué)方程。對(duì)于二維激波管問題，我們使用矩形網(wǎng)格，每個(gè)網(wǎng)格單元代表一個(gè)控制體積。9.1.4算法實(shí)現(xiàn)下面是一個(gè)使用Python實(shí)現(xiàn)的二維激波管問題的有限體積法求解器的示例代碼：importnumpyasnp

importmatplotlib.pyplotasplt

#參數(shù)設(shè)置

nx=100#x方向網(wǎng)格數(shù)

ny=100#y方向網(wǎng)格數(shù)

nt=100#時(shí)間步數(shù)

dx=1.0/(nx-1)#x方向網(wǎng)格間距

dy=1.0/(ny-1)#y方向網(wǎng)格間距

c=1.0#聲速

gamma=1.4#比熱比

#初始條件

rho_L=1.0#左側(cè)密度

rho_R=0.125#右側(cè)密度

p_L=1.0#左側(cè)壓力

p_R=0.1#右側(cè)壓力

u_L=0.0#左側(cè)x方向速度

u_R=0.0#右側(cè)x方向速度

v_L=0.0#左側(cè)y方向速度

v_R=0.0#右側(cè)y方向速度

#定義網(wǎng)格

x=np.linspace(0,1,nx)

y=np.linspace(0,1,ny)

X,Y=np.meshgrid(x,y)

#初始化狀態(tài)

rho=np.ones((nx,ny))*rho_R

p=np.ones((nx,ny))*p_R

u=np.ones((nx,ny))*u_R

v=np.ones((nx,ny))*v_R

#設(shè)置左側(cè)初始條件

rho[:,:nx//2]=rho_L

p[:,:nx//2]=p_L

u[:,:nx//2]=u_L

v[:,:nx//2]=v_L

#定義狀態(tài)向量和通量向量

defconservative_variables(rho,u,v,p):

e=p/((gamma-1)*rho)+0.5*(u**2+v**2)*rho

returnnp.array([rho,rho*u,rho*v,e])

deffluxes(F,G,u,v,p,rho):

a=np.sqrt(gamma*p/rho)

F[1]+=u*F[0]

F[2]+=v*F[0]

F[3]+=(u*F[1]+v*F[2]+p*u)

G[1]+=u*G[0]

G[2]+=v*G[0]

G[3]+=(u*G[1]+v*G[2]+p*v)

returnF,G

#主循環(huán)

forninrange(nt):

U=conservative_variables(rho,u,v,p)

F=np.zeros_like(U)

G=np.zeros_like(U)

F,G=fluxes(F,G,u,v,p,rho)

#更新狀態(tài)

rho-=(F[0,:,1:]-F[0,:,:-1])*dt/dx

u-=(F[1,:,1:]-F[1,:,:-1])*dt/dx/rho

v-=(G[2,1:,:]-G[2,:-1,:])*dt/dy/rho

p-=(F[3,:,1:]-F[3,:,:-1])*dt/dx-(G[3,1:,:]-G[3,:-1,:])*