方程：二维可压缩层流 Navier–Stokes 方程

The equations: two-dimensional compressible laminar Navier–Stokes

从守恒律到本页离散的方程

From conservation laws to the equations discretised here

二维可压缩 Navier–Stokes 方程由流体的三条守恒律组成：质量守恒；动量守恒，即流体微团的 Newton 第二定律，力来自压力和黏性应力；总能量守恒，即热力学第一定律，热量由导热带入，功由压力和黏性应力做出。平面流动中动量有 $x$、$y$ 两个分量，所以共有四个方程、四个守恒量。本页研究层流，不加湍流模型。把单位体积内的四个守恒量排成一列：

The two-dimensional compressible Navier–Stokes equations are the three conservation laws of a fluid: conservation of mass; of momentum, which is Newton's second law for a fluid parcel, with forces from pressure and from viscous stress; and of total energy, the first law of thermodynamics, with heat brought in by conduction and work done by pressure and viscous stress. In plane flow the momentum has an $x$ and a $y$ component, so there are four equations for four conserved quantities. This page studies laminar flow and adds no turbulence model. Stack the four conserved quantities per unit volume into one column:

守恒形式写成“守恒量的变化率加通量的散度等于源项”。通量分成两部分：上标 c 的对流（无黏）部分只依赖状态本身，与 Euler 方程相同；上标 v 的黏性部分由黏性应力和导热组成，还依赖速度和温度的梯度：

The conservation form reads: the rate of change of the conserved quantities plus the divergence of the flux equals the source. The flux has two parts: the convective (inviscid) part, superscript c, depends on the state alone and is the same as in the Euler equations; the viscous part, superscript v, comes from viscous stress and heat conduction and also depends on the gradients of velocity and temperature:

这里 $F^c$、$F^v$ 各有 $x$、$y$ 两列，乘以单位法向 $n=(n_x,n_y)$ 就得到穿过法向为 $n$ 的面的通量 $F_n^c$、$F_n^v$，下面三小节逐行写出。在任意一块固定区域 $V$ 上积分，散度定理把散度的面积分变成边界上的法向通量。本页离散的是这个积分形式的定常情形：时间导数为零，每个控制体上净流出的通量等于体内的源项。

Here $F^c$ and $F^v$ each have an $x$ and a $y$ column; multiplying by a unit normal $n=(n_x,n_y)$ gives the fluxes $F_n^c$ and $F_n^v$ through a face of normal $n$, written row by row in the next three subsections. Integrated over any fixed region $V$, the divergence theorem turns the area integral of the divergence into the normal flux through the boundary. This page discretises the steady case of this integral form: the time derivative is zero, and on every control volume the net outward flux equals the source inside.

物理绕流中源项 $S$ 只有一项可选的附加热源：能量方程中的 $a_5\sin^2(0.7x+0.5y)$，其幅值 $a_5$ 是后面的设计参数之一，默认为 0；制造解算例中 $S$ 另含一项解析源（见本节“本页的问题”）。第 2 节选定控制体，第 3 节在每个面上计算法向通量，第 5 节把它们组装成残差。

In the physical flow the source $S$ is only an optional added heating, the term $a_5\sin^2(0.7x+0.5y)$ of the energy equation, whose amplitude $a_5$ is one of the design parameters below (0 by default); in the manufactured case $S$ also contains an analytic source (see “The problem on this page” below). Section 2 chooses the control volumes, Section 3 computes the normal flux on each face and Section 5 assembles the fluxes into the residual.

图 1.1 画出本页的物理绕流问题：计算区域是物体（琥珀色，静止固壁）与截断的外边界之间的一圈；外边界左侧的蓝色实线是入口，右侧的绿色虚线是出口。箭头是本页用自己的代码在默认设置下解出的速度：流体从左侧流入，在物体前方减速，从两侧绕过物体，速度在壁面上为零。右侧列出设计变量 $D=(a,X)$ 的当前值和每种边界给定的数据，下文逐项解释。

Figure 1.1 draws the physical-flow problem of this page. The domain is the ring between the body (amber, a wall at rest) and the truncated outer boundary; on the outer boundary the solid blue line on the left is the inlet and the green dashes on the right are the outlet. The arrows are the velocity this page's own code computes at the default settings: the fluid enters on the left, slows down in front of the body and passes it on both sides, and the velocity is zero at the wall. On the right are the current values of the design variables $D=(a,X)$ and the data each kind of boundary takes, explained below.

状态、状态方程与物性

The state, the equation of state and the transport properties

先把四个量的物理含义分开。$\rho$ 是质量密度，$m_x=\rho u,m_y=\rho v$ 是两个动量密度，$e=\rho E$ 是总能量密度；这里 $e$ 不是比内能。理想气体的比总能量由比内能与比动能组成。先减去动能，再除以 $c_v$，才得到温度：

Separate the four physical quantities first. $\rho$ is mass density, $m_x=\rho u,m_y=\rho v$ are momentum densities, and $e=\rho E$ is total energy per volume; $e$ is not specific internal energy. Subtract specific kinetic energy from specific total energy, then divide by $c_v$ to recover temperature:

采用无量纲气体常数 1，$\gamma=1.4,c_v=2.5,c_p=3.5$。压力 $p=\rho T$，这就是无量纲的理想气体状态方程，且 $\gamma-1=1/c_v$。声速为 $c=\sqrt{\gamma p/\rho}=\sqrt{\gamma T}$，温度为 1 的气体声速是 $\sqrt{1.4}$。

Use nondimensional gas constant one and $\gamma=1.4,c_v=2.5,c_p=3.5$, with $p=\rho T$: the ideal-gas equation of state in nondimensional form, with $\gamma-1=1/c_v$. The speed of sound is $c=\sqrt{\gamma p/\rho}=\sqrt{\gamma T}$, so gas at temperature 1 has the speed of sound $\sqrt{1.4}$.

黏度 $\mu$ 与导热系数 $\kappa$ 是气体的物性。本页可以取常物性，也可以取随温度变化的 Sutherland 律（参考温度 1、常数 0.4，并归一化使 $T=1$ 时 $\mu=\mu_{\rm ref}$）；Prandtl 数 $Pr$ 把导热系数与黏度联系起来：

The viscosity $\mu$ and the thermal conductivity $\kappa$ are properties of the gas. The page takes either constant transport or Sutherland's law, which varies with temperature (reference temperature 1, constant 0.4, normalised so that $\mu=\mu_{\rm ref}$ at $T=1$); the Prandtl number $Pr$ ties the conductivity to the viscosity:

第 3 节的逐步推导对 Sutherland 律求导，并画出它在本页温度范围内的形状。

A step-by-step derivation in Section 3 differentiates Sutherland's law and plots it over this page's range of temperatures.

对流通量：质量、力和功率

The convective flux: mass, force and power

质量穿过单位面积面的速率是 $\rho u_n$，其中 $u_n$ 是法向速度；它携带的 x、y 动量分别乘以 $u,v$。压力又在面上施加法向力 $pn$，压力功率为 $p u_n$。因此能量对流项不是单独的 $e u_n$，而是 $(e+p)u_n$：

Mass crosses unit face area at rate $\rho u_n$, where $u_n$ is the normal velocity, carrying x/y momentum at that rate times $u,v$. Pressure contributes normal force $pn$ and power $pu_n$, giving energy convection $(e+p)u_n$:

对流通量只含状态本身，不含导数；这四行正是二维 Euler 方程的通量。上标 c、v 分别代表对流和黏性。

The convective flux contains the state alone and no derivatives; these four rows are exactly the flux of the two-dimensional Euler equations. Superscripts c and v denote convection and viscosity.

黏性应力与热流

Viscous stress and heat flux

相邻流体层以不同速度运动时，黏性应力阻碍它们相对滑动；温度不均匀时，热量从热处流向冷处。Newton 流体的黏性应力与速度梯度成正比，平面运动中为

When neighbouring layers of fluid move at different speeds, viscous stress resists their sliding past each other; where the temperature is uneven, heat flows from hot to cold. In a Newtonian fluid the viscous stress is proportional to the velocity gradient; in plane motion

这里的 $u_x$ 等表示空间导数；第 3 节在每个面上用重构得到的离散梯度 $G_x^u$ 等替代它们。平面运动仍是三维流体本构的限制：$w=0,\partial_z=0$，但通常 $\tau_{zz}=-(2/3)\mu(u_x+v_y)$ 不为零。这解释了系数为何是 $2/3$，而不是把空间维数改成二就改成另一系数。$\tau n$ 是面上的牵引（Cauchy 公式），既有法向分量，也有沿面的剪切分量。

Here $u_x$ and the like are spatial derivatives; Section 3 replaces them on every face by reconstructed discrete gradients $G_x^u$ and the like. Plane motion restricts a three-dimensional constitutive law: $w=0,\partial_z=0$, but generally $\tau_{zz}=-(2/3)\mu(u_x+v_y)\ne0$. This explains the factor $2/3$: it does not change to another coefficient just because the space has two dimensions. $\tau n$ is the traction on the face (Cauchy's formula), with a normal part and a shear part along the face.

导热服从 Fourier 定律 $q=-\kappa\nabla T$，穿过面的法向热流为 $q_n=-\kappa\nabla T\cdot n$。黏性应力做功为 $(u,v)\cdot\tau n$，导热携带能量而不携带质量。这给出黏性通量：

Heat conduction follows Fourier's law $q=-\kappa\nabla T$, and the normal heat flux through a face is $q_n=-\kappa\nabla T\cdot n$. Viscous traction does work $(u,v)\cdot\tau n$; heat transports energy but no mass. This gives the viscous flux:

黏性功与热流都进入能量方程。面上的总通量为对流减黏性，热流在总能量通量中以正号出现。不能把黏性应力当成四个互不相关的标量扩散项：同一个张量 $\tau$ 同时出现在两个动量行和能量行中。

Both viscous work and heat enter energy. Net flux is convection minus viscosity, so heat enters net energy with a plus sign. Viscous stress is not four independent scalar diffusion terms: the one tensor $\tau$ appears in both momentum rows and in the energy row.

图 1.2 在法向 $n=(1,0)$ 的面上画出这些量的方向：牵引 $\tau n=(\tau_{xx},\tau_{xy})^{\top}$ 的法向分量与剪切分量，以及能量行中的黏性功与热流。

Figure 1.2 draws their directions on a face of normal $n=(1,0)$: the normal and shear parts of the traction $\tau n=(\tau_{xx},\tau_{xy})^{\top}$, and the viscous work and heat flux of the energy row.

图 1.3 在本页物理绕流基准的一个实际面点上，把这两部分通量逐行逐项算出来。质量行没有黏性部分；两个动量行中压力项最大，这是低 Mach 数流动的特点；黏性部分要用到速度和温度的梯度，所以第 3 节必须先从存储的量重构出梯度。

Figure 1.3 evaluates both parts of the flux, row by row and term by term, at a real face point of this page's physical-flow baseline. The mass row has no viscous part; the pressure term is the largest in both momentum rows, as is typical at low Mach number; and the viscous part needs the gradients of velocity and temperature, which is why Section 3 must first reconstruct gradients from the stored values.

黏性改变了什么：无滑移、热壁与边界条件的个数

What viscosity changes: no slip, the thermal wall and the number of boundary conditions

没有黏性项时，方程只含一阶导数，就是 Euler 方程组。它是双曲型的：小扰动沿特征传播，法向速度为 $u_n$ 时四个特征速度是 $u_n-c$、$u_n$、$u_n$、$u_n+c$。一段边界要给几个条件，等于从外部进入区域的特征有几条：亚声速入口 3 个，亚声速出口 1 个，固壁 1 个（只要求 $u_n=0$，流体可以沿壁面滑动）。

Without the viscous terms the equations contain only first derivatives: they are the Euler equations. They are hyperbolic: small disturbances travel along characteristics, and with normal velocity $u_n$ the four characteristic speeds are $u_n-c$, $u_n$, $u_n$ and $u_n+c$. A boundary takes as many conditions as there are characteristics entering the domain through it: 3 at a subsonic inlet, 1 at a subsonic outlet and 1 at a wall (only $u_n=0$; the fluid may slide along the wall).

黏性项和导热项带来速度与温度的二阶导数，动量方程和能量方程于是像扩散方程一样：在每一段边界上，$u$、$v$、$T$ 都要给一个条件，给它本身，或给与之对应的通量（牵引、热流），每段边界 3 个。质量方程没有扩散项，仍是一阶的：密度只随流动输运，只在流体流入的边界上给定。所以二维可压缩 Navier–Stokes 方程是混合的双曲–抛物型方程组。本页的边界代码按这个个数闭合：

The viscous and conduction terms bring second derivatives of velocity and temperature, so the momentum and energy equations behave like diffusion equations: on every part of the boundary each of $u$, $v$ and $T$ needs one condition, on the value itself or on the matching flux (the traction, the heat flux), three per boundary. The mass equation has no diffusion term and stays first order: density is only carried by the flow and is prescribed only where fluid enters. The two-dimensional compressible Navier–Stokes equations are therefore a mixed hyperbolic–parabolic system. This page's boundary code closes every boundary with that many conditions:

入口给 4 个：$\rho$、$u$、$v$、$T$ 全部给定。出口给 3 个：总牵引 $t_b=pn-\tau n$ 的两个分量和法向热流 $q_b$，密度和速度取内部的值。固壁给 3 个：无滑移 $u=v=0$，以及等温 $T=T_w$ 或绝热 $q_n=0$ 二者之一。图 1.4 用本页基准解的边界状态比较这两种计数。

The inlet takes 4: $\rho$, $u$, $v$ and $T$ are all prescribed. The outlet takes 3: the two components of the total traction $t_b=pn-\tau n$ and the normal heat flux $q_b$, while the density and the velocity are taken from the interior. The wall takes 3: no slip, $u=v=0$, and either an isothermal $T=T_w$ or an adiabatic $q_n=0$. Figure 1.4 compares the two counts with the boundary states of this page's baseline.

无滑移条件要求流体在壁面上与壁面同速（本页固壁静止）：切向速度也为零，这是它与 Euler 滑移壁的区别。壁面附近的速度因此从零逐渐增加，形成速度边界层；速度梯度在壁面上产生剪切应力，也就是摩擦阻力。能量方程在壁面上同样要一个条件：等温壁给定温度（本页默认 $T_w=1.05$，比来流的 $T=1$ 略热，热量从壁面传入气体），绝热壁给定零热流。本页的“固壁热学”选项在两者之间切换；第 4 节把这些条件写进边界面的通量。

The no-slip condition makes the fluid at the wall move with the wall (at rest on this page): the tangential velocity vanishes too, which is what distinguishes it from the slip wall of the Euler equations. Near the wall the velocity therefore rises from zero across a velocity boundary layer, and the velocity gradient at the wall produces a shear stress, the skin friction. The energy equation also takes one condition at the wall: an isothermal wall prescribes the temperature (by default $T_w=1.05$, slightly warmer than the incoming $T=1$, so heat flows from the wall into the gas), an adiabatic wall prescribes zero heat flux. The page's “Wall thermal” control switches between the two; Section 4 writes these conditions into the boundary-face fluxes.

本页的问题：区域、边界数据、设计变量与基准解

The problem on this page: domain, boundary data, design variables and baseline

计算区域是物体与外边界之间的一圈（图 1.1）。物体是一个 8 边形，顶点到原点的距离在 0.28 到 0.80 之间，直径（最远两顶点的距离）为 1.23；外边界是半径 1.9 的正 8 边形。真实物体周围的流场延伸到无穷远，这里在离物体不远处截断，截断处用入口和出口条件代替。网格与二维标量页相同（第 2 节）。

The domain is the ring between the body and the outer boundary (Figure 1.1). The body is a polygon with 8 vertices at distances 0.28 to 0.80 from the origin, of diameter 1.23 (the largest distance between two vertices); the outer boundary is a regular polygon with 8 vertices at radius 1.9. The real flow around a body extends to infinity; here it is cut off not far from the body, and inlet and outlet conditions stand in for the far field. The mesh is that of the 2-D scalar pages (Section 2).

本页有两个算例。它们用同一组方程和同一张网格，却不是同一个边值问题。物理绕流（控件“层流绕流”）：入口给定 $\rho=1$、$T=1$ 和速度 $a_0\sqrt{\gamma}\,(\cos\,a_1,\ \sin\,a_1)$，其大小是 $a_0$ 倍声速，所以 $a_0$ 就是来流 Mach 数；出口给定 $t_b=n$、$q_b=0$；内圈是静止固壁，$u=v=0$，等温壁温 $T_w=a_4$（或绝热）。

The page has two cases. They share the equations and the mesh but are different boundary-value problems. Physical flow (the control “Laminar flow”): the inlet prescribes $\rho=1$, $T=1$ and the velocity $a_0\sqrt{\gamma}\,(\cos\,a_1,\ \sin\,a_1)$, whose size is $a_0$ times the speed of sound, so $a_0$ is the inflow Mach number; the outlet prescribes $t_b=n$ and $q_b=0$; the inner ring is a wall at rest, $u=v=0$, with the wall temperature $T_w=a_4$ (or adiabatic).

制造解（默认算例）：选一个光滑的解析场 $U_*$（第 15 节给出公式），把它代入方程算出源项 $S_*=\nabla\cdot(F^c-F^v)(U_*)$，于是 $U_*$ 就是加了 $S_*$ 的方程的精确解。两圈边界都按一个固定方向分成入口和出口，数据取自 $U_*$；内圈不是固壁。这个算例检查离散方程、导数和精度（第 15 节测量精度阶）；其中 $a_0$、$a_1$ 扰动入口速度，$a_4$ 是内边界入口温度的偏移。

Manufactured solution (the default case): choose a smooth analytic field $U_*$ (Section 15 gives the formulas), substitute it into the equations to get the source $S_*=\nabla\cdot(F^c-F^v)(U_*)$, and $U_*$ is then the exact solution of the equations with $S_*$ added. Both rings are split into inlet and outlet by one fixed direction, with data from $U_*$; the inner ring is not a wall. This case checks the discrete equations, the derivatives and the accuracy (Section 15 measures the order of accuracy); in it $a_0$ and $a_1$ perturb the inlet velocity and $a_4$ offsets the inlet temperature on the inner boundary.

按伴随优化文献的通常写法，设计变量记作 $D$，本页 $D=(a,X)$。$a=(a_0,\dots,a_5)$ 是六个物理参数，顺序与代码相同：入口幅度（绕流中为 Mach 数，默认 0.12）、入口方向（默认 10°）、参考黏度 $\mu_{\rm ref}$（0.2）、$Pr$（0.72）、壁温 $T_w$（绕流默认 1.05；制造解中默认为 1，内边界的入口温度取 $T_*+(a_4-1)$）和附加热源的幅度 $a_5$（0）；$X$ 是全部 24 个网格顶点的 48 个坐标。合计 54 个设计分量，页首的全页路线图就从这些分量出发；第 12 节推导对它们的导数。

As is usual in the adjoint-optimization literature, the design variables are written $D$; here $D=(a,X)$. $a=(a_0,\dots,a_5)$ are six physical parameters, in the code's order: the inlet amplitude (the Mach number in the flow case, 0.12 by default), the inlet direction (10°), the reference viscosity $\mu_{\rm ref}$ (0.2), $Pr$ (0.72), the wall temperature $T_w$ (1.05 in the flow case; 1 in the manufactured case, where the inner-boundary inlet temperature is $T_*+(a_4-1)$) and the amplitude $a_5$ of the added heating (0); $X$ holds the 48 coordinates of all 24 mesh vertices. That makes 54 design components, and the page roadmap at the top starts from them; Section 12 derives the derivatives with respect to them.

在默认设置下（二阶对流、常物性、等温壁、$\mu_{\rm ref}=0.2$），本页用 Newton 法解出基准解（第 6 节）；图 1.1 中的箭头就是物理绕流基准的速度。以来流速率 0.142、物体直径 1.23 和参考黏度计算，Reynolds 数约为 0.87：惯性力与黏性力同一量级，流动是定常层流。这张 24 个控制体的基网格很粗，壁面附近的梯度分辨不足；第 7 节给出加密网格上的同一流动。

At the default settings (second-order convection, constant transport, isothermal wall, $\mu_{\rm ref}=0.2$) the page computes the baseline solution by Newton's method (Section 6); the arrows of Figure 1.1 are the velocity of the physical-flow baseline. With the inflow speed 0.142, the body diameter 1.23 and the reference viscosity, the Reynolds number is about 0.87: inertial and viscous forces are of the same order, and the flow is steady and laminar. This base mesh of 24 control volumes is coarse and does not resolve the gradients near the wall; Section 7 shows the same flow on a refined mesh.

第 7 节定义目标函数：默认是温度跟踪 $J_{T,h}$，它是体积分，衡量温度与目标 $1+0.02\cos\,x$ 的差；也可以选物体所受合力的 $x$ 或 $y$ 分量（压力加黏性牵引），或壁面热流 $Q_{w,h}$。本页要求的就是这些目标对 $D$ 的导数。

Section 7 defines the objectives: by default the temperature tracking $J_{T,h}$, a volume integral that measures the difference between the temperature and the target $1+0.02\cos\,x$; one may also choose the $x$ or $y$ component of the force on the body (pressure plus viscous traction) or the wall heat flux $Q_{w,h}$. What this page computes is the derivative of these objectives with respect to $D$.

这组方程的作用

What these equations are for

在流体力学中，可压缩 Navier–Stokes 方程描述黏性与导热不能忽略的气体流动：物体表面的边界层与摩擦阻力，壁面与气体之间的换热，流动分离，以及尺度小或速度低、Reynolds 数不高的流动。它们是连续介质流动的基本方程；Euler 方程是略去黏性和导热后的特例。Reynolds 数很高时流动会变成湍流，直接解出所有尺度的涡代价太高，工程计算常改解时间平均的 RANS 方程：它的形式与本页的方程相同，只是黏度加上一个湍流黏度，另由湍流模型的输运方程决定这个黏度。

In fluid dynamics the compressible Navier–Stokes equations describe gas flows in which viscosity and heat conduction cannot be neglected: boundary layers and skin friction on bodies, heat transfer between a wall and the gas, flow separation, and small or slow flows at moderate Reynolds number. They are the basic equations of continuum flow; the Euler equations are the special case without viscosity and conduction. At high Reynolds number the flow becomes turbulent, resolving eddies of every size costs too much, and engineering computations usually solve the time-averaged RANS equations instead: they have the form of this page's equations with an eddy viscosity added to the viscosity, and a transport equation of a turbulence model determines that eddy viscosity.

在本课程中，这一页是从二维 Euler 到二维 RANS–SA 的中间一步。与二维 Euler 页相比，新的是黏性通量：面上除了两侧的状态，还要速度和温度的梯度，所以第 3 节用二次重构得到梯度，前向与伴随也要沿这条梯度链求导；固壁变成无滑移和热壁条件；目标函数多了黏性力和壁面热流。与一维层流 N–S 页相比，物理相同，但梯度成了二维向量，应力成了张量，面有任意方向的法向，网格坐标也成了设计变量。二维 RANS–SA 页在这些之上再加一个湍流模型方程。

In this course the page is the step between 2-D Euler and 2-D RANS–SA. Compared with the 2-D Euler page, what is new is the viscous flux: besides the states on its two sides, a face needs the gradients of velocity and temperature, so Section 3 obtains gradients by quadratic reconstruction, and the tangent and the adjoint must differentiate through this gradient chain too; the wall becomes a no-slip, thermal wall; and the objectives gain the viscous force and the wall heat flux. Compared with the 1-D laminar N–S page, the physics is the same, but the gradient becomes a two-dimensional vector, the stress a tensor, faces have normals in any direction, and the mesh coordinates become design variables. The 2-D RANS–SA page adds one turbulence-model equation on top of all this.

图 1.5 把本页放在网站的路线上：每往右一步主要增加一类新内容，每点的未知量从 1 个增加到 5 个；每一页的残差都由面循环组装，前向与伴随也沿同一个循环求导。点击方框进入该页的第 1 节。

Figure 1.5 places this page on the route through the site: each step to the right mainly adds one kind of ingredient, and the unknowns per point grow from 1 to 5; every page assembles its residual with a face loop, and the tangent and the adjoint differentiate the same loop. Click a box to open Section 1 of that page.

本页要回答的问题

The question this page answers

本页只回答一个问题：设计 $D$ 有一个小变化时，目标 $J_h$ 变化多少？这需要 54 个导数。前向（tangent）方法每个设计分量解一次线性方程，共 54 次；伴随方法只解一次转置方程。页首的全页路线图把后面各节串成这两条路线。

This page answers one question: how much does the objective $J_h$ change when the design $D$ changes a little? That takes 54 derivatives. The tangent (forward) method solves one linear system per design component, 54 in all; the adjoint method solves a single transposed system. The page roadmap at the top strings the later sections along these two routes.

阅读本页需要微积分、线性代数和数值方法（Newton 法与迭代线性求解）。自动微分与最优化不是前提：前向与伴随计算都从链式法则和矩阵转置讲起。

The page assumes calculus, linear algebra and numerical methods (Newton's method and iterative linear solvers). Automatic differentiation and optimisation are not prerequisites: the tangent and adjoint calculations start from the chain rule and matrix transposition.

要在计算机上回答这个问题，先要把连续方程变成有限个未知量的代数方程。第 2 节沿用二维标量页的网格，在每个控制体上存四个守恒量；第 3 节说明怎样在每个面上计算通量：对流部分用两侧的状态，黏性部分还要速度和温度的梯度。

To answer it on a computer, the continuous equations must first become algebraic equations in finitely many unknowns. Section 2 takes the mesh of the 2-D scalar pages and stores four conserved quantities in every control volume; Section 3 shows how the flux is computed on every face: the convective part from the states on the two sides, the viscous part also from the gradients of velocity and temperature.

网格：守恒量存在哪里

The mesh: where the conserved quantities live

第 1 节的方程在区域内每一点成立，计算机却只能存储有限个数。本节选定网格，决定未知量放在哪里（本页放在单元上：格心格式），并说明第 1 节的四个守恒量怎样存成数组、每个控制体怎样与相邻控制体共享面。

The equations of Section 1 hold at every point of the domain, but a computer stores finitely many numbers. This section chooses the mesh, decides where the unknowns live (on this page in the cells: the cell-centred scheme), and shows how the four conserved quantities of Section 1 are stored as arrays and how each control volume shares faces with its neighbours.

从已经认识的二维网格出发

Start with the familiar two-dimensional mesh

本页沿用二维标量的混合 O 型网格：24 个原始节点、8 个四边形和 16 个三角形。一个标量现在换成第 1 节的四个守恒量；第 1 节的两个算例都用这张网格。

The hybrid O-grid is the scalar lesson's mesh: 24 original vertices, eight quadrilaterals and sixteen triangles. Replace one scalar by the four conserved quantities of Section 1; both cases of Section 1 use this mesh.

本页把 $U_i$ 定义为守恒量单元平均。基网格有 24 个控制体、48 个面、96 个状态自由度。数值重构中心使用面积质心，不是顶点的算术平均；二次重构还使用单元中心二阶矩。

Store conservative cell averages $U_i$. The base mesh has 24 volumes, 48 faces and 96 unknowns. Reconstruction uses area centroids rather than vertex averages, and includes central second moments for quadratic reconstruction.

面法向从左控制体指向右控制体，边界指向流体域外。面长 $s_f$ 与单位法向分开存储，求积中只乘一次面长。点击或选择面可查看左右编号、梯度邻域及实际状态；颜色和数值均来自当前求解。数组和编号从 0 开始。

Normals point from left to right volumes, or out of the fluid domain at a boundary. Store length $s_f$ separately from the unit normal and include length exactly once in quadrature. Select a face to inspect indices, reconstruction stencils and actual state values; colours and values come from the current solve. All arrays and indices start at zero.

从“每单位体积”到“一个控制体内有多少”

From density to the amount inside a control volume

第 1 节的四个守恒量都是单位体积内的量。本页采用单位厚度二维模型，因此控制体“体积”数值上等于平面面积。格心存储平均值时，控制体内四种总量就是 $V_iU_i$。格点版存储节点点值，围绕节点积分控制方程；不能把这个点值直接说成对偶体的精确平均值。

The four conserved quantities of Section 1 are amounts per unit volume. The planar model has unit depth, so volume equals area numerically. For stored cell averages, the four amounts in a cell are $V_iU_i$. The nodal scheme stores point values and integrates equations around a vertex; a nodal value is not the exact average over its dual volume.

阅读后文时区分三种符号：$U(x,y)$ 是介绍物理模型时的连续函数；$U_i$ 是四分量存储向量；$\mathbf U_h$ 是长度 $4N$ 的有限维列向量。本页每个符号只有一种含义；下表列出最容易混淆的几个。

Distinguish the continuous field $U(x,y)$, the four-component stored vector $U_i$, and the length-$4N$ column $\mathbf U_h$. Each symbol on this page has one meaning; the table lists those most easily confused.

网格、控制体和存储的数组定好以后，第 3 节在每个面上计算第 1 节的通量：对流部分只要两侧的状态，黏性部分还要速度和温度的梯度，它们都由存储的量重构得到。

With the mesh, the control volumes and the stored arrays in place, Section 3 computes the flux of Section 1 on every face: the convective part needs only the states on the two sides, the viscous part also the gradients of velocity and temperature, all reconstructed from the stored values.

面上的通量：重构与黏性梯度

The flux on a face: reconstruction and viscous gradients

第 2 节把四个守恒量存在每个控制体上，通量却要在控制体之间的面上计算。本节从存储的量出发，得到面上的状态和梯度：先写出面上的黏性通量要用的离散梯度，再用加权最小二乘重构得到两侧的迹和梯度，最后得到一个面的净通量。

Section 2 stores the four conserved quantities in every control volume, but the flux is computed on the faces between control volumes. This section starts from the stored values and obtains the states and gradients at a face: first the discrete gradients the viscous flux on a face uses, then the weighted least-squares reconstruction that gives the traces and gradients on the two sides, and finally the net flux of a face.

面上的黏性通量用离散梯度计算

The viscous flux on a face uses discrete gradients

有了网格和存储的量，下面在每个面的求积点上计算第 1 节的通量 $F_n^c-F_n^v$。对流部分只要两侧的状态；黏性部分中的 $u_x$ 等空间导数，代码用本节后面由重构得到的离散梯度 $G_x^u$ 等替代：

With the mesh and the stored values in place, the flux $F_n^c-F_n^v$ of Section 1 is evaluated at the quadrature points of every face. The convective part needs only the states on the two sides; in the viscous part the code replaces the spatial derivatives $u_x$ and the like by the discrete gradients $G_x^u$ and the like, obtained later in this section by reconstruction:

这里 $G$ 专指下面定义的离散梯度（上标是变量名，下标是方向）。它们与第 1 节的 $\tau$、$q_n$ 是同一组公式，只是导数换成了离散梯度；物性 $\mu$、$\kappa$ 在面状态的温度上取值。

Here $G$ denotes the discrete gradients defined below (superscript: variable; subscript: direction). These are the formulas for $\tau$ and $q_n$ of Section 1 with the derivatives replaced by discrete gradients; the properties $\mu$ and $\kappa$ are evaluated at the temperature of the face state.

同时保留平均值与几何矩

Retain averages and geometric moments

格心中心矩为 $M_{xx}=\langle\Delta x^2\rangle_i$ 等，保证 $\langle P_i\rangle_i=U_i$。用邻居的几何矩计算 $\langle P_i\rangle_j$，而不是把平均值硬当成中心点值。

For cells the moments are $M_{xx}=\langle\Delta x^2\rangle_i$ and similarly for the others, ensuring $\langle P_i\rangle_i=U_i$. Neighbour equations use their geometric moments, not cell averages misinterpreted as centre point values.

邻域规则是连通图上的逐圈扩展：从 $i$ 出发，每次把当前一圈的全部面邻居加进来，直到连同自身至少有 12 个控制体（最多 5 圈）才停止。因为整圈加入，总数可能超过 12。基网格上，格心 0 用到其余 23 个格心中的 15 个；全部格心的邻居数在 11 到 15 之间。因此至少有 11 个邻居差值，多于待定的 5 个系数，加权最小二乘是超定的；四个守恒分量共享同一几何矩阵。

The stencil grows ring by ring on the face-neighbour graph: starting from $i$, add every neighbour of the current ring until, counting $i$ itself, at least 12 control volumes are collected (at most five rings). Whole rings are added, so the count can exceed 12. On the base mesh, cell 0 uses 15 of the other 23; every stencil has between 11 and 15 neighbours. Hence at least 11 neighbour differences constrain the five coefficients, and the weighted least-squares fit is overdetermined; the four conserved components share one geometric matrix.

二次多项式用于黏性迹和梯度。对流一阶使用存储状态，二阶使用 $U_i+c_{x,i}\Delta x+c_{y,i}\Delta y$；无 limiter。

Quadratic polynomials provide viscous traces and gradients. First-order convection uses stored states, while second order uses $U_i+c_{x,i}\Delta x+c_{y,i}\Delta y$. No limiter is applied.

这一对流通量是 Rusanov（局部 Lax–Friedrichs）型：两侧物理通量的平均，减去“耗散速度 × 跳跃”的一半；这里耗散速度取常数 3，所以系数是 $3/2$。它不逐面取 $\max(|u_n|+c)$，因此残差处处可微，复数步长也保持解析；代价是只在 $|u_n|+c<3$ 的范围内保证足够的耗散，第 6 节的诊断 maxWave 必须小于 3。默认基准的最大特征速度为 $1.459$，约为 3 的 49%。同一系列的取值不同：二维 Euler 页固定取 2；一维 N–S 页取两侧平均 $\tfrac12(u^-+u^++c^-+c^+)$，因为一维管道中处处 $u>0$，这个状态相关的平均本身光滑。二维面法向可指向任意方向，$u_n$ 会变号，状态相关的 $|u_n|$ 在 $u_n=0$ 处不可微，所以本页与二维 Euler 页都取常数。上面的基准只用到 3 的约一半，常数 2 也能覆盖它，但余量更小；本页取 3。

This convective flux is of Rusanov (local Lax–Friedrichs) type: the average of the two physical fluxes minus half of “dissipation speed × jump”. The dissipation speed is the constant 3, hence the coefficient $3/2$. It does not take $\max(|u_n|+c)$ face by face, so the residual stays differentiable everywhere and complex steps stay analytic; in exchange, enough dissipation is guaranteed only while $|u_n|+c<3$, and the Section 6 diagnostic maxWave must stay below 3. The default baseline's largest characteristic speed is $1.459$, about 49% of 3. The series uses different values: the 2D Euler page fixes 2; the 1D N–S page uses the two-sided average $\tfrac12(u^-+u^++c^-+c^+)$, which is smooth because $u>0$ throughout the 1D duct. In 2D a face normal can point anywhere, $u_n$ changes sign and a state-dependent $|u_n|$ is not differentiable at $u_n=0$; this page and the 2D Euler page therefore use constants. The baseline above uses about half of 3; the constant 2 would also cover it, with less margin, and this page takes 3.

最小二乘正规方程怎样得到

Deriving the least-squares normal equations

先只看一个守恒分量 $z$，例如密度。令五个未知系数排成列 $c_i$，邻居中心相对位置为 $(\Delta x,\Delta y)$。因为 $z_j$ 是平均值，格心的拟合行还必须包含邻居和本单元的中心矩：

Consider one conserved component $z$, such as density. Stack its five coefficients into $c_i$ and let $(\Delta x,\Delta y)$ be the neighbour-centre offset. Cell-average fitting must include both cells’ central moments:

格点版本令全部中心矩为零。把所有 $b_{ij}$ 叠成矩阵 $B_i$，把 $z_j-z_i$ 叠成 $d_i$，权重构成对角矩阵 $W_i$。最小化的是“邻居数据与多项式预测不一致”的加权平方和：

For nodes all central moments are zero. Stack rows $b_{ij}$ into $B_i$, differences $z_j-z_i$ into $d_i$, and weights into diagonal $W_i$. Minimize the weighted squared mismatch between data and polynomial predictions:

这是一个 $5\times5$ 系统，每个守恒分量各解一次，几何矩阵可以共享。邻居数量多于五仍不保证可解，还要有足够独立的几何方向。实现检查权重的有限性和二次多项式再现误差；这里写逆矩阵是为了说明关系，不要求读者把“大矩阵显式求逆”当作求解原则。拟合矩阵 $H_i$ 与面通量 $H_f$ 由下标和维度区分。

This is a $5\times5$ system for each conserved component, sharing geometry. More than five neighbours is insufficient without independent geometric directions. The implementation checks finite weights and quadratic reproduction. Inverse notation expresses the relation, not a recommendation to explicitly invert large matrices. Fit matrix $H_i$ and face-flux vector $H_f$ differ in indices and dimensions.

黏性迹、梯度与跳跃修正

Viscous traces, gradients and jump correction

$\theta\in\{u,v,T\}$。先对守恒多项式求空间梯度，再用原始变量转换的 Jacobian 得到各侧原始量梯度。跳跃修正使用同一面点的两侧迹，系数为 2；它不是一维的中心差分。精确二次守恒多项式的左右迹一致，修正为零。非线性转换、边界和通量装配后的整体阶数仍需制造解测量。

For $\theta\in\{u,v,T\}$, differentiate the conservative polynomial spatially, then apply the primitive-conversion Jacobian. The penalty uses both traces at the same face point with coefficient two; it is not the 1D centre difference. Exact quadratic conservative traces agree and have zero penalty. Nonlinear conversion, boundaries and flux assembly still require measured global convergence.

为什么需要跳跃项？只平均两侧梯度时，中心差分形式的单元梯度看不到一维奇偶（棋盘）模式 $\theta_j=(-1)^j$：每个单元的梯度都是零，黏性通量也为零，这个模式既不受阻尼，也不被察觉。把跳跃项系数记作 $\alpha_{\rm J}$（本页 $\alpha_{\rm J}=2$），在均匀一维网格、线性重构且面位于中点时，面梯度化为 $G_f=\alpha_{\rm J}(\theta_R-\theta_L)/h+(1-\alpha_{\rm J})(g_L+g_R)/2$，其中 $g_L,g_R$ 是两侧单元梯度。对 $\theta_j=e^{\mathrm ij\phi}$，扩散算子的符号乘以 $h^2$ 为

Why add a jump term? Averaging the two cell gradients alone cannot see the 1D odd-even (checkerboard) mode $\theta_j=(-1)^j$: every centred cell gradient vanishes, so the viscous flux vanishes and the mode is neither damped nor detected. Write the jump coefficient as $\alpha_{\rm J}$ ($\alpha_{\rm J}=2$ here). On a uniform 1D grid with linear reconstruction and the face at the midpoint, the face gradient becomes $G_f=\alpha_{\rm J}(\theta_R-\theta_L)/h+(1-\alpha_{\rm J})(g_L+g_R)/2$, with $g_L,g_R$ the cell gradients. For $\theta_j=e^{\mathrm ij\phi}$, $h^2$ times the diffusion symbol is

最短波 $\phi=\pi$ 的值依次为：$\alpha_{\rm J}=0$ 时为 $0$（奇偶解耦）；$\alpha_{\rm J}=1$ 时为 $-4$（紧致两点差分）；$\alpha_{\rm J}=4/3$ 时为 $-\tfrac{16}{3}$，并使 $\phi^2$ 修正项的系数 $\tfrac{3\alpha_{\rm J}-4}{12}$ 等于 $0$，在均匀网格上达到四阶；$\alpha_{\rm J}=2$ 时为 $-8$，最短波阻尼是紧致差分的两倍，光滑波的修正系数为 $\tfrac{1}{6}$，仍是二阶相容。可见任何 $\alpha_{\rm J}>0$ 都相容，差别只在有限网格的阻尼与截断误差；本页取 2，是偏向强阻尼的选择，不是由方程唯一决定的常数。

At the shortest wave $\phi=\pi$ the values are: $\alpha_{\rm J}=0$ gives $0$ (odd-even decoupling); $\alpha_{\rm J}=1$ gives $-4$ (compact two-point difference); $\alpha_{\rm J}=4/3$ gives $-\tfrac{16}{3}$ and makes the $\phi^2$ coefficient $\tfrac{3\alpha_{\rm J}-4}{12}$ equal $0$, fourth order on a uniform grid; $\alpha_{\rm J}=2$ gives $-8$, twice the damping of the compact difference, while the smooth-wave coefficient $\tfrac{1}{6}$ keeps second-order consistency. Any $\alpha_{\rm J}>0$ is consistent; the choice only changes finite-mesh damping and truncation error. This page uses 2, favouring strong damping; it is not a constant fixed by the equations.

守恒量梯度怎样转换成速度、温度梯度

Converting conserved gradients to velocity and temperature gradients

对 $u=m_x/\rho$ 使用商法则，再对总能量关系使用链式法则。若 $\partial_kP^z$ 表示重构守恒分量 $z$ 的多项式在方向 $k=x,y$ 上的斜率，实际计算为：

Apply the quotient rule to $u=m_x/\rho$ and the chain rule to total energy. With $\partial_kP^z$ the slope, in direction $k=x,y$, of the reconstructed polynomial of conserved component $z$, code evaluates:

不要用“分别拟合速度和温度”替代这条链：它会定义另一种离散残差。本页先拟合守恒量，再进行上述非线性转换。二阶对流仅取多项式的常数与线性部分，黏性仍使用完整二次迹及其梯度；空间精度按钮不是将所有多项式统一切成同一次数。

Fitting velocity and temperature independently would define a different residual. This implementation reconstructs conserved variables and then converts them nonlinearly. Second-order convection keeps only the constant and linear terms, while viscosity retains quadratic traces and gradients; the spatial-order control does not cut every polynomial to one degree.

逐步推导：原始变量、Jacobian 行与梯度转换

Step by step: primitive variables, Jacobian rows and the gradient conversion

面上的每个物理量都从守恒量 $U=(\rho,m_x,m_y,e)$ 的迹算出。前向、反向都要反复用到它们对 $U$ 的导数，所以先把这组导数一次推导完整，后面各节直接引用。下面的行向量都是 $1\times4$，列的顺序与代码的四个分量槽一致。

Every physical quantity on a face is computed from traces of the conserved state $U=(\rho,m_x,m_y,e)$. The tangent and the reverse pass use their derivatives with respect to $U$ again and again, so derive them once, completely; later sections refer back here. All rows below are $1\times4$, ordered like the four component slots in the code.

第一步，写出原始变量（physics.js 的 prim()）。$E=e/\rho$ 是单位质量总能：

Step 1, write the primitive variables (prim() in physics.js). $E=e/\rho$ is total energy per unit mass:

第二步，逐个求方向导数（点号）。商法则给出 $\dot u=(\dot m_x-u\,\dot\rho)/\rho$、$\dot E=(\dot e-E\,\dot\rho)/\rho$；再把它们代入 $\dot T=(\dot E-u\,\dot u-v\,\dot v)/c_v$ 并按 $\dot\rho,\dot m_x,\dot m_y,\dot e$ 归并系数：