如何将NACA翼型集成到格子玻尔兹曼流体模拟的障碍物网格？

核心问题分析 #

需要将NACA翼型的浮点坐标转换为整数网格索引，标记为LBM模拟中的障碍物（grid数组设为True），同时修复原代码语法错误，确保翼型尺寸适配模拟网格。

关键解决步骤 #

1. 生成适配网格的翼型坐标：将翼型弦长设为x方向网格尺寸，缩放并居中y坐标以适配y方向网格。
2. 创建障碍物掩码：使用多边形路径判断网格点是否位于翼型内部，生成布尔型障碍物网格。
3. 修复语法错误：移除循环中的冗余字符，避免代码运行报错。

修改后的完整代码 #

import numpy as np
from matplotlib import pyplot
from matplotlib.path import Path

# 提取NACA翼型参数
digits = input("enter airfoil digits")                                      # 从输入提取参数
m = float(str(digits)[0]) / 100
p = float(str(digits)[1]) / 10
xx = float(str(digits)[2:4]) / 100

nx = 300  # x方向网格数
ny = 75   # y方向网格数
plot_every = 100

def camber_line(x, m, p, c):
    return np.where((x >= 0) & (x <= (c*p)),
                    m * (x / np.power(p,2)) * (2.0 * p - (x / c)),
                    m * ((c - x) / np.power(1-p,2)) * (1.0 + (x / c) - 2.0 * p ))

def dyc_over_dx(x, m, p, c):
    return np.where((x >= 0) & (x <= (c*p)),
                    ((2.0 * m) / np.power(p,2)) * (p - x / c),
                    ((2.0 * m ) / np.power(1-p,2)) * (p - x / c ))

def thickness(x, xx, c):                                                  # 使用NACA标准公式
    term1 =  0.2969 * (np.sqrt(x/c))
    term2 = -0.1260 * (x/c)
    term3 = -0.3516 * np.power(x/c,2)
    term4 =  0.2843 * np.power(x/c,3)
    term5 = -0.1015 * np.power(x/c,4)
    return 5 * xx * c * (term1 + term2 + term3 + term4 + term5)

def naca4(x, m, p, xx, c=1):
    dyc_dx = dyc_over_dx(x, m, p, c)
    th = np.arctan(dyc_dx)
    yt = thickness(x, xx, c)
    yc = camber_line(x, m, p, c)  
    return ((x - yt*np.sin(th), yc + yt*np.cos(th)), 
            (x + yt*np.sin(th), yc - yt*np.cos(th)))

def main():
    tau = 0.53        # 运动粘度相关参数
    nt = 30000        # 时间步数

    # 格子玻尔兹曼节点速度与权重
    nl = 9                                           # 节点数
    cxs = np.array([0, 0, 1, 1, 1, 0, -1,-1,-1])     # x方向离散速度
    cys = np.array([0, 1, 1, 0,-1,-1, -1, 0, 1])     # y方向离散速度
    weights = np.array([4/9, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36])  # 节点权重

    # 初始化分布函数
    F = np.ones((ny, nx, nl)) + 0.01 * np.random.randn(ny, nx, nl)   
    F[:, :, 3] = 2.3  # 初始向右流动

    # --------------------------
    # 生成翼型障碍物网格
    # --------------------------
    c = nx  # 翼型弦长匹配x网格尺寸
    x_airfoil = np.linspace(0, c, 1000)  # 生成足够多的翼型点保证平滑
    (upper_x, upper_y), (lower_x, lower_y) = naca4(x_airfoil, m, p, xx, c)

    # 缩放并居中y坐标以适配ny网格
    max_y = max(np.max(upper_y), np.abs(np.min(lower_y)))
    scale_factor = (ny * 0.8) / (2 * max_y)  # 上下留10%余量
    upper_y_scaled = upper_y * scale_factor + ny/2
    lower_y_scaled = lower_y * scale_factor + ny/2

    # 创建闭合翼型多边形
    poly_x = np.concatenate([upper_x, lower_x[::-1]])
    poly_y = np.concatenate([upper_y_scaled, lower_y_scaled[::-1]])

    # 生成网格坐标点
    grid_x, grid_y = np.meshgrid(np.arange(nx), np.arange(ny))
    points = np.column_stack((grid_x.flatten(), grid_y.flatten()))

    # 判断网格点是否在翼型内部
    path = Path(np.column_stack((poly_x, poly_y)))
    inside = path.contains_points(points)

    # 转换为障碍物网格
    grid = inside.reshape((ny, nx))

    # --------------------------
    # 主模拟循环
    # --------------------------
    for it in range(nt):
        print(it)

        # 边界条件：壁面反弹
        F[:, -1, [6, 7, 8]] = F[:, -2, [6, 7, 8]]  # 右壁
        F[:, 0, [2, 3, 4]] = F[:, 1, [2, 3, 4]]    # 左壁
        F[-1, :, [1, 2, 8]] = F[-2, :, [1, 2, 8]]  # 顶壁
        F[0, :, [4, 5, 6]] = F[1, :, [4, 5, 6]]    # 底壁

        # 分布函数迁移
        for i, cx, cy in zip(range(nl), cxs, cys):
            F[:, :, i] = np.roll(F[:, :, i], cx, axis=1)
            F[:, :, i] = np.roll(F[:, :, i], cy, axis=0)

        # 障碍物边界处理：速度反转
        bndryF = F[grid, :]
        bndryF = bndryF[:, [0, 5, 6, 7, 8, 1, 2, 3, 4]]

        # 计算流体变量
        rho = np.sum(F, 2)  # 密度
        ux = np.sum(F * cxs, 2) / rho  # x方向速度
        uy = np.sum(F * cys, 2) / rho  # y方向速度

        # 应用障碍物边界条件
        F[grid, :] = bndryF
        ux[grid] = 0
        uy[grid] = 0

        # 碰撞步骤：计算平衡分布并更新
        Feq = np.zeros(F.shape)
        for i, cx, cy, w in zip(range(nl), cxs, cys, weights):
            Feq[:, :, i] = rho * w * (
                1 + 3 * (cx*ux + cy*uy) + 9 * (cx*ux + cy*uy)**2 / 2 - 3 * (ux**2 + uy**2) / 2)
        F = F + -(1/tau) * (F - Feq)

        # 定期绘制速度场
        if it % plot_every == 0:
            pyplot.imshow(np.sqrt(ux**2 + uy**2), cmap='viridis')
            pyplot.pause(0.01)
            pyplot.cla()

if __name__ == "__main__":
    main()

代码说明 #

1. 翼型坐标生成：使用1000个点生成平滑翼型轮廓，确保障碍物边缘连续。
2. 坐标缩放：将翼型y坐标缩放至网格高度的80%并居中，避免触及上下壁面。
3. 障碍物掩码：利用matplotlib.path.Path判断网格点是否在翼型内部，生成布尔型grid数组。
4. 语法修复：移除原代码循环行的冗余字符m，确保代码正常运行。

内容的提问来源于stack exchange，提问作者harryb