第4章基于Python的网格生成实战

第4章：基于Python的网格生成实战

网格生成这件事，说实话，是CFD分析里最磨人的环节之一。我见过太多人花了两周调求解器参数，结果网格质量不行，算出来的东西根本没法看。我自己刚入行那会儿也吃过这个亏——辛辛苦苦算了一周，后处理一看，升力系数曲线跟锯齿似的，最后发现是网格在翼型尾缘那里扭曲得太厉害。

这一章咱们就聊聊怎么用Python搞定翼型网格。我会用Gmsh这个工具，它开源、轻量、Python接口也还算顺手。你想想看，如果能用脚本一键生成网格，那可比在GUI里点鼠标高效多了。

核心要点：网格质量直接决定CFD结果的可靠性。边界层网格尤其关键——它捕捉的是壁面附近的速度梯度，搞不好整个阻力计算就废了。

4.1 准备工作：安装与导入

先装好Gmsh的Python接口。我个人习惯用pip安装：

pip install gmsh

导入的时候注意，Gmsh的API设计得有点特别——它用面向过程的方式操作几何和网格。我刚开始用的时候也觉得别扭，但用顺手了其实挺直观。

import gmsh
import numpy as np

# 初始化
gmsh.initialize()
gmsh.model.add("airfoil_mesh")

小技巧：每次运行完记得调用 gmsh.finalize() 释放资源。我刚开始经常忘，结果内存泄漏搞得电脑卡死。

4.2 翼型几何定义

咱们用NACA0012对称翼型做例子。它的坐标可以用解析公式生成：

def naca0012_coords(n_points=200):
    """生成NACA0012翼型坐标"""
    beta = np.linspace(0, np.pi, n_points)
    x = 0.5 * (1 - np.cos(beta))  # 余弦分布，前缘后缘加密
    
    # NACA 4位数字翼型厚度分布
    t = 0.12  # 最大厚度12%
    yt = 5 * t * (0.2969*np.sqrt(x) - 0.1260*x - 0.3516*x**2 + 0.2843*x**3 - 0.1015*x**4)
    
    # 对称翼型，上下表面
    x_upper = x
    y_upper = yt
    x_lower = x
    y_lower = -yt
    
    return x_upper, y_upper, x_lower, y_lower

这里有个坑——为什么用余弦分布而不是均匀分布？因为翼型前缘曲率大，需要更多网格点来捕捉几何特征。均匀分布的话，前缘会变得很粗糙，算出来的压力分布会失真。我在做某型风机叶片优化时吃过这个亏，后来才改过来的。

4.3 在Gmsh中创建几何

有了坐标点，接下来把它们导入Gmsh。这里要注意顺序——必须按逆时针方向闭合曲线：

def create_airfoil_geometry(x_upper, y_upper, x_lower, y_lower):
    """在Gmsh中创建翼型几何"""
    # 添加上表面点
    upper_points = []
    for i in range(len(x_upper)):
        pid = gmsh.model.geo.addPoint(x_upper[i], y_upper[i], 0, 0.001)
        upper_points.append(pid)
    
    # 添加下表面点（注意顺序：从尾缘到前缘）
    lower_points = []
    for i in range(len(x_lower)-1, -1, -1):
        pid = gmsh.model.geo.addPoint(x_lower[i], y_lower[i], 0, 0.001)
        lower_points.append(pid)
    
    # 创建样条曲线
    upper_spline = gmsh.model.geo.addSpline(upper_points)
    lower_spline = gmsh.model.geo.addSpline(lower_points)
    
    # 创建闭合曲线环
    curve_loop = gmsh.model.geo.addCurveLoop([upper_spline, lower_spline])
    
    return curve_loop, upper_spline, lower_spline

注意：点坐标的容差设置很关键。我设的是0.001，如果翼型尺寸很小（比如弦长0.1m），这个值可能太大，会导致几何失真。建议根据弦长动态调整。

4.4 远场边界设定

远场边界一般取翼型弦长的10-20倍。太小了会干扰流场，太大了浪费计算资源。我一般取15倍弦长：

def create_farfield(chord_length=1.0, farfield_ratio=15):
    """创建远场边界"""
    radius = chord_length * farfield_ratio
    
    # 创建远场圆
    center = gmsh.model.geo.addPoint(0, 0, 0)
    p1 = gmsh.model.geo.addPoint(radius, 0, 0)
    p2 = gmsh.model.geo.addPoint(0, radius, 0)
    p3 = gmsh.model.geo.addPoint(-radius, 0, 0)
    p4 = gmsh.model.geo.addPoint(0, -radius, 0)
    
    # 创建圆弧
    arc1 = gmsh.model.geo.addCircleArc(p1, center, p2)
    arc2 = gmsh.model.geo.addCircleArc(p2, center, p3)
    arc3 = gmsh.model.geo.addCircleArc(p3, center, p4)
    arc4 = gmsh.model.geo.addCircleArc(p4, center, p1)
    
    farfield_loop = gmsh.model.geo.addCurveLoop([arc1, arc2, arc3, arc4])
    
    return farfield_loop

为什么用圆形而不是矩形？因为圆形远场边界对来流方向不敏感，不管攻角怎么变，边界条件施加起来都方便。矩形的话，你得根据来流方向调整边界位置，麻烦得很。

4.5 边界层网格设置

这才是重头戏。边界层网格的核心思路是：壁面附近用很薄的棱柱层，然后逐渐过渡到三角形网格。我一般用以下参数：

参数	推荐值	说明
第一层高度	1e-5 * 弦长	保证y+ ≈ 1
增长因子	1.15 - 1.2	太大会导致网格过渡突兀
边界层层数	15 - 25	根据雷诺数调整
总厚度	0.05 * 弦长	覆盖边界层厚度

在Gmsh中设置边界层网格，需要用field来控制网格尺寸：

def setup_boundary_layer():
    """设置边界层网格"""
    # 创建边界层场
    bl_field = gmsh.model.mesh.field.add("BoundaryLayer")
    gmsh.model.mesh.field.setNumber(bl_field, "Size", 1e-5)  # 第一层高度
    gmsh.model.mesh.field.setNumber(bl_field, "Ratio", 1.15)  # 增长因子
    gmsh.model.mesh.field.setNumber(bl_field, "Thickness", 0.05)  # 总厚度
    gmsh.model.mesh.field.setNumber(bl_field, "Quads", 1)  # 使用四边形
    
    # 指定边界层附着在翼型表面
    gmsh.model.mesh.field.setNumbers(bl_field, "CurvesList", [upper_spline, lower_spline])
    
    return bl_field

经验之谈：第一层高度怎么算？用这个公式：y_wall = y+ * mu / (rho * u_tau)。其中u_tau是摩擦速度，大概估算为来流速度的5%。如果雷诺数1e6，弦长1m，第一层高度大概在1e-5量级。

4.6 网格生成与质量控制

所有设置完成后，执行网格生成：

def generate_mesh():
    """生成网格"""
    # 设置全局网格尺寸
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), 0.01)
    
    # 生成二维网格
    gmsh.model.mesh.generate(2)
    
    # 优化网格质量
    gmsh.model.mesh.optimize("Netgen")
    
    # 输出网格
    gmsh.write("airfoil.msh")
    
    # 可视化
    gmsh.fltk.run()

网格生成后一定要检查质量。我常用的指标：

最小角度：三角形单元应大于20°，四边形大于45°
长宽比：边界层内可以到1000，但主流区最好小于10
正交性：壁面附近网格线应尽量垂直于壁面

避坑指南：我曾经遇到过一个案例，网格生成时没报错，但算到一半就发散。后来发现是翼型尾缘处有几个单元扭曲得厉害，长宽比超过5000。从那以后我养成了习惯——每次生成网格后先跑一遍质量检查脚本。

4.7 完整流程整合

把上面的函数串起来，就是一个完整的网格生成脚本。我个人习惯把它封装成一个类，方便复用：

class AirfoilMeshGenerator:
    def __init__(self, chord=1.0, farfield_ratio=15):
        self.chord = chord
        self.farfield_ratio = farfield_ratio
        gmsh.initialize()
        
    def generate(self):
        # 1. 生成翼型坐标
        x_up, y_up, x_low, y_low = naca0012_coords()
        
        # 2. 创建几何
        curve_loop, up, low = create_airfoil_geometry(x_up, y_up, x_low, y_low)
        ff_loop = create_farfield(self.chord, self.farfield_ratio)
        
        # 3. 创建面
        surface = gmsh.model.geo.addPlaneSurface([ff_loop, curve_loop])
        
        # 4. 设置边界层
        bl_field = setup_boundary_layer()
        
        # 5. 生成网格
        gmsh.model.mesh.generate(2)
        
        return gmsh.model

你看，整个流程其实就五步：几何定义、远场设定、边界层配置、网格生成、质量检查。每一步都有讲究，但用脚本实现后，改参数、跑批量都非常方便。

嗯，网格生成这块儿就聊到这儿。记住一个原则：网格质量永远比数量重要。一万个高质量的网格，比十万个烂网格强得多。你花在网格上的时间，最终都会在计算结果里体现出来。

📚 课程目录