方程：可压缩 RANS 方程与 Spalart–Allmaras 模型

The equations: compressible RANS with the Spalart–Allmaras model

从守恒律到定常 RANS–SA 方程

From conservation laws to the steady RANS–SA equations

本页的方程描述可压缩气体绕物体的湍流。质量、动量、能量守恒给出 Navier–Stokes 方程；湍流在很小的尺度上快速脉动，工程计算只求平均流（雷诺平均，RANS）。平均后多出脉动造成的动量与热量输运（雷诺应力），方程不再封闭。涡粘性假设把它写成附加的粘度 μ_t，加到分子粘度上，并按湍流普朗特数 Pr_t 加到导热上；Spalart–Allmaras（SA）模型用工作变量 ν̃ 的一个输运方程给出 μ_t。于是每个点有 5 个未知量 (ρ, ρu, ρv, ρE, ρν̃)，比二维层流 N–S 多一个。

The equations of this page describe turbulent flow of a compressible gas around a body. Conservation of mass, momentum and energy gives the Navier–Stokes equations; turbulence fluctuates rapidly on small scales, and engineering calculations solve for the mean flow only (Reynolds averaging, RANS). Averaging leaves an extra transport of momentum and heat by the fluctuations (the Reynolds stresses), and the equations are no longer closed. The eddy-viscosity hypothesis writes it as an added viscosity μ_t, added to the molecular one and, through a turbulent Prandtl number Pr_t, to the heat conduction; the Spalart–Allmaras (SA) model supplies μ_t through one transport equation for a working variable ν̃. Every point then carries 5 unknowns (ρ, ρu, ρv, ρE, ρν̃), one more than two-dimensional laminar N–S.

守恒形式的 RANS–SA 方程，带时间项：

RANS–SA in conservation form, with the time term:

本页要的是定常解：∂U/∂t = 0，于是 ∇·(F^c − F^v) = S。时间项仍然要写出来，因为定常解正是靠伪时间推进得到的：从初始状态出发沿 t 推进，直到时间导数消失（第 6 节）。这里的 t 是伪时间，每个控制体取自己的步长，所以推进的中间状态不是真实的时间历程，只有推进到底得到的定常解才有物理意义。

This page wants the steady solution: ∂U/∂t = 0, so that ∇·(F^c − F^v) = S. The time term is still written, because the steady solution is reached by pseudo-time marching: from an initial state, U is advanced in t until the time derivative vanishes (Section 6). Here t is a pseudo time and every control volume takes its own step, so the intermediate states are not a physical time history; only the steady solution the marching ends at has a physical meaning.

对流通量是 Euler 通量加上随质量通量走的 ρν̃；粘性通量是应力、热流与 SA 的扩散：

The convective flux is the Euler flux plus ρν̃ carried by the mass flux; the viscous flux holds the stress, the heat flux and the SA diffusion:

其中 $p=\rho T$（取 $R=1$、$c_v=2.5$、$c_p=3.5$，所以 $\gamma =1.4$），$\mu (T)={\mu }_{\mathrm{r}\mathrm{e}\mathrm{f}}1.4T^{1.5}\mathrm{/}(T+0.4)$（Sutherland），湍流粘度 ${\mu }_t=\rho \tilde{\nu }{f}_{v1}(\chi )$，$\chi =\rho \tilde{\nu }\mathrm{/}\mu$，

with $p=\rho T$ ($R=1$, $c_v=2.5$, $c_p=3.5$, so $\gamma =1.4$), $\mu (T)={\mu }_{\mathrm{r}\mathrm{e}\mathrm{f}}1.4T^{1.5}\mathrm{/}(T+0.4)$ (Sutherland), the eddy viscosity ${\mu }_t=\rho \tilde{\nu }{f}_{v1}(\chi )$ with $\chi =\rho \tilde{\nu }\mathrm{/}\mu$, and

$Pr=0.72$、$P{r}_t=0.9$、$\sigma =2\mathrm{/}3$。源项只在第 5 个方程里，是 NASA TMR 给出的 SA 守恒形式（密度梯度项不能丢）：

$Pr=0.72$, $P{r}_t=0.9$, $\sigma =2\mathrm{/}3$. The source sits in the fifth equation only; it is the conservative form of the SA model given by NASA's TMR (the density-gradient term must not be dropped):

模型函数取标准（非负）版本：

The model functions are the standard (not the negative) version:

Ω = |∂v/∂x − ∂u/∂y| 是涡量的大小；常数 c_b1 = 0.1355、c_b2 = 0.622、κ = 0.41、c_w2 = 0.3、c_w3 = 2、c_v1 = 7.1、c_t3 = 1.2、c_t4 = 0.5、c₂ = 0.7、c₃ = 0.9，c_w1 = c_b1/κ² + (1 + c_b2)/σ。

Ω = |∂v/∂x − ∂u/∂y| is the size of the vorticity; the constants are c_b1 = 0.1355, c_b2 = 0.622, κ = 0.41, c_w2 = 0.3, c_w3 = 2, c_v1 = 7.1, c_t3 = 1.2, c_t4 = 0.5, c₂ = 0.7, c₃ = 0.9, and c_w1 = c_b1/κ² + (1 + c_b2)/σ.

模型、常数与守恒形式取自 Spalart 与 Allmaras（1994）、Allmaras、Johnson 与 Spalart（2012）以及 NASA 的湍流模型资源网站 TMR；第 3 节的 Roe 通量见 Roe（1981），粘性通量的修正平均梯度与格点页的边式格式见 SU2（Economon 等，2016）。书目列在附录 A 末尾。

The model, its constants and the conservative form follow Spalart and Allmaras (1994), Allmaras, Johnson and Spalart (2012) and NASA's Turbulence Modeling Resource (TMR); the Roe flux of Section 3 is Roe's (1981), and the corrected average gradient of the viscous flux and the node page's edge-based scheme are SU2's (Economon et al., 2016). The references are listed at the end of Appendix A.

这些函数各有来由。ν̃ 是 SA 的工作变量：离开壁面它就是涡粘性 ν_t = μ_t/ρ（f_v1 → 1）；靠近壁面 f_v1 把它压低，ν_t = ν̃f_v1(χ)，于是 ν̃ 一直到壁面都随 d 线性变化，比 ν_t 容易算准。方程右端是生成减耗散：湍流由剪切产生，所以生成项与涡量 S̃ 成正比；壁面限制涡的尺度，所以耗散项按 (ν̃/d)² 增长。S̃ = Ω + ν̃f_v2/(κ²d²) 在对数层里等于 u_τ/(κd)，f_v2 让这一关系一直保持到壁面。r = ν̃/(S̃κ²d²) 比较湍流的长度尺度与壁面距离：对数层里 r = 1、f_w = 1；离壁面越远 r 越大，f_w 趋于常数 (1 + c_w3⁶)^1/6 ≈ 2.0，所以把 r 截断在 10 并不改变模型，只免得 r⁶ 溢出。f_t2 的两项让 ν̃ = 0 成为一个稳定的解，层流区因此不会自己转捩。χ 略大于 1 时 f_v2 变负，B 随之为负，Ω + B 可能小于零；2012 年的修正在 B < −c₂Ω 时换成第二个式子，它与第一支在 B = −c₂Ω 处值与一阶导数都相接，并让 S̃ 始终大于 0.1Ω。

Each of these functions has a purpose. ν̃ is SA's working variable: away from walls it is the eddy viscosity ν_t = μ_t/ρ (f_v1 → 1); near a wall f_v1 damps it, ν_t = ν̃f_v1(χ), so that ν̃ varies linearly with d all the way to the wall and is easier to compute accurately than ν_t. The right-hand side is production minus destruction: turbulence is produced by shear, so production is proportional to the vorticity S̃; the wall limits the size of the eddies, so destruction grows as (ν̃/d)². S̃ = Ω + ν̃f_v2/(κ²d²) equals u_τ/(κd) in the log layer, and f_v2 keeps that relation all the way to the wall. r = ν̃/(S̃κ²d²) compares the turbulent length scale with the wall distance: r = 1 and f_w = 1 in the log layer; further from the wall r grows and f_w tends to the constant (1 + c_w3⁶)^1/6 ≈ 2.0, so capping r at 10 does not change the model and only keeps r⁶ from overflowing. The two f_t2 terms make ν̃ = 0 a stable solution, so a laminar region does not turn turbulent by itself. For χ slightly above 1, f_v2 turns negative, B with it, and Ω + B can drop below zero; the 2012 modification switches to the second expression for B < −c₂Ω, which meets the first branch at B = −c₂Ω in value and first derivative and keeps S̃ above 0.1Ω.

波速与边界条件

Wave speeds and boundary conditions

对流部分决定信息怎样传播。沿单位法向 n，对流通量 Jacobian 的五个特征值是两道声波 q ∓ c（c 是声速），以及都以 q = u·n 移动的熵波、剪切波和 ν̃ 波；本页 M = 0.2，处处亚声速。n 指向区域外时，特征值为负的波从外面进来；无粘极限下，边界上能给定几个值就看进来几道波：来流（q < 0）四道，出口（q > 0）一道。粘性项还要一些条件：来流给定 ρ、u、v、ν̃，T 自由、热流为零；出口给定牵引力 p n − τ·n = p∞n，热流与 ν̃ 的扩散通量为零；物面无滑移、绝热，ν̃ = 0。第 4 节的图 4.2 逐个边界数这些条件，同一节写出它们的离散。

The convective part decides how information travels. Along a unit normal n the five eigenvalues of the convective flux Jacobian are the acoustic waves q ∓ c (c the speed of sound) and the entropy, shear and ν̃ waves, all three at q = u·n; at M = 0.2 the flow is subsonic everywhere. With n pointing out of the domain, a wave with a negative eigenvalue comes in, and in the inviscid limit a boundary takes as many values as waves come in: four at the inflow (q < 0), one at the outlet (q > 0). The viscous terms need a few more: the inflow gives ρ, u, v and ν̃, with T free and zero heat flux; the outlet gives the traction p n − τ·n = p∞n and zero heat and ν̃ diffusive fluxes; the body is no-slip and adiabatic, with ν̃ = 0. Figure 4.2 counts these conditions boundary by boundary, and Section 4 writes out their discretisation.

四个 RANS–SA 页面解同一个连续问题：格心、格点两种控制体，各有两种来流边界的离散；同一种控制体的两个页面只在来流边界的离散上不同。本页沿用原来的处理：来流面用 Roe 通量，外侧放一个取来流 ρ、u、v、ν̃ 的虚拟状态；伴随一致的版本里，来流面取边界状态的物理通量，质量通量就是给定的 ρ∞u∞·n，不加 Roe 耗散。两者都给出各自离散问题的精确导数；区别在于加密时离散伴随在来流附近是否趋于连续伴随（伴随一致性）。做法见第 4 节，原因与逐级比较见附录 B。

The four RANS–SA pages solve the same continuous problem: two kinds of control volume, each with two discretisations of the inflow boundary; the two pages of one kind differ only in the inflow discretisation. This page keeps the original treatment: the inflow faces take the Roe flux against a ghost state holding the free stream's ρ, u, v and ν̃; in the adjoint-consistent version the inflow faces take the physical flux of the boundary state, whose mass flux is the given ρ∞u∞·n, with no Roe dissipation. Both give the exact derivatives of their own discrete problems; they differ in whether, under refinement, the discrete adjoint approaches the continuous one near the inflow (adjoint consistency). Section 4 gives the two treatments, Appendix B the reason and a level-by-level comparison.

本页的问题

The problem on this page

区域是物体与远场之间的环：物体是凸八边形（弦长 1.1），远场是半径 1.9 的正八边形，第 2 节的网格填满两者之间。来流取 ρ∞ = T∞ = p∞ = 1，于是声速 c∞ = √γ，U∞ = M√γ。按伴随优化文献的通常写法，设计变量记作 D = (α, M, μ_ref, X)，共 51 个，X 是 24 个节点的 48 个坐标；目标函数是第 7 节的阻力系数 C_D（下标 D 指阻力，不是设计变量）。基线解是页首默认设置（二阶）下 R(U) = 0 的收敛解（第 6 节），C_D = 6.2300，这张粗网格上的阻力有多准见第 7 节。图 1.1 把这些画在一起。

The domain is the ring between the body, a convex octagon (chord 1.1), and the far field, a regular octagon of radius 1.9; the mesh of Section 2 fills it. The free stream has ρ∞ = T∞ = p∞ = 1, so c∞ = √γ and U∞ = M√γ. As is usual in the adjoint-optimization literature, the design variables are written D = (α, M, μ_ref, X), 51 in all, with X the 48 coordinates of the 24 nodes; the objective is the drag coefficient C_D of Section 7 (the subscript D stands for drag, not for the design). The baseline is the converged solution of R(U) = 0 at the header's default settings (second order; Section 6), with C_D = 6.2300; Section 7 discusses how accurate the drag of this coarse mesh is. Figure 1.1 draws all of this.

这组方程的作用

What these equations are for

在流体力学里，RANS 加一个湍流模型是计算飞机、叶轮机械等高雷诺数定常流动的常用办法：不解析每一个涡，只算平均流，一次计算就给出升力、阻力与压力分布。SA 模型最初为航空外流设计，只多一个输运方程，工程上用得很广。与层流 N–S 相比，湍流模型带来三样东西：多一个未知量和方程；涡粘性 μ_t = ρν̃f_v1(χ) 进入动量与能量方程的粘性通量；依赖壁面距离 d 的源项，而 d 是网格坐标的函数。图 1.2 画出 f_v1 怎样把 ν̃ 变成涡粘性。

In fluid dynamics, RANS with a turbulence model is the usual way to compute steady high-Reynolds-number flows around aircraft, turbomachinery and the like: it resolves no individual eddy, only the mean flow, and one calculation gives the lift, the drag and the pressure distribution. The SA model was designed for aerodynamic external flows; it adds one transport equation and is widely used in engineering. Compared with laminar N–S the turbulence model brings three things: one more unknown and equation; the eddy viscosity μ_t = ρν̃f_v1(χ) in the viscous fluxes of momentum and energy; and source terms depending on the wall distance d, a function of the mesh coordinates. Figure 1.2 shows how f_v1 turns ν̃ into an eddy viscosity.

在本课程里，二维非结构网格上 RANS–SA 的前向与伴随是整个系列的目标。用于气动外形优化的工业伴随求解器求的正是这类方程组：一次伴随求解给出阻力这样一个目标对成百上千个外形参数的导数。前面各页一次加一样东西，本页把它们合在一起（图 1.3）。本页的前向与伴随对五个方程一起求导，涡粘性对 ν̃ 的依赖、SA 源项与壁面距离都在其中；有些工程代码冻结湍流量、只对平均流方程求伴随，本页不做这种近似。

In this course, the tangent and the adjoint of RANS–SA on a two-dimensional unstructured mesh are the goal of the whole series. Industrial adjoint solvers for aerodynamic shape optimisation solve exactly this kind of system: one adjoint solve gives the derivative of one objective, such as the drag, with respect to hundreds or thousands of shape parameters. The earlier pages add one ingredient at a time; this page puts them together (Figure 1.3). Its tangent and adjoint differentiate all five equations together, including the eddy viscosity's dependence on ν̃, the SA source and the wall distance; some production codes freeze the turbulence quantities and take the adjoint of the mean-flow equations only, and this page makes no such approximation.

本页要回答的问题

The question this page answers

本页只回答一个问题：51 个设计变量 D 各改变一点，C_D 改变多少？页首的全页路线图把后面各节串成一条线：前向每个设计变量解一次线性方程，共 51 次；伴随只解一次转置方程。

This page answers one question: how much does C_D change when each of the 51 design variables D changes a little? The page roadmap at the top strings the later sections together: the tangent solves one linear system per design variable, 51 in all, the adjoint a single transposed one.

要在计算机上回答它，先要把上面的偏微分方程变成有限个未知量的代数方程 R(U) = 0。第 2 节放下网格与控制体，第 3 节再把每个控制体上的守恒写成面通量之和：Roe 通量、最小二乘重构、粘性通量与 SA 源项。

To answer it on a computer, the partial differential equations above must first become algebraic equations R(U) = 0 in finitely many unknowns. Section 2 lays out the mesh and the control volumes, and Section 3 writes conservation on each as a sum of face fluxes: the Roe flux, the least-squares reconstruction, the viscous flux and the SA source.

网格：单元、面与壁面距离

Mesh: cells, faces and wall distance

第 1 节的方程在区域里每一点成立，计算机却只能存有限个数。本节决定未知量放在哪里：网格、控制体与它们的几何量，包括 SA 源项要用的壁面距离。

The equations of Section 1 hold at every point of the domain, but a computer stores finitely many numbers. This section decides where the unknowns live: the mesh, the control volumes and their geometry, including the wall distance the SA source needs.

网格与二维对流页面完全相同，编号也相同：节点逐圈从壁面向外，单元逐圈编号，面按“内部、壁面、远场”排序。未知量放在单元里：24 个单元就是 24 个控制体。区别在于每个控制体现在有 5 个未知量，而不是 1 个。

The mesh is exactly that of the two-dimensional advection page, with the same numbering: nodes ring by ring from the wall outwards, cells ring by ring, faces sorted interior, wall, far field. The unknowns live in the cells: the 24 cells are the 24 control volumes. The difference is that every control volume now carries 5 unknowns instead of 1.

一个面循环要读的几何量都画在下图里。和标量页面相比多了两样：从控制体中心 x_i 到面点 x_f 的重构方向 x_f − x_i（二阶格式沿它把中心值外推到面上），与两侧中心的连线 Δx = x_R − x_L（粘性通量沿它取两点差分）。

The geometric quantities a face loop reads are all in the figure below. Compared with the scalar pages there are two more: the reconstruction direction x_f − x_i from a control-volume centre x_i to the face point x_f (the second-order scheme extrapolates the centre value along it), and the line Δx = x_R − x_L between the two centres (the viscous flux takes its two-point difference along it).

壁面距离

The wall distance

SA 模型的耗散项与 (ν̃/d)² 成正比，d 是到壁面的距离。本页的 d_i 是控制体中心到离散壁面多边形的精确欧氏距离：对每条壁面边求点到线段的距离，取最小。它是坐标的函数，所以在第 12 节的几何导数里也要求导。

The SA model's destruction term goes as (ν̃/d)², with d the distance to the wall. On this page d_i is the exact Euclidean distance from the control-volume centre to the discrete wall polygon: the distance from the point to each wall edge, and the smallest. It is a function of the coordinates, so the geometric derivatives of Section 12 differentiate it too.

逐步推导：壁面距离与它对坐标的导数

Step by step: the wall distance and its derivative with respect to the coordinates

第 12 节的几何导数要用到 ∂d_i/∂X。d_i 是分段定义的：最近点可能落在一条壁面边的内部，也可能落在顶点上。先把它写成显式公式，再逐项求导。记控制体中心为 x_i，一条壁面边从壁面节点 X_a 到 X_b，e = X_b − X_a。

The geometric derivatives of Section 12 need ∂d_i/∂X. d_i is defined piecewise: the nearest point may lie inside a wall edge or at a vertex. Write it as an explicit formula first, then differentiate term by term. Let x_i be the control-volume centre and let a wall edge run from the wall node X_a to X_b, with e = X_b − X_a.

第一步，直线上的最近点。直线上的点是 X_a + ξe；到 x_i 的距离平方对 ξ 求导并令其为零：

Step 1, the nearest point on the line. The points of the line are X_a + ξe; differentiate the squared distance to x_i with respect to ξ and set it to zero:

第二步，截到线段上。ξ* ≤ 0 时最近点是顶点 X_a，ξ* ≥ 1 时是 X_b，否则是垂足。记截断后的参数 ξ = min(max(ξ*, 0), 1)，最近点 x_w = X_a + ξe，距离与它的单位方向是

Step 2, clamp to the segment. For ξ* ≤ 0 the nearest point is the vertex X_a, for ξ* ≥ 1 it is X_b, and otherwise the foot of the perpendicular. With the clamped parameter ξ = min(max(ξ*, 0), 1) the nearest point is x_w = X_a + ξe, and the distance and its unit direction are

d_i 取所有壁面边中最小的那个（mesh.js 的 segDist 与 geometry 正是这样算的，分支都按实数值选）。

d_i is the smallest over all wall edges (mesh.js segDist and geometry compute exactly this, every branch chosen on real values).

第三步，垂足情形求导。d² = |x_i − X_a − ξe|²，而 ξ = ξ*(x_i, X_a, X_b) 本身也随坐标变。但 ξ* 使这个距离平方取极小，∂(d²)/∂ξ = 0：ξ 的变化对 d 没有一阶影响（包络定理）。所以只对显式出现的三个点求前向导数（图 2.4 (a)）：

Step 3, differentiate in the foot case. d² = |x_i − X_a − ξe|², and ξ = ξ*(x_i, X_a, X_b) moves with the coordinates as well. But ξ* minimises this squared distance, so ∂(d²)/∂ξ = 0: a change of ξ has no first-order effect on d (the envelope theorem). Only the three points that appear explicitly are differentiated, in forward mode (Figure 2.4 (a)):

n_w 与 e 垂直，所以沿壁面边移动 X_a 或 X_b 不改变 d；三个导数之和为零，整体平移不改变 d。

n_w is perpendicular to e, so moving X_a or X_b along the edge leaves d unchanged; the three derivatives sum to zero, and a rigid translation does not change d.

第四步，顶点情形。ξ* ≤ 0 时 d = |x_i − X_a|，直接求导得 ∂d/∂x_i = n_w、∂d/∂X_a = −n_w、∂d/∂X_b = 0：正是第三步的公式取 ξ = 0（ξ* ≥ 1 时取 ξ = 1）。一个公式覆盖两种情形；在 ξ* = 0 处两边的值与导数都相同，所以最近点从顶点滑进边内时，d 的导数不跳。

Step 4, the vertex case. For ξ* ≤ 0, d = |x_i − X_a|, and direct differentiation gives ∂d/∂x_i = n_w, ∂d/∂X_a = −n_w, ∂d/∂X_b = 0: Step 3's formula with ξ = 0 (and ξ = 1 for ξ* ≥ 1). One formula covers both cases; at ξ* = 0 both sides have the same value and the same derivative, so the derivative of d does not jump when the nearest point slides from the vertex into the edge.

第五步，回到节点坐标。格心页的中心是单元顶点的平均，x_i = (1/k)Σ_nX_n，所以 ∂d_i/∂x_i = n_w 平分给单元的 k 个顶点，每个 n_w/k；X_a、X_b 是壁面节点，各得到第三步的一项。

Step 5, back to the node coordinates. On the cell page the centre is the average of the cell's vertices, x_i = (1/k)Σ_nX_n, so ∂d_i/∂x_i = n_w is shared equally by the cell's k vertices, n_w/k each; X_a and X_b are wall nodes and receive the terms of Step 3.

核对：本页 24 个中心里，16 个的最近点在边内，8 个在顶点上；对每个 d_i，第三至第五步的公式与 geometry() 的磁带反向扫描给出的 ∂d_i/∂X 最大相对差 6.2×10⁻¹⁶（本页代码在构建时计算）。

Check: of this page's 24 centres, 16 have their nearest point inside an edge and 8 at a vertex; for every d_i, the formulas of Steps 3 to 5 and a reverse sweep of a tape of geometry() give ∂d_i/∂X with a largest relative difference of 6.2×10⁻¹⁶ (computed by this page's code at build time).

网格与几何量到这里就齐了。第 3 节把第 1 节的守恒律对每个控制体积分，写成穿过这些面的通量之和。

The mesh and its geometry are now complete. Section 3 integrates the conservation laws of Section 1 over each control volume and writes them as sums of fluxes through these faces.

面上的通量

The flux on a face

第 2 节给出了控制体和它们的面。本节把第 1 节的方程对每个控制体积分，写出每个面上的通量：Roe 通量、最小二乘重构、粘性通量，以及每个控制体的 SA 源项。

Section 2 gave the control volumes and their faces. This section integrates the equations of Section 1 over each control volume and writes the flux through each face: the Roe flux, the least-squares reconstruction, the viscous flux, and each control volume's SA source.

有限体积形式

The finite-volume form

对每个控制体积分（U_i 取控制体上的平均值），散度变成面上的通量和，每个控制体得到一个关于 t 的常微分方程（半离散形式）：

Integrated over each control volume (U_i is the average over it), the divergence becomes a sum of face fluxes, and each control volume gets an ordinary differential equation in t (the semi-discrete form):

roe 是对流部分，visc 是粘性部分。面通量 Φ_f 从 L 流向 R：对 L 加、对 R 减，保证守恒。W = (ρ, u, v, T, ν̃) 是原始变量，G 是它们的最小二乘梯度，W⁺ 是重构到面上的状态。

roe is the convective part and visc the viscous one. The face flux Φ_f runs from L to R: added to L and subtracted from R, which makes the scheme conservative. W = (ρ, u, v, T, ν̃) are the primitive variables, G their least-squares gradients, and W⁺ the state reconstructed to the face.

定常解满足 R(U) = 0：这就是本页要解、并在之后各节求导的离散方程；时间项只在求解时出现（第 6 节）。

The steady solution satisfies R(U) = 0: that is the discrete equation this page solves and, in the later sections, differentiates; the time term appears in the solve only (Section 6).

每个面仍是一个一维问题：Roe 通量

Every face is still a one-dimensional problem: the Roe flux

面的两侧各有一个状态，Roe 通量把它们的差分解成五道波：两道声波（速度 q̃ ∓ c̃）、熵波、剪切波与 ν̃ 波（速度都是 q̃），每道波按它速度的绝对值加耗散：

Each side of a face carries a state, and the Roe flux splits their difference into five waves: two acoustic waves (speeds q̃ ∓ c̃), and the entropy, shear and ν̃ waves (all at speed q̃); each wave is dissipated in proportion to the size of its speed:

~ 表示 Roe 平均（以 √ρ 加权：ρ̃ = √(ρ_Lρ_R)，ũ、H̃、ν̃~ 是两侧的加权平均），c̃² = (γ − 1)(H̃ − ½|ũ|²)，q̃ = ũ·n̂。Δ 表示 R 侧减 L 侧的跳跃，Δq = Δu·n̂。Λ_k 是第 k 道波的波速，r_k 是 Roe 平均状态处通量 Jacobian 的右特征向量（分量按 ρ、ρu、ρv、ρE、ρν̃ 排列），Δ_k 是波强：把跳跃 ΔU 写成 Σ_kΔ_kr_k 时的系数。

The ~ marks a Roe average (weighted by √ρ: ρ̃ = √(ρ_Lρ_R), and ũ, H̃, ν̃~ are the weighted averages of the two sides), c̃² = (γ − 1)(H̃ − ½|ũ|²) and q̃ = ũ·n̂. Δ is a jump, the R side minus the L side, and Δq = Δu·n̂. Λ_k is the speed of wave k, r_k the right eigenvector of the flux Jacobian at the Roe-averaged state (components ordered ρ, ρu, ρv, ρE, ρν̃), and Δ_k the wave strength: the coefficient of r_k when the jump ΔU is written as Σ_kΔ_kr_k.

两道声波（k = 1、3）、熵波（k = 2）、剪切波（k = 4）与 ν̃ 波（k = 5）；后三道都以 q̃ 移动。剪切波与 ν̃ 波直接写成乘积 Δ_kr_k。ρν̃ 作为被动标量跟着质量通量走：r_1,2,3 的最后一个分量是 ν̃~，所以它的耗散是 ν̃~ 乘质量方程的耗散，再加 ν̃ 波的 |q̃|ρ̃Δν̃。残差里的特征值不做熵修正：熵修正是为了在声速点（|q̃ ± c̃| → 0）附近防止出现非物理的膨胀激波，本页 M = 0.2，处处亚声速，用不到它。代价是 |Λ_k|（k = 2, 4, 5）在 q̃ = 0 处不可导（第 4 节列出它离当前解多远）。

The two acoustic waves (k = 1, 3), the entropy wave (k = 2), the shear wave (k = 4) and the ν̃ wave (k = 5); the last three all move at q̃. The shear and ν̃ waves are written directly as the products Δ_kr_k. ρν̃ rides on the mass flux as a passive scalar: the last component of r_1,2,3 is ν̃~, so its dissipation is ν̃~ times the mass equation's, plus the ν̃ wave's |q̃|ρ̃Δν̃. The residual's eigenvalues carry no entropy fix: an entropy fix keeps non-physical expansion shocks from forming near sonic points (|q̃ ± c̃| → 0), and this page's flow at M = 0.2 is subsonic everywhere. The price is that |Λ_k| (k = 2, 4, 5) is not differentiable at q̃ = 0 (Section 4 lists how far the solution is from it).

逐步推导：Roe 平均、波分解与 Roe 通量的线性化

Step by step: Roe averages, the wave decomposition and the linearisation of the Roe flux

第一步，Roe 平均的乘积法则。任一量 φ 的 Roe 平均是 φ̃ = (√ρ_Lφ_L + √ρ_Rφ_R)/(√ρ_L + √ρ_R)，ρ̃ = √(ρ_Lρ_R)。直接展开（用 ρ_R − ρ_L = (√ρ_R − √ρ_L)(√ρ_R + √ρ_L)）：

Step 1, the product rule of Roe averages. The Roe average of any quantity φ is φ̃ = (√ρ_Lφ_L + √ρ_Rφ_R)/(√ρ_L + √ρ_R), and ρ̃ = √(ρ_Lρ_R). Expanding directly (with ρ_R − ρ_L = (√ρ_R − √ρ_L)(√ρ_R + √ρ_L)):

所以对 φ = u、v、H、ν̃，跳跃 Δ(ρφ) 精确地等于 φ̃Δρ + ρ̃Δφ：在 Roe 平均处把跳跃线性化没有误差。

So for φ = u, v, H, ν̃ the jump Δ(ρφ) is exactly φ̃Δρ + ρ̃Δφ: linearising the jump at the Roe average makes no error.

第二步，压力的跳跃。本页 p = ρT、ρE = ρ(c_vT + ½|u|²)，所以 p = (γ − 1)(ρE − ½ρ|u|²)。写 ρ|u|² = |√ρ u|²，并用 √ρ_Lu_L + √ρ_Ru_R = (√ρ_L + √ρ_R)ũ，得 ½Δ(ρ|u|²) = ũ·Δ(ρu) − ½|ũ|²Δρ，于是

Step 2, the pressure jump. On this page p = ρT and ρE = ρ(c_vT + ½|u|²), so p = (γ − 1)(ρE − ½ρ|u|²). Writing ρ|u|² = |√ρ u|² and using √ρ_Lu_L + √ρ_Ru_R = (√ρ_L + √ρ_R)ũ gives ½Δ(ρ|u|²) = ũ·Δ(ρu) − ½|ũ|²Δρ, hence

第三步，求波强。把跳跃写成 ΔU = Σ_kΔ_kr_k，逐个分量比较（r_k 见上式）。质量分量：Δ₁ + Δ₂ + Δ₃ = Δρ。法向动量 n̂·Δ(ρu) = q̃Δρ + ρ̃Δq：c̃(Δ₃ − Δ₁) = ρ̃Δq。把第二步 Δp 的式子（γ − 1 乘方括号）作用在每个 r_k 上，r₁、r₃ 给 c̃²，r₂ 与剪切、ν̃ 两个向量给 0，所以 c̃²(Δ₁ + Δ₃) = Δp。三式联立：

Step 3, the wave strengths. Write the jump as ΔU = Σ_kΔ_kr_k and compare component by component (the r_k as above). The mass component: Δ₁ + Δ₂ + Δ₃ = Δρ. The normal momentum n̂·Δ(ρu) = q̃Δρ + ρ̃Δq: c̃(Δ₃ − Δ₁) = ρ̃Δq. Step 2's formula for Δp ((γ − 1) times the bracket) applied to each r_k gives c̃² for r₁ and r₃ and 0 for r₂ and the shear and ν̃ vectors, so c̃²(Δ₁ + Δ₃) = Δp. Together:

动量里剩下的 ρ̃(Δu − Δq n̂) 是剪切波；ρν̃ 里剩下的 ρ̃Δν̃（第一步取 φ = ν̃）是 ν̃ 波：正是上面的 Δ₄r₄、Δ₅r₅。

What remains in the momentum, ρ̃(Δu − Δq n̂), is the shear wave; what remains in ρν̃, ρ̃Δν̃ (Step 1 with φ = ν̃), is the ν̃ wave: the Δ₄r₄ and Δ₅r₅ above.

第四步，耗散矩阵。r_k 是 Roe 平均状态处通量 Jacobian 的特征向量，特征值 Λ_k。把跳跃按特征向量展开、逐波乘 |Λ_k|，就是 |J̃^c|ΔU。记 𝓡 = (r₁, …, r₅)：

Step 4, the dissipation matrix. The r_k are the eigenvectors of the flux Jacobian at the Roe-averaged state, with eigenvalues Λ_k. Expanding the jump in eigenvectors and multiplying each wave by |Λ_k| gives |J̃^c|ΔU. With 𝓡 = (r₁, …, r₅):

第二式来自 r₁、r₂、r₃ 的最后一个分量都是 ν̃~、剪切波的是 0（model.js 的 d4 = st·d0 + sh·ds）。所以 |J̃^c| 的 ρν̃ 列只有对角元 |q̃|：ν̃ 不进入平均流方程的耗散，而 ρν̃ 的耗散是质量耗散乘 ν̃~ 再加自己的一项。本页一个面上的 |J̃^c| 见图 3.2。

The second formula holds because the last component of r₁, r₂ and r₃ is ν̃~ and that of the shear wave is 0 (model.js's d4 = st·d0 + sh·ds). So the ρν̃ column of |J̃^c| holds only the diagonal |q̃|: ν̃ does not enter the mean-flow equations' dissipation, and ρν̃'s dissipation is the mass dissipation times ν̃~ plus a term of its own. Figure 3.2 shows |J̃^c| on one face of this page.

第五步，线性化。一阶格式 Φ = |n|[½(F(U_L) + F(U_R))·n̂ − ½|J̃^c|ΔU]，而 J̃^c 同时依赖两侧状态，所以

Step 5, the linearisation. At first order Φ = |n|[½(F(U_L) + F(U_R))·n̂ − ½|J̃^c|ΔU], and J̃^c depends on both states, so

最后一项与跳跃 ΔU 成正比，含 |Λ_k| 的导数 sgn Λ_k（Λ_k = 0 处不可导，第 4 节列出了余量）。局部磁带对核的每个运算求导，三项全部精确算到；第 6 节的近似 Jacobian 只留前两项。二阶格式的面状态是 W_L⁺ = W_L + G_L·(x_f − x_L)，所以 Roe 部分对 W_L 的导数等于它对 W_L⁺ 的导数，对 G_L 的导数是它再乘 (x_f − x_L)^T；粘性通量另外直接读 W_L、W_R 与 G_L、G_R，所以面核在两种格式里都有这 30 个输入。

The last term is proportional to the jump ΔU and contains the derivative of |Λ_k|, sgn Λ_k (not differentiable at Λ_k = 0; Section 4 lists the margins). The local tape differentiates every operation of the kernel, so all three terms come out exactly; the approximate Jacobian of Section 6 keeps only the first two. At second order the face state is W_L⁺ = W_L + G_L·(x_f − x_L), so the Roe part's derivative with respect to W_L equals its derivative with respect to W_L⁺, and with respect to G_L it is that times (x_f − x_L)^T; the viscous flux reads W_L, W_R and G_L, G_R directly as well, which is why the face kernel has these 30 inputs at either order.

核对：本页全部 32 个内部面上（U* 处的二阶面状态），J^c(W̃)r_k = Λ_kr_k 的相对残差最大 2.7×10⁻¹⁶；Σ_kΔ_kr_k 重建跳跃 ΔU 到 9.5×10⁻¹⁴；𝓡|Λ|𝓡⁻¹ΔU 与 model.js 的耗散相差 3.3×10⁻¹⁴；Roe 性质 F(U_R) − F(U_L) = J^c(W̃)ΔU 成立到 4.5×10⁻¹⁴；|J̃^c| 的 ρν̃ 行与列符合第四步到 4.1×10⁻¹⁸。第五步被丢掉的最后一项（一阶、中心状态）占磁带给出的 ∂Φ/∂U_L 的比例（Frobenius 范数）：最大 5.7%，中位数 1.4%（本页代码在构建时计算）。

Check: on all 32 interior faces of this page (the second-order face states at U*), J^c(W̃)r_k = Λ_kr_k holds to a relative residual of at most 2.7×10⁻¹⁶; Σ_kΔ_kr_k rebuilds the jump ΔU to 9.5×10⁻¹⁴; 𝓡|Λ|𝓡⁻¹ΔU agrees with model.js's dissipation to 3.3×10⁻¹⁴; Roe's property F(U_R) − F(U_L) = J^c(W̃)ΔU holds to 4.5×10⁻¹⁴; and the ρν̃ row and column of |J̃^c| follow Step 4 to 4.1×10⁻¹⁸. The last term dropped in Step 5 (first order, centre states) is at most 5.7% of the taped ∂Φ/∂U_L (Frobenius norm), 1.4% at the median (computed by this page's code at build time).

第二步的恒等式与第三步的分解另由 sympy 符号验证（check_enrich_sympy.py）。

The identity of Step 2 and the decomposition of Step 3 are also verified symbolically by sympy (check_enrich_sympy.py).

最小二乘重构：一阶与二阶

Least-squares reconstruction: first and second order

每个控制体的梯度来自无权最小二乘：

Each control volume's gradient comes from unweighted least squares:

模板是面相邻的单元。贴壁单元还有一个“壁面点”（壁面面的中点），那里 u = v = ν̃ = 0 已知；ρ 与 T 在绝热壁上不知道，所以贴壁单元用两组权重：u、v、ν̃ 用含壁面点的一组，ρ、T 用只含邻居的一组。两组都对线性场精确。远场面不进入模板。二阶格式把中心值沿 x_f − x_i 外推：W⁺ = W_i + G_i·(x_f − x_i)，不加限制器（流动无激波）；一阶格式直接取 W⁺ = W_i。粘性通量与 SA 源项在两种格式里都用这些梯度，所以一阶／二阶开关只改变对流通量；壁面面与来流面例外：它们的粘性系数 μ(T) 用重构到面上的温度 T_w、T_f，一阶时取单元值。

The stencil is the face neighbours. A wall cell also has a “wall point” (its wall face's midpoint), where u = v = ν̃ = 0 are known; ρ and T are not known on an adiabatic wall, so a wall cell has two weight sets: u, v and ν̃ use the one with the wall point, ρ and T the one with the neighbours only. Both are exact for linear fields. Far-field faces take no part in any stencil. The second-order scheme extrapolates the centre value along x_f − x_i: W⁺ = W_i + G_i·(x_f − x_i), without a limiter (the flow has no shocks); the first-order scheme takes W⁺ = W_i. The viscous flux and the SA source use these gradients in both schemes, so the first/second-order switch changes the convective flux only; the wall and inflow faces are the exception: their viscous coefficients μ(T) use the reconstructed face temperature T_w, T_f (the cell's at first order).

逐步推导：最小二乘梯度、它的权重与转置

Step by step: the least-squares gradient, its weights and its transpose

第一步，拟合问题。在中心 x_i 附近用线性函数 W_i + G·(x − x_i) 拟合模板点上的值，取使误差平方和最小的 G（W 的每个分量分别做）：

Step 1, the fitting problem. Near the centre x_i, fit the values at the stencil points by the linear function W_i + G·(x − x_i), and take the G that minimises the sum of squared misfits (for each component of W separately):

第二步，法方程。对 G 求梯度并令其为零：

Step 2, the normal equations. Set the gradient with respect to G to zero:

𝒩_i 是 2×2 对称矩阵；只要模板点不全在一条过 x_i 的直线上，它就正定。

𝒩_i is a symmetric 2×2 matrix, positive definite as long as the stencil points do not all lie on one line through x_i.

第三步，显式的权重。2×2 矩阵的逆直接写出（mesh.js 的 weights 正是这几行）：

Step 3, explicit weights. The inverse of a 2×2 matrix written out (mesh.js's weights are exactly these lines):

权重 Γ_ij 只依赖几何，与 W 无关：梯度那一遍是线性算子 Γ，Γ_ij 就是它的元素。

The weights Γ_ij depend on the geometry only, not on W: the gradient pass is the linear operator Γ, and the Γ_ij are its entries.

第四步，线性场精确。若 W 是线性场，W_j − W_i = Δx_ij·∇W 对所有 j 成立，则

Step 4, exact for linear fields. If W is linear, W_j − W_i = Δx_ij·∇W for every j, and then

贴壁单元的两组权重各自满足这个恒等式：u、v、ν̃ 用的一组含壁面点（那里 W_j = 0 已知），ρ、T 用的一组只含邻居。

A wall cell's two weight sets satisfy it each: the one for u, v, ν̃ includes the wall point (where W_j = 0 is known), the one for ρ and T only the neighbours.

第五步，前向与转置。Γ 是线性的，前向就是同一组权重作用在 $\dot{W}$ 上；伴随把它转置，梯度 G_i 的伴随 ${\bar{G}}_i$ 按累加送回整个模板：

Step 5, the tangent and the transpose. Γ is linear, so the tangent applies the same weights to $\dot{W}$; the adjoint transposes it, and the adjoint ${\bar{G}}_i$ of the gradient G_i goes back over the whole stencil as accumulations:

两者是同一个双线性型 ${\sum }_i{\bar{G}}_i\cdot {\sum }_j{\mathrm{\Gamma }}_{ij}({\dot{W}}_j-{\dot{W}}_i)$ 按两种次序求和（core.js applyAT 的 “pass 1 reversed”）；两个方向画在图 3.3 (b) 里。

Both are the same bilinear form ${\sum }_i{\bar{G}}_i\cdot {\sum }_j{\mathrm{\Gamma }}_{ij}({\dot{W}}_j-{\dot{W}}_i)$, summed in two orders (core.js applyAT's “pass 1 reversed”); both directions are drawn in Figure 3.3 (b).

第六步，对坐标的前向导数（第 12 节用）。Δx_ij 随中心移动，$\mathrm{\Delta }{\dot{x}}_{ij}={\dot{x}}_j-{\dot{x}}_i$；由 $({𝒩}^{-1}{)}^{\mathit{\cdot }}=-{𝒩}^{-1}\dot{𝒩}{𝒩}^{-1}$，

Step 6, the forward derivative with respect to the coordinates (used in Section 12). Δx_ij moves with the centres, $\mathrm{\Delta }{\dot{x}}_{ij}={\dot{x}}_j-{\dot{x}}_i$; with $({𝒩}^{-1}{)}^{\mathit{\cdot }}=-{𝒩}^{-1}\dot{𝒩}{𝒩}^{-1}$,

核对：本页 32 组权重（邻居个数：8 个控制体有 2 个、16 个控制体有 3 个，贴壁单元另有壁面点）：第三步的 Γ_ij 与 mesh.js 算出的权重最大相对差 0；Σ_jΓ_ijΔx_ij^T 与 I 相差 5.6×10⁻¹⁶；一个线性场的梯度误差 6.7×10⁻¹⁶；沿一个随机坐标方向，第六步的 ${\dot{\mathrm{\Gamma }}}_{ij}$ 与几何磁带的前向扫描相差 3.5×10⁻¹⁵；随机向量上 $⟨\bar{G},\mathrm{\Gamma }\dot{W}⟩$ 与 $⟨{\mathrm{\Gamma }}^T\bar{G},\dot{W}⟩$ 相差 2.2×10⁻¹⁶（本页代码在构建时计算）。

Check: this page's 32 weight sets (neighbours: 8 control volumes have 2, 16 control volumes have 3, and the wall cells also have the wall point): the Γ_ij of Step 3 and mesh.js's weights differ by at most 0; Σ_jΓ_ijΔx_ij^T differs from I by 5.6×10⁻¹⁶; a linear field's gradient is off by 6.7×10⁻¹⁶; along a random coordinate direction, Step 6's ${\dot{\mathrm{\Gamma }}}_{ij}$ and a forward sweep of the geometry's tape differ by 3.5×10⁻¹⁵; and on random vectors $⟨\bar{G},\mathrm{\Gamma }\dot{W}⟩$ and $⟨{\mathrm{\Gamma }}^T\bar{G},\dot{W}⟩$ differ by 2.2×10⁻¹⁶ (computed by this page's code at build time).

粘性通量：边上的修正平均梯度

The viscous flux: the corrected average gradient on the edge

逐步推导：修正平均梯度与粘性通量的线性化

Step by step: the corrected average gradient and the linearisation of the viscous flux

第一步，拆开平均梯度。对 u、v、T、ν̃ 中的某个分量 φ，记两侧最小二乘梯度的平均为 ∇φ^avg = ½(G_L + G_R)，Δx = x_R − x_L，单位向量 Δx̂ = Δx/|Δx|。把平均梯度拆成沿 Δx 与垂直于 Δx 的两部分：

Step 1, split the average gradient. For one component φ of u, v, T, ν̃, let ∇φ^avg = ½(G_L + G_R) be the mean of the two least-squares gradients, Δx = x_R − x_L, and Δx̂ = Δx/|Δx| its unit vector. Split the mean into its parts along and across Δx:

第二步，替换沿 Δx 的部分。沿 Δx 的方向导数改用两个中心值的差 (φ_R − φ_L)/|Δx|，得到面上的梯度 ∇_fφ（model.js 的 corrected()）：

Step 2, replace the part along Δx. The directional derivative along Δx is replaced by the difference of the two centre values, (φ_R − φ_L)/|Δx|, which gives the face gradient ∇_fφ (model.js corrected()):

第三步，面梯度的导数。∇_fφ 对它的四个输入都是线性的：

Step 3, the derivatives of the face gradient. ∇_fφ is linear in each of its four inputs:

两侧的梯度只经过“垂直于 Δx 的投影”进入；沿 Δx 只剩两点差，所以不会出现棋盘模式。

The two gradients enter only through the projection across Δx; along Δx only the two-point difference remains, so no checkerboard mode appears.

第四步，通量与系数。面上的粘性通量（第 1 节的 F^v·n）读四个面梯度与两侧平均的系数：

Step 4, the flux and its coefficients. The face's viscous flux (Section 1's F^v·n) reads the four face gradients and coefficients averaged over the two sides:

这里 ∇_fu 是以 ∇_fu、∇_fv 为行的 2×2 矩阵，u^avg、μ^avg、μ_t^avg 是两侧的平均（与第一步的 ∇φ^avg 同一个写法）。

Here ∇_fu is the 2×2 matrix with rows ∇_fu and ∇_fv, and u^avg, μ^avg, μ_t^avg are the averages of the two sides (written like Step 1's ∇φ^avg).

第五步，系数的导数。Sutherland 律 μ = μ_ref·1.4T^1.5/(T + 0.4) 给 μ_T = dμ/dT = μ(1.5/T − 1/(T + 0.4))；μ_t = ρν̃f_v1(χ)，χ = ρν̃/μ，f_v1′ = 3c_v1³χ²/(χ³ + c_v1³)²，于是（其余原始变量固定）

Step 5, the derivatives of the coefficients. Sutherland's law μ = μ_ref·1.4T^1.5/(T + 0.4) gives μ_T = dμ/dT = μ(1.5/T − 1/(T + 0.4)); μ_t = ρν̃f_v1(χ) with χ = ρν̃/μ and f_v1′ = 3c_v1³χ²/(χ³ + c_v1³)², so (the other primitives held)

第六步，链式法则（前向模式）。τn 对 μ_e 与四个面梯度都是线性的；面梯度 ∇_f 对它的四个输入也是线性的，作用在带点的输入上就是它的前向导数（第三步），法向 n 固定。所以

Step 6, the chain rule (forward mode). τn is linear in μ_e and in the four face gradients; the face gradient ∇_f is linear in its four inputs as well, so applied to the dotted inputs it is its own forward derivative (Step 3), with the normal n held. So

能量与 ν̃ 两行同理，按乘积法则逐项加点：

The energy and ν̃ rows follow the same way, a dot on each factor in turn by the product rule:

局部磁带记下的正是这些乘积：面核对 W_L、W_R 的导数含系数与两点差两部分，对 G_L、G_R 的导数只含第三步的投影。

These products are exactly what the local tape records: the face kernel's derivative with respect to W_L, W_R has a coefficient part and a two-point part, and with respect to G_L, G_R only Step 3's projection.

核对：本页每个内部面上，corrected() 的磁带导数与第三步的公式相差 2.0×10⁻¹⁶；每个控制体上，viscosity() 与 eddy() 的磁带给出的 μ_T 与第五步三个偏导数与公式相差 6.9×10⁻¹⁶（本页代码在构建时计算）。

Check: on every interior face of this page, the taped derivatives of corrected() and Step 3's formulas differ by 2.0×10⁻¹⁶; in every control volume, μ_T and the three partials of Step 5 from the tapes of viscosity() and eddy() differ from the formulas by 6.9×10⁻¹⁶ (computed by this page's code at build time).

SA 源项

The SA source

源项是每个控制体一个点的积分：V_i 乘以中心处的 S₅。它读三样东西：中心的 W、最小二乘梯度 G（涡量、∇ν̃、∇ρ）、壁面距离 d。

The source is a one-point integral per control volume: V_i times S₅ at the centre. It reads three things: the centre's W, the least-squares gradients G (the vorticity, ∇ν̃, ∇ρ), and the wall distance d.

逐步推导：SA 源项对 W、G 与 d 的导数

Step by step: the derivatives of the SA source with respect to W, G and d

源项核的局部 Jacobian J_s（15 个输入、1 个输出；第 12 节还要它对壁面距离 d 的导数）由局部磁带求出。这里把磁带做的事按前向模式一层一层写出来：每一层给出它的点。控制体 i 的源项进入残差的 ρν̃ 行：

The source kernel's local Jacobian J_s (15 inputs and 1 output; Section 12 also needs its derivative with respect to the wall distance d) comes from the local tape. Here is what the tape does in forward mode, layer by layer: each layer gives its dot. The source of control volume i enters the residual's ρν̃ row:

Q 是第 1 节 S₅ 方括号里的净源项，s = −S₅。输入是中心值 W_i = (ρ, u, v, T, ν̃)、梯度 G_i（10 个分量）与 d_i。

Q is the net source inside the bracket of Section 1's S₅, and s = −S₅. The inputs are the centre values W_i = (ρ, u, v, T, ν̃), the gradients G_i (10 entries) and d_i.

第一步，分子粘性。ν = μ/ρ，μ(T) 是 Sutherland 律：

Step 1, the molecular viscosity. ν = μ/ρ, and μ(T) is Sutherland's law:

第二步，χ = ν̃/ν = ρν̃/μ：

Step 2, χ = ν̃/ν = ρν̃/μ:

第三步，三个阻尼函数只依赖 χ（f_v2′ 用商的求导法则，1 + χf_v1 − χ(f_v1 + χf_v1′) = 1 − χ²f_v1′）：

Step 3, the three damping functions depend on χ only (f_v2′ by the quotient rule, with 1 + χf_v1 − χ(f_v1 + χf_v1′) = 1 − χ²f_v1′):

第四步，B = ν̃f_v2/(κ²d²)。ν̃ 既直接出现，又经 χ 出现，而 ν̃·∂χ/∂ν̃ = χ；ρ 与 T 只经 χ 进入：

Step 4, B = ν̃f_v2/(κ²d²). ν̃ appears directly and through χ, with ν̃·∂χ/∂ν̃ = χ; ρ and T enter through χ only:

第五步，涡量。ω = ∂v/∂x − ∂u/∂y 取自最小二乘梯度 G_i，Ω = |ω|：

Step 5, the vorticity. ω = ∂v/∂x − ∂u/∂y comes from the least-squares gradients G_i, and Ω = |ω|:

第六步，S̃ 的两支。第一支 S̃ = Ω + B：F_B = ∂S̃/∂B = 1，F_Ω = ∂S̃/∂Ω = 1。第二支 S̃ = Ω + Ω(c₂²Ω + c₃B)/Ξ，Ξ = (c₃ − 2c₂)Ω − B，用商的求导法则：

Step 6, the two branches of S̃. First branch, S̃ = Ω + B: F_B = ∂S̃/∂B = 1 and F_Ω = ∂S̃/∂Ω = 1. Second branch, S̃ = Ω + Ω(c₂²Ω + c₃B)/Ξ with Ξ = (c₃ − 2c₂)Ω − B, by the quotient rule:

在接缝 B = −c₂Ω 处 Ξ = (c₃ − c₂)Ω，第二支的分式等于 −c₂Ω：两支的值都是 (1 − c₂)Ω，F_B 与 F_Ω 也都等于 1，所以 S̃ 在接缝处一阶连续。于是 $\dot{\tilde{S}}=F_B\dot{B}+F_{\mathrm{\Omega }}\dot{\mathrm{\Omega }}$。本页求源项的 24 个控制体里，第二支在 8 个上起作用。

At the seam B = −c₂Ω, Ξ = (c₃ − c₂)Ω and the second branch's fraction equals −c₂Ω: both branches give (1 − c₂)Ω, and F_B and F_Ω both equal 1, so S̃ is continuously differentiable across the seam. Hence $\dot{\tilde{S}}=F_B\dot{B}+F_{\mathrm{\Omega }}\dot{\mathrm{\Omega }}$. Of this page's 24 control volumes, the second branch is active in 8.

第七步，r 与截断。r₀ = ν̃/(κ²d²S̃)；未截断时

Step 7, r and its cap. r₀ = ν̃/(κ²d²S̃); uncapped,

截断时 r = 10 是常数，$\dot{r}=0$：那里 f_w 对任何输入都没有导数。本页 7 个控制体被截断。

When capped, r = 10 is a constant and $\dot{r}=0$: there f_w has no derivative with respect to any input. On this page 7 control volumes are capped.

第八步，g 与 f_w。对 f_w = g(1 + c_w3⁶)^1/6(g⁶ + c_w3⁶)^−1/6 求导，(g⁶ + c_w3⁶)^−1/6 − g⁶(g⁶ + c_w3⁶)^−7/6 = c_w3⁶(g⁶ + c_w3⁶)^−7/6：

Step 8, g and f_w. Differentiating f_w = g(1 + c_w3⁶)^1/6(g⁶ + c_w3⁶)^−1/6, with (g⁶ + c_w3⁶)^−1/6 − g⁶(g⁶ + c_w3⁶)^−7/6 = c_w3⁶(g⁶ + c_w3⁶)^−7/6:

第九步，净源项。把以上各层代进 Q（${\dot{f}}_{t2}=f_{t2}^{\mathrm{\prime }}\dot{\chi }$）：

Step 9, the net source. Put the layers above into Q (${\dot{f}}_{t2}=f_{t2}^{\mathrm{\prime }}\dot{\chi }$):

第十步，守恒形式的两项。密度梯度项 a∇ρ·∇ν̃ 与 c_b2 项 ρb|∇ν̃|² 的导数直接写出：

Step 10, the two terms of the conservative form. The derivatives of the density-gradient term a∇ρ·∇ν̃ and the c_b2 term ρb|∇ν̃|² written directly:

第十一步，装配 J_s。s 对每个输入的导数是以上各项之和，其中 −ρQ 对 ρ 的导数是 −Q − ρ∂Q/∂ρ；乘以 V_i 就是 J_s 的一行。16 个输入里有 10 个出现：ρ、T、ν̃、∇ρ 的两个分量、∂u/∂y、∂v/∂x、∇ν̃ 的两个分量与 d；u、v、∂u/∂x、∂v/∂y 与 ∇T 的导数恒为零。各部分在本页的大小见图 3.6。

Step 11, assemble J_s. The derivative of s with respect to each input is the sum of the terms above, where −ρQ contributes −Q − ρ∂Q/∂ρ to the ρ entry; times V_i it is the row J_s. 10 of the 16 inputs appear: ρ, T, ν̃, the two components of ∇ρ, ∂u/∂y, ∂v/∂x, the two components of ∇ν̃, and d; the derivatives with respect to u, v, ∂u/∂x, ∂v/∂y and ∇T are identically zero. Figure 3.6 shows how large each part is on this page.

核对：把上面十一步照写成代码（enrich_data.js 的 saHand，不调用 model.js 的源项），在 U* 的每个求源项的控制体（24 个）上与 model.js source() 的磁带比较：16 个导数的最大相对差 4.6×10⁻¹⁴，源项值 5.6×10⁻¹⁶。其中第二支 8 个、r 截断 7 个、生成项为负（f_t2 > 1）9 个（本页代码在构建时计算）。

Check: the eleven steps written out as code (enrich_data.js saHand, which does not call model.js's source) and the tape of model.js source() compared in every control volume where the source is evaluated (24 of them) at U*: the 16 derivatives differ by at most 4.6×10⁻¹⁴ (relative), the source's value by 5.6×10⁻¹⁶. Among them 8 are on the second branch, 7 capped, and 9 have negative production (f_t2 > 1) (computed by this page's code at build time).

第三、六、七、八步的公式与接缝处的一阶连续另由 sympy 符号验证（check_enrich_sympy.py）。

The formulas of Steps 3, 6, 7 and 8 and the continuity at the seam are also verified symbolically by sympy (check_enrich_sympy.py).

到这里，内部面已经完整：两侧都有控制体，面循环读两侧的状态与梯度。边界上的面只有一侧；物面、来流与出口的条件怎样写进通量，是第 4 节的内容。

Interior faces are now complete: they have a control volume on each side, and the face loop reads the states and gradients of both. A boundary face has one side only; how the conditions of the wall, the inflow and the outlet enter its flux is the subject of Section 4.

边界条件

Boundary conditions

内圈：弱条件的无滑移绝热壁

Inner ring: the no-slip adiabatic wall, weakly

壁面面的通量里没有质量、能量与 ν̃ 的对流（u = 0）；动量通量是壁面压力 p_w = ρ_wT_w（ρ_w、T_w 是单元自己重构到壁面的值；一阶格式直接取单元值）减去粘性应力；能量通量为零（绝热，且 u_w = 0 时应力不做功）；ν̃ 方程有扩散通量 ρa ∇ν̃·n 流向壁面，因为壁面上 ν̃ = 0。

The wall face's flux carries no convection of mass, energy or ν̃ (u = 0); the momentum flux is the wall pressure p_w = ρ_wT_w (ρ_w, T_w are the cell's own reconstruction to the wall; the first-order scheme takes the cell values) minus the viscous stress; the energy flux is zero (adiabatic, and the stress does no work where u_w = 0); the ν̃ equation has a diffusive flux ρa ∇ν̃·n into the wall, since ν̃ = 0 there.

外圈：远场

Outer ring: the far field

来流面的 Roe 通量右侧放虚拟状态 W_g = (ρ∞, u∞, v∞, T_f, ν̃∞)，T_f 是面状态自己的温度（T 不给定）：入射的熵波、剪切波、声波与 ν̃ 波取 W_g 的值，出射声波取单元的值，所以不强加压力。粘性通量沿中心到面点取两点差，面点上 u、v、ν̃ 是来流值；没有热流，功项用 u∞。h → 0 时粘性通量把 u、v、ν̃ 固定在来流值，Roe 通量的熵波把 ρ 也固定下来（两侧温度相同时它的强度是 (1 − 1/γ)(ρ∞ − ρ)）。出口面：面状态的对流通量，压力换成 p∞ = 1，能量通量里是 p∞u·n（牵引力做的功），没有粘性通量。

An inflow face puts the ghost state W_g = (ρ∞, u∞, v∞, T_f, ν̃∞) on the right of the Roe flux, T_f being the face state's own temperature (T is not prescribed): the incoming entropy, shear, acoustic and ν̃ waves take W_g's values and the outgoing acoustic wave the cell's, so no pressure is imposed. The viscous flux takes the two-point difference from the centre to the face point, where u, v and ν̃ are the free stream's; there is no heat flux, and the work term uses u∞. As h → 0 the viscous flux holds u, v and ν̃ at the free stream, and the entropy wave of the Roe flux holds ρ as well (with equal temperatures on the two sides its strength is (1 − 1/γ)(ρ∞ − ρ)). An outlet face: the face state's convective flux with its pressure replaced by p∞ = 1, p∞u·n in the energy flux (the work of the traction), and nothing viscous.

在解处冻结的分支

The branches, frozen at the solution

离散残差里有几处不可导点。导数（前向、伴随、复数步长）都在当前解所在的那一侧分支上取；只要解离分支点足够远，差分的步长就跨不过去。下表列出二阶格式在收敛解处的距离：

The discrete residual has a few points where it is not differentiable. Every derivative (tangent, adjoint, complex step) is taken on the branch the solution sits on; as long as the solution is far enough from a switch, a finite-difference step cannot cross it. The table lists the distances at the converged second-order solution:

修正分支在 8 个控制体上起作用，r 在 7 个控制体上被截断：两个分支都真的在用。

The modified branch is active in 8 control volumes and r is capped in 7: both branches are really in use.

下面的演示把这些分支画出来：f_v1(χ)、f_w(r) 与它在 r = 10 的截断、S̃/Ω 随 B/Ω 的两支，并把当前解（页首所选的阶）每个计算源项的控制体的值标在曲线上。点离开分支点有多远，就是上表的余量。

The demo below draws these branches: f_v1(χ), f_w(r) with its cap at r = 10, and the two branches of S̃/Ω against B/Ω, with every control volume of the current solution (at the order chosen in the header) whose source is evaluated marked on the curves. How far the points sit from a switch is the margin of the table above.

残差 = 几遍循环：原始变量、梯度、面、源项

The residual = passes: primitives, gradients, faces, source

标量页面的残差是一个面循环。这里一个面要读两侧的梯度，而梯度要读整个模板的值，所以残差分成几遍：先把所有控制体的原始变量算出来，再算所有梯度，然后才是面通量和源项。

On the scalar pages the residual is one face loop. Here a face reads the gradients on both sides, and a gradient reads the values over its whole stencil, so the residual comes in passes: first the primitive variables of every control volume, then every gradient, and only then the face fluxes and the source.

π 在每个控制体上独立地把守恒变量换成原始变量；Γ 是线性的最小二乘算子（权重只依赖几何）；F 是面通量的装配，S 是源项。下面的演示逐句执行这几遍：410 步，可以用上方的按钮直接跳到某一遍。

π turns each control volume's conserved variables into primitives, independently; Γ is the linear least-squares operator (its weights depend on the geometry only); F assembles the face fluxes and S the source. The demo below executes the passes statement by statement: 410 steps, and the buttons above it jump straight to a pass.

最后一步 max|R| 在舍入误差量级：U* 是收敛解。

At the last step max|R| is at round-off: U* is the converged solution.

求解：Jacobian、近似 Jacobian 与 Gauss–Seidel

The solve: the Jacobian, the approximate Jacobian and Gauss–Seidel

从残差到 Jacobian

From the residual to the Jacobian

A = ∂R/∂U 是 120×120 的矩阵，由 5×5 的块组成。一个面的通量读两侧的值与梯度，梯度又读整个模板，所以 A 的一行不只连到邻居，还连到邻居的邻居：200 个非零块。

A = ∂R/∂U is a 120×120 matrix of 5×5 blocks. A face's flux reads the values and gradients on both sides, and a gradient reads its whole stencil, so a row of A reaches past the neighbours to the neighbours' neighbours: 200 non-zero blocks.

近似 Jacobian：一阶 Roe 加薄剪切层

The approximate Jacobian: first-order Roe plus the thin shear layer

Gauss–Seidel 用的不是 A，而是一个只连接面两侧的近似 P：对流部分是一阶 Roe（J^c(U) = ∂(F^c(U)·n̂)/∂U 是物理通量的 Jacobian，|J̃^c| = Σ_k|Λ_k|r_kl_k^T 是它在两侧中心状态的 Roe 平均处取绝对值，l_k 是左特征向量），粘性部分是薄剪切层（只保留修正梯度沿 Δx 的两点差分，系数冻结）：

Gauss–Seidel does not use A but an approximation P that couples only the two sides of each face: first-order Roe for the convective part (J^c(U) = ∂(F^c(U)·n̂)/∂U is the Jacobian of the physical flux, and |J̃^c| = Σ_k|Λ_k|r_kl_k^T its absolute value at the Roe average of the two centre states, l_k being the left eigenvectors), and the thin shear layer for the viscous part (only the two-point difference along Δx of the corrected gradient, with frozen coefficients):

Θ_f 是薄剪切层矩阵乘上两点系数 θ_f = n_f·Δx/|Δx|²：动量行 θ_fμ_e(I + ⅓n̂n̂^T)，能量行加上功与导热 θ_fk_e，ν̃ 行 θ_fρa。(ρν̃, ρν̃) 对角元再加上 SA 源项对 ρν̃ 的导数，只在它为正（耗散占优、使对角增强）时加；对角块还加上伪时间项 V/Δt。P 有 88 个非零块。

Θ_f is the thin-shear-layer matrix times the two-point coefficient θ_f = n_f·Δx/|Δx|²: the momentum rows θ_fμ_e(I + ⅓n̂n̂^T), the energy row adds the work and the conduction θ_fk_e, and the ν̃ row θ_fρa. The (ρν̃, ρν̃) entry also receives the SA source's derivative with respect to ρν̃ where it is positive (destruction dominating, so it strengthens the diagonal), and the diagonal blocks the pseudo-time term V/Δt. P has 88 non-zero blocks.

逐步推导：薄剪切层矩阵从哪里来

Step by step: where the thin-shear-layer matrix comes from

第一步，只留两点部分。在第 3 节的面梯度里只保留沿 Δx 的两点差，∇_fφ ≈ (φ_R − φ_L)Δx/|Δx|²，并冻结所有系数。于是速度的面梯度矩阵是 ∇_fu ≈ Δu Δx^T/|Δx|²（Δu = u_R − u_L）。

Step 1, keep only the two-point part. In the face gradient of Section 3 keep only the two-point difference along Δx, ∇_fφ ≈ (φ_R − φ_L)Δx/|Δx|², and freeze every coefficient. The velocity's face-gradient matrix is then ∇_fu ≈ Δu Δx^T/|Δx|² (Δu = u_R − u_L).

第二步，代入应力：

Step 2, put it into the stress:

第三步，Δx 平行于 n 时。Δx = |Δx|n̂、n = |n|n̂，两点系数 θ_f = n·Δx/|Δx|² = |n|/|Δx|，第二步化为

Step 3, when Δx is parallel to n. With Δx = |Δx|n̂ and n = |n|n̂, the two-point coefficient is θ_f = n·Δx/|Δx|² = |n|/|Δx|, and Step 2 becomes

第四步，能量与 ν̃ 两行。功 u^avg·τn 用第三步的 τn，导热与 ν̃ 的扩散也只剩两点差（u_n^avg = u^avg·n̂）：

Step 4, the energy and ν̃ rows. The work u^avg·τn uses Step 3's τn, and the conduction and the ν̃ diffusion keep their two-point differences only (u_n^avg = u^avg·n̂):