Over the last few lectures we have shown how to model various kinds of astrophysical problems with fluid equations and finally getting out the famed Navier Stokes Equation. Sadly we have yet to tame such an equation analytically, and so until someone does that, we can instead solve it numerically and beat the fluid equations with compute alone. To get started let us introduce ourselves to how numerical simulatiors / integrators work, their draw backs and design. We start with the basics.
Aside from errors in programming or discretizing the physical model, there are 2 types of errors in computations:
Reason: computers are machines, with inherent limitations.
Under what circumstances is the following possible? $$(x+y)+z \neq x + (y + z)$$
Let's try it in Python:
x = 1e16 # this will actually work
y = -1e16 # need to decrease
z = 1.
print("(x+y)+z = {}".format((x+y)+z))
print("x+(y+z) = {}".format(x+(y+z)))
(x+y)+z = 1.0 x+(y+z) = 0.0
Lmao, this is not right, what happened? Round-off error!
Consider this system representing phasor rotation in the complex plane: $$\dot Z = i\omega Z, \quad\text{given}\quad Z_0 = Z(t=0).$$
Analytical solution is $Z(t) = Z_0 \exp(i\omega t).$
How can we solve it numerically?
Taylor expansion: $$\dot Z(t) = \frac{Z(t+\Delta t)-Z(t)}{\Delta t} + H.O.T. = i\omega Z(t).$$
Suggests algorithm:
Let's code it up and see
# Unstable solution to dZ/dt
from numpy import pi
# initialize Z0, omega and dt below;
Z = complex(1., 0.)
nsteps = 200
omega = 1.
dt = 2*pi/(nsteps*omega)
# print initial value
print('t = {0}, Z = {1}, |Z| = {2}'.format(0., Z, abs(Z)))
# rotate by 2pi
for k in range(nsteps):
Z = (1+1j*omega*dt) * Z
#print final value
print('t = {0}, Z = {1}, |Z| = {2}'.format(k*dt, Z, abs(Z)))
#print('t = {0:.3f}'.format(k*dt))
#print('Z = {0:.3f}'.format(Z))
print('|Z| = {0:.3f}'.format(abs(Z)))
t = 0.0, Z = (1+0j), |Z| = 1.0 t = 6.2517693806436885, Z = (1.1036746878104917-0.002280042726220316j), |Z| = 1.1036770429380236 |Z| = 1.104
What happened? What is the problem? Why did it happen?
Errors propagate statistically like they do in experimental physics.
These errors in differences become very important when taking numerical derivatives: $$\frac{df}{dt} \approx \frac{f_{i+1} - f_i}{\Delta t} = \text{danger zone}$$ as we will see later.
Never, ever. Never ever ever. Do something like this:
if 7./3. - 4./3. - 1. == 0:
print('7/3 - 4/3 - 1 == 0')
else:
print('7/3 - 4/3 - 1 != 0')
7/3 - 4/3 - 1 != 0
Instead:
delta = 1e-15
if abs(7./3. - 4./3. - 1.) < delta:
print('7/3 - 4/3 - 1 == 0 (or close enough anyway...)')
else:
print('7/3 - 4/3 - 1 != 0')
7/3 - 4/3 - 1 == 0 (or close enough anyway...)
Before we get to PDE's let us crack the ODE problem first. With two methods the Euler Method and the Runge-Kutta Methods of integrations.
$\quad x_{i} = x_{i-1} + hf(x_i) $
Consider the simple problem $$\frac{d^2x}{dt^2} + w^2 x = 0$$ We can solve this analytically and obtain the solution set of $x(t) = C_0e^{iw t}$ To numerically integrate this we have the simple "forward" scheme. We define a velocity $v_i$ and displacement $x_i$ and write the second order equation as 2 first order equations coupled together. $$x_{i+1} = x_i + \Delta t v_i $$ $$v_{i+1} = v_i - \Delta t w^2 x_i$$ We call this the forward Euler method becuase of the order at which we updated the equations. First position then momenta.
import numpy as np
import matplotlib.pyplot as plt
# Parameters
omega = 1.0 # Angular frequency
dt = 0.01 # Time step
T = 10 # Total time
# Time array
time = np.arange(0, T, dt)
# Initial conditions
x0 = 1.0 # Initial displacement
v0 = 0.0 # Initial velocity
# Arrays to store results
x = np.zeros(len(time))
v = np.zeros(len(time))
# Initial values
x[0] = x0
v[0] = v0
# Forward Euler loop
for i in range(len(time) - 1):
x[i + 1] = x[i] + dt * v[i]
v[i + 1] = v[i] - dt * omega**2 * x[i]
# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label="Displacement (x)")
plt.plot(time, v, label="Velocity (v)")
plt.title("Simple Harmonic Oscillator: Forward Euler Method")
plt.xlabel("Time (t)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()
The forward Euler method is NOT energy conserving. We can check the plot in the bottom for this case.
E_t = []
E_k_total = []
E_p_total = []
for i, t in enumerate(time):
E_k = 0.5 * v[i]**2
E_p = 0.5 * x[i]**2
E_k_total.append(E_k)
E_p_total.append(E_p)
E_t.append(E_k +E_p)
plt.plot(time, E_t , label='total energy')
plt.plot(time, E_k_total, label="kinetic energy")
plt.plot(time, E_p_total, label='spring potential')
plt.legend()
plt.grid()
plt.show()
We can further see that this is because the forward Euler method produces a phase space that is NOT closed. This is
plt.plot(x,v)
plt.xlabel("position")
plt.ylabel("velocity")
plt.grid()
plt.show()
We consider the Euler integration but back subsituting the values of $x_i$ after updating $v_i$. Let us first write out the forward case $$x_{i+1} = x_i + \Delta t v_{i+1}$$ $$v_{i+1} = v_i - w^2 \Delta t x_{i+1} $$ But let us say we only have access to the $v_i$ how do we update $v_{i+1}$ before we compute it? Well we can just rearrange to solve for the expressions we want (right hand side to not have any forward terms). Sub eq2 into eq1 $$x_{i+1} = x_i + \Delta t [v_i - w^2 \Delta t x_{i+1}]$$ $$x_{i+1} = x_i + \Delta t v_i - w^2 \Delta t^2 x_{i+1}$$ $$x_{i+1} + w^2 \Delta t^2 x_{i+1} = x_i + \Delta t v_i $$ $$x_{i+1} (1+ w^2 \Delta t^2 ) = x_i + \Delta t v_i $$ $$\boxed{x_{i+1} = \frac{x_i + \Delta t v_i}{1+ w^2 \Delta t^2 }} $$ $$\boxed{v_{i+1} = v_i - w^2 \Delta t x_{i+1}} $$ Now this is something we can write in a algorithm! The back subsitution makes it a "backward" euler.
import numpy as np
import matplotlib.pyplot as plt
# Parameters
omega = 1.0 # Angular frequency
dt = 0.01 # Time step
T = 10 # Total time
# Time array
time = np.arange(0, T, dt)
# Initial conditions
x0 = 1.0 # Initial displacement
v0 = 0.0 # Initial velocity
# Arrays to store results
x = np.zeros(len(time))
v = np.zeros(len(time))
# Initial values
x[0] = x0
v[0] = v0
# Backward Euler loop
for i in range(len(time) - 1):
# Solve for v_{i+1} first
v_next = (v[i] - dt * omega**2 * x[i]) / (1 + dt**2 * omega**2)
# Solve for x_{i+1} using the updated v_{i+1}
x_next = x[i] + dt * v_next
# Update the arrays
v[i + 1] = v_next
x[i + 1] = x_next
# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label="Displacement (x)")
plt.plot(time, v, label="Velocity (v)")
plt.title("Simple Harmonic Oscillator: Backward Euler Method")
plt.xlabel("Time (t)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()
E_t = []
E_k_total = []
E_p_total = []
for i, t in enumerate(time):
E_k = 0.5 * v[i]**2
E_p = 0.5 * x[i]**2
E_k_total.append(E_k)
E_p_total.append(E_p)
E_t.append(E_k +E_p)
plt.plot(time, E_t , label='total energy')
plt.plot(time, E_k_total, label="kinetic energy")
plt.plot(time, E_p_total, label='spring potential')
plt.legend()
plt.grid()
plt.show()
Energy is still not conserved !..... this is also wrong!
Okay so none of these really worked that well... We notice that one over estimates and the other underestimates the total energy. Let us combine them together through a central difference technique. Or better called the Symplectic Euler method. The key is to calculate "half step" corrections.
$$\boxed{v_{i+1/2} = v_i - w^2 \frac{\Delta t}{2} x_{i}} $$Then using the half velocity correction we update the $x_i$ $$\boxed{x_{i+1} = x_i + \frac{\Delta t}{2} v_{i+1/2}} $$ Then we fully update the velocity $$\boxed{v_{i+1} = v_{i+1/2} - w^2 \frac{\Delta t}{2} x_{i+1}} $$ This is basically the combination of both the foreward and backward euler implementation! Here the the mention of symplectic means that this method preserves the symplectic structure of the system. This just means a phase space volume is preserved in a hamiltonian system. (i think)
import numpy as np
import matplotlib.pyplot as plt
# Parameters
omega = 1.0 # Angular frequency
dt = 0.01 # Time step
T = 10 # Total time
# Time array
time = np.arange(0, T, dt)
# Initial conditions
x0 = 1.0 # Initial displacement
v0 = 0.0 # Initial velocity
# Arrays to store results
x = np.zeros(len(time))
v = np.zeros(len(time))
# Initial values
x[0] = x0
v[0] = v0 - 0.5 * dt * omega**2 * x0 # Half-step velocity update
# Central difference loop
for i in range(len(time) - 1):
# Update position
x[i + 1] = x[i] + dt * v[i]
# Update velocity at the next half-step
v[i + 1] = v[i] - dt * omega**2 * x[i + 1]
# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label="Displacement (x)")
plt.plot(time, v, label="Velocity (v)")
plt.title("Simple Harmonic Oscillator: Central Difference Method")
plt.xlabel("Time (t)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()
E_t = []
E_k_total = []
E_p_total = []
for i, t in enumerate(time):
E_k = 0.5 * v[i]**2
E_p = 0.5 * x[i]**2
E_k_total.append(E_k)
E_p_total.append(E_p)
E_t.append(E_k +E_p)
plt.plot(time, E_t , label='total energy')
plt.plot(time, E_k_total, label="kinetic energy")
plt.plot(time, E_p_total, label='spring potential')
plt.legend()
plt.grid()
plt.show()
Energy is conserved! (pretty much at least)
The RK method looks at using the mid points to evaluate Eulers method. Intuitvely the Euler method is just piecing together slopes together to approximate the general curve and trajectory given IC. Instead a better way is estimate the next timestep is the compute the slope midway between two points and use that slope to update you to the next step.
$\quad\displaystyle x\left(t + \frac{h}2\right) \approx x(t) + \frac{h}2 f[x(t), t]$
$\displaystyle \Rightarrow \boxed{x(t+h) = x(t) + h f\! \left[x(t) + \frac{h}2 f[x(t), t], t+\frac{h}2\right]}$
Let us use a simple toy example $$\frac{dx}{dt} = -kx$$ And let us try to integrate out the solution
import numpy as np
import matplotlib.pyplot as plt
# Parameters
k = 1.0 # Decay constant
dt = 0.3 # Time step
T = 5.0 # Total time
# Time array
time = np.arange(0, T + dt, dt)
# Initial condition
y0 = 1.0 # Initial value of y
# Array to store results
y = np.zeros(len(time))
y[0] = y0
# RK2 loop
for i in range(len(time) - 1):
t_i = time[i]
y_i = y[i]
# Compute k1 and k2
k1 = -k * y_i
k2 = -k * (y_i + k1 * dt)
# Update y
y[i + 1] = y_i + dt * (k1 + k2) / 2 # RK2 update formula
# Analytical solution for comparison
y_analytical = y0 * np.exp(-k * time)
# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, y, label="RK2 Numerical Solution")
plt.plot(time, y_analytical, label="Analytical Solution", linestyle="dashed")
plt.title("Runge-Kutta 2nd Order Method for ODE: dy/dt = -ky")
plt.xlabel("Time (t)")
plt.ylabel("y(t)")
plt.legend()
plt.grid()
plt.show()
For even higher accuracy, consider using Runge-Kutta 4th Order (RK4), which is a widely used numerical integration method.
To get the 4th order expansion you can implement by hand but its annoying but the framework for doing so is the following: This is the common go-to integration scheme used!
Firstly, the classification of partial differential equations (PDEs) is critical for numerical solvers because different types of PDEs exhibit different stable conditions which is needed for accurate solutions. Misclassifying a PDE can lead to instability, inefficiency, or incorrect results. Some examples.
Solutions may develop discontinuities (shocks) even if initial conditions are smooth. Initial value and boundary value problems. Methods: Explicit methods like finite difference, finite volume, or finite element methods.
We will see what makes the Navier Stokes Eq solvers to be so difficult is because it combines the challenges of multiple classes of PDEs. This requires hybrid solvers that adapt to each behavior.
We will not have time to explore each of these. But we will explore 2 example to understand the general framework of a PDE solver. Before we approach a full NSE problem.
Example problem: Modeling Electric potential. We want to model the electric potential for an empty 2D box, 10cm x 10cm in size, where the top wall is held at $V$ = 1.0V and the other walls at 0V.
$$0 = \nabla^2 \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2},$$We have boundary conditions. We note that this is time independent. $$ \phi(y = 10) = 1.0 V$$ $$ \phi(y = 0) = \phi(x = 0) = \phi(x = 10) = 0$$
To solve this we will use a method called Jacobi Relaxation the update scheme is derived as follows:
Recall that finte element difference to computing derivative is
$$ \frac{\partial \phi(x,y)}{\partial x} \approx \frac{ \phi(x+\Delta x,y) - \phi(x,y)}{\Delta x}$$We can compute the approximate second order version of this $$ \frac{\partial^2 \phi(x,y)}{\partial x^2} \approx \frac{\frac{\partial}{\partial x} \phi(x+\Delta x,y) - \frac{\partial}{\partial x}\phi(x,y)}{\Delta x}$$ $$ \frac{\partial^2 \phi(x,y)}{\partial x^2} \approx \frac{\frac{ \phi(x+2\Delta x,y) - \phi(x+\Delta x,y)}{\Delta x} - \frac{ \phi(x+\Delta x,y) - \phi(x,y)}{\Delta x}}{\Delta x}$$ $$ \frac{\partial^2 \phi(x,y)}{\partial x^2} \approx \frac{\phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y)}{\Delta x^2}$$ We can repeat for the y version $$ \frac{\partial^2 \phi(x,y)}{\partial y^2} \approx \frac{\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y)}{\Delta y^2}$$ Together we have the following according the Laplace equation $$ \frac{\phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y)}{\Delta x^2}+ \frac{\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y)}{\Delta y^2} = 0$$ We assume that the spacing is the same $\Delta x = \Delta y$ same size grids. $$ \frac{\phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y) +\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y)}{\Delta x^2} = 0$$ $$ \phi(x+2\Delta x,y) - 2\phi(x+\Delta x,y) + \phi(x,y) +\phi(x ,y+2\Delta y) - 2\phi(x,y+\Delta y) + \phi(x,y) = 0$$
$$\phi(x+a, y) + \phi(x-a, y) + \phi(x, y+a) + \phi(x, y-a) - 4\phi(x,y) = 0$$$\phi_{new}(x, y) = \frac14\left[\phi(x+a, y) + \phi(x-a, y) + \phi(x, y+a) + \phi(x, y-a)\right].$
import numpy as np
from tqdm import tqdm
def initPhi(M):
phi = np.zeros((M+1,M+1), np.float32)
phi[0,:] = 1.0
return phi
def jacobi(target, M, phi):
delta = 1
phinew = np.empty([M+1,M+1], np.float32)
iterations = 0
while delta > target:
for i in range(M+1):
for j in range(M+1):
if (i==0 or i ==M or j==0 or j==M):
phinew[i,j] = phi[i,j]
else:
phinew[i,j] = (phi[i+1,j] + phi[i-1,j] + phi[i,j+1]+ phi[i,j-1])/4
delta = np.amax(abs(phi-phinew))
phi, phinew = phinew, phi
iterations += 1
return iterations, phi
numIterations, phiJacobi = jacobi(1e-4, 100, initPhi(100))
print(numIterations)
print(phiJacobi)
1909 [[1.0000000e+00 1.0000000e+00 1.0000000e+00 ... 1.0000000e+00 1.0000000e+00 1.0000000e+00] [0.0000000e+00 4.9966612e-01 6.9698507e-01 ... 6.9698513e-01 4.9966609e-01 0.0000000e+00] [0.0000000e+00 3.0168000e-01 4.9866584e-01 ... 4.9866581e-01 3.0168000e-01 0.0000000e+00] ... [0.0000000e+00 1.5559126e-05 3.1097243e-05 ... 3.1097246e-05 1.5559128e-05 0.0000000e+00] [0.0000000e+00 7.7394761e-06 1.5468486e-05 ... 1.5468488e-05 7.7394761e-06 0.0000000e+00] [0.0000000e+00 0.0000000e+00 0.0000000e+00 ... 0.0000000e+00 0.0000000e+00 0.0000000e+00]]
import matplotlib.pyplot as plt
plt.imshow(phiJacobi, cmap='hot')
<matplotlib.image.AxesImage at 0x10f590a50>
Recall the 1D wave equation: $$\frac{\partial^2\phi}{\partial t^2}=v^2\frac{\partial^2\phi}{\partial x^2}$$
Consider a piano string of length $L$, initially at rest. At time $t=0$ the string is struck by the piano hammer a distance $d$ from the end of from the string. The string vibrates as a result of being struck, except at the ends, $x=0$, and $x=L$, where it is held fixed.
Consider the case $v=100\mathrm{ms^{-1}}$, with the initial condition that $\phi(x)=0$ everywhere but the velocity $\psi(x)$ is nonzero, with profile \begin{equation} \psi(x) = C\frac{x(L-x)}{L^2}\exp\left[-\frac{(x-d)^2}{2\sigma^2}\right], \end{equation} where $L=1$m, $d=10$cm, $C=1\mathrm{ms^{-1}}$, and $\sigma=0.3$m.
Approach using the FTCS method
import numpy as np
import matplotlib.pyplot as plt
# Parameters
L = 1.0 # Length of the domain
T = 0.11 # Total simulation time
c = 100
# Discretization
h = 1e-6
Nt = int(T/h)
Nx = 2**7
dx = L/Nx
x = np.linspace(0, L, Nx)
u = np.zeros((Nx, Nt + 1))
u_new = np.zeros(Nx)
# Initial condition
d = 0.1
C = 0.001
sigma = 0.3
u0 =0
u1 = C*(x*(L-x))/(L**2) * np.exp( - ((x-d)**2/(2*sigma**2)))
u[:, 0] = u0
u[:, 1] = u0 + 0.5 * c * h * u1
plt.figure(figsize=(8, 6))
for n in range(1, Nt):
for i in range(1, Nx - 1):
u_new[i] = 2 * u[i, n] - u[i, n - 1] + (c**2) * (h**2) / (dx**2) * (u[i + 1, n] - 2 * u[i, n] + u[i - 1, n])
u_new[0] = 0
u_new[-1] = 0
u[:, n + 1] = u_new
if n*h in [0.006, 0.004, 0.002, 0.012, 0.100002] :
plt.plot(u_new)
plt.xlabel('Time')
plt.ylabel('Position')
plt.title('Solution of 1D Wave Equation using FTCS')
Text(0.5, 1.0, 'Solution of 1D Wave Equation using FTCS')
We recall that the full NSE looks like the following $$\frac{\partial \rho}{\partial t} + \nabla \rho v = 0 \tag{mass cons.}$$ $$\frac{\partial v}{\partial t} + v \cdot \nabla v = \frac{-1}{\rho}\nabla P - \nabla \Phi + \nu[\nabla^2 v + \frac{1}{3}\nabla(\nabla \cdot v)] \tag{mom. cons}$$
In our lid driven example, we assume a few things, namely the fluid is incompressible and so $\nabla \cdot u = 0 $ we also assume this is 2-D problem.
In a cavity the top of the fluid has some velocity but everything inside is initially static. The top moving fluid drags the fluid internally in the cavity which creates swirly patterns (intuitively speaking). Initially velocity and pressure have zero initial condition. At the walls of the cavity we enforce a homogeneous dirichlet boundary conditions everywhere except for velocity at top. It is driven by an external flow.
The primary challenge with these simulations is enforcing the incompressibility condition, which requires that the velocity field has zero divergence. Chorin’s Projection solves this by "projecting" the velocity field onto a divergence-free space. The method operates by splitting the computation of the velocity field into intermediate steps, bypassing the need to solve for pressure and velocity simultaneously loosely speaking.
Below is an implementation following code derived from video here. Credits to machine-learning-and-simulation youtube channel
N_POINTS = 41
DOMAIN_SIZE = 1.0
N_ITERATIONS = 500
TIME_STEP_LENGTH = 0.001
KINEMATIC_VISCOSITY = 0.1
DENSITY = 1.0
HORIZONTAL_VELOCITY_TOP = 1.0
N_PRESSURE_POISSON_ITERATIONS = 50
STABILITY_SAFETY_FACTOR = 0.5
element_length = DOMAIN_SIZE / (N_POINTS - 1)
x = np.linspace(0.0, DOMAIN_SIZE, N_POINTS)
y = np.linspace(0.0, DOMAIN_SIZE, N_POINTS)
X, Y = np.meshgrid(x, y)
u_prev = np.zeros_like(X)
v_prev = np.zeros_like(X)
p_prev = np.zeros_like(X)
def central_difference_x(f):
diff = np.zeros_like(f)
diff[1:-1, 1:-1] = (
f[1:-1, 2: ]
-
f[1:-1, 0:-2]
) / (
2 * element_length
)
return diff
def central_difference_y(f):
diff = np.zeros_like(f)
diff[1:-1, 1:-1] = (
f[2: , 1:-1]
-
f[0:-2, 1:-1]
) / (
2 * element_length
)
return diff
def laplace(f):
diff = np.zeros_like(f)
diff[1:-1, 1:-1] = (
f[1:-1, 0:-2]
+
f[0:-2, 1:-1]
-
4
*
f[1:-1, 1:-1]
+
f[1:-1, 2: ]
+
f[2: , 1:-1]
) / (
element_length**2
)
return diff
maximum_possible_time_step_length = (
0.5 * element_length**2 / KINEMATIC_VISCOSITY
)
if TIME_STEP_LENGTH > STABILITY_SAFETY_FACTOR * maximum_possible_time_step_length:
raise RuntimeError("Stability is not guarenteed")
for _ in tqdm(range(N_ITERATIONS)):
d_u_prev__d_x = central_difference_x(u_prev)
d_u_prev__d_y = central_difference_y(u_prev)
d_v_prev__d_x = central_difference_x(v_prev)
d_v_prev__d_y = central_difference_y(v_prev)
laplace__u_prev = laplace(u_prev)
laplace__v_prev = laplace(v_prev)
# Perform a tentative step by solving the momentum equation without the
# pressure gradient
u_tent = (
u_prev
+
TIME_STEP_LENGTH * (
-
(
u_prev * d_u_prev__d_x
+
v_prev * d_u_prev__d_y
)
+
KINEMATIC_VISCOSITY * laplace__u_prev
)
)
v_tent = (
v_prev
+
TIME_STEP_LENGTH * (
-
(
u_prev * d_v_prev__d_x
+
v_prev * d_v_prev__d_y
)
+
KINEMATIC_VISCOSITY * laplace__v_prev
)
)
# Velocity Boundary Conditions: Homogeneous Dirichlet BC everywhere
# except for the horizontal velocity at the top, which is prescribed
u_tent[0, :] = 0.0
u_tent[:, 0] = 0.0
u_tent[:, -1] = 0.0
u_tent[-1, :] = HORIZONTAL_VELOCITY_TOP
v_tent[0, :] = 0.0
v_tent[:, 0] = 0.0
v_tent[:, -1] = 0.0
v_tent[-1, :] = 0.0
d_u_tent__d_x = central_difference_x(u_tent)
d_v_tent__d_y = central_difference_y(v_tent)
# Compute a pressure correction by solving the pressure-poisson equation
rhs = (
DENSITY / TIME_STEP_LENGTH
*
(
d_u_tent__d_x
+
d_v_tent__d_y
)
)
for _ in range(N_PRESSURE_POISSON_ITERATIONS):
p_next = np.zeros_like(p_prev)
p_next[1:-1, 1:-1] = 1/4 * (
+
p_prev[1:-1, 0:-2]
+
p_prev[0:-2, 1:-1]
+
p_prev[1:-1, 2: ]
+
p_prev[2: , 1:-1]
-
element_length**2
*
rhs[1:-1, 1:-1]
)
# Pressure Boundary Conditions: Homogeneous Neumann Boundary
# Conditions everywhere except for the top, where it is a
# homogeneous Dirichlet BC
p_next[:, -1] = p_next[:, -2]
p_next[0, :] = p_next[1, :]
p_next[:, 0] = p_next[:, 1]
p_next[-1, :] = 0.0
p_prev = p_next
d_p_next__d_x = central_difference_x(p_next)
d_p_next__d_y = central_difference_y(p_next)
# Correct the velocities such that the fluid stays incompressible
u_next = (
u_tent
-
TIME_STEP_LENGTH / DENSITY
*
d_p_next__d_x
)
v_next = (
v_tent
-
TIME_STEP_LENGTH / DENSITY
*
d_p_next__d_y
)
# Velocity Boundary Conditions: Homogeneous Dirichlet BC everywhere
# except for the horizontal velocity at the top, which is prescribed
u_next[0, :] = 0.0
u_next[:, 0] = 0.0
u_next[:, -1] = 0.0
u_next[-1, :] = HORIZONTAL_VELOCITY_TOP
v_next[0, :] = 0.0
v_next[:, 0] = 0.0
v_next[:, -1] = 0.0
v_next[-1, :] = 0.0
# Advance in time
u_prev = u_next
v_prev = v_next
p_prev = p_next
# The [::2, ::2] selects only every second entry (less cluttering plot)
plt.figure()
plt.contourf(X[::2, ::2], Y[::2, ::2], p_next[::2, ::2], cmap="coolwarm")
plt.colorbar()
plt.quiver(X[::2, ::2], Y[::2, ::2], u_next[::2, ::2], v_next[::2, ::2])
# plt.streamplot(X[::2, ::2], Y[::2, ::2], u_next[::2, ::2], v_next[::2, ::2], color="black")
plt.xlim((0, 1))
plt.ylim((0, 1))
plt.show()
100%|███████████████████████████████████████| 500/500 [00:00<00:00, 1392.49it/s]