This repository was archived by the owner on Jul 17, 2023. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathchannelflow_with_turbulence.py
175 lines (143 loc) · 7.13 KB
/
channelflow_with_turbulence.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
DYN_VISCOSITY_MOL = 1e-3 # Dynamic viscosity
RHO = 1000
KIN_VISCOSITY_MOL = DYN_VISCOSITY_MOL / RHO
VON_KARMAN = 0.41
# Spatial constants and variables
aspect_ratio = 25 # Aspect ratio between y and x direction
n_y = 30 # Points in y direction
n_x = (n_y - 1) * aspect_ratio + 1 # Points in x direction
height = 1
length = height * aspect_ratio
delta_x = length / (n_x - 1)
delta_y = height / (n_y - 1)
x_range = np.linspace(0, length, n_x)
y_range = np.linspace(0, height, n_y)
coord_x, coord_y = np.meshgrid(x_range, y_range)
y = np.linspace(0, height, n_y - 1)
# Temporal constants and variables
delta_t = 5e-5 # time step size
n_t = int(6e5) # number times steps
plot_frequency = 5000
# X velocity
u_inlet = 1
# Initial Conditions
u_prev = np.ones((n_y + 1, n_x)) * u_inlet
# u_prev = np.zeros_like(u_prev)
u_prev[0, :] = - u_prev[1, :]
u_prev[1:-1, 0] = u_inlet
u_prev[-1, :] = - u_prev[-2, :]
# Y velocity
# Initial Conditions
v_prev = np.zeros((n_y, n_x + 1))
v_prev[0, :] = - v_prev[1, :]
v_prev[-1, :] = - v_prev[-2, :]
# Pressure initialization
pressure_prev = np.zeros((n_y + 1, n_x + 1))
n_poisson_pressure = 20 # Pressure Poisson iterations
# Pre-allocate predictor velocity arrays
u_star = np.zeros_like(u_prev)
u_next = np.zeros_like(u_prev)
v_star = np.zeros_like(v_prev)
v_next = np.zeros_like(v_prev)
# Define function that check if y is at <10% or >90% of the width of the pipe,
# then linear function is used, otherwise a constant value is used.
region_function = np.ones(n_y + 1)
region_function[:int(.1 * n_y) + 1] = np.array([((n - 1) * height * delta_y)**2 for n in range(1, int(.1 * n_y) + 2)])
region_function[int(.9 * n_y):] = np.array([((n_y - n) * height * delta_y)**2 for n in range(int(.9 * n_y), int(n_y) + 1)])
region_function[int(.1 * n_y):int(.9 * n_y) + 1] = (height * 0.1) ** 2
region_function = region_function[:, np.newaxis] * VON_KARMAN**2
for iter in tqdm(range(n_t)):
eddy_viscosity = np.abs(u_prev[2:, 1:-1] - u_prev[:-2, 1:-1]) * region_function[1:-1, :] / delta_y
# x velocity (u)
diffusive_u = (KIN_VISCOSITY_MOL + eddy_viscosity) * \
(u_prev[1:-1, 2:] + \
u_prev[1:-1, :-2] + \
u_prev[2:, 1:-1] + \
u_prev[:-2, 1:-1] - \
4 * u_prev[1:-1, 1:-1]) / delta_y ** 2
convective_u = (u_prev[1:-1, 2:]**2 - u_prev[1:-1, :-2]**2) / 2 / delta_x + \
(v_prev[1:, 1:-2] + v_prev[1:, 2:-1] + v_prev[:-1, 1:-2] + v_prev[:-1, 2:-1]) / 4 * \
(u_prev[2:, 1:-1] - u_prev[:-2, 1:-1]) / 2 / delta_x
pressure_gradient_x = (pressure_prev[1:-1, 2:-1] - pressure_prev[1:-1, 1:-2]) / delta_x
u_star[1:-1, 1:-1] = u_prev[1:-1, 1:-1] + delta_t * (-pressure_gradient_x + diffusive_u - convective_u)
# Boundary conditions
u_star[1:-1, 0] = u_inlet
u_star[1:-1, -1] = u_star[1:-1, -2]
u_star[0, :] = - u_star[1, :]
u_star[-1, :] = - u_star[-2, :]
# y velocity (v)
diffusive_v = KIN_VISCOSITY_MOL * (v_prev[1:-1, 2:] + \
v_prev[1:-1, :-2] + \
v_prev[2:, 1:-1] + \
v_prev[:-2, 1:-1] - \
4 * v_prev[1:-1,1:-1]) / delta_y ** 2
convective_v = (v_prev[2:, 1:-1]**2 - v_prev[:-2, 1:-1]**2) / 2 / delta_x + \
(u_prev[2:-1, 1:] + u_prev[2:-1, :-1] + u_prev[1:-2, 1:] + u_prev[1:-2, :-1]) / 4 * \
(v_prev[1:-1, 2:] - v_prev[1:-1, :-2]) / 2 / delta_x
pressure_gradient_y = (pressure_prev[2:-1, 1:-1] - pressure_prev[1:-2, 1:-1]) / delta_x
v_star[1:-1, 1:-1] = v_prev[1:-1, 1:-1] + delta_t * (-pressure_gradient_y + diffusive_v - convective_v)
# Boundary Conditions
v_star[1:-1, 0] = - v_star[1:-1, 1]
v_star[1:-1, -1] = v_star[1:-1, -2]
v_star[0, :] = 0.0
v_star[-1, :] = 0.0
# Righthandside poisson pressure
poisson_pressure_rhs = (u_star[1:-1, 1:] - u_star[1:-1, :-1] + v_star[1:, 1:-1] - v_star[:-1, 1:-1]) / delta_x / delta_t
# Pressure correction
pressure_correction_prev = np.zeros_like(pressure_prev)
for _ in range(n_poisson_pressure):
pressure_correction_next = np.zeros_like(pressure_correction_prev)
pressure_correction_next[1:-1, 1:-1] = (pressure_correction_prev[1:-1, 2:] + \
pressure_correction_prev[1:-1, :-2] + \
pressure_correction_prev[2:,1:-1] + \
pressure_correction_prev[-2,1:-1] - \
poisson_pressure_rhs * delta_x**2) / 4
# Boundary Conditions, use Neumann every expect for outlet where we have Dirichlet
pressure_correction_next[1:-1, 0] = pressure_correction_next[1:-1, 1]
pressure_correction_next[1:-1, -1] = pressure_correction_next[1:-1, -2]
pressure_correction_next[0, :] = -pressure_correction_next[1, :]
pressure_correction_next[-1, :] = pressure_correction_next[-2, :]
# Advance
pressure_correction_prev = pressure_correction_next
# Update Pressure
pressure_next = pressure_prev + pressure_correction_next
# Incompresibility
pressure_correction_gradient_x = (pressure_correction_next[1:-1, 2:-1] - pressure_correction_next[1:-1, 1:-2]) / delta_x
pressure_correction_gradient_y = (pressure_correction_next[2:-1, 1:-1] - pressure_correction_next[1:-2, 1:-1]) / delta_x
u_next[1:-1, 1:-1] = u_star[1:-1, 1:-1] - delta_t * pressure_correction_gradient_x
v_next[1:-1, 1:-1] = v_star[1:-1, 1:-1] - delta_t * pressure_correction_gradient_y
# BC again
u_next[1:-1, 0] = u_inlet
inflow_flux = np.sum(u_next[1:-1, 0])
outflow_flux = np.sum(u_next[1:-1, -2])
u_next[1:-1, -1] = u_next[1:-1, -2] * inflow_flux / outflow_flux
u_next[0, :] = - u_next[1, :]
u_next[-1, :] = - u_next[-2, :]
v_next[1:-1, 0] = - v_next[1:-1, 1]
v_next[1:-1, -1] = v_next[1:-1, -2]
v_next[0, :] = 0.0
v_next[-1, :] = 0.0
# Advance
u_prev = u_next
v_prev = v_next
pressure_prev = pressure_next
# Visualize simulation
if iter % plot_frequency == 0:
u_center = (u_next[1:, :] + u_next[:-1, :]) / 2
v_center = (v_next[:, 1:] + v_next[:, :-1]) / 2
plt.figure(dpi=200)
plt.contourf(coord_x, coord_y, u_center, levels=10)
plt.colorbar()
plt.quiver(coord_x[:, ::6], coord_y[:, ::6], u_center[:, ::6], v_center[:, ::6], alpha=0.4)
# plt.plot(5 * delta_x + u_center[:, 5], coord_y[:, 5], linewidth=3)
# plt.plot(20 * delta_x + u_center[:, 5], coord_y[:, 20], linewidth=3)
# plt.plot(u_center[:, int(0.8 * n_x)], coord_y[:, int(0.8 * n_x)], linewidth=3)
plt.draw()
plt.pause(0.05)
plt.clf()
plt.show()
plt.title("U velocity profile at exit of channel.")
plt.plot(coord_y[:, -3], u_center[:, -3])