Функция завихренности, дающая необычный результат ⇐ Python

[Anonymous]

Я работаю над задачей потока в полости, управляемой крышкой, используя дискретные версии уравнений Навье-Стокса для завихренности:
Уравнения Навье-Стокса в дискретной форме завихренности
По какой-то причине моя функция потока работает (по крайней мере, я так считаю), но мои графики завихренности — полная чушь. Я пробовал что-то менять, но безрезультатно.
Полный код:

Код: Выделить всё

import numpy as np
import matplotlib.pyplot as plt

#Initialising Parameters:
N = 30
L = 1.0
Re = 5
dt = 0.0001
Vwall = 1.0
h = L/(N-1)

#Initialising Arrays:
u = np.zeros((N, N))
w = np.zeros((N,N))
Vx = np.zeros((N, N))
Vy = np.zeros((N, N))
P = np.zeros((N, N))

#Update stream function using 1a (Using SOR):
def updateStreamfn(u, w, h, tolerance=1e-6, maxIter=1000, omega=1.5):
#SOR for solving Poisson Equation (iterating using while & for loops):
maxDiff = tolerance + 1
iter = 0

#New while loop replacing for loop:
while maxDiff > tolerance and iter <  maxIter:
uOld = np.copy(u)

for i in range(1, u.shape[0] - 1):
for j in range(1, u.shape[1] - 1):

#Update Equation (using discretised 1a):
u[i, j] = (1-omega) * u[i, j] + (omega * 0.25) *(
u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1] + (h**2 * w[i, j])
)

#Break at convergence after maximum difference value reached:
maxDiff = np.max(np.abs(u - uOld))
iter += 1

if iter >= maxIter:
print("Reached max iterations without convergence.")
break

return u

#Update vorticity on wall with eqn 10 for all boundaries:
def updateWallVorticity(w, u, Vwall, h):
N = w.shape[0]

#Top wall (y=max at moving lid)
for i in range(N):
w[i, N-1] = ((2 / h**2) * (u[i, N-2] - u[i, N-1]) + (2 * Vwall / h))

#Bottom wall (y=0)
for i in range(N):
w[i, 0] = (2 / h**2) * (u[i, 1] - u[i, 0])

#Left wall (x=0)
for j in range(N):
w[0, j] = (2 / h**2) * (u[1, j] - u[0, j])

#Right wall (x=max)
for j in range(N):
w[N-1, j] = (2 / h**2) * (u[N-2, j] - u[N-1, j])

#Update vorticity in interior with 1b:
def updateInteriorVorticity(w, u, dt, Re, h):
wNew = np.copy(w)

#For loop to update through iterations:
for i in range(1, w.shape[0] - 1):
for j in range(1, w.shape[1] - 1):

#Advection and diffusion terms:
advectionx = ((u[i, j+1] - u[i, j-1]) / (2 * h) * ((w[i+1, j] - w[i-1, j]) / (2 * h)))
advectiony = ((u[i+1, j] - u[i-1, j]) / (2 * h)) * ((w[i, j+1] - w[i, j-1]) / (2 * h))

diffusion = (1 / Re) * ((w[i+1, j] - 2 * w[i, j] + w[i-1, j]) / h**2 + (
w[i, j+1] - 2 * w[i, j] + w[i, j-1]) / h**2)

#Update formula:
wNew[i, j] = w[i, j] + dt * (-advectionx + advectiony + diffusion)

return wNew

#Task 4 plots:
def plotStreamfn(u, Vx, Vy, title="Stream Function (u)", time_step=None):
plt.figure(figsize=(8, 6))
plt.contourf(u, levels=50, cmap='viridis')
plt.colorbar(label='Stream Function Value')

#Plotting time steps:
if time_step is not None:
plt.title(f"{title} at Time Step {time_step}")
else:
plt.title(title)
plt.xlabel("Grid X")
plt.ylabel("Grid Y")

#Streamlines to indicate flow direction:
Y, X = np.mgrid[0:u.shape[0], 0:u.shape[1]]
plt.streamplot(X, Y, Vx, Vy, color='white', linewidth=0.5)
plt.show()

#Vorticity plot:
def plotVorticity(w, title="Vorticity (w)", time_step=None):
plt.figure(figsize=(8, 6))
contour = plt.contourf(w, levels=50, cmap='viridis')
plt.colorbar(contour, label='Vorticity Value')

#Overlay contour lines:
plt.contour(w, levels=10, colors='white', linewidths=0.5, linestyles='solid')

#Title and labels:
if time_step is not None:
plt.title(f"{title} at Time Step {time_step}")
else:
plt.title(title)
plt.xlabel("Grid X")
plt.ylabel("Grid Y")
plt.show()

#Main time-stepping loop:
numTimesteps = 100
visualInter = 10

for n in range(numTimesteps):
#Update stream function:
u = updateStreamfn(u, w, h)

#Update vorticity on boundaries:
updateWallVorticity(w, u, Vwall, h)

#Update vorticity in interior:
w = updateInteriorVorticity(w, u, dt, Re, h)

if n % visualInter == 0:
print(f"Time Step:  {n}")

#Calculating Vx and Vy:
Vx[1:-1, 1:-1] = (u[1:-1, 2:] - u[1:-1, :-2]) / (2 * h)
Vy[1:-1, 1:-1] = -(u[2:, 1:-1] - u[:-2, 1:-1]) / (2 * h)

#Plotting stream function with streamlines:
plotStreamfn(u, Vx, Vy, title=f"Stream Function (u)", time_step=n)

#Plotting the vorticity:
plotVorticity(w, title=f"Vorticity (w)", time_step=n)

Я думаю, проблема в следующем:

Код: Выделить всё

#Update vorticity on wall with eqn 10 for all boundaries:
def updateWallVorticity(w, u, Vwall, h):
N = w.shape[0]

#Top wall (y = max at moving lid)
for i in range(N):
w[i, N-1] = ((2 / h**2) * (u[i, N-2] - u[i, N-1]) + (2 * Vwall / h))

#Bottom wall (y = 0)
for i in range(N):
w[i, 0] = (2 / h**2) * (u[i, 1] - u[i, 0])

#Left wall (x = 0)
for j in range(N):
w[0, j] = (2 / h**2) * (u[1, j] - u[0, j])

#Right wall (x = max)
for j in range(N):
w[N-1, j] = (2 / h**2) * (u[N-2, j] - u[N-1, j])

#Update vorticity in interior with 1b:
def updateInteriorVorticity(w, u, dt, Re, h):
wNew = np.copy(w)

#For loop to update through iterations:
for i in range(1, w.shape[0] - 1):
for j in range(1, w.shape[1] - 1):

# Corrected advection and diffusion terms
advectionx = ((u[i, j+1] - u[i, j-1]) / (2 * h)) * ((w[i+1, j] - w[i-1, j]) / (2 * h))

advectiony = ((u[i+1, j] - u[i-1, j]) / (2 * h)) * ((w[i, j+1] - w[i, j-1]) / (2 * h))

diffusion = (1 / Re) * ((w[i+1, j] - 2 * w[i, j] + w[i-1, j]) / h**2 + (
w[i, j+1] - 2 * w[i, j] + w[i, j-1]) / h**2)

# Update formula
wNew[i, j] = w[i, j] + dt * (-advectionx + advectiony + diffusion)

return wNew

и я получаю:


то есть бессмысленный результат. Я также получаю пустые графики позже при большем количестве итераций.
Есть ли какие-либо ошибки, которые я по какой-то причине не вижу, мелкие детали или что-то серьезное, о чем я не знаю?

Подробнее здесь: https://stackoverflow.com/questions/791 ... iar-result