空氣動(dòng)力學(xué)數(shù)值方法：有限差分法(FDM)：二維穩(wěn)態(tài)對(duì)流方程的有限差分解法

空氣動(dòng)力學(xué)數(shù)值方法：有限差分法(FDM)：二維穩(wěn)態(tài)對(duì)流方程的有限差分解法1空氣動(dòng)力學(xué)數(shù)值方法：有限差分法(FDM)：二維穩(wěn)態(tài)對(duì)流方程的有限差分解法1.1緒論1.1.1有限差分法的基本概念有限差分法(FiniteDifferenceMethod,FDM)是一種廣泛應(yīng)用于偏微分方程數(shù)值求解的方法，尤其在空氣動(dòng)力學(xué)領(lǐng)域中，用于模擬流體流動(dòng)、熱傳導(dǎo)等物理現(xiàn)象。FDM的基本思想是將連續(xù)的物理空間離散化，即將連續(xù)的偏微分方程在離散的網(wǎng)格點(diǎn)上用差分近似代替導(dǎo)數(shù)，從而將偏微分方程轉(zhuǎn)化為代數(shù)方程組，通過(guò)求解這些方程組來(lái)得到原問(wèn)題的近似解。1.1.1.1示例：一維熱傳導(dǎo)方程的有限差分解考慮一維熱傳導(dǎo)方程：?其中，u是溫度，α是熱擴(kuò)散率。在有限差分法中，我們首先定義一個(gè)離散網(wǎng)格，假設(shè)網(wǎng)格間距為Δx，時(shí)間步長(zhǎng)為Δt。在網(wǎng)格點(diǎn)xi和時(shí)間tu這里，uin表示在網(wǎng)格點(diǎn)xi1.1.2對(duì)流方程的物理意義對(duì)流方程描述了流體中物質(zhì)或能量的對(duì)流過(guò)程，即物質(zhì)或能量隨流體的運(yùn)動(dòng)而轉(zhuǎn)移。在空氣動(dòng)力學(xué)中，對(duì)流方程常用于描述空氣或氣體的流動(dòng)，特別是當(dāng)流體的速度遠(yuǎn)大于擴(kuò)散速度時(shí)，對(duì)流效應(yīng)更為顯著。1.1.2.1維穩(wěn)態(tài)對(duì)流方程二維穩(wěn)態(tài)對(duì)流方程可以表示為：u其中，u和v分別是流體在x和y方向的速度分量，?是隨流體一起移動(dòng)的標(biāo)量量，如溫度或濃度。這個(gè)方程表明，在穩(wěn)態(tài)條件下，?的梯度與流體速度的方向和大小相乘，其結(jié)果為零，意味著?在流體中的分布不會(huì)隨時(shí)間改變。1.1.2.2有限差分解在二維穩(wěn)態(tài)對(duì)流方程的有限差分解中，我們同樣需要定義一個(gè)離散網(wǎng)格，假設(shè)網(wǎng)格在x和y方向的間距分別為Δx和Δy。在網(wǎng)格點(diǎn)xiu這里，?i,j表示在網(wǎng)格點(diǎn)xi,1.1.3代碼示例：二維穩(wěn)態(tài)對(duì)流方程的有限差分解importnumpyasnp

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

nx,ny=100,100

dx,dy=1.0/(nx-1),1.0/(ny-1)

u,v=1.0,1.0#流體速度

#初始化\phi值

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

#邊界條件

phi[0,:]=1.0#左邊界\phi=1

phi[-1,:]=0.0#右邊界\phi=0

#內(nèi)部點(diǎn)的有限差分解

foriinrange(1,nx-1):

forjinrange(1,ny-1):

phi[i,j]=(u*(dx/2)*(phi[i+1,j]-phi[i-1,j])+

v*(dy/2)*(phi[i,j+1]-phi[i,j-1]))/(u*dx/2+v*dy/2)

#輸出結(jié)果

print(phi)這段代碼展示了如何使用有限差分法求解二維穩(wěn)態(tài)對(duì)流方程。首先，定義了網(wǎng)格參數(shù)和流體速度，然后初始化了?值并設(shè)置了邊界條件。通過(guò)迭代求解內(nèi)部點(diǎn)的差分方程，最終得到了?在二維空間中的分布。1.1.4結(jié)論有限差分法是空氣動(dòng)力學(xué)數(shù)值模擬中的一種重要工具，它通過(guò)將連續(xù)的物理空間離散化，將偏微分方程轉(zhuǎn)化為代數(shù)方程組，從而實(shí)現(xiàn)對(duì)復(fù)雜物理現(xiàn)象的數(shù)值求解。對(duì)流方程的有限差分解是理解流體流動(dòng)和物質(zhì)傳輸?shù)年P(guān)鍵，通過(guò)適當(dāng)?shù)木W(wǎng)格劃分和差分格式選擇，可以有效地模擬空氣動(dòng)力學(xué)中的對(duì)流過(guò)程。2空氣動(dòng)力學(xué)數(shù)值方法：有限差分法(FDM)：二維穩(wěn)態(tài)對(duì)流方程的有限差分解法2.1有限差分法的數(shù)學(xué)基礎(chǔ)2.1.1泰勒級(jí)數(shù)展開(kāi)泰勒級(jí)數(shù)展開(kāi)是有限差分法的核心數(shù)學(xué)工具，它允許我們將連續(xù)函數(shù)在某一點(diǎn)的值用該點(diǎn)及其鄰域的導(dǎo)數(shù)的無(wú)窮級(jí)數(shù)來(lái)表示。對(duì)于一個(gè)在點(diǎn)x處的函數(shù)fxf其中，h是離散步長(zhǎng)，f′x、f″x、2.1.1.1示例假設(shè)我們有一個(gè)函數(shù)fx=ex，我們想要在importmath

#定義函數(shù)f(x)=e^x

deff(x):

returnmath.exp(x)

#定義泰勒級(jí)數(shù)展開(kāi)的前幾項(xiàng)

deftaylor_expansion(x,h,n_terms=3):

result=f(x)

forninrange(1,n_terms):

result+=(h**n)/math.factorial(n)*f(x)

returnresult

#在x=0處計(jì)算f(0.1)的泰勒級(jí)數(shù)近似值

x=0

h=0.1

approx_value=taylor_expansion(x,h)

#輸出結(jié)果

print("f(0.1)的泰勒級(jí)數(shù)近似值為:",approx_value)

print("f(0.1)的實(shí)際值為:",f(0.1))2.1.2差分逼近的構(gòu)造差分逼近是通過(guò)在離散點(diǎn)上使用函數(shù)值來(lái)估計(jì)導(dǎo)數(shù)的方法。常見(jiàn)的差分逼近有向前差分、向后差分和中心差分。向前差分：f向后差分：f中心差分：f2.1.2.1示例假設(shè)我們有函數(shù)fx=x#定義函數(shù)f(x)=x^2

deff(x):

returnx**2

#定義向前差分逼近

defforward_difference(x,h):

return(f(x+h)-f(x))/h

#定義向后差分逼近

defbackward_difference(x,h):

return(f(x)-f(x-h))/h

#定義中心差分逼近

defcentral_difference(x,h):

return(f(x+h)-f(x-h))/(2*h)

#在x=1處計(jì)算導(dǎo)數(shù)的差分逼近

x=1

h=0.001

forward_approx=forward_difference(x,h)

backward_approx=backward_difference(x,h)

central_approx=central_difference(x,h)

#輸出結(jié)果

print("向前差分逼近的導(dǎo)數(shù)值為:",forward_approx)

print("向后差分逼近的導(dǎo)數(shù)值為:",backward_approx)

print("中心差分逼近的導(dǎo)數(shù)值為:",central_approx)2.1.3差分格式的穩(wěn)定性分析差分格式的穩(wěn)定性是有限差分法成功的關(guān)鍵。穩(wěn)定性分析通常涉及確定差分格式的穩(wěn)定性條件，即步長(zhǎng)h和時(shí)間步長(zhǎng)Δt2.1.3.1示例考慮一維對(duì)流方程ut+aux=0，其中au其中，uin表示在時(shí)間nΔt和空間位置ih處的uimportnumpyasnp

#定義差分格式的穩(wěn)定性條件

defstability_condition(a,h,dt):

returna*dt/(2*h)

#定義參數(shù)

a=1

h=0.1

dt=0.01

#計(jì)算穩(wěn)定性條件

stability=stability_condition(a,h,dt)

#輸出結(jié)果

print("穩(wěn)定性條件為:",stability)在本例中，我們計(jì)算了穩(wěn)定性條件aΔt2h，如果該值小于等于1，則差分格式是穩(wěn)定的。這為選擇合適的3維穩(wěn)態(tài)對(duì)流方程的離散化3.1方程的有限差分形式在空氣動(dòng)力學(xué)中，二維穩(wěn)態(tài)對(duì)流方程描述了流體在固定坐標(biāo)系中隨時(shí)間不變化的流動(dòng)特性。該方程的一般形式可以表示為：u其中，u和v分別是流體在x和y方向的速度分量，?是流體的物理量（如溫度、濃度等），S是源項(xiàng)。3.1.1離散化過(guò)程為了使用有限差分法求解上述方程，我們首先需要將其離散化。離散化是將連續(xù)的偏微分方程轉(zhuǎn)換為離散的代數(shù)方程組的過(guò)程，以便在計(jì)算機(jī)上進(jìn)行數(shù)值求解。離散化通常涉及在空間上將域劃分為網(wǎng)格，并在網(wǎng)格點(diǎn)上用差商近似偏導(dǎo)數(shù)。3.1.1.1階上風(fēng)差分在一維情況下，對(duì)流項(xiàng)的離散化可以使用一階上風(fēng)差分。在二維中，我們同樣可以應(yīng)用這一技術(shù)。例如，對(duì)于x方向的對(duì)流項(xiàng)u???xuu對(duì)于y方向的對(duì)流項(xiàng)v?vv3.1.2代碼示例假設(shè)我們有一個(gè)二維網(wǎng)格，其中u和v已知，我們想要離散化對(duì)流方程。以下是一個(gè)使用Python的示例代碼：importnumpyasnp

#定義網(wǎng)格尺寸

nx,ny=100,100

dx,dy=1.0,1.0

#初始化速度場(chǎng)和物理量場(chǎng)

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

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

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

#假設(shè)速度場(chǎng)和物理量場(chǎng)的初始值

u[50:75,50:75]=1.0

v[25:50,25:50]=-1.0

phi[30:80,30:80]=100.0

#離散化對(duì)流項(xiàng)

defconvective_flux(u,v,phi,dx,dy):

phi_x=np.zeros_like(phi)

phi_y=np.zeros_like(phi)

#x方向的上風(fēng)差分

foriinrange(1,nx):

forjinrange(ny):

ifu[i,j]>0:

phi_x[i,j]=u[i,j]*(phi[i+1,j]-phi[i,j])/dx

else:

phi_x[i,j]=u[i,j]*(phi[i,j]-phi[i-1,j])/dx

#y方向的上風(fēng)差分

foriinrange(nx):

forjinrange(1,ny):

ifv[i,j]>0:

phi_y[i,j]=v[i,j]*(phi[i,j+1]-phi[i,j])/dy

else:

phi_y[i,j]=v[i,j]*(phi[i,j]-phi[i,j-1])/dy

returnphi_x,phi_y

#計(jì)算對(duì)流通量

phi_x,phi_y=convective_flux(u,v,phi,dx,dy)3.2網(wǎng)格的生成與選擇網(wǎng)格的生成和選擇是有限差分法中的關(guān)鍵步驟。網(wǎng)格的選擇直接影響到數(shù)值解的準(zhǔn)確性和計(jì)算效率。在二維穩(wěn)態(tài)對(duì)流方程的求解中，網(wǎng)格的密度、形狀和分布都至關(guān)重要。3.2.1網(wǎng)格類型常見(jiàn)的網(wǎng)格類型包括：結(jié)構(gòu)網(wǎng)格：網(wǎng)格點(diǎn)按照規(guī)則的模式排列，如矩形網(wǎng)格。非結(jié)構(gòu)網(wǎng)格：網(wǎng)格點(diǎn)的排列沒(méi)有固定模式，適用于復(fù)雜幾何形狀的流場(chǎng)。3.2.2網(wǎng)格生成網(wǎng)格生成可以通過(guò)多種軟件工具完成，如Gmsh、Gridgen等。在Python中，可以使用numpy和matplotlib庫(kù)來(lái)生成和可視化簡(jiǎn)單的結(jié)構(gòu)網(wǎng)格。3.2.2.1代碼示例以下是一個(gè)使用Python生成和可視化二維結(jié)構(gòu)網(wǎng)格的示例：importnumpyasnp

importmatplotlib.pyplotasplt

#定義網(wǎng)格尺寸

nx,ny=100,100

#生成網(wǎng)格

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

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

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

#可視化網(wǎng)格

plt.figure()

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('2DStructuredGrid')

plt.show()3.2.3網(wǎng)格選擇網(wǎng)格的選擇應(yīng)基于流場(chǎng)的特性。在對(duì)流主導(dǎo)的區(qū)域，網(wǎng)格應(yīng)更密集以捕捉對(duì)流邊界層的細(xì)節(jié)。在源項(xiàng)S集中的區(qū)域，網(wǎng)格也應(yīng)更密集以準(zhǔn)確表示源項(xiàng)的影響。3.2.3.1網(wǎng)格適應(yīng)性網(wǎng)格適應(yīng)性（gridadaptivity）是一種技術(shù)，它允許在計(jì)算過(guò)程中動(dòng)態(tài)調(diào)整網(wǎng)格的密度，以提高計(jì)算效率和準(zhǔn)確性。這通常通過(guò)監(jiān)測(cè)解的梯度或誤差來(lái)實(shí)現(xiàn)，從而在需要更高分辨率的區(qū)域增加網(wǎng)格點(diǎn)。3.2.3.2代碼示例網(wǎng)格適應(yīng)性通常涉及到復(fù)雜的算法，這里提供一個(gè)簡(jiǎn)化的示例，說(shuō)明如何根據(jù)物理量?的梯度調(diào)整網(wǎng)格密度：importnumpyasnp

#定義初始網(wǎng)格尺寸

nx,ny=100,100

#初始化物理量場(chǎng)

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

#假設(shè)物理量場(chǎng)的初始值

phi[30:80,30:80]=100.0

#計(jì)算物理量場(chǎng)的梯度

phi_grad_x=np.gradient(phi,axis=0)

phi_grad_y=np.gradient(phi,axis=1)

#根據(jù)梯度調(diào)整網(wǎng)格密度

dx_adaptive=1.0/(np.abs(phi_grad_x)+1)

dy_adaptive=1.0/(np.abs(phi_grad_y)+1)

#注意：實(shí)際應(yīng)用中，dx_adaptive和dy_adaptive的值將用于重新生成網(wǎng)格在實(shí)際應(yīng)用中，dx_adaptive和dy_adaptive的值將用于重新生成網(wǎng)格，以適應(yīng)物理量場(chǎng)的局部變化。然而，上述代碼僅用于說(shuō)明目的，實(shí)際的網(wǎng)格適應(yīng)性算法會(huì)更復(fù)雜，可能需要重新生成網(wǎng)格并重新計(jì)算所有相關(guān)的物理量和通量。通過(guò)上述原理和代碼示例，我們可以看到二維穩(wěn)態(tài)對(duì)流方程的有限差分解法涉及對(duì)流項(xiàng)的離散化和網(wǎng)格的生成與選擇。這些步驟對(duì)于準(zhǔn)確求解空氣動(dòng)力學(xué)問(wèn)題至關(guān)重要。4邊界條件的處理在空氣動(dòng)力學(xué)數(shù)值模擬中，邊界條件的正確處理對(duì)于獲得準(zhǔn)確的解至關(guān)重要。邊界條件反映了流體與邊界之間的相互作用，包括周期性邊界條件、壁面邊界條件和遠(yuǎn)場(chǎng)邊界條件。下面將詳細(xì)討論這三種邊界條件的原理和處理方法。4.1周期性邊界條件周期性邊界條件通常用于模擬具有周期性特征的流場(chǎng)，如湍流通道或管道流動(dòng)。在二維穩(wěn)態(tài)對(duì)流方程的有限差分解法中，周期性邊界條件意味著流場(chǎng)在邊界上的物理量（如速度、壓力）在周期性邊界上是相等的。4.1.1實(shí)現(xiàn)示例假設(shè)我們有一個(gè)二維流場(chǎng)，其中x方向和y方向的速度分別為u和v，壓力為p。在周期性邊界上，我們可以通過(guò)以下方式更新邊界上的速度和壓力：#周期性邊界條件處理

defapply_periodic_boundary_conditions(u,v,p):

"""

應(yīng)用周期性邊界條件

:paramu:二維數(shù)組，x方向速度

:paramv:二維數(shù)組，y方向速度

:paramp:二維數(shù)組，壓力

"""

#假設(shè)邊界在x=0和x=Lx，y=0和y=Ly

Lx,Ly=u.shape

#更新x方向的周期性邊界

u[:,0]=u[:,Lx-1]

u[:,Ly]=u[:,1]

#更新y方向的周期性邊界

v[0,:]=v[Lx-1,:]

v[Ly,:]=v[1,:]

#更新壓力的周期性邊界

p[:,0]=p[:,Lx-1]

p[:,Ly]=p[:,1]

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

u=np.zeros((10,10))

v=np.zeros((10,10))

p=np.zeros((10,10))

#應(yīng)用周期性邊界條件

apply_periodic_boundary_conditions(u,v,p)4.2壁面邊界條件壁面邊界條件通常用于模擬流體與固體壁面的相互作用。在壁面上，流體的速度通常為零（無(wú)滑移條件），并且流體與壁面之間的剪應(yīng)力與速度梯度成正比（牛頓內(nèi)摩擦定律）。4.2.1實(shí)現(xiàn)示例在有限差分解法中，壁面邊界條件可以通過(guò)設(shè)置邊界上的速度為零來(lái)實(shí)現(xiàn)。假設(shè)我們有一個(gè)二維流場(chǎng)，其中x方向和y方向的速度分別為u和v，并且壁面位于y=#壁面邊界條件處理

defapply_wall_boundary_condition(u,v):

"""

應(yīng)用壁面邊界條件

:paramu:二維數(shù)組，x方向速度

:paramv:二維數(shù)組，y方向速度

"""

#假設(shè)壁面在y=0的位置

v[:,0]=0.0#y方向速度為零

u[:,0]=u[:,1]#x方向速度使用中心差分法更新

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

u=np.zeros((10,10))

v=np.zeros((10,10))

#應(yīng)用壁面邊界條件

apply_wall_boundary_condition(u,v)4.3遠(yuǎn)場(chǎng)邊界條件遠(yuǎn)場(chǎng)邊界條件用于模擬遠(yuǎn)離流體區(qū)域的邊界，通常假設(shè)流體在這些邊界上不受內(nèi)部流動(dòng)的影響。在二維穩(wěn)態(tài)對(duì)流方程的有限差分解法中，遠(yuǎn)場(chǎng)邊界條件可以通過(guò)設(shè)置邊界上的速度和壓力為自由流條件來(lái)實(shí)現(xiàn)。4.3.1實(shí)現(xiàn)示例假設(shè)我們有一個(gè)二維流場(chǎng)，其中x方向和y方向的速度分別為u和v，壓力為p，并且遠(yuǎn)場(chǎng)邊界位于x=Lx的位置。自由流條件為u#遠(yuǎn)場(chǎng)邊界條件處理

defapply_far_field_boundary_condition(u,v,p,u_inf,p_inf):

"""

應(yīng)用遠(yuǎn)場(chǎng)邊界條件

:paramu:二維數(shù)組，x方向速度

:paramv:二維數(shù)組，y方向速度

:paramp:二維數(shù)組，壓力

:paramu_inf:浮點(diǎn)數(shù)，自由流速度

:paramp_inf:浮點(diǎn)數(shù)，自由流壓力

"""

#假設(shè)遠(yuǎn)場(chǎng)邊界在x=Lx的位置

Lx,Ly=u.shape

u[Lx-1,:]=u_inf#x方向速度等于自由流速度

p[Lx-1,:]=p_inf#壓力等于自由流壓力

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

u=np.zeros((10,10))

v=np.zeros((10,10))

p=np.zeros((10,10))

u_inf=1.0#自由流速度

p_inf=101325.0#自由流壓力

#應(yīng)用遠(yuǎn)場(chǎng)邊界條件

apply_far_field_boundary_condition(u,v,p,u_inf,p_inf)以上示例展示了如何在二維穩(wěn)態(tài)對(duì)流方程的有限差分解法中處理周期性邊界條件、壁面邊界條件和遠(yuǎn)場(chǎng)邊界條件。通過(guò)這些邊界條件的正確應(yīng)用，可以確保數(shù)值模擬的準(zhǔn)確性和穩(wěn)定性。5數(shù)值求解方法在空氣動(dòng)力學(xué)中，數(shù)值方法是解決復(fù)雜流體動(dòng)力學(xué)問(wèn)題的關(guān)鍵工具。有限差分法（FDM）作為其中一種主要的數(shù)值求解技術(shù)，被廣泛應(yīng)用于求解二維穩(wěn)態(tài)對(duì)流方程。本教程將深入探討迭代求解技術(shù)和直接求解方法在二維穩(wěn)態(tài)對(duì)流方程中的應(yīng)用。5.1迭代求解技術(shù)迭代求解技術(shù)是通過(guò)一系列逐步逼近的過(guò)程來(lái)求解線性方程組的方法。在有限差分法中，迭代技術(shù)特別適用于大型稀疏矩陣的求解，這類矩陣在空氣動(dòng)力學(xué)問(wèn)題中很常見(jiàn)。5.1.1原理迭代求解技術(shù)基于一個(gè)初始猜測(cè)值，通過(guò)重復(fù)應(yīng)用一個(gè)迭代公式來(lái)逐步改進(jìn)解的精度。常見(jiàn)的迭代方法包括Jacobi迭代法、Gauss-Seidel迭代法和SOR（SuccessiveOver-Relaxation）迭代法。5.1.1.1Jacobi迭代法Jacobi迭代法是最簡(jiǎn)單的迭代求解技術(shù)之一。它將線性方程組的系數(shù)矩陣分解為對(duì)角矩陣、上三角矩陣和下三角矩陣，然后基于對(duì)角矩陣的逆來(lái)更新解向量。importnumpyasnp

defjacobi(A,b,x0,tol,max_iter):

"""

Jacobi迭代法求解線性方程組Ax=b。

參數(shù):

A:系數(shù)矩陣

b:右側(cè)向量

x0:初始猜測(cè)向量

tol:容忍誤差

max_iter:最大迭代次數(shù)

返回:

x:迭代解向量

"""

D=np.diag(np.diag(A))#對(duì)角矩陣

R=A-D#剩余矩陣

x=x0.copy()

for_inrange(max_iter):

x_new=np.dot(np.linalg.inv(D),b-np.dot(R,x))

ifnp.linalg.norm(x_new-x)<tol:

returnx_new

x=x_new

returnx

#示例：求解二維穩(wěn)態(tài)對(duì)流方程的線性方程組

A=np.array([[4,-1,0,-1],

[-1,4,-1,0],

[0,-1,4,-1],

[-1,0,-1,4]])

b=np.array([1,2,3,4])

x0=np.zeros(4)

tol=1e-6

max_iter=1000

x=jacobi(A,b,x0,tol,max_iter)

print("Jacobi迭代解：",x)5.1.1.2Gauss-Seidel迭代法Gauss-Seidel迭代法是Jacobi迭代法的改進(jìn)版，它在每次迭代中使用最新的解來(lái)更新剩余的解向量元素，從而通常能更快收斂。defgauss_seidel(A,b,x0,tol,max_iter):

"""

Gauss-Seidel迭代法求解線性方程組Ax=b。

參數(shù):

A:系數(shù)矩陣

b:右側(cè)向量

x0:初始猜測(cè)向量

tol:容忍誤差

max_iter:最大迭代次數(shù)

返回:

x:迭代解向量

"""

x=x0.copy()

for_inrange(max_iter):

x_new=x.copy()

foriinrange(len(x)):

x_new[i]=(b[i]-np.dot(A[i,:i],x_new[:i])-np.dot(A[i,i+1:],x[i+1:]))/A[i,i]

ifnp.linalg.norm(x_new-x)<tol:

returnx_new

x=x_new

returnx

#使用Gauss-Seidel迭代法求解同一線性方程組

x=gauss_seidel(A,b,x0,tol,max_iter)

print("Gauss-Seidel迭代解：",x)5.2直接求解方法直接求解方法通過(guò)矩陣分解或消元等技術(shù)直接求得線性方程組的解，而不需要迭代過(guò)程。LU分解和Cholesky分解是兩種常見(jiàn)的直接求解方法。5.2.1原理直接求解方法基于矩陣的分解，將原矩陣分解為更簡(jiǎn)單的矩陣形式，然后通過(guò)前向和后向代換來(lái)求解未知向量。5.2.1.1LU分解LU分解將系數(shù)矩陣分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣U的乘積。一旦分解完成，線性方程組的求解就轉(zhuǎn)化為兩個(gè)三角矩陣的求解問(wèn)題。deflu_decomposition(A,b):

"""

使用LU分解求解線性方程組Ax=b。

參數(shù):

A:系數(shù)矩陣

b:右側(cè)向量

返回:

x:解向量

"""

P,L,U=scipy.linalg.lu(A)

Pb=np.dot(P,b)

y=scipy.linalg.solve_triangular(L,Pb,lower=True)

x=scipy.linalg.solve_triangular(U,y)

returnx

#使用LU分解求解線性方程組

x=lu_decomposition(A,b)

print("LU分解解：",x)5.2.1.2Cholesky分解Cholesky分解適用于對(duì)稱正定矩陣，將矩陣分解為一個(gè)下三角矩陣和其轉(zhuǎn)置的乘積。這種方法在求解空氣動(dòng)力學(xué)中的某些特定問(wèn)題時(shí)非常有效。defcholesky_decomposition(A,b):

"""

使用Cholesky分解求解線性方程組Ax=b。

參數(shù):

A:系數(shù)矩陣（對(duì)稱正定）

b:右側(cè)向量

返回:

x:解向量

"""

L=scipy.linalg.cholesky(A,lower=True)

y=scipy.linalg.solve_triangular(L,b,lower=True)

x=scipy.linalg.solve_triangular(L.T,y)

returnx

#示例：使用Cholesky分解求解對(duì)稱正定矩陣的線性方程組

A_sym=np.array([[4,-1,0,-1],

[-1,4,-1,0],

[0,-1,4,-1],

[-1,0,-1,4]])

b_sym=np.array([1,2,3,4])

x=cholesky_decomposition(A_sym,b_sym)

print("Cholesky分解解：",x)5.3結(jié)論迭代求解技術(shù)和直接求解方法各有優(yōu)勢(shì)，選擇哪種方法取決于具體問(wèn)題的性質(zhì)和求解需求。在空氣動(dòng)力學(xué)數(shù)值模擬中，迭代方法通常用于大型問(wèn)題，而直接方法則適用于較小或特定類型的矩陣問(wèn)題。理解這些方法的原理和應(yīng)用，對(duì)于高效求解空氣動(dòng)力學(xué)中的數(shù)值問(wèn)題至關(guān)重要。6案例分析與結(jié)果驗(yàn)證6.1簡(jiǎn)單二維對(duì)流問(wèn)題的求解在空氣動(dòng)力學(xué)中，二維穩(wěn)態(tài)對(duì)流方程描述了流體在固定幾何結(jié)構(gòu)中流動(dòng)時(shí)，其物理量（如速度、溫度、濃度）隨空間變化的規(guī)律。該方程可以表示為：?其中，u和v分別是流體在x和y方向的速度分量。為了求解這個(gè)方程，我們可以使用有限差分法（FDM），將連續(xù)的偏微分方程轉(zhuǎn)換為離散的代數(shù)方程組。6.1.1示例代碼假設(shè)我們有一個(gè)簡(jiǎn)單的二維矩形區(qū)域，其中流體以恒定速度u=1,v=0沿x方向流動(dòng)。邊界條件為：左側(cè)邊界u=1importnumpyasnp

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

nx,ny=101,101

dx,dy=0.1,0.1

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

u[:,0]=1#左側(cè)邊界條件

#定義迭代參數(shù)

tolerance=1e-6

max_iterations=5000

iteration=0

residual=1

#迭代求解

whileresidual>toleranceanditeration<max_iterations:

un=u.copy()

u[1:-1,1:-1]=(un[1:-1,2:]+un[1:-1,:-2])/2#中心差分

residual=np.max(np.abs(u-un))

iteration+=1

#輸出結(jié)果

print("迭代次數(shù):",iteration)

print("殘差:",residual)6.1.2解釋上述代碼中，我們首先定義了網(wǎng)格的大小和邊界條件。然后，通過(guò)迭代更新內(nèi)部網(wǎng)格點(diǎn)的速度值，使用中心差分近似對(duì)流方程的導(dǎo)數(shù)項(xiàng)。迭代直到殘差小于給定的容差或達(dá)到最大迭代次數(shù)為止。6.2結(jié)果的收斂性檢查收斂性檢查是驗(yàn)證數(shù)值解是否接近真實(shí)解的關(guān)鍵步驟。在FDM中，收斂性通常通過(guò)檢查迭代過(guò)程中殘差的變化來(lái)評(píng)估。殘差是當(dāng)前迭代解與前一次迭代解之間的差異，當(dāng)殘差足夠小，可以認(rèn)為解已經(jīng)收斂。6.2.1示例代碼我們可以修改上述代碼，添加一個(gè)繪制殘差隨迭代次數(shù)變化的圖，以直觀地檢查收斂性。importmatplotlib.pyplotasplt

#...上述代碼...

#存儲(chǔ)殘差歷史

residual_history=[]

#迭代求解

whileresidual>toleranceanditeration<max_iterations:

un=u.copy()

u[1:-1,1:-1]=(un[1:-1,2:]+un[1:-1,:-2])/2

residual=np.max(np.abs(u-un))

residual_history.append(residual)

iteration+=1

#繪制殘差變化圖

plt.figure()

plt.semilogy(residual_history)

plt.xlabel('迭代次數(shù)')

plt.ylabel('殘差')

plt.title('收斂性檢查')

plt.grid(True)

plt.show()6.2.2解釋通過(guò)繪制殘差的對(duì)數(shù)圖，我們可以清晰地看到殘差隨迭代次數(shù)的下降趨勢(shì)，從而判斷解是否收斂。6.3與解析解的比較在可能的情況下，將數(shù)值解與解析解進(jìn)行比較是驗(yàn)證數(shù)值方法準(zhǔn)確性的有效方式。對(duì)于簡(jiǎn)單的二維對(duì)流問(wèn)題，解析解可能不易獲得，但可以通過(guò)理論分析或數(shù)值模擬的高精度解來(lái)近似。6.3.1示例代碼假設(shè)我們有一個(gè)解析解，可以表示為ux,y#...上述代碼...

#定義解析解

u_analytical=1-np.linspace(0,1,nx)

#計(jì)算數(shù)值解在邊界上的值

u_numerical=u[:,0]

#計(jì)算誤差

error=np.abs(u_analytical-u_numerical)

#輸出誤差

print("最大誤差:",np.max(error))6.3.2解釋在這個(gè)例子中，我們計(jì)算了數(shù)值解與解析解在左側(cè)邊界上的差異，以此來(lái)評(píng)估FDM的準(zhǔn)確性。通過(guò)比較最大誤差，我們可以判斷數(shù)值方法是否能夠準(zhǔn)確地模擬物理現(xiàn)象。以上示例展示了如何使用有限差分法求解二維穩(wěn)態(tài)對(duì)流問(wèn)題，以及如何檢查結(jié)果的收斂性和與解析解的準(zhǔn)確性。這些步驟對(duì)于任何使用FDM的空氣動(dòng)力學(xué)數(shù)值模擬都是基本的。7空氣動(dòng)力學(xué)數(shù)值方法：有限差分法(FDM)：高級(jí)主題7.1非結(jié)構(gòu)化網(wǎng)格上的有限差分法在空氣動(dòng)力學(xué)中，非結(jié)構(gòu)化網(wǎng)格允許更靈活地適應(yīng)復(fù)雜幾何形狀，尤其是在處理具有不規(guī)則邊界或需要局部細(xì)化的流場(chǎng)時(shí)。非結(jié)構(gòu)化網(wǎng)格可以是三角形、四邊形、或更高維度的多邊形，它們的節(jié)點(diǎn)和邊的連接沒(méi)有固定模式，這與結(jié)構(gòu)化網(wǎng)格形成對(duì)比，后者通常由規(guī)則排列的矩形或六面體組成。7.1.1原理在非結(jié)構(gòu)化網(wǎng)格上應(yīng)用有限差分法，首先需要確定每個(gè)網(wǎng)格單元的節(jié)點(diǎn)和邊。然后，對(duì)于每個(gè)節(jié)點(diǎn)，我們構(gòu)建一個(gè)基于周?chē)?jié)點(diǎn)的差分方程。這通常涉及到計(jì)算節(jié)點(diǎn)的局部坐標(biāo)系，以及確定與該節(jié)點(diǎn)相鄰的網(wǎng)格單元。對(duì)于對(duì)流方程，我們可能使用中心差分或上風(fēng)差分格式，具體取決于流體的速度方向和網(wǎng)格的局部結(jié)構(gòu)。7.1.2內(nèi)容在非結(jié)構(gòu)化網(wǎng)格上解二維穩(wěn)態(tài)對(duì)流方程，我們首先定義網(wǎng)格和流體速度場(chǎng)。然后，對(duì)于每個(gè)節(jié)點(diǎn)，我們基于其周?chē)?jié)點(diǎn)的值來(lái)構(gòu)建差分方程。例如，對(duì)于節(jié)點(diǎn)i，其差分方程可以表示為：a其中，ai是節(jié)點(diǎn)i的系數(shù)，aj是與節(jié)點(diǎn)i相鄰的節(jié)點(diǎn)j的系數(shù)，N是與節(jié)點(diǎn)i相鄰的節(jié)點(diǎn)數(shù)，?是求解的變量，7.1.3示例假設(shè)我們有一個(gè)非結(jié)構(gòu)化網(wǎng)格，其中包含三角形單元。我們將使用Python和SciPy庫(kù)來(lái)構(gòu)建和求解差分方程。importnumpyasnp

fromscipy.sparseimportlil_matrix

fromscipy.sparse.linalgimportspsolve

#假設(shè)網(wǎng)格節(jié)點(diǎn)和連接信息

nodes=np.array([[0,0],[1,0],[1,1],[0,1]])

elements=np.array([[0,1,2],[0,2,3]])

#流體速度場(chǎng)

velocity=np.array([1,0])

#定義差分方程的系數(shù)矩陣和源項(xiàng)向量

A=lil_matrix((len(nodes),len(nodes)),dtype=np.float64)

b=np.zeros(len(nodes),dtype=np.float64)

#構(gòu)建差分方程

forelementinelements:

#計(jì)算每個(gè)三角形單元的面積

area=0.5*np.abs(np.cross(nodes[element[1]]-nodes[element[0]],nodes[element[2]]-nodes[element[0]]))

#計(jì)算每個(gè)節(jié)點(diǎn)的系數(shù)

foriinrange(3):

forjinrange(3):

ifi!=j:

#使用上風(fēng)差分格式

ifnp.dot(velocity,nodes[element[j]]-nodes[element[i]])>0:

A[element[i],element[i]]+=np.dot(velocity,nodes[element[j]]-nodes[element[i]])/area

A[element[i],element[j]]-=np.dot(velocity,nodes[element[j]]-nodes[element[i]])/area

else:

A[element[i],element[i]]+=0

A[element[i],element[j]]-=0

#假設(shè)邊界條件和初始條件

b[0]=1#設(shè)定邊界條件

#求解差分方程

phi=spsolve(A.tocsr(),b)

#輸出結(jié)果

print("Solution:",phi)在這個(gè)例子中，我們首先定義了一個(gè)非結(jié)構(gòu)化網(wǎng)格，由四個(gè)節(jié)點(diǎn)和兩個(gè)三角形元素組成。然后，我們構(gòu)建了一個(gè)稀疏矩陣A來(lái)表示差分方程的系數(shù)，以及一個(gè)向量b來(lái)表示源項(xiàng)。我們使用上風(fēng)差分格式來(lái)處理對(duì)流項(xiàng)，這確保了數(shù)值穩(wěn)定性。最后，我們使用SciPy的spsolve函數(shù)來(lái)求解線性方程組，得到每個(gè)節(jié)點(diǎn)的解?。7.2高階差分格式的應(yīng)用高階差分格式可以提高有限差分法的精度，尤其是在處理高梯度區(qū)域或需要高分辨率的流場(chǎng)時(shí)。高階格式通常涉及更多的節(jié)點(diǎn)，以更準(zhǔn)確地近似導(dǎo)數(shù)，從而減少數(shù)值誤差。7.2.1原理高階差分格式通過(guò)使用更多的節(jié)點(diǎn)來(lái)構(gòu)建差分方程，從而提高導(dǎo)數(shù)的近似精度。例如，二階中心差分格式使用三個(gè)節(jié)點(diǎn)來(lái)近似一階導(dǎo)數(shù)，而四階中心差分格式則使用五個(gè)節(jié)點(diǎn)。高階格式可以減少數(shù)值擴(kuò)散和振蕩，但可能需要更復(fù)雜的穩(wěn)定性分析和更大的計(jì)算資源。7.2.2內(nèi)容在二維穩(wěn)態(tài)對(duì)流方程中應(yīng)用高階差分格式，我們首先需要確定每個(gè)節(jié)點(diǎn)的局部差分方程。對(duì)于中心差分格式，我們可能使用以下公式來(lái)近似一階導(dǎo)數(shù)：??其中，?是求解的變量，Δx和Δ7.2.3示例我們將使用Python和NumPy庫(kù)來(lái)構(gòu)建和求解使用四階中心差分格式的二維穩(wěn)態(tài)對(duì)流方程。importnumpyasnp

#定義網(wǎng)格和流體速度場(chǎng)

nx,ny=10,10

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

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

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

velocity=np.array([1,0])

#定義差分方程的系數(shù)矩陣和源項(xiàng)向量

A=np.zeros((nx*ny,nx*ny),dtype=np.float64)

b=np.zeros(nx*ny,dtype=np.float64)

#構(gòu)建差分方程

foriinrange(1,nx-1):

forjinrange(1,ny-1):

index=i*ny+j

#使用四階中心差分格式

A[index,index-2*ny]=-1/12

A[index,index-ny]=8/12

A[index,index+ny]=-8/12

A[index,index+2*ny]=1/12

A[index,index-2]=-1/12

A[index,index-1]=8/12

A[index,index+1]=-8/12

A[index,index+2]=1/12

A[index,index]=-20/12+velocity[0]*(1/12+8/12)/