我正在尝试使用该方法模拟特斯拉阀门，从 此代码 开始。 问题在于遮罩、边界或反弹未正确应用。 如果我反转创建障碍物遮罩的条件，则流量似乎更多地在阀门内部流动。
图像：
模拟输出
使用相反遮罩的模拟|| |完整代码：
我不太理解所有数学，所以我不知道该尝试什么。
我尝试询问 ChatGPT，但它只是告诉我确保正确应用掩模，未正确应用边界条件，但它不知道如何修复它们，并检查初始条件。
import jax
from PIL import Image
import jax.numpy as jnp
import matplotlib.pyplot as plt
import cmasher as cmr
from tqdm import tqdm
N_ITERATIONS = 30_000
REYNOLDS_NUMBER = 80
N_POINTS_X = 1024
N_POINTS_Y = 301
MAX_HORIZONTAL_INFLOW_VELOCITY = 0.04
VISUALIZE = True
PLOT_EVERY_N_STEPS = 100
SKIP_FIRST_N_ITERATIONS = 0
r"""
LBM Grid: D2Q9
6 2 5
\ | /
3 - 0 - 1
/ | \
7 4 8
"""
N_DISCRETE_VELOCITIES = 9
LATTICE_VELOCITIES = jnp.array(
[
[0, 1, 0, -1, 0, 1, -1, -1, 1],
[0, 0, 1, 0, -1, 1, 1, -1, -1],
]
)
LATTICE_INDICES = jnp.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
OPPOSITE_LATTICE_INDICES = jnp.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
# fmt: off
LATTICE_WEIGHTS = jnp.array(
[
4/9, # Center Velocity [0,]
1/9, 1/9, 1/9, 1/9, # Axis-Aligned Velocities [1, 2, 3, 4]
1/36, 1/36, 1/36, 1/36 # 45 ° Velocities [5, 6, 7, 8]
]
)
# fmt: on
RIGHT_VELOCITIES = jnp.array([1, 5, 8])
UP_VELOCITIES = jnp.array([2, 5, 6])
LEFT_VELOCITIES = jnp.array([3, 6, 7])
DOWN_VELOCITIES = jnp.array([4, 7, 8])
PURE_VERTICAL_VELOCITIES = jnp.array([0, 2, 4])
PURE_HORIZONTAL_VELOCITIES = jnp.array([0, 1, 3])
def create_tesla_valve_mask_from_image(image_path, lattice_size):
image = Image.open(image_path).convert("L")
image_resized = image.resize(lattice_size, Image.BILINEAR)
image_array = jnp.array(image_resized)
mask = image_array > 128
return mask
image_path = "valve-image-mask.png"
tesla_valve_mask = create_tesla_valve_mask_from_image(
image_path, (N_POINTS_X, N_POINTS_Y)
)
def get_density(discrete_velocities):
"""
For a single cell, if the velocities are [f0, f1, ..., f8],
the density is equal to f0 + f1 + ... + f8
"""
density = jnp.sum(discrete_velocities, axis=-1)
return density
def get_macroscopic_velocities(discrete_velocities, density):
macroscopic_velocities = (
jnp.einsum(
"NMQ,dQ->NMd",
discrete_velocities,
LATTICE_VELOCITIES,
)
/ density[..., jnp.newaxis]
)
return macroscopic_velocities
def get_equilibrium_discrete_velocities(macroscopic_velocities, density):
projected_discrete_velocities = jnp.einsum(
"dQ,NMd->NMQ",
LATTICE_VELOCITIES,
macroscopic_velocities,
)
macroscopic_velocity_magnitude = jnp.linalg.norm(
macroscopic_velocities,
axis=-1,
ord=2,
)
equilibrium_discrete_velocities = (
density[..., jnp.newaxis]
* LATTICE_WEIGHTS[jnp.newaxis, jnp.newaxis, :]
* (
1
+ 3 * projected_discrete_velocities
+ 9 / 2 * projected_discrete_velocities**2
- 3 / 2 * macroscopic_velocity_magnitude[..., jnp.newaxis] ** 2
)
)
return equilibrium_discrete_velocities
def main():
jax.config.update("jax_enable_x64", True)
kinematic_viscosity = (
MAX_HORIZONTAL_INFLOW_VELOCITY * N_POINTS_Y
) / REYNOLDS_NUMBER
relaxation_omega = (1.0) / (3.0 * kinematic_viscosity + 0.5)
# Define a mesh
x = jnp.arange(N_POINTS_X)
y = jnp.arange(N_POINTS_Y)
X, Y = jnp.meshgrid(x, y, indexing="ij")
# Obstacle Mask: An array of the shape like X or Y, but contains True if the
# point belongs to the obstacle and False if not
# (calculates the distance between a point and the center)
obstacle_mask = tesla_valve_mask.T
velocity_profile = jnp.zeros((N_POINTS_X, N_POINTS_Y, 2))
velocity_profile = velocity_profile.at[:, :, 0].set(MAX_HORIZONTAL_INFLOW_VELOCITY)
@jax.jit
def update(discrete_velocities_prev):
# (1) Prescribe the outflow BC on the right boundary
discrete_velocities_prev = discrete_velocities_prev.at[
-1, :, LEFT_VELOCITIES
].set(discrete_velocities_prev[-2, :, LEFT_VELOCITIES])
# (2) Macroscopic Velocities
density_prev = get_density(discrete_velocities_prev)
macroscopic_velocities_prev = get_macroscopic_velocities(
discrete_velocities_prev,
density_prev,
)
# (3) Prescribe Inflow Dirichlet BC using Zou/He scheme
macroscopic_velocities_prev = macroscopic_velocities_prev.at[0, 1:-1, :].set(
velocity_profile[0, 1:-1, :]
)
density_prev = density_prev.at[0, :].set(
(
get_density(discrete_velocities_prev[0, :, PURE_VERTICAL_VELOCITIES].T)
+ 2 * get_density(discrete_velocities_prev[0, :, LEFT_VELOCITIES].T)
)
/ (1 - macroscopic_velocities_prev[0, :, 0])
)
# (4) Compute discrete Equilibria velocities
equilibrium_discrete_velocities = get_equilibrium_discrete_velocities(
macroscopic_velocities_prev,
density_prev,
)
# (3) Belongs to the Zou/He scheme
discrete_velocities_prev = discrete_velocities_prev.at[
0, :, RIGHT_VELOCITIES
].set(equilibrium_discrete_velocities[0, :, RIGHT_VELOCITIES])
# (5) Collide according to BGK
discrete_velocities_post_collision = (
discrete_velocities_prev
- relaxation_omega
* (discrete_velocities_prev - equilibrium_discrete_velocities)
)
# (6) Bounce-Back Boundary Conditions to enfore the no-slip
# if in the direction of the velocity th obstacle mask is positive, then reverse the velocity
for i in range(N_DISCRETE_VELOCITIES):
discrete_velocities_post_collision = discrete_velocities_post_collision.at[
obstacle_mask, LATTICE_INDICES[i]
].set(discrete_velocities_prev[obstacle_mask, OPPOSITE_LATTICE_INDICES[i]])
# (7) Stream alongside lattice velocities
# for each cell, moves the velocity to a neighbor cell in its direction
discrete_velocities_streamed = discrete_velocities_post_collision
for i in range(N_DISCRETE_VELOCITIES):
discrete_velocities_streamed = discrete_velocities_streamed.at[:, :, i].set(
jnp.roll(
jnp.roll(
discrete_velocities_post_collision[:, :, i],
LATTICE_VELOCITIES[0, i],
axis=0,
),
LATTICE_VELOCITIES[1, i],
axis=1,
)
)
return discrete_velocities_streamed
discrete_velocities_prev = get_equilibrium_discrete_velocities(
velocity_profile,
jnp.ones((N_POINTS_X, N_POINTS_Y)),
)
plt.style.use("dark_background")
plt.figure(figsize=(15, 6), dpi=100)
for iteration_index in tqdm(range(N_ITERATIONS)):
discrete_velocities_next = update(discrete_velocities_prev)
discrete_velocities_prev = discrete_velocities_next
if (
iteration_index % PLOT_EVERY_N_STEPS == 0
and VISUALIZE
and iteration_index > SKIP_FIRST_N_ITERATIONS
):
density = get_density(discrete_velocities_next)
macroscopic_velocities = get_macroscopic_velocities(
discrete_velocities_next,
density,
)
velocity_magnitude = jnp.linalg.norm(
macroscopic_velocities,
axis=-1,
ord=2,
)
d_u__d_x, d_u__d_y = jnp.gradient(macroscopic_velocities[..., 0])
d_v__d_x, d_v__d_y = jnp.gradient(macroscopic_velocities[..., 1])
curl = d_u__d_y - d_v__d_x
"""
# Velocity Magnitude Contour Plot in the top
plt.subplot(211)
plt.contourf(
X,
Y,
velocity_magnitude,
levels=50,
cmap=cmr.amber,
)
plt.colorbar().set_label("Velocity Magnitude")
plt.imshow(tesla_valve_mask, cmap="gray", alpha=0.5, origin="lower")
# Vorticity Magnitude Contour PLot in the bottom
plt.subplot(212)
plt.contourf(
X,
Y,
curl,
levels=50,
cmap=cmr.redshift,
vmin=-0.02,
vmax=0.02,
)
plt.colorbar().set_label("Vorticity Magnitude")
plt.imshow(tesla_valve_mask, cmap="gray", alpha=0.5, origin="lower")
plt.draw()
plt.pause(0.0001)
plt.clf()
"""
plt.subplot(1, 3, 1)
plt.title("Obstacle Mask")
plt.imshow(obstacle_mask.T, cmap="gray", origin="upper")
plt.subplot(1, 3, 2)
plt.title("Velocity Magnitude")
plt.imshow(velocity_magnitude, cmap="viridis", origin="upper")
plt.subplot(1, 3, 3)
plt.title("Vorticity")
plt.imshow(curl, cmap="RdBu", origin="upper")
plt.draw()
plt.pause(0.0001)
plt.clf()
# """
if VISUALIZE:
plt.show()
if __name__ == "__main__":
main()
I don't really understand all the math so I don't know what to try. I tried asking ChatGPT but it just tells me to ensure that the mask is applied correctly, the boundary conditions are not applied correctly but it doesn't know how to fix them, and to check the initial conditions.
The problem with your Lattice Boltzmann Method (LBM) simulation of a Tesla valve lies in the implementation of the boundary conditions, specifically the bounce-back rule at the obstacle interface. Let's break down the issue and how to address it:
Understanding the Issue
- Bounce-Back Rule Misinterpretation: The current code seems to misinterpret the bounce-back rule. It flips the velocity direction for all velocities at the obstacle boundary if the mask is positive in the velocity's direction. This is not entirely accurate.
- Incorrect Implementation: The bounce-back rule should apply to fluid particles that would cross the boundary into the obstacle after streaming.
Solution
-
Identify Boundary Cells: Before the streaming step, you need to identify the fluid cells that are adjacent to the obstacle boundary. These are cells where at least one of their neighboring cells in the direction of a lattice velocity is an obstacle cell.
-
Apply Bounce-Back After Streaming:
- Instead of applying the bounce-back rule before streaming, apply it after streaming to the velocities that have streamed from a fluid cell into an obstacle cell.
- For each velocity direction, find the cells where the fluid has streamed into an obstacle.
- Swap the distribution functions of these velocities with their opposite counterparts in the corresponding fluid cells.
Revised Code (Conceptual)
# ... (previous code) ...
# (7) Stream alongside lattice velocities
discrete_velocities_streamed = discrete_velocities_post_collision
for i in range(N_DISCRETE_VELOCITIES):
discrete_velocities_streamed = discrete_velocities_streamed.at[:, :, i].set(
jnp.roll(
jnp.roll(
discrete_velocities_post_collision[:, :, i],
LATTICE_VELOCITIES[0, i],
axis=0,
),
LATTICE_VELOCITIES[1, i],
axis=1,
)
)
# (8) Apply Bounce-Back after streaming
for i in range(N_DISCRETE_VELOCITIES):
# Find cells where fluid streamed into obstacle
boundary_cells = jnp.logical_and(
jnp.roll(obstacle_mask, -LATTICE_VELOCITIES[0, i], axis=0), # Shifted obstacle mask
jnp.roll(obstacle_mask, -LATTICE_VELOCITIES[1, i], axis=1) # to check streaming direction
)
# Apply bounce-back
discrete_velocities_streamed = discrete_velocities_streamed.at[boundary_cells, LATTICE_INDICES[i]].set(
discrete_velocities_streamed[boundary_cells, OPPOSITE_LATTICE_INDICES[i]]
)
return discrete_velocities_streamed
# ... (rest of the code) ...
Explanation
- We stream the velocities first.
-
We then iterate through each lattice velocity and identify the
boundary_cells. We do this by shifting the obstacle mask in the opposite direction of the lattice velocity. This shifted mask, when combined with the original mask, highlights the cells where fluid has flowed into an obstacle. -
Finally, we apply the bounce-back rule only on those
boundary_cellsby swapping the distribution functions.
Additional Notes
- Accuracy: This implementation assumes a full bounce-back scheme, which is first-order accurate. Consider using a higher-order boundary condition for improved accuracy.
-
Optimization:
The current implementation loops through all velocities. Optimizations using JAX's
jnp.whereor masking techniques can significantly improve performance.
This revised approach ensures that the bounce-back boundary condition is applied correctly at the fluid-obstacle interface, leading to a more accurate simulation of the Tesla valve.
标签：python,numpy,matrix,jax,fluid-dynamics From： 78770964