中科院计算流体力学讲义CFD117讲有限体积法

有限体积法中科院计算流体力学讲义·第117讲主讲人：李新亮研究员·中科院力学所Contents目录中科院计算流体力学讲义·有限体积法专题01有限体积法基本概念与核心原理02重构与反演：格式构造的技术核心03迎风型有限体积法与Riemann求解器04多维迎风型有限体积法与坐标旋转CHAPTER01有限体积法基本概念与核心原理从守恒律的积分形式出发，建立离散求解的基本框架ComputationalFluidDynamics有限体积法的基本思想有限体积法以守恒律的积分形式为出发点，将计算域划分为互不重叠的控制体积，通过对每个控制体积积分建立离散方程。其核心优势在于天然保证守恒性，且对包含间断的流场具有坚实的数学基础。01控制体积划分：将计算域划分为一系列互不重叠的控制体积（cell），每个控制体积围绕一个网格节点或位于单元中心02积分方程转化：对每个控制体积上的守恒方程进行体积分，将微分方程转化为积分方程，使间断解在弱解意义下仍然成立03离散守恒律：控制体积内物理量的时间变化率等于通过其所有界面流入流出的净通量，离散层面严格满足守恒律04与有限差分的本质区别：有限差分基于微分形式的点态离散，有限体积基于积分形式的区域平均离散有限体积法控制体积划分示意CFD·NumericalMethods有限体积法与有限差分法的本质对比有限体积法与有限差分法在离散基础、守恒性保证、网格适应性和物理直观性四个维度存在根本差异。有限差分法FiniteDifferenceMethodTaylor展开近似—基于微分形式在网格节点处展开，要求网格光滑且拓扑规则DifferentialForm守恒性需特殊设计—普通中心差分在间断附近可能产生非物理解Non-Conservative高精度格式成熟—WENO、紧致差分等，适合基础研究与规则几何域WENO/Compact有限体积法FiniteVolumeMethod天然离散守恒—基于积分形式对控制体积内平均值演化，无需额外设计IntegralForm非结构网格适应—三角形、四面体等，对复杂几何外形适应力优异UnstructuredMesh物理图像直观—界面通量对应真实输运过程，便于理解和实现边界条件PhysicalFluxFINITEVOLUMEMETHOD一维守恒律的有限体积离散推导一维守恒律方程的有限体积离散将微分方程转化为控制体积平均值的时间演化方程。半离散格式的核心在于界面数值通量的计算，这正是重构与反演技术要解决的问题。01定义控制体积与单元平均值——在计算域上划分控制体积[xj-½,xj+½]，引入ūj=(1/Δx)∫udx，对守恒律方程在整个控制体积上积分。02Newton-Leibniz转化为边界通量——dūj/dt=−(1/Δx)[fj+½−fj-½]，体积积分化为界面通量差，即半离散有限体积格式。03核心信息缺口——已知的是单元平均值ūj，但格式需要的是界面处的数值通量fj+½，从平均值到界面值之间存在信息缺口。04时间与空间离散——时间推进可采用多种Runge-Kutta方法，空间精度完全取决于界面通量fj+½的重构精度。FINITEVOLUMEMETHOD控制体积的定义方式与选择有限体积法的控制体积有节点中心型和单元中心型两种基本定义方式。在一维均匀网格上两者等价，但在多维非结构网格上各有优劣，选择需综合考虑精度、稳定性和计算效率。节点中心型Vertex-Centered控制体积围绕网格节点构建，相邻单元中点连线形成对偶网格作为实际控制体积存储量较少（节点数通常少于单元数），适合大规模三维计算，但非结构网格上重构较复杂对偶网格单元中心型Cell-Centered未知量存储在网格单元中心，控制体积即为网格单元本身，概念更直观、实现更简单在激波捕捉和间断分辨率方面表现更稳健，高阶重构（如MUSCL、WENO）更容易在单元中心框架下实现高阶重构FINITEVOLUMEMETHOD离散守恒性：有限体积法的理论基石有限体积法在离散层面严格保持守恒性，相邻控制体积共享界面通量使得内部通量精确抵消。这一性质保证了激波的Rankine-Hugoniot条件在离散层面自动满足，是激波捕捉格式的数学基础。01相邻控制体积共享界面通量：fj+1/2对第j个单元是流出，对第j+1个单元是流入，符号相反、数值相同界面通量共享02全域求和时内部通量精确抵消（telescopingproperty），总物理量变化仅由边界通量决定，离散守恒严格成立TelescopingProperty03Rankine-Hugoniot跳跃条件在离散层面自动满足，保证激波传播速度的正确性，无需人工干预激波位置R-H条件04守恒性是获得正确弱解的必要条件：Lax-Wendroff定理表明，守恒型格式的收敛极限必为原方程的弱解Lax-Wendroff定理CHAPTER02重构与反演：格式构造的技术核心从单元平均值到界面数值通量的两步关键技术FINITEVOLUMEMETHOD重构与反演：有限体积法的两步核心技术有限体积法的空间离散分为重构和反演两个串联步骤。重构从单元平均值恢复界面处的物理量分布，决定空间精度；反演从界面两侧状态计算数值通量，决定波传播和间断处理能力。重构Reconstruction空间分布构造从单元平均值ūj构造单元内物理量的空间分布uj(x)，可采用分段常数、分段线性或高阶多项式逼近，是有限体积法实现高阶精度的基础环节精度决定阶数k阶多项式重构配合适当的数值通量可获得k+1阶空间精度，重构阶数直接决定了格式的理论收敛阶间断处理高阶重构需引入限制器或WENO权值，防止间断附近产生虚假振荡，保证解的TVD/TVB有界性反演FluxEvaluation数值通量计算利用界面左右状态uL和uR计算数值通量fj+½，本质为求解局部Riemann问题，是有限体积法的核心计算步骤耗散与分辨率从简单迎风格式到精确Riemann求解器，方法选择决定格式的数值耗散特性与激波分辨率，影响计算结果的锐度相容性条件当uL=uR=u时，数值通量应退化为物理通量f(u)，确保格式在光滑区域具有一致性FINITEVOLUMEMETHOD重构阶数：从分段常数到分段线性重构的阶数直接决定有限体积格式的空间精度。分段常数重构仅一阶精度但绝对稳定，分段线性重构（MUSCL）达到二阶精度但需配合限制器防止振荡。重构阶数的选择需要在精度和鲁棒性之间权衡。分段常数重构一阶精度u_j(x)=ū_j，界面左右值即为相邻单元平均值。实现简单、计算量小，耗散较大但绝对稳定，适用于对稳定性要求极高的场景。绝对稳定·无振荡分段线性重构MUSCL·二阶u_j(x)=ū_j+σ_j(x−x_j)，斜率由相邻单元差分计算。界面值精度显著提升，是工程计算中最常用的重构策略。需配合限制器使用斜率计算方式精度⇄稳定中心差分精度最高但可能引入振荡，单侧差分更稳定但精度降低。实际应用中需根据流场特征和网格质量进行权衡选择。Minmod·Superbee·MC高阶重构≥三阶·WENO三阶及以上需使用二次或更高次多项式，配合WENO等自适应权值技术。在光滑区实现高分辨率，间断处自动降阶保持稳定性。高分辨率·计算成本高LIMITERFUNCTIONS斜率限制器：高阶重构的振荡控制斜率限制器是保证高阶重构在间断附近不产生虚假振荡的关键技术。通过在光滑区保留完整斜率、在间断附近自动压缩斜率，实现TVD性质。经典限制器函数MinMod限制器φ(r)=max(0,min(1,r))，取左右差分中绝对值较小者，最保守但保证无振荡，激波附近耗散偏大VanLeer限制器φ(r)=(r+|r|)/(1+r)，用连续函数实现平滑过渡，兼顾稳定性和分辨率，工程中应用最广Superbee限制器φ(r)=max(0,min(2r,1),min(r,2))，在接触间断附近允许更陡斜率，分辨率高但对强激波可能产生轻微振荡理论基础Godunov定理线性单调格式最高一阶精度，高阶无振荡格式必须是非线性的——限制器正是引入非线性的机制TVD条件TotalVariationDiminishing要求限制器函数满足特定约束区域，保证总变差不增大从而杜绝虚假振荡FINITEVOLUMEMETHODWENO重构：自适应高阶精度方法WENO重构通过多候选模板的自适应加权组合，在光滑区域实现最优阶数精度，在间断附近自动降阶并保持无振荡特性。相比传统限制器，WENO在复杂波系结构中展现出更优异的分辨率和鲁棒性。构造r个候选子模板，每个子模板生成一个k阶重构多项式，整体最优模板可达2k−1阶精度2k−1阶光滑度指示子βk衡量各子模板上解的光滑程度，含间断的模板权重被自动压低至接近零βk最终重构值为各子模板多项式的非线性加权平均：光滑区恢复最优精度，间断区自适应降阶保证无振荡自适应降阶WENO-Z改进版通过全局光滑度指标τ重新定义权重，在临界点处恢复设计精度，减少经典WENO的过度耗散τFINITEVOLUMEMETHOD数值通量计算：基础反演方法数值通量的计算（反演）是有限体积法的第二步核心技术。从中心通量到Lax-Friedrichs再到Rusanov通量，通过引入与波速相关的数值耗散实现格式稳定。01中心通量f=[f(uL)+f(uR)]/2，无耗散但绝对不稳定，不能单独用于双曲型方程的数值求解。不稳定02Lax-Friedrichs通量加入全局最大波速α控制的耗散项−α(uR−uL)/2实现稳定，但耗散过大。全局α03Rusanov通量使用界面局部最大特征速度αj+½代替全局α，显著减少光滑区域的数值耗散。局部α04共同特点不需要求解Riemann问题，实现简单且通用，适合多物理场耦合中的初步计算。简单通用CHAPTER03迎风型有限体积法与Riemann求解器基于波传播物理图像的数值通量构造方法ComputationalFluidDynamicsRiemann问题：迎风型方法的数学基础Riemann问题是初始条件为分段常数双曲型方程的初值问题。间断分解产生的波系结构构成了迎风型数值通量的物理基础。01Riemann问题定义：双曲守恒律在初始时刻具有分段常数初值的初值问题，界面处02解的波系结构：初始间断分解为若干基本波——激波、稀疏波、接触间断，波之间由恒定03Euler方程的Riemann解：通常包含三个波——左行波（激波或膨胀波）、接触间断、右行波，中间有两个星区状态u*L和u*R04面通量，保证格式的物理一致性激波管实验设备——Riemann问题的物理实现FINITEVOLUMEMETHOD·CFD117Godunov格式：迎风型方法的奠基之作Godunov格式（1959）通过在每个界面求解精确Riemann问题来计算数值通量，是迎风型有限体积法的开创性工作。格式天然满足熵条件但仅一阶精度，高精度扩展需结合MUSCL重构，计算效率受限于精确Riemann求解器的迭代开销。01核心步骤将单元平均值视为分段常数初值→在每个界面构造局部Riemann问题→求精确解在x/t=0处的值→计算界面通量02天然满足熵条件精确Riemann求解器自动排除膨胀激波等非物理解，无需额外添加熵修正03一阶精度限制空间精度仅为一阶（分段常数假设），需结合MUSCL重构升级为二阶Godunov型格式（即MUSCL-Hancock方法）04计算开销精确Riemann求解器需要对非线性方程组进行Newton迭代，计算量大，尤其对含复杂状态方程的实际气体问题NumericalFlux·FiniteVolumeMethodRoe近似Riemann求解器Roe求解器（1981）通过构造满足特定性质的Roe平均矩阵，将非线性Riemann问题线性化后用特征分解精确求解。格式具有高分辨率和良好的计算效率，是工程CFD中应用最广泛的数值通量方法之一。Roe性质—构造Roe平均矩阵Ã(uL,uR)，满足相容性（uL=uR时退化为Jacobian）、守恒性（通量差等于矩阵乘状态差）、双曲性三个条件通量公式—F=½[f(uL)+f(uR)]−½|Ã|(uR−uL)，第一项为中心通量，第二项为基于特征分解的迎风耗散项特征分解—|Ã|=R|Λ|R⁻¹，R为右特征向量矩阵，Λ为特征值对角阵，耗散按各特征波分量独立施加分辨率优势—对接触间断的分辨率显著优于L-F和Rusanov格式，因接触间断对应的特征值耗散项较小CFD·有限体积法·LECTURE117Roe格式的熵修正与鲁棒性分析Roe格式在跨声速膨胀区可能因特征值过零而产生不满足熵条件的膨胀激波。通过Harten-Yee等熵修正方法在特征值接近零时引入适当耗散，可有效消除这一缺陷。极端条件下的正定性保持是保证格式鲁棒性的关键。熵修正EntropyFix01Roe格式中特征值|λ|在跨声速点附近趋近于零，耗散项消失，可能产生违反熵条件的膨胀激波02Harten-Yee修正：当|λ|<ε时用ε²+λ²/(2ε)代替|λ|，在零点附近保持最小耗散量ε，消除膨胀激波03修正参数ε选取需谨慎：过小无法完全消除膨胀激波，过大则增加不必要的数值耗散降低接触间断分辨率鲁棒性增强策略01极端流动条件下（高马赫数、强激波相互作用）Roe平均可能产生负压力/密度，需配合positivity-preserving限制器02HLLC等替代求解器在保持接触间断分辨率的同时天然保证正定性，成为近年工程应用的热门选择Harten-YeeEntropyFixλ̃=|λ||λ|≥ε/ε²+λ²─────────2ε|λ|<εMinimumDissipationatSonicPointεWaveDecompositionRoe格式的新理解：特征波分解视角从特征波分解的角度重新理解Roe格式：通量差被分解为各特征波分量的叠加，耗散沿每个波的传播方向独立施加。01通量差的特征分解f(uR)−f(uL)=RΛR⁻¹(uR−uL)=Σkαkλkrk每个分量对应一个特征波的贡献，将非线性通量差线性化为特征模态的叠加02波强度与迎风耗散αk=lk·(uR−uL)，Roe耗散项对每个波分量乘以|λk|实现迎风处理，保证数值耗散的方向与信息传播方向一致03物理意义：波传播方向正特征值（右行波）耗散来自左侧状态，负特征值（左行波）耗散来自右侧状态符合双曲型方程的信息传播特性，确保数值解满足熵条件04统一框架：近似Riemann求解器不同格式的区别在于如何处理特征波分解和耗散施加方式为设计新型求解器提供理论基础，可扩展至高分辨率格式和隐式方法APPROXIMATERIEMANNSOLVERSHLL与HLLC近似Riemann求解器HLL求解器通过双波假设实现极强的鲁棒性但牺牲了接触间断分辨率；HLLC在此基础上增加中间波恢复接触间断分辨率并保持正定性，成为现代工程CFD中最受欢迎的数值通量方法之一。HLL求解器01假设Riemann解由两个波（速度SL和SR）分隔为三个区域，中间仅一个恒定状态，极大简化了求解过程02通量公式简洁且无条件保证正定性，适合高马赫数流动，但接触间断和剪切波被完全耗散，分辨率低于Roe双波结构HLLC求解器01在HLL的双波结构中增加一个中间波S*（接触间断速度），恢复三个中间状态，显著提升接触间断分辨率02天然保证正定性（密度和压力非负），在强激波相互作用等极端条件下比Roe更鲁棒，已成为主流工程CFD软件的默认选择三波恢复FiniteVolumeMethod·FluxSchemes数值通量方法特性系统对比不同数值通量方法在精度、分辨率、鲁棒性和计算效率四个维度各有优劣。工程实践中需根据流动特征选择合适的方法：低速不可压缩流适合中心格式+人工耗散，高速可压缩流推荐HLLC或Roe+MUSCL组合。常见数值通量方法特性对比通量方法基础精度接触间断分辨率鲁棒性熵条件计算开销Lax-Friedrichs一阶差（过度耗散）极强自动满足最低Rusanov一阶较差极强自动满足低Roe一阶*优中等（需熵修正）需修正中等HLL一阶差（单中间态）极强自动满足低HLLC一阶*优强自动满足中等ExactRiemann一阶*最优强自动满足最高*配合MUSCL/WENO重构可提升至二阶或更高阶精度；HLLC在鲁棒性和分辨率之间取得了最佳平衡FluxSplittingEuler方程的通量分裂方法通量分裂通过将Euler方程通量分解为正负两部分实现迎风离散，是构造迎风格式的另一重要途径。01Steger-Warming分裂基于Jacobian矩阵特征值符号将f=A⁺u+A⁻u，A⁺对应正特征值（右行波），A⁻对应负特征值（左行波）A⁺u+A⁻u02vanLeer分裂构造关于马赫数的光滑分裂函数，在声速点（M=1）处连续可微，消除了Steger-Warming分裂在声速点的非光滑缺陷M=103Roe特征投影分裂基于Roe平均矩阵的R|Λ|R⁻¹分解实现通量差分裂，物理意义最清晰，推荐使用Roe平均R|Λ|R⁻¹04与Riemann求解器通量分裂本质上是Riemann求解器的另一种表述形式，两者在数学上可以互相转化等价转化CFD·有限体积法隐式时间推进：LU-SGS方法LU-SGS方法通过近似LU分解将隐式方程组的求解转化为前后两次扫描推进，在保持隐式方法大时间步长优势的同时避免了求解大型联立方程组。LU-SGS方法原理线性化：引入差量q=Δu将隐式方程线性化为(I+Δt·∂A/∂x)q=−Δt·∂f/∂x≡RHS，避免直接求解非线性方程组近似分解：对系数矩阵进行LU分解(L+D+U)≈(L+D)D⁻¹(D+U)，L步前向扫描q*=D⁻¹(RHS−Lq*)，U步后向修正效率提升：单侧差分使每次扫描只需代入已知值即可推进，无需解联立方程组，计算量从O(N³)降至O(N)工程应用要点大CFL收敛：CFL数可取到10–100量级，相比显式方法（CFL<1）定常问题收敛速度提升一个数量级以上通量分裂：A矩阵有正有负时强行单侧差分会不稳定，需利用通量分裂A=A⁺+A⁻分别进行迎风处理CHAPTER04多维迎风型有限体积法与坐标旋转从一维到多维的关键技术跨越与工程实现FINITEVOLUMEMETHOD多维有限体积法的基本框架多维有限体积法将一维方法自然推广到二维和三维：对每个控制体积的所有界面分别计算法向数值通量。界面通量计算可转化为沿法线方向的一维Riemann问题，使得一维Riemann求解器可直接应用于多维场景。01二维积分形式—dū/dt=-(1/|Ω|)ΣkFk·nk·lk，对控制体积Ω的每条边k计算法向通量F·n并乘以边长lk。该形式将面积分转化为沿各边的线积分，是有限体积法的核心离散框架。02旋转不变性—Euler方程在坐标旋转下保持不变，界面法向通量可转化为沿法线方向的一维Riemann问题求解。这一性质保证了多维问题求解的数学严谨性。03局部坐标系分解—每个界面定义法向n和切向τ，将速度分解为vn和vτ，法向分量参与一维Riemann求解。通过局部坐标变换，将多维问题降维处理。04切向速度被动输运—切向速度在Riemann求解中作为被动标量输运，不影响法向波系结构，求解完毕后再转回全局坐标系。这种处理方式大幅简化了计算复杂度。FINITEVOLUMEMETHOD坐标旋转技术的数学实现坐标旋转通过正交变换矩阵将多维界面通量问题转化为法向一维Riemann问题。旋转矩阵的正交性保证了变换的高效性和数值精确性。01旋转矩阵与速度分解旋转矩阵T将全局速度(u,v)分解为法向分量vₙ=unₓ+vnᵧ和切向分量vτ=−unᵧ+vnₓ，密度和压力保持不变。vₙ·vτ02通量计算三步流程U′L=TUL，U′R=TUR→F′½=Riemann(U′L,U′R)→F·n=T⁻¹F′½，三步完成通量求解。3Steps03正交矩阵高效性T为正交矩阵（T⁻¹=Tᵀ），逆变换仅需转置运算，计算开销极低，且保证变换过程的数值精确性。T⁻¹=Tᵀ04三维空间推广需要两个切向分量，旋转矩阵扩展为3×3速度子块，但坐标旋转与法向通量求解的基本思路完全一致。3×3CFD·有限体积法多维网格类型与工程实践多维有限体积法支持结构化、非结构和混合三类网格。结构化网格实现简单但几何适应性差，非结构网格适应复杂外形但数据结构复杂，混合网格结合两者优势。网格类型选择01结构化网格（四边形/六面体）：拓扑关系明确、内存访问连续、计算效率高，但复杂几何外形下网格生成困难HighEfficiency02非结构网格（三角形/四面体）：自动生成容易、几何适应性强，但邻居关系需额外存储、内存访问不连续FlexibleGeometry03混合网格：边界层用棱柱层/六面体保证壁面法向分辨率，远场用四面体适应外形，工程中应用最广BestPractice梯度计算与网格质量01非结构网格上梯度计算采用最小二乘法或Green-Gauss方法，精度和鲁棒性受网格质量直接影响LSQ/Green-Gauss02网格质量关键指标：纵横比（aspectratio）、偏斜度（skewness）、最小角度，差质量网格导致重构精度严重退化Aspect·Skewness·AngleFINITEVOLUMEMETHOD高阶多维有限体积方法高阶多维有限体积法通过k-exact重构或DG方法实现三阶及以上空间精度。k-exact重构基于扩展模板的最小二乘拟合，DG方法在单元内独立存储高阶多项式。两者在精度、效率和实现复杂度上各有取舍。01k-exact重构构造k阶多项式p(x,y)使得对每个邻居单元的积分平均值等于已知值，超定系统用最小二乘求解。最小二乘02三阶重构（k=2）需要确定6个系数，至少6个约束条件（自身+5个邻居），模板选取和条件数控制是关键难点。6个系数03DG方法每个单元独立存储k阶多项式的(k+1)(k+2)/2个系数，无需跨单元重构，但存储量和计算量成倍增加。(k+1)(k+2)/204工程现状二阶MUSCL+限制器仍是主流选择，高阶方法主要用于气动声学、大涡模拟等对精度要求极高的领域。二阶MUSCLBOUNDARYCONDITIONS边界条件处理与GhostCell方法有限体积法通过GhostCell（虚