方程：可压缩 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 页面解同一个连续问题：格心、格点两种控制体，各有两种来流边界的离散；同一种控制体的两个页面只在来流边界的离散上不同。本页沿用原来的处理：来流节点的 ρ、ρu、ρv、ρν̃ 四行换成 U − U∞；伴随一致的版本里，来流节点保留质量行，只以节点自己的 ρ 固定 u、v、ν̃。两者都给出各自离散问题的精确导数；区别在于加密时离散伴随在来流附近是否趋于连续伴随（伴随一致性）。做法见第 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 ρ, ρu, ρv and ρν̃ rows of an inflow node are replaced by U − U∞; in the adjoint-consistent version an inflow node keeps its mass row and fixes u, v and ν̃ at its own ρ. 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 = 54.2003，这张粗网格上的阻力有多准见第 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 = 54.2003; 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: dual control volumes, dual 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.

网格与二维对流页面完全相同，编号也相同：节点逐圈从壁面向外，单元逐圈编号，面按“内部、壁面、远场”排序。未知量放在节点上：每个节点的控制体是中位对偶，由相邻单元的边中点与单元中心围成。区别在于每个控制体现在有 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 at the nodes: each node's control volume is its median dual, bounded by the edge midpoints and the centres of the cells around it. 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 = X_i，所以 ∂d_i/∂X_i = n_w；X_a、X_b 是壁面节点，各得到第三步的一项。壁面节点自己的 d = 0，不求源项，也不求导。

Step 5, back to the node coordinates. On the node page the centre is the node itself, x_i = X_i, so ∂d_i/∂X_i = n_w; X_a and X_b are wall nodes and receive the terms of Step 3. A wall node's own d = 0: it has no source and no derivative.

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

Check: of this page's 16 centres, 2 have their nearest point inside an edge and 14 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 4.6×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 节）。格点页壁面节点的 ρu、ρv、ρν̃ 三行与来流节点的 ρ、ρu、ρv、ρν̃ 四行不是上面的常微分方程：它们被替换成边界条件本身（第 4 节），推进中始终精确满足 u = v = ν̃ = 0 与 ρ、u、v、ν̃ 的来流值。

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). The ρu, ρv and ρν̃ rows of a wall node and the ρ, ρu, ρv and ρν̃ rows of an inflow node are not ODEs like the one above: they are replaced by the boundary conditions themselves (Section 4), and u = v = ν̃ = 0 and the inflow values of ρ, u, v, ν̃ hold exactly throughout the marching.

每个面仍是一个一维问题：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.

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

Check: on all 48 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 3.2×10⁻¹⁶; Σ_kΔ_kr_k rebuilds the jump ΔU to 1.2×10⁻¹⁴; 𝓡|Λ|𝓡⁻¹ΔU agrees with model.js's dissipation to 7.2×10⁻¹⁵; Roe's property F(U_R) − F(U_L) = J^c(W̃)ΔU holds to 1.3×10⁻¹⁴; and the ρν̃ row and column of |J̃^c| follow Step 4 to 3.3×10⁻¹⁸. The last term dropped in Step 5 (first order, centre states) is at most 32.3% of the taped ∂Φ/∂U_L (Frobenius norm), 8.9% 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），所以不需要额外的壁面点。二阶格式把中心值沿 x_f − x_i 外推：W⁺ = W_i + G_i·(x_f − x_i)，不加限制器（流动无激波）；一阶格式直接取 W⁺ = W_i。粘性通量与 SA 源项在两种格式里都用这些梯度，所以一阶／二阶开关只改变对流通量。

The stencil is the edge neighbours. A wall node is itself an unknown (u = v = ν̃ = 0 under the strong condition), so no extra wall point is needed. 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.

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

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、ν̃ 本身就是未知量）。

Every node's weight set satisfies it, the wall nodes' too (their u, v, ν̃ are unknowns themselves).

第五步，前向与转置。Γ 是线性的，前向就是同一组权重作用在 $\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}$,

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

Check: this page's 24 weight sets (neighbours: 8 control volumes have 3, 8 control volumes have 4, 8 control volumes have 5): the Γ_ij of Step 3 and mesh.js's weights differ by at most 0; Σ_jΓ_ijΔx_ij^T differs from I by 2.2×10⁻¹⁶; a linear field's gradient is off by 3.3×10⁻¹⁶; along a random coordinate direction, Step 6's ${\dot{\mathrm{\Gamma }}}_{ij}$ and a forward sweep of the geometry's tape differ by 5.2×10⁻¹⁵; and on random vectors $⟨\bar{G},\mathrm{\Gamma }\dot{W}⟩$ and $⟨{\mathrm{\Gamma }}^T\bar{G},\dot{W}⟩$ differ by 3.8×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.2×10⁻¹⁶；每个控制体上，viscosity() 与 eddy() 的磁带给出的 μ_T 与第五步三个偏导数与公式相差 7.2×10⁻¹⁶（本页代码在构建时计算）。

Check: on every interior face of this page, the taped derivatives of corrected() and Step 3's formulas differ by 2.2×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 7.2×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 }}$。本页求源项的 11 个控制体里，第二支在 0 个上起作用。

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 11 control volumes, the second branch is active in 0.

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

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

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

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 1 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* 的每个求源项的控制体（11 个）上与 model.js source() 的磁带比较：16 个导数的最大相对差 4.9×10⁻¹⁵，源项值 2.8×10⁻¹⁵。其中第二支 0 个、r 截断 1 个、生成项为负（f_t2 > 1）4 个（本页代码在构建时计算）。

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 (11 of them) at U*: the 16 derivatives differ by at most 4.9×10⁻¹⁵ (relative), the source's value by 2.8×10⁻¹⁵. Among them 0 are on the second branch, 1 capped, and 4 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, and which rows they replace, is the subject of Section 4.

边界条件

Boundary conditions

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

Inner ring: the no-slip adiabatic wall, strongly

壁面节点的 ρu、ρv、ρν̃ 三行在所有循环之后被替换成 U 本身，于是解满足 u = v = ν̃ = 0。连续方程与能量方程不能替换，它们仍是通量平衡。壁面半面对它们没有贡献：质量通量 ρu·n = 0，能量通量为零（绝热，且应力不做功），所以面循环跳过壁面半面。

The ρu, ρv and ρν̃ rows of a wall node are replaced by U itself after all the loops, so the solution satisfies u = v = ν̃ = 0. The continuity and energy equations cannot be replaced; they remain flux balances. The wall half-faces add nothing to them: the mass flux ρu·n = 0 and the energy flux is zero (adiabatic, and the stress does no work), so the face loop skips the wall half-faces.

外圈：远场

Outer ring: the far field

来流节点的 ρ、ρu、ρv、ρν̃ 四行在循环之后替换成 U − U∞（强条件）；能量行保留，来流半面上的通量是本节点状态的对流通量减去本节点梯度给出的粘性通量，没有热流。出口半面：本节点的状态，压力换成 p∞ = 1，没有粘性通量。一个节点只要有一个来流半面就是来流节点，所以两个交界节点也是来流节点，它们的出口半面只进入能量行。

The ρ, ρu, ρv and ρν̃ rows of an inflow node are replaced by U − U∞ after the loops (a strong condition); the energy row stays, and the flux through an inflow half-face is the node state's convective flux minus the viscous flux of the node's gradients, with no heat flux. An outlet half-face: the node's own state with its pressure replaced by p∞ = 1, and nothing viscous. A node with one inflow half-face is an inflow node, so the two junction nodes are inflow nodes too, and their outlet half-faces enter their energy rows only.

在解处冻结的分支

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:

r 在 1 个控制体上被截断。修正分支只出现在 5 个来流节点上，而这些节点的 ν̃ 行被替换、源项不计算，所以本页二阶格式的残差用不到它（上表的余量也只在计算源项的控制体上取）。格点页一阶格式有一处例外：两个壁面节点之间的面，两侧速度都正好是 0，所以 q̃ = 0，正好落在 |q̃| 的不可导点上。导数取 |q̃| 的正分支；这些导数项只与壁面节点的速度相乘，而任何前向解 v 在壁面节点的速度分量上恒为零（它们的行被替换），所以设计导数与分支的选法无关，只有壁面节点动量分量的伴随值依赖它。

r is capped in 1 control volume. The modified branch occurs only at 5 inflow nodes, whose ν̃ rows are replaced and whose source is not evaluated, so this page's second-order residual does not use it (the margins of the table are taken over the control volumes whose source is evaluated). The node page at first order has one exception: on a face between two wall nodes both velocities are exactly 0, so q̃ = 0 and the face sits exactly on the kink of |q̃|. The derivative takes the positive branch of |q̃|; these terms multiply only the wall nodes' velocities, and every tangent solution v is zero in the wall nodes' velocity components (their rows are replaced), so the design derivatives do not depend on the choice; only the adjoint of the wall nodes' momentum components does.

下面的演示把这些分支画出来：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 (the wall and inflow nodes, whose ν̃ rows are replaced, are left out). How far the points sit from a switch is the margin of the table above.

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

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

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

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 是源项；格点页最后还有第 4 遍，替换壁面节点与来流节点的强边界行。下面的演示逐句执行这几遍：538 步，可以用上方的按钮直接跳到某一遍。

π 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: 538 steps, and the buttons above it jump straight to a pass.

最后一步 max|R| 在舍入误差量级：U* 是收敛解。注意第 4 遍：壁面节点的三行在所有循环之后才被替换，面循环写进这些行的值全部被覆盖；来流节点的四行也在这一遍替换。第 9 节的伴随要先处理这件事。

At the last step max|R| is at round-off: U* is the converged solution. Note pass 4: the three rows of each wall node are replaced after all the loops, overwriting whatever the face loop wrote there; the four rows of each inflow node are replaced in this pass too. The adjoint of Section 9 has to deal with this first.

求解：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 的一行不只连到邻居，还连到邻居的邻居：264 个非零块。

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: 264 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 有 120 个非零块。

Θ_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. The three rows of a wall node and the four rows of an inflow node are unit rows. P has 120 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: