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

[Anonymous]

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

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

import numpy as np
import matplotlib.pyplot as plt

#Initialising Parameters:
N = 50
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))

#Setting Boundary Conditions (Velocity):
Vx[:, 0] = 0  #Left wall
Vx[:, -1] = 0  #Right wall
Vx[0, :] = 0  #Bottom wall
Vx[-1, :] = 1  #Top wall

Vy[:, 0] = 0  #Left wall
Vy[:, -1] = 0  #Right wall
Vy[0, :] = 0  #Bottom wall
Vy[-1, :] = 0  #Top wall

#Update stream function using 1a (Using SOR):
def updateStreamfn(u, w, h, tolerance=1e-6, maxIter=1000, omega=1.3):
#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, 0] = -((2 / h**2) * (u[i, 1] - u[i, 0]) + (2 * Vwall / h))

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

#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, tolerance=1e-6, maxIter=1000, omega=1.3):
wNew = np.copy(w)
maxDiff = tolerance + 1
iter_count = 0

while maxDiff > tolerance and iter_count <  maxIter:
wOld = np.copy(wNew)

#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] = (1 - omega) * w[i, j] + omega * (w[i, j] + dt * (-advectionx + advectiony + diffusion))

#Convergence check
maxDiff = np.max(np.abs(wNew - wOld))
iter_count += 1

if iter_count >= maxIter:
print("Vorticity update reached max iterations without full convergence.")
break

return wNew

#Task 4 plots:
def plotStreamfnFin(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()

def plotStreamfnInt(u, 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")

#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}")

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

#Plotting the vorticity:
plotVorticity(w, title=f"Vorticity (w)", 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 the results at the final timestep
plotStreamfnFin(u, Vx, Vy, title="Final Stream Function (u)")
plotVorticity(w, title="Final Vorticity (w)")

Теперь мои окончательные графики (поток и завихренность) выглядят так:
Функция потока с линиями тока
График завихренности
Честно говоря, я не уверен, что происходит с графиком завихренности. Он все еще не выглядит... как ожидалось по завихренности? Я ошибаюсь?

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