简介库埃特流动是流体力学中经典的剪切流动问题描述两平行平板相对运动时流体层的速度分布。压缩包面向流体力学初学者、相关课程学生以及希望用Matlab开展数值模拟的开发者提供一个可直接运行的库埃特流动求解示例。包内共3个文件包括Matlab脚本.m、说明文档.txt和附加资源压缩包.zip整个资源仅1KB结构简洁。已有520人学习过该资源。其中couete.m用于实现控制方程的离散与迭代求解帮助理解纳维-斯托克斯方程的数值处理方式文本文件说明运行参数与边界条件设置zip内为补充数据或辅助文件。通过实际运行脚本读者可以掌握有限差分法等基础数值方法了解粘度、密度等物性参数及无滑移边界条件对结果的影响并学会根据稳定性和精度要求调整空间网格与时间步长。这个低门槛实践入口适合在课程实验或自主学习中搭配教材同步使用为进一步研究更复杂的剪切流动和流体动力学问题打下基础。1. 库埃特流不是凑数的算例它是检验 CFD 代码是否撒谎的试金石库埃特流Couette flow的解析解只有一条直线两平行平板相对运动中间黏性流体被带动速度剖面呈线性壁面剪切应力恒定。正因为简单它成了验证任何黏性流求解器的第一道关卡。但你要是真的动手跑一遍就会发现解析解只给一个方程数值解却要面对网格、边界、压力参考点、松弛因子和上千次迭代。很多工程师在库埃特流上翻车不是因为不懂物理而是因为求解器和边界条件在撒谎。这份资源把解析解、Python 有限差分代码、OpenFOAM 案例和常见的排错路径放在一起适合刚接触 CFD 的从业者也适合手头有几个求解器但不确定它们还能不能信的人。它能解决的就是一件事让算出来的库埃特流对得上教科书。2. 从解析解到离散格式为什么库埃特流是最难翻车的验证算例2.1 解析解线性速度剖面与剪切应力的逻辑闭环稳定的库埃特流中速度只沿垂直板面方向变化。设下板静止上板以速度 U 沿 x 方向运动板间距为 h则充分发展后的流向速度 u 满足u(y) U * y / h这里 y 从下板 0 到上板 h。物理含义很直接黏性剪切是唯一驱动力压力梯度为零对流项也为零所以速度剖面是纯线性插值。壁面剪切应力 tau_w mu * U / h恒定不变。这个解之所以是“闭环”是因为它把几何参数、运动参数和物性参数全部收敛到一个简洁表达式里任何一个变量改错最终数值解都会偏离这条直线。但这里有一个反直觉的点库埃特流在流向x方向没有梯度所以计算域沿流向做多长都不会改变结果。很多新手习惯性地把管道拉得很长担心入口段效应。实际上库埃特流的关键是上下壁面的相对运动而不是入口动量。你完全可以把流向计算域设成极薄的一片甚至用一个网格单元只要流向边界用周期性条件或对称条件解就不会变。解析解的另一层意义是给验证设定了一个绝对标准数值解必须逼近 u U*y/h壁面剪切应力误差必须在可接受范围内否则就是数值方法本身出了问题。2.2 控制方程与无量纲参数Re、压力梯度与边界谁在起作用从不可压缩 Navier-Stokes 方程出发x 方向动量方程是rho * (du/dt udu/dx vdu/dy) -dp/dx mu * (d²u/dx² d²u/dy²)库埃特流稳态下 v 0du/dx 0对流项整体为零。若压力梯度也为零方程直接退化成 mu * d²u/dy² 0积分两次就能得到线性速度剖面。这个推导过程在数值实现里对应的就是你必须在求解器中确保对流项不产生伪贡献压力梯度也不能被人为构造出来。雷诺数 Re U*h / nu 在这个问题里地位特殊。稳态线性解对 Re 完全不敏感因为对流项为零惯性项不起作用。但在瞬态启动过程中Re 决定了从上板启动到速度剖面完全发展的时间尺度也决定了如果初始场给得不好会不会在启动瞬间出现数值振荡。做瞬态模拟时Re 还影响最大允许时间步长因为显式格式的稳定性约束本质上和扩散系数 nu 以及网格尺寸 dy 绑定。无量纲化之后压力场不是完全无足轻重。虽然速度剖面线性但压力场可以整体浮动。如果求解器没有指定压力参考点压力矩阵奇异收敛残差会卡住。所以实际模拟中压力需要设定一个参考值比如 pRefCell 和 pRefValue。物理上压力绝对值没有意义但数值上必须钉住它才能让线性系统有唯一解。这也是很多人在库埃特流里第一次接触压力参考点的原因。下面是一个常用参数组合的参考表参数数值示例说明板间距 h0.01 m控制网格分辨率的关键几何尺寸上板速度 U0.1 m/s驱动速度决定剪切应力量级运动黏度 nu1e-5 m²/s常见于水或轻质油Re 约 100Re100层流范围无湍流问题计算域流向长度 Lx0.001 m可取得极短库埃特流无流向梯度2.3 离散格式选型有限差分、有限体积还是 LBM实现库埃特流并不是只有一条路。有限差分实现最简单因为问题退化成 y 方向的一维二阶常微分方程中心差分配合 Thomas 算法可以直接求解不需要迭代。有限体积是 OpenFOAM 的默认底盘好处是几何适应性更强但同样的物理问题要写进控制字典边界条件类型的选择直接影响求解成败。LBM格子玻尔兹曼方法也可以做边界处理用的是反弹格式速度剖面依然能恢复但网格步长和松弛时间 tau 的换算关系容易出错。我的建议是第一阶段用有限差分在 Python 里把物理逻辑跑通第二阶段再进 OpenFOAM 验证代码和边界条件设置。这样你手里有两个独立实现可以互相印证出问题时能快速区分是物理建模错误还是求解器配置错误。本资源的两条路线都给了完整案例下面先走有限差分这条路。3. 用 Python 从零实现库埃特流求解器三步逼近稳态3.1 构建网格与边界条件上下板速度怎么给一维有限差分模型不需要搞懂 OpenFOAM 那些边界字典只需要把网格点和边界条件写清楚。首先构建 y 方向均匀网格将下板放在索引 0上板放在索引 N-1。import numpy as np N 101 # 网格点数量包含上下两个边界点 h 0.01 # 板间距单位 m U 0.1 # 上板运动速度单位 m/s nu 1e-5 # 运动黏度单位 m^2/s Re U * h / nu # 雷诺数约 100 y np.linspace(0, h, N) dy y[1] - y[0] # 初始场直接用线性分布比全零场更快逼近稳态 u U * y / h # 边界条件无滑移 u[0] 0.0 u[-1] U这段代码的要点在初始场上。有人习惯把全场初始化成零然后靠时间推进慢慢发展出速度剖面。但库埃特流的稳态解就是线性分布用线性分布作为初始场等于把真解的一半提前交给了求解器能显著减少迭代次数。边界条件必须放在数组两端并且在每次迭代之后重新强制赋值否则内部的离散方程会覆盖边界值。参数上N 决定 dy 也就决定了数值误差。对于一维线性剖面101 个点已经足够输出非常漂亮的直线如果你想做网格收敛性分析N 可以从 21、41、81、161 这样依次翻倍。注意 Re 在这里只做物理判断不会改变稳态解但会影响显式时间推进的稳定性。3.2 时间推进与松弛显式、隐式还是定常迭代稳态库埃特流满足 d²u/dy² 0这可以直接写成三对角线性方程组。用 numpy 的linalg.solve一步到位比时间推进更干净也完全不涉及稳定性条件。# 构造三对角矩阵内部点用中心差分 A np.zeros((N, N)) b np.zeros(N) for i in range(1, N - 1): A[i, i - 1] 1.0 A[i, i] -2.0 A[i, i 1] 1.0 # 边界行直接赋值 A[0, 0] 1.0 b[0] 0.0 A[-1, -1] 1.0 b[-1] U u_soln np.linalg.solve(A, b)这里没有除以 dy²因为方程右侧为零除以一个常数不改变解。如果未来扩展成带源项的泊松方程就一定要补上 dy² 系数。边界行的处理很关键把对角线设为 1其他位置设为 0右端项直接放边界速度这样求解器不会把边界点当未知量处理。如果你想模拟上板突然启动的瞬态过程就不能再用这个定常矩阵需要时间推进。显式格式最简单但稳定性条件要求 nu * dt / dy² 0.5其中 dy 是网格间距dt 是时间步长。以 N101、h0.01 为例dy1e-4nu1e-5那么 dt 最多只能取约 0.5 秒实际取 0.1 秒才安全。隐式格式无条件稳定代价是要在每一时间步解一次三对角方程代码量稍大但不用担心步长上限。3.3 收敛判据和结果输出什么时候说算完了如果用直接解法不存在收敛问题。但如果你用迭代法比如 SOR 或者简单的迭代松弛就必须定义一个收敛判据。常见做法是用相邻两次迭代解之间的 L2 相对残差。def l2_residual(u_new, u_old): return np.linalg.norm(u_new - u_old) / np.linalg.norm(u_new) # 迭代示例 u_old np.zeros_like(u) residual 1.0 it 0 while residual 1e-10 and it 10000: # 这里替换成你的迭代更新逻辑 # u_new 迭代更新(u_old) residual l2_residual(u_new, u_old) u_old u_new it 1 if it % 100 0: print(fstep {it}, residual {residual:.2e})收敛阈值取多小取决于你的目标。如果只是看速度剖面形状1e-6 就够如果要提取壁面剪切应力做对比建议放到 1e-10。我习惯把阈值设得比实际需求高两个数量级因为迭代求解器的残差下降曲线经常会出现“伪收敛”也就是残差卡在某个平台期不再变化这时候改小阈值也没用必须回头检查边界条件是否被破坏。算完之后还要和解析解做一次误差分析u_analytic U * y / h error u_soln - u_analytic l2_error np.linalg.norm(error) / np.linalg.norm(u_analytic) print(fL2 error: {l2_error:.2e})这个 L2 误差不是随便看看的后续做网格收敛性时它会成为核心指标。建议把速度和误差一起画成两张图速度图看趋势误差图看分布。如果误差在近壁面明显增大说明边界处理有瑕疵如果误差在中间区域波动说明网格不够密或离散格式有问题。4. 把库埃特流跑进 OpenFOAM运动壁面、周期边界与层流求解器4.1 用 blockMesh 生成一个轻薄网格有限差分跑通以后OpenFOAM 的作用是验证边界条件和求解器在真实工程环境中的表现。首先要生成网格。由于库埃特流在流向无梯度计算域可以做得非常薄但为了让周期性边界正常建立联系x 方向至少要有 1 个网格。下面是 blockMeshDict 的核心设置。// blockMeshDict convertToMeters 1; vertices ( (0 0 0) (0.001 0 0) (0.001 0.01 0) (0 0.01 0) (0 0 0.001) (0.001 0 0.001) (0.001 0.01 0.001) (0 0.01 0.001) ); blocks ( hex (0 1 2 3 4 5 6 7) (1 40 1) simpleGrading (1 1 1) ); boundary ( lowerWall { type wall; faces ( (0 1 5 4) ); } upperWall { type wall; faces ( (3 2 6 7) ); } cyclicX { type cyclic; faces ( (0 3 7 4) (1 2 6 5) ); } frontBack { type empty; faces ( (0 4 7 3) (1 5 6 2) ); } );网格在 y 方向分成 40 份x 方向和 z 方向各 1 份这样整体只有 40 个六面体单元。这里的关键是cyclicX边界它把 x 方向的两个面变成一对周期面第一个 face 列表属于 x0 的面第二个 face 列表属于 x0.001 的面。OpenFOAM 要求两个 cycle face 完全匹配blockMesh 会自动建立变换关系。frontBack是 z 方向的面因为二维模拟只需要 x-y 平面所以设为empty。这个网格没有任何正交问题simpleGrading (1 1 1)表示均匀网格。如果你想测试近壁面分辨率对剪切应力的影响可以把(1 40 1)改为(1 40 1)加上simpleGrading (1 10 1)让网格向上下壁面加密。加密后壁面剪切应力的输出值会更接近解析解但要注意 y 方向网格比不能太极端否则非正交误差反而增大。4.2 初始场与运动壁面速度指定OpenFOAM 的初始场文件在0/目录下。速度场需要三个文件0/U、0/p有时候还要0/epsilon和0/k如果跑湍流模型。库埃特流是层流只需要 U 和 p。0/U文件内容/*--------------------------------*- C -*----------------------------------*/ FoamFile { version 2.0; format ascii; class volVectorField; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { lowerWall { type fixedValue; value uniform (0 0 0); } upperWall { type fixedValue; value uniform (0.1 0 0); } cyclicX { type cyclic; } frontBack { type empty; } }internalField给全零场可以但和 Python 脚本一样用线性分布初始化会更快收敛。OpenFOAM 里写非均匀初始场比较麻烦常见做法是用setFields或者写一个简单的解析函数文件但对于这个算例全零场下 PISO 求解器多迭代几步也就到了稳态不影响最终结果。上壁面的速度边界必须用fixedValue不能想当然地写movingWallVelocity。movingWallVelocity是给压力边界条件用的本质上是让壁面运动带来的通量在压力方程里被正确处理。如果你把它用在0/U的 upperWall 上OpenFOAM 运行几秒后往往会出现压力振荡或速度越界。压力边界0/p中 upperWall 和 lowerWall 都用zeroGradientcyclicX 用 cyclicfrontBack 用 emptyinternalField设 uniform 0。4.3 求解器与湍流模型为什么 icoFoam 比 simpleFoam 更合适OpenFOAM 里跑不可压缩层流找稳态解很多人第一反应用 simpleFoam。实际上 simpleFoam 是为稳态设计的但它在入口段问题上的收敛路径不如瞬态 PISO 光滑。库埃特流没有入口内部也没有回流区用 icoFoam 跑几百个时间步残差降到 1e-6 以下是很容易的。icoFoam 使用的是 PISO 算法对压力速度偶合的处理在低雷诺数下很稳。system/fvSolution中关键的 PISO 参数solvers { p { solver PCG; preconditioner DIC; tolerance 1e-8; relTol 0.01; } U { solver PBiCGStab; preconditioner DILU; tolerance 1e-8; relTol 0.01; } } PISO { nCorrectors 2; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; }nCorrectors设 2 是标准做法。设为 1 也能跑但压力场收敛速度会慢。nNonOrthogonalCorrectors在这种全正交六面体网格上完全不需要设 0 就对了。pRefCell指定压力参考点这里选最左下角的体心pRefValue设为 0。这一行是很多人忽视的缺了它压力矩阵奇异残差会永远卡在一个固定值上。system/fvScheme中速度离散格式建议用Gauss linear这是中心差分格式对库埃特流这种光滑问题精度足够。对流项如果非要用带限制器的格式反而会把曲线修圆。时间项用Euler足够CrankNicolson 虽然二阶但在这个问题上收益不大反而可能因为初始振荡引入额外误差。5. 避坑指南库埃特流模拟里的五个经典翻车现场5.1 速度剖面是二次曲线而不是直线对流项和压力梯度没被消掉现象输出的速度剖面在 y 方向明显弯曲看起来像抛物线和线性解析解完全对不上。原因有两种一是计算域流向太短入口和出口效应穿透整个计算域破坏了“流向无梯度”的前提二是周期性边界没有真正生效数值上在流向建立了一个虚假压力梯度。解决先把计算域沿流向压缩到几十个网格以内并检查cyclicX是否正确配对。在 Python 里对应的问题则是你写的离散方程混入了对流项比如给u*du/dx保留了上游迎风格式。解除对流项之后剖面立刻回到直线。5.2 壁面剪切力对不上解析解粗网格和梯度插值现象数值算出的壁面剪切应力 tau_w mu * du/dy 比解析解低 20% 以上且网格加密后误差显著缩小。原因均匀粗网格下中心差分算出的壁面梯度是从第一层网格中心到边界点的差商它只具有一阶精度即使整体格式是二阶。解决在 OpenFOAM 中给 blockMesh 加 near-wall 加密或者在 Python 中用一次性单边差分公式提取壁面梯度。我一般会在后处理脚本里专门写一个函数用二阶单边差分算壁面速度梯度而不是直接用numpy.gradient因为numpy.gradient在边界处默认只用一阶精度。5.3 残差卡在 1e-4 不再下降压力参考点悬空现象icoFoam 的残差曲线刚开始正常下降到 1e-4 附近后变成一条水平线不管把tolerance调多小都没反应。原因不可压缩流体的压力场整体可以通过加一个常数而不改变速度场。如果fvSolution里没写pRefCell和pRefValue压力矩阵奇异迭代到一定精度后残差不再下降。解决在PISO子字典里加pRefCell 0;和pRefValue 0;并重新运行。这个坑等你在三维复杂流域里再遇到时就晚了库埃特流是提前踩雷最好的地方。5.4 OpenFOAM 周期性边界报错cyclic 面的法向和相邻配置不匹配现象checkMesh报错信息类似face is not consistent with cyclic或cyclic patch has no transform。原因blockMeshDict 中 cyclic 面的 face 列表没有成对给出或者两个面的顶点顺序不符合 OpenFOAM 的变换要求。解决确认 cyclic 边界块里第一个列表是 min 方向的面第二个列表是 max 方向的面并且两个面上的顶点按投影坐标系排序。如果你用的是导入网格再用createPatch创建周期面要确保主从面网格密度和节点位置严格一致否则就得用cyclicAMI。普通 blockMesh 网格不存在这个问题所以严格按 4.1 节的顺序写就行。5.5 时间步长调小还是发散非线性稳定的隐性要求现象显式格式已经把 dt 缩到极小u 仍然在壁面附近振荡OpenFOAM 里即便把maxCo调到 0.1照样在启动阶段发散。原因不只是线性 CFL 条件库埃特流中上下壁面速度突跳导致初始压力场不协调PISO 前几步会有强烈的压力修正波。如果nCorrectors设得过小压力修正跟不上速度变化就会发散。解决在 OpenFOAM 中用adjustTimeStep把maxCo设为 0.5 以下同时把nCorrectors调到 2。Python 的一维模型里则要对上板速度做线性斜坡启动不要让它从零瞬间跳到 0.1 m/s否则显式时间推进的前几个时间步会产生很大的梯度项。6. 用库埃特流给你的求解器做一次“体检”误差分析与网格收敛性验证到这你已经能算出正确的库埃特流了但“正确”到底是多少一个更严苛的检验方式是用网格收敛性分析来量化误差下降速率。对有限差分格式真解已知我们可以直接比较不同网格数下的 L2 误差。以 Python 求解器为例取 N 21, 41, 81, 161分别跑一遍记录 L2 误差。二阶中心差分格式下当网格数翻倍、dy 减半时L2 误差应该大约缩小到原来的 1/4。如果误差缩小比例接近 1/2说明格式只有一阶精度某个地方退化了如果误差不降反升说明边界条件写错了或是离散方程里混入了不守恒的项。N_list [21, 41, 81, 161] for N in N_list: h 0.01 U 0.1 y np.linspace(0, h, N) dy y[1] - y[0] A np.zeros((N, N)) b np.zeros(N) for i in range(1, N - 1): A[i, i-1] 1.0 A[i, i] -2.0 A[i, i1] 1.0 A[0,0] 1.0 b[0] 0.0 A[-1,-1] 1.0 b[-1] U u_soln np.linalg.solve(A, b) u_exact U * y / h l2_err np.linalg.norm(u_soln - u_exact) / np.linalg.norm(u_exact) print(fN{N:3d} L2 error{l2_err:.6e})输出的误差序列如果呈现 1/4 比例说明离散格式和边界处理都符合理论预期。这个方法可以延伸到你自己的求解器不管做的是管道流还是圆柱绕流只要有一个基准解就能用同一套误差分析流程检测空间离散的阶数。如果误差收敛速度不对首先怀疑壁面梯度计算和边界赋值的位置其次是离散方程的系数是否少乘了 dy。另一个实用技巧是用 Richardson 外推估计真解。取三套网格 N1, N2, N3利用误差随 dy 以 p 阶收敛的性质把数值解外推到 dy 趋近零的极限。这个极限和解析解之间的差值就是你的求解器在这个问题上的系统误差。我习惯每次跑完库埃特流都把这两个数打进后处理脚本里L2 误差和 Richardson 外推的真解估计。如果某个版本改了离散格式后误差阶掉到 1.5 以下那基本可以判定这改动不受控制。从那以后我每换一次求解器版本、每改一次边界条件代码都会先跑一遍库埃特流网格收敛测试确认误差阶数没掉才敢继续往后做工程算例。这个方法能挡掉很多看似正常实际上已经悄悄退化的代码改动。希望这份库埃特流资源里的 Python 脚本和 OpenFOAM 案例也能给你留下一个同样可靠的基准测试流程。本文还有配套的精品资源点击获取