简介资源为西北工业大学计算流体力学课程大作业的完整 Python 实现源码面向航空宇航、力学等专业本科生与研究生系统覆盖 O 型贴体网格生成、拉瓦尔喷管流动计算与 Burgers 方程数值求解三大经典算例可直接应对课程提交与实验报告撰写需求。压缩包共 53 个文件以三个 py 主程序为核心配套 dat 结果数据、lay 流场图层、nb 计算笔记与二十余张可视化 png 图像整体仅 1.5MB轻量易用。已有 468 人学习下载。各算例均附带输入参数与说明文档图像涵盖速度场、压力系数、马赫数分布以及不同粘性系数下的三维对比便于逐项复现、核对结果代码按算例分目录组织、主程序独立清晰适合快速理解有限差分、贴体网格与喷管流场的数值求解思路也是期末冲刺、备战课程答辩的高质量参考。1. 三个Python算例拼起的一门CFD课程西北工业大学的计算流体力学大作业直接给了三个独立的Python求解器burgers.py负责非线性模型方程ogrid.py做圆柱绕流的O型贴体网格与流场求解laval.py算拉瓦尔喷管的可压缩流动。每个算例都有input.txt作为参数入口输出统一为Tecplot的.dat和.lay格式工程里也预先渲染好了多张结果图包括不同粘性系数下的速度场、压力系数分布和激波位置。对正在补CFD数值格式或者需要快速复现出图的人来说这三个脚本比课堂讲义直观得多而且从Burgers方程到变截面超声速流动的难度梯度也正好对应课程考核的覆盖范围。2. Burgers方程先搞清数值耗散与粘性的区别2.1 方程与离散格式Burgers方程是CFD里最常用的非线性标量模型方程u_t u u_x μ u_xx。左边是对流项右边是粘性项。当μ0时方程退化为纯双曲型解在有限时间内可能发展出间断当μ0时方程属于抛物型间断被物理粘性抹平成有限宽度的过渡层。这个算例的价值在于把同一套显式格式放在不同μ下运行能直观区分物理粘性和数值耗散的区别。burgers.py采用的是时间前向欧拉、空间对流通量一阶迎风、扩散项中心差分的组合。这样一个一阶格式在CFL数小于1时能稳定运行且代码逻辑简单适合作为课程作业的起步算例。真正需要理解的是迎风差分的隐含耗散即使μ设为0.0数值解仍然有类似粘性的抹平效果这也是为什么mu0.0.dat中的间断不是一条真正的直线而是一条略有坡度的过渡带。2.2 从mu0.0到0.05能看出什么burgers/目录下三个数据文件mu0.0.dat、mu0.01.dat、mu0.05.dat分别对应三种粘性设置。第一种没有物理粘性波形保持陡峭后两种粘性依次增大波形明显变缓。同时目录下的2D与3D PNG图显示了不同时刻和不同μ值下的曲面形态方便直观对比。μ行为特征对应数据文件0.0纯对流间断被格式自带数值耗散抹平mu0.0.dat0.01物理粘性强于数值耗散过渡层较窄mu0.01.dat0.05粘性主导波形接近高斯型光滑曲线mu0.05.dat初场在input.txt里配置通常是方波、正弦波或双曲正切分布。如果改成更陡的初场μ0.0时的数值解会出现明显的非物理振荡因为一阶迎风在间断附近仍然会引入色散误差。把这些图的纵轴统一成无量纲速度后还能计算出过渡层宽度与μ的平方根关系这是一个很好的从作业转向研究的小实验。2.3 burgers.py的迭代骨架# burgers.py 核心迭代循环简化版 import numpy as np nx 201 # 网格点数 dx 2.0 / (nx - 1) # 计算域 [-1, 1]空间步长 dt 0.001 # 时间步长需满足 CFL 条件 mu 0.01 # 粘性系数实际从 input.txt 读取 u np.zeros(nx) # 数值解数组初场已在前面初始化 for n in range(nt): un u.copy() for i in range(1, nx - 1): # 对流项按当地速度方向取迎风差分 if un[i] 0: conv un[i] * (un[i] - un[i-1]) / dx else: conv un[i] * (un[i1] - un[i]) / dx # 扩散项二阶中心差分 diff mu * (un[i1] - 2 * un[i] un[i-1]) / dx**2 # 前向欧拉时间推进 u[i] un[i] - dt * conv dt * diff对流项的方向由当地速度u[i]的符号决定速度为正时从上游取点速度为负时从下游取点这是保证单调性的关键。扩散项采用中心差分系数就是物理粘性μ直接调成0.0、0.01或0.05会得到三种不同形态的解。dt设置过大会出现数值发散表现为解中出现NaN或非物理正负交替振荡因此input.txt里一般会把CFL数控制在0.5以下。如果最终目标是验证格式精度还应该增加网格数nx并同时减小dt观察误差随网格步长的变化斜率。提示如果发现mu0.0时解仍然被抹平不要怀疑程序错误这是迎风差分数值耗散造成的。想进一步削弱耗散可以改用MUSCL线性重构或通量限制器Burgers方程是测试这些高阶格式最廉价的平台。3. O-Grid网格生成与圆柱绕流流场求解3.1 为什么用O型贴体网格对圆柱绕流这类外流场问题直接用笛卡尔网格会在物面附近产生阶梯状边界很难保证法向分辨率。O型网格的思路是围绕圆柱生成一个环形计算域内边界是圆柱表面外边界是远场圆。网格沿周向均匀分布沿径向采用指数递增方式加密近壁区域这样既避免了极坐标原点处的奇异性又让圆柱表面法向网格获得足够高的分辨率。ogrid.py先生成网格坐标X、Y再在贴体网格上求解流函数或涡量流函数方程。3.2 网格坐标生成的数学形式一个可复现的O型网格生成方式是把圆柱半径r0作为起始半径远场半径rmax作为终止半径沿径向按等比数列分布网格点。周向则直接对角度做等距切分。下面的代码逻辑与ogrid.py中常见做法一致# ogrid.py 中 O 型网格生成的简化实现 import numpy as np r0 0.5 # 圆柱半径 rmax 10.0 # 远场半径 nr 120 # 径向网格数 ntheta 160 # 周向网格数 # 径向按照等比数列加密近壁网格 r r0 * (rmax / r0) ** (np.arange(nr) / (nr - 1)) theta np.linspace(0.0, 2 * np.pi, ntheta, endpointFalse) R, Theta np.meshgrid(r, theta, indexingij) X R * np.cos(Theta) Y R * np.sin(Theta)径向网格点分布在指数增长序列上近壁处r变化小边界层内能布置足够多的点远场网格变粗减少总计算量。周向点数直接决定圆柱表面压力系数Cp的光滑程度太密会增加求解时间太疏则会在驻点附近产生波动。对于课设级别nr取100到150、ntheta取120到200已经足够复现出光滑的Cp曲线。3.3 流函数与速度场计算在不可压理想流假设下圆柱绕流满足拉普拉斯方程工程上常用流函数ψ作为求解变量。ogrid.py内部通常使用逐点超松弛迭代SOR求解流函数得到ψ后再用数值差分计算两个速度分量u ∂ψ/∂y v -∂ψ/∂x然后由无量纲伯努利方程计算表面压力系数Cp 1 - (u² v²) / V∞²。这里的V∞是来流速度会在input.txt里设置。下面给出SOR迭代的骨架# 流函数逐点松弛迭代伪代码 psi np.zeros_like(X) for it in range(max_iter): err 0.0 for i in range(1, nr - 1): for j in range(1, ntheta - 1): # 实际代码中需要按曲线坐标系数加权 psi_new 0.25 * (psi[i1,j] psi[i-1,j] psi[i,j1] psi[i,j-1]) err max(err, abs(psi_new - psi[i,j])) psi[i,j] psi_new omega * (psi_new - psi[i,j])SOR松弛因子omega取值一般在1.5到1.8之间过小收敛慢过大会出现迭代抖动。err是残差通常以小于1e-6为收敛标准。这个简化伪代码没有包含贴体坐标变换中的雅可比系数实际程序里会把这些系数作为数组提前计算好否则光滑圆柱表面的解会叠加周期性误差。网格参数对结果的影响常见取值范围r0圆柱半径决定表面曲率0.5rmax远场半径越大越接近无界流8~15nr径向网格数决定边界层分辨率100~200ntheta周向网格数决定Cp曲线光滑度120~2003.4 不同mu下的尾涡结构目录下的2Dmu0.0.png和2Dmu0.05.png直观显示了粘性系数的影响μ0.0时流场前后对称圆柱上下表面速度分布一致没有尾迹μ增大后下游出现非对称分离区圆柱背风面的负压区变宽低速尾流区域明显。u.png和v.png分别单独展示了x方向与y方向速度分量用于判断边界层分离点和涡的位置。如果进一步增大μ分离区会继续扩大阻力系数也随之上升这个趋势与实验定性一致。该算例把网格生成、迭代求解、后处理三个环节串成一条完整的CFD流程很多同学在答辩时会被问到“Cp曲线为什么在0°和180°附近有波动”答案就藏在周向网格数选得过少上。4. 拉瓦尔喷管可压缩流动的守恒型求解4.1 控制方程与面积变化源项拉瓦尔喷管是典型的一维变截面可压缩流动问题。控制方程是带面积源项的一维欧拉方程守恒形式为∂Q/∂t ∂F/∂x S其中守恒变量 Q [ρ, ρu, e]通量 F [ρu, ρu²p, u(ep)]源项与截面积A的空间导数有关。laval.py采用有限体积法把喷管沿轴向划分成若干等距单元每个单元的体积和界面面积都由几何关系预先算好。面积变化率越快流动参数梯度越陡因此网格在喉道附近需要加密。4.2 数值通量的离散有限体积格式下最关键的是计算界面上的数值通量。常见做法是采用Roe格式或Lax-Friedrichs格式。这里以Roe格式作为示意# laval.py 中 Roe 界面通量的简化版本 def primitives(Q, gamma): rho Q[0] u Q[1] / rho p (gamma - 1.0) * (Q[2] - 0.5 * rho * u * u) return rho, u, p def compute_flux_roe(QL, QR, gamma): rhoL, uL, pL primitives(QL, gamma) rhoR, uR, pR primitives(QR, gamma) # Roe 平均 rho 0.5 * (rhoL rhoR) u 0.5 * (uL uR) p 0.5 * (pL pR) c np.sqrt(gamma * p / rho) # 界面通量为左右通量平均加上耗散项 FL np.array([rhoL * uL, rhoL*uL*uL pL, uL * (rhoL * uL**2/2 gamma*pL/(gamma-1))]) FR np.array([rhoR * uR, rhoR*uR*uR pR, uR * (rhoR * uR**2/2 gamma*pR/(gamma-1))]) F 0.5 * (FL FR) - 0.5 * abs(c) * (QR - QL) return F注意这里为了可读性省略了Roe平均的完整表达式和波强分裂只保留耗散项的核心结构。实际laval.py中通量计算会更完整包含激波和膨胀波捕捉能力。界面通量计算完后再用前向欧拉或Runge-Kutta时间推进更新单元平均值。对于拉瓦尔喷管喉道附近会从亚声速连续加速到超声速如果耗散不足喉道附近会出现非物理振荡如果耗散过大激波会被抹成多层的过渡带。4.3 边界条件与截面输出入口为亚声速需要给定滞止压力和滞止温度出口如果流动已经转为超声速则可以全部外推。程序内部会根据当地马赫数自动切换边界处理。计算结果按不同流向截面输出到文件中Cx0.05.dat、Cx0.15.dat、Cx0.2.dat分别是三个截面上的流动参数分布其中Cx0.15通常接近喉道位置压力梯度最大。输出文件截面位置典型观察点Cx0.05.dat收缩段亚声速加速过程Cx0.15.dat喉道附近马赫数接近1压力最低Cx0.2.dat扩张段超声速段激波是否出现precise.dat提供了该工况下的参考解可以直接与数值解对比。用Python读取这两个数据文件后计算最大误差能定量评估程序的三阶精度还是二阶精度。laval.py里还保存了.lay布局文件打开Tecplot后能一键复现压力云图网格precise.nb则是一个Mathematica笔记本里面一般是解析解的符号推导或一维等熵关系式。对照查看时要注意解析解假设等熵流动若计算中出现激波解析解与数值解的偏差会集中在激波前后几个网格点内。5. 复现、验证与后处理技巧5.1 统一读取input.txt参数整个zip包三个算例都带input.txt但键名不完全一样。一个通用做法是写一个不区分类型的加载函数def load_input(path): params {} with open(path, r, encodingutf-8) as f: for line in f: line line.strip() if not line or line.startswith(#): continue key, value line.split(, 1) params[key.strip()] float(value.strip()) return params这种解析方式能忽略注释行和空行键值对以分隔。load_input返回dict后用params.get(mu, 0.01)这样的方式读取具体参数避免不同算例键名差异导致的KeyError。实际作业中把三个算例的input.txt都跑一遍确认图谱与自带PNG一致后再逐个修改参数就能快速完成参数敏感性分析。5.2 从.dat到可视化的路径.dat文件是Tecplot或者通用文本格式通常第一行写标题第二行是变量名之后是数据列。用matplotlib直接读的时候跳过前两行即可def read_tecplot(path): with open(path, r) as f: lines f.readlines() rows [] for line in lines[2:]: if line.strip() : continue rows.append([float(x) for x in line.split()]) return np.array(rows)如果你只想快速复现论文级的图优先用Tecplot打开对应.lay文件布局文件里已经绑定好数据范围、配色和视角。如果不想装商业软件就把上面读出的数据用matplotlib的contourf画云图。需要注意坐标系顺序Tecplot保存数据的顺序可能是I循环在外也可能是J循环在外画图前先打印数组维度确定reshape的参数否则云图会出现全部颠倒的问题。5.3 用precise.nb或解析解做定量校验laval/下的precise.nb可以通过Wolfram Engine的命令行接口批量运算也可以导出成pdf后直接对照。选定一个截面把数值解与解析解的距离定义为L2范数观察网格加密一倍后误差是否下降约4倍如果是说明格式达到二阶精度。对于burgers.py还可以调整mu0.01时用精确解Crank-Nicolson作对照能发现显式欧拉格式在同样网格下的相位误差更大。最后一个小技巧在迭代过程中每100步输出一次当前残差观察残差是否单调下降如果残差进入平台期说明网格数或松弛因子需要调整而不是程序出错。本文还有配套的精品资源点击获取