空氣動力學(xué)方程：能量方程：能量方程的數(shù)值解法

空氣動力學(xué)方程：能量方程：能量方程的數(shù)值解法1緒論1.1空氣動力學(xué)與能量方程的重要性空氣動力學(xué)是研究氣體與物體相互作用的科學(xué)，尤其關(guān)注氣體流動對物體的力和能量交換。在航空、汽車設(shè)計、風(fēng)力發(fā)電等領(lǐng)域，理解氣體流動的特性對于設(shè)計高效、安全的系統(tǒng)至關(guān)重要。能量方程，作為空氣動力學(xué)中的一個核心方程，描述了流體流動過程中能量的守恒，包括動能、位能和內(nèi)能的變化。它在分析飛行器的熱力學(xué)性能、預(yù)測發(fā)動機(jī)效率、優(yōu)化風(fēng)力渦輪機(jī)設(shè)計等方面發(fā)揮著關(guān)鍵作用。1.2能量方程的基本概念能量方程基于能量守恒定律，適用于不可壓縮和可壓縮流體。對于不可壓縮流體，能量方程可以簡化為伯努利方程，表達(dá)為：p其中：-p是流體的壓力，-ρ是流體的密度，-v是流體的速度，-g是重力加速度，-z是流體所在位置的高度。對于可壓縮流體，能量方程更為復(fù)雜，需要考慮流體的溫度、壓力和密度的變化。一般形式的能量方程可以表示為：ρ其中：-e是單位質(zhì)量的總能量，-v是流體的速度向量，-q是熱傳導(dǎo)通量，-f是單位質(zhì)量的外力。1.2.1示例：使用Python求解不可壓縮流體的能量方程假設(shè)我們有一個簡單的不可壓縮流體流動問題，流體在管道中流動，管道的直徑在不同位置發(fā)生變化。我們將使用Python來求解在不同位置的壓力分布，假設(shè)流體的速度和高度已知。importnumpyasnp

#定義流體的密度和重力加速度

rho=1.225#流體密度，單位：kg/m^3

g=9.81#重力加速度，單位：m/s^2

#定義管道中不同位置的速度和高度

v=np.array([10,15,20,25])#速度，單位：m/s

z=np.array([0,5,10,15])#高度，單位：m

#定義管道中不同位置的初始壓力

p=np.array([101325,101325,101325,101325])#壓力，單位：Pa

#求解伯努利方程

foriinrange(1,len(v)):

p[i]=p[0]-0.5*rho*(v[i]**2-v[0]**2)-rho*g*(z[i]-z[0])

#打印結(jié)果

print("在不同位置的壓力分布為：")

print(p)在這個例子中，我們首先定義了流體的密度和重力加速度，然后給出了管道中不同位置的速度和高度。我們假設(shè)管道的初始壓力在所有位置都是相同的。通過伯努利方程，我們計算了在不同位置的壓力分布。這個簡單的例子展示了如何使用Python來求解不可壓縮流體的能量方程。1.2.2解釋上述代碼中，我們使用了伯努利方程來計算壓力分布。伯努利方程表明，在流體流動過程中，流體的靜壓、動壓和位壓之和保持不變。通過給定的速度和高度，我們可以計算出在不同位置的壓力變化。這個例子雖然簡單，但它展示了數(shù)值解法的基本思想：將復(fù)雜的物理問題轉(zhuǎn)化為數(shù)學(xué)方程，然后使用計算機(jī)進(jìn)行求解。在實際應(yīng)用中，能量方程的求解可能涉及到更復(fù)雜的數(shù)學(xué)模型和算法，例如有限差分法、有限元法或有限體積法。這些方法通常需要更高級的數(shù)學(xué)知識和編程技巧，但它們能夠處理更復(fù)雜、更實際的流體動力學(xué)問題。在后續(xù)的教程中，我們將深入探討這些數(shù)值解法的原理和應(yīng)用。2能量方程的數(shù)學(xué)基礎(chǔ)2.1連續(xù)性方程簡介在空氣動力學(xué)中，連續(xù)性方程描述了流體在流動過程中質(zhì)量守恒的原理。對于不可壓縮流體，連續(xù)性方程可以簡化為：?其中，ρ是流體的密度，u是流體的速度向量，t是時間。這個方程表明，在任意封閉體積內(nèi)，流體的質(zhì)量不會隨時間改變，除非有流體流入或流出這個體積。2.2動量方程與能量方程的關(guān)系動量方程描述了流體流動時受到的力與流體動量變化之間的關(guān)系。在理想流體中，動量方程可以表示為：ρ其中，p是流體的壓力，g是重力加速度。能量方程則描述了流體流動時能量的守恒，包括動能、位能和內(nèi)能。動量方程與能量方程緊密相關(guān)，因為流體的運動（動量變化）直接影響其能量分布。2.3能量方程的推導(dǎo)能量方程可以從第一定律熱力學(xué)出發(fā)推導(dǎo)。對于單位質(zhì)量的流體，第一定律熱力學(xué)可以表示為：D其中，e是單位質(zhì)量的內(nèi)能，?是單位質(zhì)量的位能，q是單位質(zhì)量的熱源。上式左邊描述了單位質(zhì)量流體的總能量隨時間的變化率，右邊則描述了能量的損失和增加。2.3.1示例：能量方程的數(shù)值解法假設(shè)我們有一個簡單的二維不可壓縮流體流動問題，其中流體的密度和粘度是常數(shù)。我們將使用有限差分法來求解能量方程。首先，我們定義網(wǎng)格和初始條件：importnumpyasnp

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

nx,ny=100,100

dx,dy=1.0,1.0

dt=0.01

#初始化速度和溫度場

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

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

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

#設(shè)置邊界條件

T[0,:]=100.0#下邊界溫度為100

T[-1,:]=50.0#上邊界溫度為50接下來，我們定義能量方程的離散形式，并在時間步長上迭代求解：#定義能量方程的離散形式

defenergy_equation(T,u,v,dt,dx,dy):

T_new=np.zeros_like(T)

foriinrange(1,nx-1):

forjinrange(1,ny-1):

T_new[i,j]=T[i,j]-dt*(u[i,j]*(T[i,j]-T[i-1,j])/dx+

v[i,j]*(T[i,j]-T[i,j-1])/dy)

returnT_new

#迭代求解

fortinrange(1000):

T=energy_equation(T,u,v,dt,dx,dy)2.3.2解釋在上述代碼中，我們首先初始化了速度場和溫度場，并設(shè)置了邊界條件。然后，我們定義了一個函數(shù)energy_equation，它使用有限差分法來離散化能量方程。在每個時間步長，我們更新溫度場T，通過計算流體在網(wǎng)格點上的速度和溫度梯度，來預(yù)測下一個時間步長的溫度分布。這個簡單的例子展示了如何使用數(shù)值方法求解能量方程，但在實際應(yīng)用中，還需要考慮更多的物理效應(yīng)，如熱傳導(dǎo)、粘性耗散等，并且可能需要更復(fù)雜的數(shù)值方法，如有限體積法或有限元法，來確保解的準(zhǔn)確性和穩(wěn)定性。以上內(nèi)容詳細(xì)介紹了能量方程的數(shù)學(xué)基礎(chǔ)，包括連續(xù)性方程、動量方程與能量方程的關(guān)系，以及能量方程的推導(dǎo)。通過一個具體的數(shù)值解法示例，我們展示了如何在計算機(jī)上求解能量方程，盡管這個示例非常簡化，但它為理解更復(fù)雜問題的數(shù)值求解提供了一個基礎(chǔ)框架。3能量方程的數(shù)值方法3.1有限差分法基礎(chǔ)3.1.1原理有限差分法是求解偏微分方程的一種數(shù)值方法，它通過將連續(xù)的偏微分方程離散化為一系列離散的代數(shù)方程來近似求解。在空氣動力學(xué)中，能量方程描述了流體的能量守恒，包括動能、位能和內(nèi)能的變化。有限差分法通過在空間和時間上對能量方程進(jìn)行差分，將偏導(dǎo)數(shù)轉(zhuǎn)換為差商，從而將能量方程轉(zhuǎn)化為可以數(shù)值求解的形式。3.1.2內(nèi)容3.1.2.1空間差分考慮一維能量方程：?其中，E是總能量，u是流速，p是壓力，μ是動力粘度，Q是熱源項。在空間上，我們可以通過中心差分來近似?Eu???3.1.2.2時間差分時間上，我們可以使用顯式歐拉法或隱式歐拉法來近似?E?其中，Δt是時間步長，n3.1.2.3數(shù)值解法將上述差分公式代入能量方程，可以得到離散化的能量方程。然后，通過迭代求解這些離散方程，可以得到能量E在每個網(wǎng)格點上的數(shù)值解。3.1.3代碼示例importnumpyasnp

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

L=1.0#空間長度

N=100#空間網(wǎng)格數(shù)

T=1.0#時間長度

M=1000#時間步數(shù)

dx=L/N

dt=T/M

mu=0.1#動力粘度

Q=0.0#熱源項

#初始條件和邊界條件

E=np.zeros(N+1)

E[0]=1.0#左邊界條件

E[N]=0.0#右邊界條件

#迭代求解

forninrange(M):

foriinrange(1,N):

E[i]=E[i]-dt/dx*(E[i]*u[i]-E[i-1]*u[i-1])+dt/dx*(mu*(u[i+1]-u[i-1])/(2*dx))+dt*Q

#輸出結(jié)果

print(E)3.2有限體積法簡介3.2.1原理有限體積法是另一種求解偏微分方程的數(shù)值方法，它基于守恒定律，將計算域劃分為一系列控制體積，然后在每個控制體積上應(yīng)用守恒定律。在空氣動力學(xué)中，能量方程的有限體積法求解，就是將能量守恒定律應(yīng)用于每個控制體積，得到能量在每個控制體積內(nèi)的平均值。3.2.2內(nèi)容3.2.2.1控制體積考慮一維能量方程，我們可以在空間上劃分一系列控制體積，每個控制體積的中心點為網(wǎng)格點。在每個控制體積上，能量方程可以寫作：d其中，Eui+1/3.2.2.2數(shù)值通量數(shù)值通量可以通過不同的數(shù)值格式來計算，例如，一階迎風(fēng)格式或二階中心格式。這里以一階迎風(fēng)格式為例：E3.2.2.3數(shù)值解法將上述控制體積方程和數(shù)值通量公式代入，可以得到離散化的能量方程。然后，通過迭代求解這些離散方程，可以得到能量E在每個控制體積內(nèi)的平均值。3.2.3代碼示例importnumpyasnp

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

L=1.0#空間長度

N=100#空間網(wǎng)格數(shù)

T=1.0#時間長度

M=1000#時間步數(shù)

dx=L/N

dt=T/M

mu=0.1#動力粘度

Q=0.0#熱源項

#初始條件和邊界條件

E=np.zeros(N+1)

E[0]=1.0#左邊界條件

E[N]=0.0#右邊界條件

#迭代求解

forninrange(M):

E_flux=np.zeros(N+1)

foriinrange(1,N):

ifu[i]>0:

E_flux[i]=E[i]*u[i]

else:

E_flux[i]=E[i+1]*u[i+1]

foriinrange(1,N):

E[i]=E[i]-dt/dx*(E_flux[i+1]-E_flux[i])+dt/dx*(mu*(u[i+1]-u[i-1])/(2*dx))+dt*Q

#輸出結(jié)果

print(E)3.3有限元法概述3.3.1原理有限元法是一種基于變分原理的數(shù)值方法，它將連續(xù)的偏微分方程轉(zhuǎn)化為離散的代數(shù)方程組。在空氣動力學(xué)中，能量方程的有限元法求解，就是將能量方程轉(zhuǎn)化為一個能量泛函的極小化問題，然后通過有限元方法求解這個泛函的極小值。3.3.2內(nèi)容3.3.2.1變分原理考慮一維能量方程，我們可以將其轉(zhuǎn)化為一個能量泛函的極小化問題：min其中，ρ是密度，cp是比熱容，T3.3.2.2有限元離散在有限元方法中，我們首先將計算域劃分為一系列有限元，然后在每個有限元上，使用插值函數(shù)來近似能量E。插值函數(shù)通常是一些低階多項式，例如線性插值函數(shù)或二次插值函數(shù)。3.3.2.3數(shù)值解法將上述能量泛函和插值函數(shù)代入，可以得到離散化的能量方程。然后，通過求解這個離散化的能量方程，可以得到能量E在每個有限元上的數(shù)值解。3.3.3代碼示例importnumpyasnp

fromscipy.sparseimportdiags

fromscipy.sparse.linalgimportspsolve

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

L=1.0#空間長度

N=100#空間網(wǎng)格數(shù)

T=1.0#時間長度

M=1000#時間步數(shù)

dx=L/N

dt=T/M

mu=0.1#動力粘度

Q=0.0#熱源項

#初始條件和邊界條件

E=np.zeros(N+1)

E[0]=1.0#左邊界條件

E[N]=0.0#右邊界條件

#構(gòu)建有限元矩陣

data=np.array([[-1.0],[2.0],[-1.0]])

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

A=diags(data,diags,shape=(N-1,N-1))

#迭代求解

forninrange(M):

b=np.zeros(N-1)

foriinrange(1,N):

b[i-1]=E[i]-dt/dx*(E[i]*u[i]-E[i-1]*u[i-1])+dt/dx*(mu*(u[i+1]-u[i-1])/(2*dx))+dt*Q

E[1:N]=spsolve(A,b)

#輸出結(jié)果

print(E)請注意，上述代碼示例僅為示意圖，實際應(yīng)用中需要根據(jù)具體問題和邊界條件進(jìn)行調(diào)整。4數(shù)值解法的實施步驟在解決空氣動力學(xué)方程中的能量方程時，數(shù)值解法是一種常用且有效的方法。本教程將詳細(xì)探討實施數(shù)值解法的三個關(guān)鍵步驟：網(wǎng)格生成技術(shù)、邊界條件的設(shè)定、以及初始條件的選擇。4.1網(wǎng)格生成技術(shù)網(wǎng)格生成是數(shù)值模擬的第一步，它將連續(xù)的物理域離散化為一系列有限的、非重疊的單元，以便于數(shù)值計算。網(wǎng)格的質(zhì)量直接影響到數(shù)值解的準(zhǔn)確性和計算效率。4.1.1均勻網(wǎng)格在簡單幾何形狀中，可以使用均勻網(wǎng)格。例如，考慮一個二維矩形區(qū)域，我們可以創(chuàng)建一個由等間距節(jié)點組成的網(wǎng)格。importnumpyasnp

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

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

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

Lx=1.0#x方向的總長度

Ly=0.5#y方向的總長度

#生成網(wǎng)格

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

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

X,Y=np.meshgrid(x,y)4.1.2非均勻網(wǎng)格對于復(fù)雜的幾何形狀或需要在特定區(qū)域提高分辨率的情況，非均勻網(wǎng)格更為適用。例如，可以使用非均勻網(wǎng)格來更好地捕捉流體邊界層的細(xì)節(jié)。importnumpyasnp

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

nx=100

ny=50

Lx=1.0

Ly=0.5

#使用非均勻分布生成網(wǎng)格

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

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

X,Y=np.meshgrid(x,y)4.2邊界條件的設(shè)定邊界條件描述了在計算域邊界上的物理狀態(tài)，對于能量方程而言，邊界條件通常涉及溫度、熱流或絕熱條件。4.2.1絕熱邊界在絕熱邊界上，沒有熱量交換，這意味著邊界上的熱流為零。#絕熱邊界條件

defset_adiabatic_boundary(T,dx,dy):

"""

設(shè)置絕熱邊界條件

T:溫度矩陣

dx,dy:網(wǎng)格間距

"""

#上下邊界

T[0,:]=T[1,:]#底部邊界

T[-1,:]=T[-2,:]#頂部邊界

#左右邊界

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

T[:,-1]=T[:,-2]#右側(cè)邊界4.2.2固定溫度邊界在固定溫度邊界上，邊界上的溫度保持不變。#固定溫度邊界條件

defset_fixed_temperature_boundary(T,T_left,T_right,T_top,T_bottom):

"""

設(shè)置固定溫度邊界條件

T:溫度矩陣

T_left,T_right,T_top,T_bottom:邊界溫度

"""

T[0,:]=T_bottom#底部邊界

T[-1,:]=T_top#頂部邊界

T[:,0]=T_left#左側(cè)邊界

T[:,-1]=T_right#右側(cè)邊界4.3初始條件的選擇初始條件是數(shù)值模擬開始時的物理狀態(tài)，對于能量方程，這通常意味著初始溫度分布。4.3.1均勻初始溫度在許多情況下，可以假設(shè)初始溫度在計算域內(nèi)是均勻的。#均勻初始溫度

defset_uniform_initial_temperature(T,T_initial):

"""

設(shè)置均勻初始溫度

T:溫度矩陣

T_initial:初始溫度

"""

T.fill(T_initial)4.3.2非均勻初始溫度在某些情況下，初始溫度分布可能不均勻，例如，當(dāng)計算域內(nèi)存在熱源時。importnumpyasnp

#非均勻初始溫度

defset_non_uniform_initial_temperature(T,T_initial,heat_source):

"""

設(shè)置非均勻初始溫度

T:溫度矩陣

T_initial:初始溫度

heat_source:熱源位置和強(qiáng)度的列表

"""

T.fill(T_initial)

forsourceinheat_source:

x,y,intensity=source

T[int(y/dy),int(x/dx)]+=intensity4.3.3示例：非均勻初始溫度和邊界條件假設(shè)我們有一個2D矩形區(qū)域，初始溫度為300K，底部邊界溫度為350K，其余邊界為絕熱。計算域內(nèi)存在一個熱源，位于(0.25,0.25)，強(qiáng)度為50K。importnumpyasnp

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

nx=100

ny=50

Lx=1.0

Ly=0.5

dx=Lx/(nx-1)

dy=Ly/(ny-1)

#溫度矩陣初始化

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

#設(shè)置初始溫度

T_initial=300

set_uniform_initial_temperature(T,T_initial)

#設(shè)置邊界條件

T_bottom=350

set_fixed_temperature_boundary(T,T_initial,T_initial,T_bottom,T_initial)

set_adiabatic_boundary(T,dx,dy)

#設(shè)置熱源

heat_source=[(0.25,0.25,50)]

set_non_uniform_initial_temperature(T,T_initial,heat_source)通過以上步驟，我們已經(jīng)為能量方程的數(shù)值解法準(zhǔn)備好了網(wǎng)格、邊界條件和初始條件。接下來，可以使用如有限差分、有限體積或有限元等方法來求解能量方程。在實際應(yīng)用中，這些步驟可能需要根據(jù)具體問題進(jìn)行調(diào)整，以確保數(shù)值解的準(zhǔn)確性和穩(wěn)定性。5能量方程的離散化5.1時間離散化方法在數(shù)值求解空氣動力學(xué)方程中的能量方程時，時間離散化是將連續(xù)的時間域轉(zhuǎn)換為離散時間步的關(guān)鍵步驟。這通常通過時間積分方法實現(xiàn)，其中最常見的是顯式和隱式方法。5.1.1顯式方法顯式方法是一種直接計算未來狀態(tài)的方法，它基于當(dāng)前狀態(tài)和時間步長。例如，一維能量方程的顯式歐拉方法可以表示為：E其中，E是能量，u是速度，λ是熱導(dǎo)率，Δt是時間步長，n5.1.2隱式方法隱式方法則需要解一個方程組來找到未來狀態(tài)，它考慮了未來狀態(tài)對當(dāng)前狀態(tài)的影響。例如，隱式歐拉方法可以表示為：E5.1.3Python代碼示例：顯式歐拉方法importnumpyasnp

defexplicit_euler(E,u,lambda_,dx,dt):

"""

使用顯式歐拉方法離散化能量方程。

參數(shù):

E(np.array):當(dāng)前能量分布。

u(np.array):當(dāng)前速度分布。

lambda_(np.array):當(dāng)前熱導(dǎo)率分布。

dx(float):空間步長。

dt(float):時間步長。

返回:

np.array:下一時間步的能量分布。

"""

dE_dx=np.gradient(E,dx)

dlambda_dE_dx=np.gradient(lambda_*dE_dx,dx)

E_next=E+dt*(-u*dE_dx+dlambda_dE_dx)

returnE_next

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

E=np.array([100,105,110,115,120])

u=np.array([0.1,0.2,0.3,0.4,0.5])

lambda_=np.array([0.01,0.02,0.03,0.04,0.05])

dx=1.0

dt=0.1

#調(diào)用函數(shù)

E_next=explicit_euler(E,u,lambda_,dx,dt)

print("下一時間步的能量分布:",E_next)5.2空間離散化技術(shù)空間離散化技術(shù)用于將連續(xù)的空間域轉(zhuǎn)換為離散的網(wǎng)格點。常見的空間離散化技術(shù)包括有限差分、有限體積和有限元方法。5.2.1有限差分方法有限差分方法通過在網(wǎng)格點上用差商代替導(dǎo)數(shù)來離散化方程。例如，能量方程的一階向前差分可以表示為：?5.2.2有限體積方法有限體積方法將空間域劃分為體積單元，并在每個單元上應(yīng)用守恒定律。能量方程的有限體積離散化可以表示為：E其中，F(xiàn)是能量通量。5.2.3Python代碼示例：有限差分方法deffinite_difference(E,dx):

"""

使用有限差分方法計算能量的一階導(dǎo)數(shù)。

參數(shù):

E(np.array):能量分布。

dx(float):空間步長。

返回:

np.array:能量的一階導(dǎo)數(shù)。

"""

dE_dx=np.gradient(E,dx)

returndE_dx

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

E=np.array([100,105,110,115,120])

dx=1.0

#調(diào)用函數(shù)

dE_dx=finite_difference(E,dx)

print("能量的一階導(dǎo)數(shù):",dE_dx)5.3離散方程的穩(wěn)定性分析穩(wěn)定性分析是確保數(shù)值解法不會產(chǎn)生無物理意義的解的關(guān)鍵步驟。常見的穩(wěn)定性分析方法包括Fourier分析和VonNeumann穩(wěn)定性分析。5.3.1Fourier分析Fourier分析通過將解表示為Fourier級數(shù)來檢查離散方程的穩(wěn)定性。如果Fourier級數(shù)的系數(shù)隨時間增加而減小，則方程是穩(wěn)定的。5.3.2VonNeumann穩(wěn)定性分析VonNeumann穩(wěn)定性分析通過假設(shè)解是空間和時間的指數(shù)函數(shù)來檢查穩(wěn)定性。如果解的模隨時間增加而減小，則方程是穩(wěn)定的。5.3.3Python代碼示例：VonNeumann穩(wěn)定性分析importnumpyasnp

defvon_neumann_stability(u,lambda_,dx,dt):

"""

使用VonNeumann穩(wěn)定性分析檢查能量方程的穩(wěn)定性。

參數(shù):

u(float):速度。

lambda_(float):熱導(dǎo)率。

dx(float):空間步長。

dt(float):時間步長。

返回:

bool:方程是否穩(wěn)定。

"""

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

stability_condition=lambda_*dt/(dx**2)

#檢查穩(wěn)定性條件是否滿足

ifstability_condition<=0.5:

returnTrue

else:

returnFalse

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

u=0.2

lambda_=0.02

dx=1.0

dt=0.1

#調(diào)用函數(shù)

is_stable=von_neumann_stability(u,lambda_,dx,dt)

print("能量方程是否穩(wěn)定:",is_stable)以上示例展示了如何使用Python實現(xiàn)能量方程的時間和空間離散化，以及如何進(jìn)行穩(wěn)定性分析。這些方法是數(shù)值求解空氣動力學(xué)方程中能量方程的基礎(chǔ)，通過調(diào)整時間步長和空間步長，可以確保數(shù)值解的準(zhǔn)確性和穩(wěn)定性。6數(shù)值求解算法在空氣動力學(xué)方程的數(shù)值解法中，能量方程的求解是關(guān)鍵步驟之一。本教程將詳細(xì)介紹三種數(shù)值求解算法：迭代求解方法、直接求解技術(shù)以及多網(wǎng)格方法。我們將通過原理闡述和具體示例來理解這些方法如何應(yīng)用于能量方程的求解。6.1迭代求解方法迭代求解方法是求解非線性方程組的常用技術(shù)，尤其適用于大型稀疏矩陣。在空氣動力學(xué)中，能量方程通常與連續(xù)性方程、動量方程等耦合，形成復(fù)雜的非線性系統(tǒng)。迭代方法通過逐步逼近的方式，最終達(dá)到方程的解。6.1.1原理迭代求解方法基于一個初始猜測值，通過反復(fù)應(yīng)用一個迭代公式來逐步改進(jìn)解的精度。常見的迭代方法包括Jacobi迭代、Gauss-Seidel迭代和SOR（SuccessiveOver-Relaxation）迭代。6.1.2示例：Gauss-Seidel迭代假設(shè)我們有以下形式的能量方程：?其中，E是能量，u和v分別是x和y方向的速度，α是熱導(dǎo)率，S是源項。在離散化后，我們得到一個線性方程組，可以使用Gauss-Seidel迭代方法求解。以下是一個使用Python實現(xiàn)的Gauss-Seidel迭代示例：importnumpyasnp

defgauss_seidel(A,b,x0,omega,max_iter=100,tol=1e-6):

"""

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

參數(shù):

A:系數(shù)矩陣

b:右側(cè)向量

x0:初始猜測向量

omega:松弛因子

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

tol:收斂容差

"""

x=x0.copy()

for_inrange(max_iter):

x_new=x.copy()

foriinrange(A.shape[0]):

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]

x_new[i]=omega*x_new[i]+(1-omega)*x[i]

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

returnx_new

x=x_new

returnx

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

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

[1,15.5,3,8],

[0,-1.3,-4,1.1],

[14,5,-2,30]])

b=np.array([1,1,1,1])

x0=np.zeros(4)

omega=1.0

#運行迭代

x=gauss_seidel(A,b,x0,omega)

print("迭代解:",x)在實際應(yīng)用中，A矩陣和b向量將由能量方程的離散化得到，而松弛因子ω的選擇對收斂速度有重要影響。6.2直接求解技術(shù)直接求解技術(shù)通過矩陣分解或消元等方法，直接求得方程的解。對于小型或中型問題，直接求解方法通常比迭代方法更有效。6.2.1原理直接求解技術(shù)包括LU分解、Cholesky分解和QR分解等。其中，LU分解是最常用的方法之一，它將系數(shù)矩陣分解為一個下三角矩陣和一個上三角矩陣的乘積，然后通過前向和后向代換求解未知向量。6.2.2示例：LU分解以下是一個使用Python和SciPy庫實現(xiàn)的LU分解示例：importnumpyasnp

fromscipy.linalgimportlu

defdirect_solve(A,b):

"""

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

參數(shù):

A:系數(shù)矩陣

b:右側(cè)向量

"""

P,L,U=lu(A)

y=np.linalg.solve(L,np.dot(P.T,b))

x=np.linalg.solve(U,y)

returnx

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

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

[1,15.5,3,8],

[0,-1.3,-4,1.1],

[14,5,-2,30]])

b=np.array([1,1,1,1])

#運行直接求解

x=direct_solve(A,b)

print("直接解:",x)6.3多網(wǎng)格方法多網(wǎng)格方法是一種加速迭代求解的高效技術(shù)，通過在不同網(wǎng)格尺度上交替求解，可以顯著提高收斂速度。6.3.1原理多網(wǎng)格方法基于一個基本思想：在粗網(wǎng)格上求解可以快速消除低頻誤差，而在細(xì)網(wǎng)格上求解可以精確處理高頻誤差。通過在不同網(wǎng)格尺度上迭代求解，可以同時處理不同頻率的誤差，從而加速收斂。6.3.2示例：多網(wǎng)格迭代多網(wǎng)格方法的實現(xiàn)較為復(fù)雜，涉及到網(wǎng)格的細(xì)化和粗化、限制和插值等操作。以下是一個簡化的多網(wǎng)格迭代示例，使用Python實現(xiàn)：importnumpyasnp

defmultigrid(A,b,x0,max_iter=100,tol=1e-6):

"""

使用多網(wǎng)格迭代方法求解線性方程組Ax=b。

參數(shù):

A:系數(shù)矩陣

b:右側(cè)向量

x0:初始猜測向量

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

tol:收斂容差

"""

x=x0.copy()

for_inrange(max_iter):

#粗化網(wǎng)格

A_coarse=A[::2,::2]

b_coarse=b[::2]

x_coarse=x[::2].copy()

#在粗網(wǎng)格上求解

x_coarse=gauss_seidel(A_coarse,b_coarse,x_coarse,1.0)

#插值回細(xì)網(wǎng)格

x[::2]=x_coarse

x[1::2]=x[::2]

#在細(xì)網(wǎng)格上求解

x=gauss_seidel(A,b,x,1.0)

ifnp.linalg.norm(A@x-b)<tol:

returnx

returnx

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

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

[1,15.5,3,8],

[0,-1.3,-4,1.1],

[14,5,-2,30]])

b=np.array([1,1,1,1])

x0=np.zeros(4)

#運行多網(wǎng)格迭代

x=multigrid(A,b,x0)

print("多網(wǎng)格迭代解:",x)請注意，上述示例僅用于說明多網(wǎng)格方法的基本思想，實際應(yīng)用中需要更復(fù)雜的網(wǎng)格管理和求解策略。通過上述介紹和示例，我們了解了迭代求解方法、直接求解技術(shù)和多網(wǎng)格方法在空氣動力學(xué)能量方程數(shù)值解法中的應(yīng)用。每種方法都有其適用場景和優(yōu)缺點，選擇合適的方法對于提高求解效率和精度至關(guān)重要。7案例分析與應(yīng)用7.1維流場的能量方程求解7.1.1原理在空氣動力學(xué)中，能量方程描述了流體流動過程中能量的守恒。對于二維流場，能量方程可以簡化為：ρ其中，ρ是流體密度，u和v分別是流體在x和y方向的速度，h是焓，k是熱導(dǎo)率，T是溫度，q是熱源項。7.1.2內(nèi)容7.1.2.1數(shù)值解法在數(shù)值模擬中，我們通常使用有限差分法或有限體積法來離散化上述方程。這里，我們以有限差分法為例，將方程在空間上離散化為：ρ7.1.2.2代碼示例importnumpyasnp

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

nx,ny=100,100

dx,dy=0.1,0.1

rho=1.225#空氣密度

k=0.026#空氣熱導(dǎo)率

q=0#無熱源

u=np.zeros((nx,ny))#x方向速度

v=np.zeros((nx,ny))#y方向速度

h=np.zeros((nx,ny))#焓

T=np.zeros((nx,ny))#溫度

#初始化邊界條件

T[:,0]=300#下邊界溫度

T[:,-1]=300#上邊界溫度

T[0,:]=300#左邊界溫度

T[-1,:]=300#右邊界溫度

#迭代求解

foriterinrange(1000):

T_new=np.copy(T)

foriinrange(1,nx-1):

forjinrange(1,ny-1):

T_new[i,j]=T[i,j]+(k/(rho*dx**2))*(T[i+1,j]-2*T[i,j]+T[i-1,j])+(k/(rho*dy**2))*(T[i,j+1]-2*T[i,j]+T[i,j-1])+(q/rho)

T=np.copy(T_new)

#輸出結(jié)果

print(T)7.1.3描述上述代碼示例展示了如何使用有限差分法求解二維流場的能量方程。首先，我們定義了網(wǎng)格參數(shù)和流體屬性，然后初始化了邊界條件。通過迭代更新內(nèi)部網(wǎng)格點的溫度，最終得到整個流場的溫度分布。7.2維流場的數(shù)值模擬7.2.1原理三維流場的能量方程更為復(fù)雜，包括了在x、y和z三個方向上的能量傳輸：ρ其中，w是流體在z方向的速度。7.2.2內(nèi)容7.2.2.1數(shù)值解法三維流場的能量方程可以使用有限體積法進(jìn)行離散化，這種方法在處理復(fù)雜幾何形狀和邊界條件時更為有效。離散化后的方程如下：V7.2.2.2代碼示例importnumpyasnp

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

nx,ny,nz=50,50,50

dx,dy,dz=0.1,0.1,0.1

rho=1.225#空氣密度

k=0.026#空氣熱導(dǎo)率

q=0#無熱源

u=np.zeros((nx,ny,nz))#x方向速度

v=np.zeros((nx,ny,nz))#y方向速度

w=np.zeros((nx,ny,nz))#z方向速度

h=np.zeros((nx,ny,nz))#焓

T=np.zeros((nx,ny,nz))#溫度

#初始化邊界條件

T[:,:,0]=300#下邊界溫度

T[:,:,-1]=300#上邊界溫度

T[:,0,:]=300#左邊界溫度

T[:,-1,:]=300#右邊界溫度

T[0,:,:]=300#前邊界溫度

T[-1,:,:]=300#后邊界溫度

#迭代求解

foriterinrange(1000):

T_new=np.copy(T)

foriinrange(1,nx-1):

forjinrange(1,ny-1):

forkinrange(1,nz-1):

T_new[i,j,k]=T[i,j,k]+(k/(rho*dx**2))*(T[i+1,j,k]-2*T[i,j,k]+T[i-1,j,k])+(k/(rho*dy**2))*(T[i,j+1,k]-2*T[i,j,k]+T[i,j-1,k])+(k/(rho*dz**2))*(T[i,j,k+1]-2*T[i,j,k]+T[i,j,k-1])+(q/rho)

T=np.copy(T_new)

#輸出結(jié)果

print(T)7.2.3描述這段代碼示例展示了三維流場能量方程的數(shù)值求解過程。與二維情況類似，我們定義了網(wǎng)格參數(shù)和流體屬性，初始化了邊界條件，并通過迭代更新內(nèi)部網(wǎng)格點的溫度，最終得到三維流場的溫度分布。7.3實際飛行器的能量方程分析7.3.1原理在實際飛行器的能量方程分析中，除了考慮流體的動能和內(nèi)能，還需要考慮飛行器表面的熱交換、飛行器內(nèi)部的熱源以及飛行器運動引起的壓縮性和粘性效應(yīng)。7.3.2內(nèi)容7.3.2.1數(shù)值解法實際飛行器的能量方程分析通常需要使用更為復(fù)雜的數(shù)值方法，如高階有限體積法或有限元法，以及非結(jié)構(gòu)化網(wǎng)格。此外，還需要考慮非線性效應(yīng)和湍流模型。7.3.2.2代碼示例importnumpyasnp

fromscipy.sparseimportdiags

fromscipy.sparse.linalgimportspsolve

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

nx,ny=100,100

dx,dy=0.1,0.1

rho=1.225#空氣密度

k=0.026#空氣熱導(dǎo)率

q=0#無熱源

u=np.zeros((nx,ny))#x方向速度

v=np.zeros((nx,ny))#y方向速度

h=np.zeros((nx,ny))#焓

T=np.zeros((nx,ny))#溫度

#初始化邊界條件

T[:,0]=300#下邊界溫度

T[:,-1]=300#上邊界溫度

T[0,:]=300#左邊界溫度

T[-1,:]=300#右邊界溫度

#構(gòu)建系數(shù)矩陣

data=[np.ones(nx-2),-2*np.ones(nx-2),np.ones(nx-2)]

diags=[0,-1,1]

A=diags(data,diags,shape=(nx-2,nx-2))

A=A/(dx**2)

data=[np.ones(ny-2),-2*np.ones(ny-2),np.ones(ny-2)]

diags=[0,-1,1]

B=diags(data,diags,shape=(ny-2,ny-2))

B=B/(dy**2)

#迭代求解

foriterinrange(1000):

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

foriinrange(1,nx-1):

forjinrange(1,ny-1):

T_new[i,j]=spsolve(A+B,(k/(rho*dx**2))*(T[i+1,j]-2*T[i,j]+T[i-1,j])+(k/(rho*dy**2))*(T[i,j+1]-2*T[i,j]+T[i,j-1])+(q/rho))

T=np.copy(T_new)

#輸出結(jié)果

print(T)7.3.3描述這段代碼示例展示了如何使用高階有限體積法求解實際飛行器的能量方程。我們使用了scipy庫中的spsolve函數(shù)來求解系數(shù)矩陣和熱源項的線性方程組，這種方法在處理大型矩陣時更為高效。通過迭代更新內(nèi)部網(wǎng)格點的溫度，最終得到飛行器周圍流場的溫度分布。請注意，上述代碼示例僅為簡化版，實際飛行器的能量方程分析需要考慮更多因素，如飛行器表面的熱交換、飛行器內(nèi)部的熱源以及飛行器運動引起的壓縮性和粘性效應(yīng)。此外，還需要使用更為復(fù)雜的數(shù)值方法和非結(jié)構(gòu)化網(wǎng)格來提高計算精度和效率。8結(jié)果驗證與后處理8.1數(shù)值解的收斂性檢查數(shù)值解法在求解空氣動力學(xué)方程時，收斂性檢查是確保解的準(zhǔn)確性和穩(wěn)定性的重要步驟。收斂性意味著隨著迭代次數(shù)的增加，解逐漸接近一個穩(wěn)定的值。檢查收斂性的方法通常包括觀察殘差、比較迭代解的差異以及使用網(wǎng)格獨立性測試。8.1.1觀察殘差殘差是迭代過程中方程左側(cè)和右側(cè)的差值，它反映了當(dāng)前解與方程的符合程度。當(dāng)殘差下降到一個預(yù)設(shè)的閾值時，可以認(rèn)為解已經(jīng)收斂。#示例代碼：使用Python檢查殘差

importnumpyasnp

#假設(shè)我們有能量方程的殘差數(shù)據(jù)

residuals=np.array([1.0,0.5,0.25,0.125,0.0625,0.03125,0.015625])

#預(yù)設(shè)的收斂閾值

convergence_threshold=0.01

#檢查殘差是否低于閾值

converged=np.all(residuals<convergence_threshold)

#輸出結(jié)果

ifconverged:

print("數(shù)值解已收斂")

else:

print("數(shù)值解未收斂")8.1.2比較迭代解的差異另一種方法是通過比較連續(xù)迭代之間的解的差異來判斷收斂性。如果差異足夠小，可以認(rèn)為解已經(jīng)穩(wěn)定。#示例代碼：使用Python比較迭代解的差異

#假設(shè)我們有連續(xù)迭代的解

solutions=np.array([[1.0,2.0,3.0],[1.01,2.01,3.01],[1.005,2.005,3.005]])

#計算連續(xù)迭代之間的差異

differences=np.abs(solutions[1:]-solutions[:-1])

#預(yù)設(shè)的差異閾值

difference_threshold=0.001

#檢查所有差異是否低于閾值

converged=np.all(differences<difference_threshold)

#輸出結(jié)果

ifconverged:

print("數(shù)值解已收斂")

else:

print("數(shù)值解未收斂")8.1.3網(wǎng)格獨立性測試網(wǎng)格獨立性測試是通過比較不同網(wǎng)格密度下的解來確保解的準(zhǔn)確性。如果在不同網(wǎng)格密度下解的差異足夠小，可以認(rèn)為解是網(wǎng)格獨立的。#示例代碼：使用Python進(jìn)行網(wǎng)格獨立性測試

#假設(shè)我們有不同網(wǎng)格密度下的解

solution_coarse=np.array([1.0,2.0,3.0])

solution_medium=np.array([1.001,2.001,3.001])

solution_fine=np.array([1.0001,2.0001,3.0001])

#計算不同網(wǎng)格密度下的解的差異

difference_coarse_medium=np.abs(solution_coarse-solution_medium)

difference_medium_fine=np.abs(solution_medium-solution_fine)

#預(yù)設(shè)的差異閾值

grid_independence_threshold=0.001

#檢查差異是否低于閾值

grid_independence=np.all(difference_coarse_medium<grid_independence_threshold)andnp.all(difference_medium_fine<grid_independence_threshold)

#輸出結(jié)果

ifgrid_independence:

print("解是網(wǎng)格獨立的")

else:

print("解不是網(wǎng)格獨立的")8.2物理量的后處理后處理階段涉及對數(shù)值解進(jìn)行分析，提取關(guān)鍵的物理量，如壓力、溫度、速度等，以便于理解和解釋結(jié)果。8.2.1提取物理量使用數(shù)值解，可以計算出流場中的各種物理量。例如，從速度場中計算出動能。#示例代碼：使用Python從速度場計算動能

importnumpyasnp

#假設(shè)我們有速度場數(shù)據(jù)

velocity_field=np.array([[1.0,2.0,3.0],[1.5,2.5,3.5],[2.0,3.0,4.0]])

#計算動能

kinetic_energy=0.5*np.sum(velocity_field**2,axis=1)

#輸出動能

print("動能:",kinetic_energy)8.2.2可視化結(jié)果可視化是理解流場行為的關(guān)鍵工具?？梢允褂肞ython的matplotlib庫來繪制速度矢量圖、壓力分布圖等。#示例代碼：使用Python的matplotlib庫繪制速度矢量圖

importmatplotlib.pyplotasplt

importnumpyasnp

#假設(shè)我們有速度場數(shù)據(jù)

velocity_field=np.array([[1.0,2.0],[1.5,2.5]])

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

x=np.linspace(0,1,velocity_field.shape[1])

y=np.linspace(0,1,velocity_field.shape[0])

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

#繪制速度矢量圖

plt.quiver(X,Y,velocity_field[:,:,0],velocity_field[:,:,1])

plt.xlabel('x')

plt.ylabel('y')

plt.title('速度矢量圖')

plt.show()8.3與實驗數(shù)據(jù)的對比分析將數(shù)值解與實驗數(shù)據(jù)進(jìn)行對比，是驗證模型準(zhǔn)確性的關(guān)鍵步驟。這通常涉及到計算誤差百分比或相關(guān)系數(shù)。8.3.1計算誤差百分比誤差百分比是數(shù)值解與實驗數(shù)據(jù)之間的差異的量化指標(biāo)。#示例代碼：使用Python計算誤差百分比

importnumpyasnp

#假設(shè)我們有實驗數(shù)據(jù)和數(shù)值解

experimental_data=np.array([1.0,2.0,3.0])

numerical_solution=np.array([1.01,2.01,3.01])

#計算誤差百分比

error_percentage=np.abs((experimental_data-numerical_solution)/experimental_data)*100

#輸出誤差百分比

print("誤差百分比:",error_percentage)8.3.2計算相關(guān)系數(shù)相關(guān)系數(shù)衡量了數(shù)值解與實驗數(shù)據(jù)之間的線性關(guān)系。一個接近1的相關(guān)系數(shù)表示兩者之間有很強(qiáng)的線性關(guān)系。#示例代碼：使用Python計算相關(guān)系數(shù)

importnumpyasnp

fromscipy.statsimportpearsonr

#假設(shè)我們有實驗數(shù)據(jù)和數(shù)值解

experimental_data=np.array([1.0,2.0,3.0])

numerical_solution=np.array([1.01,2.01,3.01])

#計算相關(guān)系數(shù)

correlation_coefficient,_=pearsonr(experimental_data,numerical_solution)

#輸出相關(guān)系數(shù)

print("相關(guān)系數(shù):",correlation_coefficient)通過上述步驟，可以確保數(shù)值解的準(zhǔn)確性和可靠性，以及與實驗數(shù)據(jù)的一致性，從而對空氣動力學(xué)方程的能量方程的數(shù)值解法有更深入的理解和應(yīng)用。9湍流能量方程的數(shù)值解法9.1湍流能量方程簡介湍流能量方程是描述湍流中能量傳輸和轉(zhuǎn)換的數(shù)學(xué)表達(dá)式。在湍流中，能量不僅在平均流場中傳輸，還通過湍流脈動在不同尺度上進(jìn)行交換。能量方程通常包括動能和內(nèi)能的傳輸，其中內(nèi)能傳輸方程尤為重要，因為它直接關(guān)聯(lián)到溫度分布和熱傳遞。9.1.1內(nèi)能傳輸方程內(nèi)能傳輸方程可以表示為：ρ其中，ρ是流體密度，Cp是定壓比熱，T是溫度，u是流體速度向量，k是熱導(dǎo)率，?9.2數(shù)值解法9.2.1離散化在數(shù)值模擬中，連續(xù)的能量方程需要被離散化，以便在計算機(jī)上求解。常用的離散化方法包括有限差分法、有限體積法和有限元法。這里，我們使用有限體積法進(jìn)行說明。9.2.1.1有限體積法有限體積法將計算域劃分為一系列控制體積，然后在每個控制體積上應(yīng)用守恒定律。對于能量方程，這意味著在每個控制體積上積分能量方程，得到離散形式的能量守恒方程。9.2.2時間離散化時間離散化通常采用顯式或隱式方法。顯式方法簡單直觀，但可能受到時間步長的限制；隱式方法則更為穩(wěn)定，但需要求解線性方程組。9.2.2.1顯式歐拉法顯式歐拉法是一種簡單的時間離散化方法，適用于初學(xué)者理解。它基于當(dāng)前時間步的信息預(yù)測下一個時間步的狀態(tài)。T9.2.2.2隱式歐拉法隱式歐拉法考慮了未來時間步的信息，因此更為穩(wěn)定。T9.2.3空間離散化空間離散化通常采用中心差分或上風(fēng)差分等方法。中心差分在平滑流場中效果較好，而上風(fēng)差分則適用于處理對流主導(dǎo)的流場。9.2.3.1中心差分中心差分是一種二階精度的空間離散化方法，適用于內(nèi)部點的計算。?9.2.3.2上風(fēng)差分上風(fēng)差分是一種一階精度的方法，適用于處理對流項。u9.2.4數(shù)值求解流程初始化：設(shè)定初始條件和邊界條件。離散化：將能量方程離散化。迭代求解：使用迭代方法求解離散方程，直到收斂。后處理：分析和可視化結(jié)果。9.2.5Python代碼示例下面是一個使用Python和NumPy庫求解一維穩(wěn)態(tài)能量方程的簡單示例。假設(shè)流體在管道中流動，管道兩端有溫度差，且沒有熱源。importnumpyasnp

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

L=1.0#管道長度

N=100#網(wǎng)格點數(shù)

rho=1.0#密度

Cp=1.0#定壓比熱

k=1.0#熱導(dǎo)率

T_left=100.0#左端溫度

T_right=50.0#右端溫度

#網(wǎng)格生成

dx=L/(N-1)

x=np.linspace(0,L,N)

T=np.zeros(N)#溫度分布

#離散化

A=np.zeros((N,N))

b=np.zeros(N)

#內(nèi)部點

foriinrange(1,N-1):

A[i,i-1]=-k/dx**2

A[i,i]=2*k/dx**2

A[i,i+1]=-k/dx**2

#邊界條件

A[0,0]=1

b[0]=T_left

A[N-1,N-1]=1

b[N-1]=T_right

#求解

T=np.linalg.solve(A,b)

#可視化結(jié)果

importmatplotlib.pyplotasplt

plt.plot(x,T)

plt.xlabel('位置(m)')

plt.ylabel('溫度(°C)')

plt.title('一維穩(wěn)態(tài)能量方程的數(shù)值解')

plt.show()9.3非定常流場的能量方程求解非定常流場的能量方程求解需要考慮時間變化，通常采用時間步進(jìn)方法，如上述的顯式或隱式歐拉法。9.3.1時間步進(jìn)時間步進(jìn)是通過迭代更新流場狀態(tài)來模擬非定常流場的演化。在每個時間步，需要求解離散化后的能量方程，更新溫度分布。9.3.2Python代碼示例下面是一個使用Python求解一維非定常能量方程的示例，采用顯式歐拉法進(jìn)行時間離散化。importnumpyasnp

importmatplotlib.pyplotasplt

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

L=1.0#管道長度

N=100#網(wǎng)格點數(shù)

rho=1.0#密度

Cp=1.0#定壓比熱

k=0.1#熱導(dǎo)率

T_left=100.0#左端溫度

T_right=50.0#右端溫度

dt=0.01#時間步長

t_end=1.0#模擬結(jié)束時間

#網(wǎng)格生成

dx=L/(N-1)

x=np.linspace(0,L,N)

T=np.zeros(N)#溫度分布

#時間步進(jìn)

t=0.0

whilet<t_end:

T_new=np.zeros(N)

foriinrange(1,N-1):

T_new[i]=T[i]+dt*(k/(rho*Cp))*(T[i+1]-2*T[i]+T[i-1])/dx**2

T_new[0]=T_left

T_new[N-1]=T_right

T=T_new

t+=dt

#可視化結(jié)果

plt.plot(x,T)

plt.xlabel('位置(m)')

plt.ylabel('溫度(°C)')

plt.