用SIMPLE算法求解自然对流方程

问：

在随附的代码中，我正在尝试使用 SIMPLE 算法模拟自然对流，并使用 Upwind 等收敛方法。我有两个问题：

该代码不会生成有意义的绘图，因为我注意到在生成的文本文件中，我得到了 NaN（非数字）值。

我想改进图，因为我想观察的格拉斯霍夫数变化是流线末端的对流单元。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# Parametri della griglia e delle simulazioni
nx, ny = 64, 64
lx, ly = 1.0, 1.0
Gr, Re, Pr = 10000.0, 1000.0, 0.71

# Parametri di simulazione
dx, dy, dt = lx / nx, ly / ny, 0.00005
rel, omega = 1.97, 1.7

# Inizializzazione dei campi
u, v = np.zeros((ny + 2, nx + 2)), np.zeros((ny + 2, nx + 2))
p = np.ones((ny + 2, nx + 2)) * 0.5
T = np.ones((ny + 2, nx + 2)) * 0.5
p1, p2, T1 = np.zeros((ny + 2, nx + 2)), np.zeros((ny + 2, nx + 2)), np.zeros((ny + 2, nx + 2))
px, py = np.zeros((ny + 2, nx + 2)), np.zeros((ny + 2, nx + 2))

# Inizializzazione delle posizioni iniziali delle particelle per le streamline
start_x, start_y = np.linspace(0.1, 0.9, 10), np.linspace(0.1, 0.9, 10)
x_coords, y_coords = start_x.copy(), start_y.copy()

# Imposta le condizioni iniziali della temperatura
T_initial = 293.0  # Temperatura iniziale uniforme (in Kelvin)
T_hot_wall = 298.0  # Temperatura della parete superiore (in Kelvin)
T_cold_wall = 288.0  # Temperatura della parete inferiore (in Kelvin)

# Inizializzazione della temperatura
T[:, :] = T_initial  # Inizializza la temperatura uniformemente

# Imposta le condizioni di temperatura alle pareti
T[0, :] = T_cold_wall  # Parete inferiore
T[-1, :] = T_hot_wall  # Parete superiore

# Inizializzazione delle velocità (vicine a zero)
u[:, :] = 1.0
v[:, :] = 0.5

# Altri parametri
i_dx, i_dy = 1.0 / dx, 1.0 / dy
i_dx2, i_dy2 = 1.0 / (dx * dx), 1.0 / (dy * dy)
i_pr, i_gr = 1.0 / Pr, 1.0 / Gr

def calculate_streamlines(u, v):
ny, nx = u.shape
X, Y = np.meshgrid(np.linspace(0, lx, nx), np.linspace(0, ly, ny))

stream_x, stream_y = np.meshgrid(np.linspace(0, lx, 100), np.linspace(0, ly, 100))
u_interp = griddata((X.flatten(), Y.flatten()), u.flatten(), (stream_x, stream_y), method='linear')
v_interp = griddata((X.flatten(), Y.flatten()), v.flatten(), (stream_x, stream_y), method='linear')

return u_interp, v_interp, stream_x, stream_y

# Crea una griglia regolare per le coordinate (X, Y) per calcolare le velocità 
interpolate
X_reg, Y_reg = np.meshgrid(np.linspace(0, lx, 100), np.linspace(0, ly, 100))

# Monitoraggio della convergenza
max_iterations = 100
convergence_threshold = 1e-4
results_interval = 25

# Inizializza le variabili precedenti per il monitoraggio della convergenza
p_previous = np.copy(p)
T_previous = np.copy(T)

# Iterazioni principali
for it in range(max_iterations):
print(f"Iterazione n. {it + 1}")

# Calcolo delle componenti temporanee delle velocità con Upwind
for i in range(1, nx):
    for j in range(1, ny + 1):
        vt = 0.5 * (v[j - 1, i] + v[j - 1, i + 1])
        vb = 0.5 * (v[j, i] + v[j, i + 1])
        ut = 0.5 * (u[j - 1, i] + u[j, i])
        ub = 0.5 * (u[j + 1, i] + u[j, i])
        ur = 0.5 * (u[j, i] + u[j, i + 1])
        ul = 0.5 * (u[j, i] + u[j, i - 1])

        a = -(ur**2 - ul**2) * i_dx - (ub * vb - ut * vt) * i_dy
        b = (u[j, i + 1] - 2 * u[j, i] + u[j, i - 1]) * i_dx2
        c = (u[j + 1, i] - 2 * u[j, i] + u[j - 1, i]) * i_dy2
        AA = -a + np.sqrt(i_gr) * (b + c)
        u[j, i] = u[j, i] + dt * (AA - (p[j, i + 1] - p[j, i]) * i_dx)

for i in range(1, nx + 1):
    for j in range(1, ny):
        ur = 0.5 * (u[j, i] + u[j + 1, i])
        ul = 0.5 * (u[j, i - 1] + u[j + 1, i - 1])
        vr = 0.5 * (v[j, i + 1] + v[j, i])
        vl = 0.5 * (v[j, i - 1] + v[j, i])
        vt = 0.5 * (v[j - 1, i] + v[j, i])
        vb = 0.5 * (v[j + 1, i] + v[j, i])

        a = -(vb**2 - vt**2) * i_dy - (vr * ur - vl * ul) * i_dx
        b = (v[j, i + 1] - 2 * v[j, i] + v[j, i - 1]) * i_dx2
        c = (v[j + 1, i] - 2 * v[j, i] + v[j - 1, i]) * i_dy2
        BB = -a + np.sqrt(i_gr) * (b + c)
        v[j, i] = v[j, i] + dt * (BB - (p[j + 1, i] - p[j, i]) * i_dy)

# Calcola le differenze
err_x = np.max(np.abs(u - p))
err_y = np.max(np.abs(v - p))

# Calcola la pressione su griglia sfalsata
for i in range(1, nx + 1):
    for j in range(1, ny + 1):
        px[j, i] = (p[j, i + 1] + p[j, i - 1]) * i_dx2
        py[j, i] = (p[j + 1, i] + p[j - 1, i]) * i_dy2
        p2[j, i] = (px[j, i] + py[j, i] - rel * (p1[j, i + 1] + p1[j, i - 1] + p1[j + 1, i] + p1[j - 1, i]) / (
                    i_dx2 + i_dy2)) / (2 * (i_dx2 + i_dy2))

# Aggiorna la pressione utilizzando il metodo SOR
for i in range(1, nx + 1):
    for j in range(1, ny + 1):
        p[j, i] = (1.0 - omega) * p1[j, i] + omega * p2[j, i]

# Calcola l'errore della pressione
error_pressure = np.max(np.abs(p - p_previous))

# Calcola la differenza di temperatura su griglia sfalsata
for i in range(1, nx + 1):
    for j in range(1, ny + 1):
        if i_dx != 0 and i_dy != 0:
            Tx = (T[j, i + 1] - T[j, i - 1]) * 0.5 * i_dx
            Ty = (T[j + 1, i] - T[j - 1, i]) * 0.5 * i_dy
        else:
            Tx = 0.0
            Ty = 0.0
            
        T2 = Tx * Tx + Ty * Ty
        T1[j, i] = T[j, i] + dt * (i_pr * (T[j, i + 1] - 2 * T[j, i] + T[j, i - 1]) * i_dx2 +
                                   i_pr * (T[j + 1, i] - 2 * T[j, i] + T[j - 1, i]) * i_dy2 -
                                   0.5 * (u[j, i] * Tx + v[j, i] * Ty) -
                                   0.5 * (Tx + Ty) * T2)

# Aggiorna la temperatura
T[:, :] = T1[:, :]

# Calcola l'errore della temperatura
error_temperature = np.max(np.abs(T - T_previous))

if error_pressure < convergence_threshold and error_temperature < convergence_threshold:
    print(f"Convergenza raggiunta alla iterazione {it + 1}")
    break

# Aggiorna le variabili precedenti per il monitoraggio della convergenza
p_previous = np.copy(p)
T_previous = np.copy(T)

# Visualizza i risultati ogni 'results_interval' iterazioni
if it % results_interval == 0:
    plt.figure(figsize=(12, 6))

    # Campo vettoriale delle velocità
    plt.subplot(2, 2, 1)
    X, Y = np.meshgrid(np.linspace(0, lx, nx + 2), np.linspace(0, ly, ny + 2))
    plt.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2], scale=20)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title(f'Campo di velocità (Iterazione {it + 1})')
    
    # Mappa di calore della pressione
    plt.subplot(2, 2, 2)
    plt.imshow(p, cmap='jet', extent=(0, lx, 0, ly), origin='lower', aspect='equal')
    plt.colorbar(label='Pressione')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title(f'Campo di Pressione (Iterazione {it + 1})')

    # Mappa di calore della temperatura
    plt.subplot(2, 2, 3)
    plt.imshow(T, cmap='inferno', extent=(0, lx, 0, ly), origin='lower', aspect='equal')
    plt.colorbar(label='Temperatura')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title(f'Campo di Temperatura (Iterazione {it + 1})')
    
    # Streamlines
    plt.subplot(2, 2, 4)
    u_interp, v_interp, stream_x, stream_y = calculate_streamlines(u, v)
    plt.streamplot(stream_x, stream_y, u_interp, v_interp, density=2, linewidth=1, color='black')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title(f'Streamlines (Iterazione {it + 1})')

    plt.tight_layout()
    plt.show()

# Salvataggio dei dati in file separati
np.savetxt('pressure.txt', p, fmt='%f')
np.savetxt('temperature.txt', T, fmt='%f')
np.savetxt('velocity_x.txt', u, fmt='%f')
np.savetxt('velocity_y.txt', v, fmt='%f')