his is the paper's DOI: 10.1017/S0022112003003835, system is on page 190 and procedure is on page 191.
import numpy as np
import matplotlib.pyplot as plt
Parameters
R = 1000 # Reynolds number
Delta = 0.3 # Oscillation amplitude
L = 1.0 # Length of the domain
N = 101 # Number of spatial grid points
eta = np.linspace(0, L, N) # Spatial grid
h = eta[1] - eta[0] # Spatial step size
T_end = 10000.0 # Total simulation time
num_time_steps = 10000 # Number of time steps (increased for finer resolution)
dt = T_end / num_time_steps # Time step size
Wall motion functions
H = lambda t: 1 + Delta * np.cos(2 * t)
H_dot = lambda t: -2 * Delta * np.sin(2 * t)
Initialize F and G
F = np.zeros(N)
G = np.zeros(N)
Helper functions for spatial derivatives
def dFdx(F):
"""Compute first derivative using central differences."""
dF = np.zeros_like(F)
dF[1:-1] = (F[2:] - F[:-2]) / (2 * h)
return dF
def d2Fdx2(F):
"""Compute second derivative using central differences."""
d2F = np.zeros_like(F)
d2F[1:-1] = (F[2:] - 2 * F[1:-1] + F[:-2]) / h**2
return d2F
def spatial_operator_G(G):
"""Apply the operator L to G."""
d2G = d2Fdx2(G)
dG = dFdx(G)
L_G = d2G + dG / eta - G / eta**2
return L_G
Right-hand side for Gt
def compute_rhs(F, G, t):
"""Compute the RHS for F and G."""
dF = dFdx(F)
d2F = d2Fdx2(F)
dG = dFdx(G)
L_G = spatial_operator_G(G)
RHS_F = d2F + dF / eta - F / eta**2 + H(t)**2 * G
RHS_G = (
-H_dot(t) / H(t) * eta * dG
- F * dG / H(t)
+ dF * G / H(t)
+ 2 * F * G / (H(t) * eta)
+ L_G / (H(t)**2 * R)
)
return RHS_F, RHS_G
RK4 integration
F_sol = [] # Store solutions of F
G_sol = [] # Store solutions of G
times = [] # Store time steps
t = 0.0
for n in range(num_time_steps):
F_sol.append(F.copy())
G_sol.append(G.copy())
times.append(t)
# RK4 steps
k1_F, k1_G = compute_rhs(F, G, t)
k2_F, k2_G = compute_rhs(F + 0.5 * dt * k1_F, G + 0.5 * dt * k1_G, t + 0.5 * dt)
k3_F, k3_G = compute_rhs(F + 0.5 * dt * k2_F, G + 0.5 * dt * k2_G, t + 0.5 * dt)
k4_F, k4_G = compute_rhs(F + dt * k3_F, G + dt * k3_G, t + dt)
F += dt * (k1_F + 2 * k2_F + 2 * k3_F + k4_F) / 6
G += dt * (k1_G + 2 * k2_G + 2 * k3_G + k4_G) / 6
# Boundary conditions
F[0], G[0] = 0, 0 # Dirichlet at eta = 0
F[-1] = -H_dot(t + dt) # Dirichlet at eta = 1
F[-2] = F[-1] - h * H(t + dt) # Neumann at eta = 1
# Increment time
t += dt
Convert solutions to arrays
F_sol = np.array(F_sol)
G_sol = np.array(G_sol)
times = np.array(times)
Plot results for the time range of t = 8100 to t = 8200
time_range = np.where((times >= 8100) & (times <= 8200))
plt.figure(figsize=(10, 6))
plt.plot(times[time_range], np.squeeze(F_sol[time_range, -1]), label='F at eta=1/4')
plt.plot(times[time_range], np.squeeze(G_sol[time_range, -1]), label='G at eta=1/4')
plt.xlabel('Time t')
plt.ylabel('Solution values')
plt.title('Time evolution of F and G from t = 8100 to t = 8200 at eta=1/4')
plt.legend()
plt.grid()
plt.show()
