Modellazione del problema della conduzione del calore bidimensionale utilizzando Python
In questo tutorial vedremo come modellare l'equazione di conduzione del calore 2D utilizzando Python. Un'equazione 2D, stabile, di conduzione del calore con generazione di calore può essere scritta in coordinate cartesiane come segue −
$$\mathrm{\triangledown^{2} T \: + \: \frac{q_{g}}{k} \:=\: \frac{\parziale^{2 }T}{\parziale x^{2}} \: + \: \frac{\parziale^{2}T}{\parziale y^{2}} \: + \: \frac{q_{g}}{k} \:=\: 0 \:\:\dotso\dotso (1)}$$
Questo deve essere discretizzato per ottenere un'equazione alle differenze finite. Consideriamo una griglia rettangolare come quella mostrata di seguito.
L'indice 𝑖 corre verticalmente, cioè per riga, mentre l'indice 𝑗 corre orizzontalmente, cioè per colonna. Qualsiasi nodo interno (i,j) è circondato dai quattro nodi pertanto le informazioni su questo nodo possono essere trasferite solo da questi quattro nodi. Per i nodi al contorno, invece, l'informazione proverrà dalle condizioni al contorno.
Per discretizzare l’Eq. Verrà utilizzato 1 schema di differenza centrale che è accurato al 2° ordine. Inoltre, la spaziatura tra i punti della griglia nella direzione orizzontale e nella direzione verticale viene mantenuta la stessa, ovvero $\mathrm{\Delta x \:=\: \Delta y \:=\: \Delta} $. La formula della differenza centrale per la derivata del 2° ordine di una funzione di singola variabile è data come −
$$\mathrm{\frac{d^{2}f(x)}{dx^{2}} \:=\: \frac{f_{i-1} \: − \ : 2f_{i} \: + \: f_{i+1}}{\Delta x^{2}}\:\:\dotso\dotso (2)}$$
Allo stesso modo, la discretizzazione dell’Eq. 1 può essere scritto come −
$$\mathrm{\frac{T_{i-1, \: j} \: − \: 2T_{i,\: j} \: + \: T_{i+1, \: j}}{\Delta x^{2}} \: + \: \frac{T_{i, \: j-1} \: − \: 2T_{i,&bsol ;: j} \: + \: T_{i, \: j+1}}{\Delta 6^{2}} \: + \: \frac{q_{g}}{ k} \:=\: 0 \:\:\dotso\dotso (3)}$$
La forma semplificata dell’Eq. 3 sarà −
$$\mathrm{T_{i,\: j} \:=\: \frac{1}{4}(T_{i-1, \: j}\: + \: T_{i+1, \: j} \: + \: T_{i,\: j-1} \: + \: T_{i, \: j+1} \: + \: \Delta^{2} \frac{q_{g}}{k}) \:\:\dotso\dotso (4)}$$
L'equazione 4 è quella che deve essere risolta per ciascun nodo interno del dominio.
Per questo tipo di problemi l'informazione viaggia dal confine verso l'interno quindi la conoscenza delle condizioni al contorno è della massima importanza. Fondamentalmente in qualsiasi analisi numerica e analitica vengono spesso utilizzati tre tipi di condizioni al contorno: −
Dirichlet: In questo il valore della variabile dipendente è noto al confine. Per esempio
$$\mathrm{x \:=\: 0 \: \rightarrow \: T \:=\: 500 \: K }$$
-
Newmann: In questo la derivata della variabile dipendente è data come zero o funzione della variabile indipendente. Per esempio
$$\mathrm{x \:=\: 0 \: \rightarrow \: \frac{\partial T}{\partial x} \:=\: 0 }$$
Robin: In questo la derivata della variabile dipendente è una funzione della variabile dipendente stessa. Per esempio
$$\mathrm{− \: k \: (\frac{\partial T}{\partial x})_{x=0} \:=\: h(T \: − \: T_{\infty}) }$$
Algoritmo
Seleziona il dominio
Selezionare il numero di punti della griglia nelle direzioni x e y
Definisci $\mathrm{\Delta x}$e $\mathrm{\Delta y}$.
Definire la griglia.
Fornire l'ipotesi iniziale della temperatura.
Condizioni al contorno dell'offerta
-
Risolvi l'equazione 4 per i nodi interni
Valutare l'errore.
Se l'errore è superiore al criterio di convergenza, impostare il valore corrente di T come una nuova ipotesi e seguire i passaggi 4-9. Altrimenti, si raggiunge la risposta e si traccia il risultato.
In questo tutorial, risolveremo due casi, vale a dire. con e senza generazione di calore. Inoltre, assumeremo condizioni al contorno di Dirichlet.
I problemi da risolvere sono i seguenti −
Problema 1
Problema 2
Per entrambi i problemi tutto rimarrà uguale tranne il parametro; generazione di calore. Questo è il motivo per cui in questo articolo viene mostrato il codice solo una volta.
Importing modules
from pylab import *
from numpy import*
# Defining thermal properties
# case 1
# qg=0
# case 2
qg= 1000000
k=100
# Define domain
ℓx=1.0
ℓy=1.0
# Number of grid points
nx=20
ny=20
# Grid spacing
Δx=ℓx/(nx-1)
Δy=ℓy/(ny-1)
Δ=Δx
# Grid generation
x=linspace(0,ℓx,nx)
y=linspace(0,ℓy,ny)
X,Y=meshgrid(x,flipud(y))
# Initial Guess
T=zeros(shape=(nx,ny))
# Boundary conditions
#left
T[:,0]=300
#right
T[:,-1]=300
#top
T[0,:]=800
#bottom
T[-1,:]=300
# Guess array for comparison
Tg=T.copy()
# Initial error or entry in loop
error=1
# Iteration counter
count=1
# Comparison loop
while error>1.E-5:
# Sweeping in the domain
for i in range(1,nx-1):
for j in range(1,ny-1):
T[i,j]=(T[i,j-1]+T[i,j+1]+T[i-1,j]+T[i+1,j]+Δ**2*qg/k)/4
# Evaluating and printing error
error=sqrt(sum(abs(T-Tg)**2))
figure(1,dpi=300)
if count%100==0:
print(error)
semilogy(count,error,'ko')
xlabel("Iteration Counter")
ylabel("Error")
savefig("Error.jpg")
# Updating guess for next iteration
Tg=T.copy()
# Incrementing counter
count+=
# Result Plotting (Temperature contour plot)
figure(3,dpi=300)
cp1=contourf(X,Y,T,10,cmap='jet')
colorbar()
cp1=contour(X,Y,T,10,colors='k')
clabel(cp1,inline=True, fontsize=10)
savefig("temp.jpg")
Grafici di errore e temperatura per il primo caso
Grafici di errore e temperatura per il secondo caso
Ora puoi osservare che l'effetto della generazione di calore ha portato ad un aumento della temperatura interna così tanto che l'impatto delle condizioni al contorno non può eguagliarlo. Nel primo problema, hai osservato che l'effetto inizia dall'alto e si sposta verso il basso e lateralmente. Nel secondo caso, tuttavia, le condizioni al contorno cercano di influenzare, ma a causa dell'elevata generazione di calore, il flusso di calore ora prende il comando dal centro e si sposta verso l'esterno. Tuttavia, poiché il problema deve corrispondere alla condizione al contorno, le temperature al contorno rimangono inalterate.
Conclusione
In questo tutorial, l'equazione della conduzione del calore 2D è stata modellata in Python. La metodologia delle differenze finite è stata presentata per risolvere due problemi; uno con e uno senza generazione di calore. È stato osservato che l'entità della generazione di calore considerata nel codice non ha un grande impatto sulla soluzione finale in stato stazionario.