方程：二维定常可压缩 Euler 方程组

The equations: steady two-dimensional compressible Euler

从守恒律到 Euler 方程组

From conservation laws to the Euler equations

在流体中任取一块固定的区域。区域里的质量、两个方向的动量和总能量，只能通过它的边界改变：流体穿过边界时把这些量带进带出，边界上的压力推动流体并对它做功。不计黏性和热传导时，这四条守恒律就是二维可压缩 Euler 方程组。

Take any fixed region of the fluid. The mass, the momentum in two directions and the total energy inside it can change only through its boundary: the fluid carries these quantities in and out as it crosses the boundary, and the pressure on the boundary pushes the fluid and does work on it. Without viscosity and heat conduction, these four conservation laws are the two-dimensional compressible Euler equations.

拟一维 Euler 的三分量状态现在多了一个横向动量。二维没有变截面积源项：几何直接进入每个面的面积向量（第 2 节）。采用理想气体 $\gamma=1.4$，$\rho E$ 是单位体积总能量，包含两个方向的动能；压力计算漏掉横向动能，就会破坏原始方程及所有导数。

The quasi-one-dimensional Euler state gains a transverse momentum component. There is no variable-area source: geometry enters through each face area vector (Section 2). Use an ideal gas with $\gamma=1.4$. The total energy per unit volume $\rho E$ includes kinetic energy in both directions; omitting the transverse kinetic energy corrupts both the primal equations and their derivatives.

第一行是质量守恒；第二、三行是两个方向的动量守恒，压力 $p$ 像通量一样出现在其中；第四行是能量守恒，能量通量 $(\rho E+p)\boldsymbol u$ 包括流体携带的总能量和压力所做的功。四个方程里有 $\rho,u,v,E,p$ 五个量，状态方程 $p=(\gamma-1)\bigl(\rho E-\tfrac12\rho(u^2+v^2)\bigr)$ 把它们闭合；声速 $c=\sqrt{\gamma p/\rho}$ 决定下面的特征速度。

The first row is conservation of mass; the second and third are conservation of momentum in the two directions, where the pressure $p$ appears like a flux; the fourth is conservation of energy, whose flux $(\rho E+p)\boldsymbol u$ holds the total energy carried by the fluid and the work done by the pressure. The four equations contain five quantities $\rho,u,v,E,p$; the equation of state $p=(\gamma-1)\bigl(\rho E-\tfrac12\rho(u^2+v^2)\bigr)$ closes them, and the speed of sound $c=\sqrt{\gamma p/\rho}$ sets the characteristic speeds below.

对任意控制体 $V$ 积分，并用散度定理把通量散度的体积分变成边界积分，就得到积分形式。记 $\hat{\boldsymbol n}=(\hat n_x,\hat n_y)$ 为 $V$ 边界的单位外法向，$u_n=\boldsymbol u\cdot\hat{\boldsymbol n}$ 为法向速度，$F_{\hat n}$ 是穿过边界的法向通量：

Integrating over any control volume $V$, and turning the volume integral of the flux divergence into a boundary integral with the divergence theorem, gives the integral form. With $\hat{\boldsymbol n}=(\hat n_x,\hat n_y)$ the outward unit normal of the boundary of $V$ and $u_n=\boldsymbol u\cdot\hat{\boldsymbol n}$ the normal velocity, $F_{\hat n}$ is the normal flux through the boundary:

本页只研究定常解。令 $\partial_tU=0$，要离散的就是定常的积分形式：对计算域 $\Omega$ 中的每个控制体，穿过它边界的净通量为零。第 2 节把控制体取成网格的单元（格心）或中位对偶体（格点），第 3 节把边界积分写成各个面上的通量之和。

This page studies steady solutions only. With $\partial_tU=0$, the equation to be discretised is the steady integral form: for every control volume in the domain $\Omega$, the net flux through its boundary vanishes. Section 2 takes the control volumes to be the cells (cell-centred) or the median-dual volumes (node-centred) of a mesh, and Section 3 writes the boundary integral as a sum of fluxes over faces.

本页的来流马赫数低，绕流光滑、没有激波：第 3 节的一阶通量不需要限制器，第 14 节的复数步长也适用于整条计算链。激波及其伴随不在本页讨论。

The free-stream Mach number of this page is low and the flow past the body is smooth, without shocks: the first-order flux of Section 3 needs no limiter, and the complex step of Section 14 applies to the whole chain of computations. Shocks and their adjoints are not treated on this page.

图 1.1 把这个问题画在一张图上（基准状态 $M_\infty=0.2$、$\theta=10^\circ$，本页代码的收敛解，单元平均值；第 2 节起讲怎样算）。左边是计算域：物面八边形（琥珀色）与远场八边形（蓝色）之间的区域；每个控制体按当地马赫数 $|\boldsymbol u|/c$ 着色，短箭头是当地速度，远场的实线面是入流面，虚线面是出流面。来流从左下方流入：流体在物体上游一侧减速，在物体上下两侧加速；当地马赫数在 0.176 到 0.236 之间，最高值在物体下侧紧贴物面的控制体里，处处远小于 1。右边是物面压力：迎风一侧高于来流压力 $p_\infty=1$，背风一侧低于它，流体在物体后面没有恢复压力。第 7 节说明这主要是本页固定数值耗散造成的离散误差，这里的阻力不为零也是因为它。

Figure 1.1 draws the problem in one picture (base state $M_\infty=0.2$, $\theta=10^\circ$; the converged solution of this page’s code, cell averages; Section 2 onwards shows how it is computed). On the left is the domain, the region between the body octagon (amber) and the far-field octagon (blue); every control volume is coloured by its local Mach number $|\boldsymbol u|/c$, the short arrows are the local velocity, and the solid far-field faces are inflow faces, the dashed ones outflow faces. The free stream comes from the lower left: the fluid slows down on the upstream side of the body and speeds up above and below it; the local Mach number lies between 0.176 and 0.236, highest in a volume next to the wall on the lower side of the body, and far below 1 everywhere. On the right is the wall pressure: above the free-stream pressure $p_\infty=1$ on the upstream side and below it on the downstream side, where the flow does not recover its pressure behind the body. Section 7 shows that this is mostly discretisation error caused by the page’s fixed numerical dissipation, which is also why the drag here is not zero.

特征速度：一个边界要给几个条件

Characteristic speeds: how many conditions a boundary takes

把法向通量 $F_{\hat n}$ 对 $U$ 求导，得到 $4\times4$ 的 Jacobian $\partial F_{\hat n}/\partial U$。它有四个实特征值：

Differentiating the normal flux $F_{\hat n}$ with respect to $U$ gives the $4\times4$ Jacobian $\partial F_{\hat n}/\partial U$. It has four real eigenvalues:

它们是小扰动沿 $\hat{\boldsymbol n}$ 传播的速度：两个声波以 $u_n\pm c$ 传播，另外两个量（熵与切向速度）以流体的法向速度 $u_n$ 随流输运。第 4 节由 Riemann 不变量得到同样的四个速度；第 8 节在代码的通量 Jacobian 上核对它们的和与积。

They are the speeds at which small disturbances travel along $\hat{\boldsymbol n}$: two sound waves move at $u_n\pm c$, and two more quantities (the entropy and the tangential velocity) are carried with the normal velocity $u_n$ of the fluid. Section 4 finds the same four speeds from the Riemann invariants; Section 8 checks their sum and product on the flux Jacobian of the code.

在边界上取 $\hat{\boldsymbol n}$ 指向域外。速度为负的波从外部进入计算域，它携带的信息只能由边界提供；速度为正的波离开计算域，它的值由域内决定，若在边界上再指定它，问题就过定了。所以边界条件的个数等于负特征值的个数：

At a boundary take $\hat{\boldsymbol n}$ pointing out of the domain. A wave with a negative speed enters the domain from outside, and the information it carries can only come from the boundary; a wave with a positive speed leaves the domain, its value is decided inside, and prescribing it at the boundary as well would over-determine the problem. The number of boundary conditions is therefore the number of negative eigenvalues:

• 超声速入流（$u_n<-c$）：四个特征值都为负，给 4 个条件，即整个状态；
• 亚声速入流（$-c<u_n<0$）：$\lambda_1,\lambda_2,\lambda_3<0$，给 3 个条件；
• 亚声速出流（$0<u_n<c$）：只有 $\lambda_1<0$，给 1 个条件；
• 超声速出流（$u_n>c$）：没有负特征值，不给条件；
• 滑移壁面（$u_n=0$）：只有 $\lambda_1=-c<0$，给 1 个条件，就是无穿透条件 $u_n=0$ 本身；$\lambda_2=\lambda_3=0$ 的两个量沿壁面输运。

• supersonic inflow ($u_n<-c$): all four eigenvalues are negative; 4 conditions, the whole state;
• subsonic inflow ($-c<u_n<0$): $\lambda_1,\lambda_2,\lambda_3<0$; 3 conditions;
• subsonic outflow ($0<u_n<c$): only $\lambda_1<0$; 1 condition;
• supersonic outflow ($u_n>c$): no negative eigenvalue; no condition;
• slip wall ($u_n=0$): only $\lambda_1=-c<0$; 1 condition, the no-penetration condition $u_n=0$ itself; the two quantities with $\lambda_2=\lambda_3=0$ travel along the wall.

图 1.2 在本页的收敛解上画出这些速度。上半部分是三个边界面：远场入流面 $f=44$ 的四个速度中有 3 个为负，出流面 $f=40$ 只有 $u_n-c$ 为负，物面 $f=32$ 上 $u_n=0$，只有 $-c$ 为负。下半部分把条件个数画成法向马赫数 $u_n/c$ 的函数：本页全部 8 个远场面的 $|u_n|/c$ 都不超过 0.18，远在亚声速带内。所以远场入流面从来流取 3 个条件，出流面取 1 个，物面只要 $u_n=0$。第 4 节用 Riemann 不变量 $r_\pm=u_n\pm2c/(\gamma-1)$、熵函数 $S=p/\rho^\gamma$ 与切向速度 $u_t$ 实现远场的这种取舍，并把 $u_n=0$ 写成只含压力的壁面通量。

Figure 1.2 draws these speeds on this page’s converged solution. The upper part shows three boundary faces: on the far-field inflow face $f=44$, 3 of the 4 speeds are negative; on the outflow face $f=40$, only $u_n-c$ is; and on the wall face $f=32$, $u_n=0$ and only $-c$ is negative. The lower part draws the number of conditions as a function of the normal Mach number $u_n/c$: all 8 far-field faces of this page have $|u_n|/c$ at most 0.18, far inside the subsonic band. A far-field inflow face therefore takes 3 conditions from the free stream, an outflow face 1, and the wall only $u_n=0$. Section 4 implements this choice at the far field with the Riemann invariants $r_\pm=u_n\pm2c/(\gamma-1)$, the entropy function $S=p/\rho^\gamma$ and the tangential velocity $u_t$, and writes $u_n=0$ as a wall flux that contains only the pressure.

本页的问题：绕物体的亚声速流动

The problem on this page: subsonic flow past a body

计算域是一个闭合物体与外边界之间的区域（图 1.1 左）。物体是一条光滑的闭曲线，外边界是半径 1.9 的圆；计算实际使用两条曲线上各 8 个网格节点连成的八边形，第 2 节给出物面曲线并构造网格。物面是滑移壁面；外边界是远场，在那里流动应当回到均匀的来流。

The domain is the region between a closed body and an outer boundary (left of Figure 1.1). The body is a smooth closed curve and the outer boundary a circle of radius 1.9; the computation uses the octagons through the 8 mesh nodes on each of the two curves, and Section 2 gives the body curve and builds the mesh. The body is a slip wall; the outer boundary is the far field, where the flow should return to the uniform free stream.

来流固定为 $\rho_\infty=1$、$p_\infty=1$，速度大小为 $M_\infty\sqrt\gamma$，方向为攻角 $\theta$。于是 $c_\infty=\sqrt\gamma$，$|\boldsymbol u_\infty|/c_\infty=M_\infty$ 正好是来流马赫数。每次原始求解都从所有控制体都等于 $U_\infty$ 的均匀流出发。第 7 节的目标值，以及第 10–14 节对 $M_\infty$、$\theta$ 的导数，都以这个无量纲化为准。

The free stream is fixed at $\rho_\infty=1$ and $p_\infty=1$, with speed $M_\infty\sqrt\gamma$ in the direction of the incidence $\theta$. Hence $c_\infty=\sqrt\gamma$ and $|\boldsymbol u_\infty|/c_\infty=M_\infty$ is exactly the free-stream Mach number. Every primal solve starts from the uniform state $U_\infty$ in all volumes. The objective values of Section 7 and the derivatives with respect to $M_\infty$ and $\theta$ in Sections 10–14 refer to this scaling.

设计变量是 $D=(M_\infty,\theta,X)$：两个来流参数，加上全部 24 个网格节点的坐标 $X$（48 个数），共 50 个分量。按伴随优化文献的通常写法，设计变量记作 $D$。移动物面节点改变物体的形状，移动远场节点改变外边界，移动其余节点只改变网格；它们都会改变离散解和目标。第 10 节给出这些分量的排列，第 12 节求对全部坐标的导数。

The design variables are $D=(M_\infty,\theta,X)$: the two free-stream parameters and the coordinates $X$ of all 24 mesh nodes (48 numbers), 50 components in all. As is usual in the adjoint-optimization literature, the design variables are written $D$. Moving a body node changes the shape of the body, moving a far-field node changes the outer boundary, and moving any other node changes only the mesh; each of them changes the discrete solution and the objective. Section 10 lists the components, and Section 12 differentiates with respect to all coordinates.

第 7 节给出设计的目的：目标 $J$ 是流体作用在物体上的压力合力的一个分量，默认是 $y$ 分量 $J_y$，页首可以改成 $x$ 分量 $J_x$。它由物面压力（图 1.1 右）逐面求和得到；第 7 节还讨论在这些粗网格上，这个力有多少是物理的。

Section 7 states the purpose of the design: the objective $J$ is one component of the pressure force that the fluid exerts on the body, by default the $y$ component $J_y$; the header can switch it to the $x$ component $J_x$. It is summed face by face from the wall pressures (right of Figure 1.1); Section 7 also discusses how much of this force is physics on these coarse meshes.

这组方程的作用

What these equations are for

在流体力学中，Euler 方程组描述不计黏性和热传导的可压缩流动。黏性只在物面附近的薄边界层和尾迹里起作用时，例如高雷诺数下的流线型物体，它给出其余区域的速度和压力；它包含声波，流速更高时还包含激波（本页不讨论）；它也是 Navier–Stokes 方程组去掉黏性应力和热传导后剩下的部分。无黏流中物体受的力全部来自物面压力；另一方面，稳态、亚声速、无激波的无黏绕流阻力为零（d'Alembert 佯谬，第 7 节），这提醒我们它没有描述什么。气动外形设计常先用 Euler 方程组计算压力载荷，再用伴随求载荷对大量外形参数的梯度，这正是本页的内容。

In fluid dynamics the Euler equations describe compressible flow without viscosity and heat conduction. Where viscosity matters only in thin boundary layers next to walls and in wakes, as for streamlined bodies at high Reynolds numbers, they give the velocity and the pressure everywhere else; they contain sound waves and, at higher speeds, shock waves (not treated here); and they are what remains of the Navier–Stokes equations without viscous stresses and heat conduction. In inviscid flow the force on a body comes entirely from the wall pressure; on the other hand, steady, subsonic, shock-free inviscid flow past a body has zero drag (d’Alembert’s paradox, Section 7), a reminder of what the equations leave out. Aerodynamic shape design often computes pressure loads with the Euler equations first and then uses an adjoint for the gradient of a load with respect to many shape parameters, which is the subject of this page.

在本课程中，这组方程把前面两类页面接在一起：

In this course the equations join the two earlier kinds of page:

• 同一张网格：二维标量页的混合网格与面积向量（第 2 节）、面循环（第 5 节）原样沿用，每个控制体的未知量从 1 个变成 4 个；
• 同样的热力学与特征边界：状态方程、声速与 Riemann 不变量来自拟一维 Euler，二维的新内容是在远场先把速度分解到面的法向与切向（第 4 节）；
• 一个方程组：Jacobian 的每个非零元变成一个 $4\times4$ 块（第 6 节），前向与伴随按块传播；
• 几何也是设计变量：节点坐标进入面积向量、面长与单位法向，第 12 节对它们求导；
• 一个光滑的小例子：96 个未知量，没有激波，每个导数都能用有限差分与复数步长核对（第 13、14 节）。

• The same mesh. The hybrid mesh and the area vectors (Section 2) and the face loop (Section 5) of the 2-D scalar pages carry over unchanged; each control volume now holds 4 unknowns instead of 1.
• The same thermodynamics and characteristic boundaries. The equation of state, the sound speed and the Riemann invariants come from quasi-1-D Euler; what is new in two dimensions is that the far field first splits the velocity into the normal and tangential components of a face (Section 4).
• A system. Every nonzero entry of the Jacobian becomes a $4\times4$ block (Section 6), and the tangent and the adjoint propagate blocks.
• Geometry as design variables. The node coordinates enter the area vectors, face lengths and unit normals, and Section 12 differentiates through them.
• A small, smooth example. 96 unknowns and no shocks, so every derivative can be checked with finite differences and complex steps (Sections 13 and 14).

图 1.3 把本页放在整个网站的路线上，每个方框都链接到那一页的第 1 节。每往右一步主要增加一类新内容：一维 Burgers 只有一个非线性对流方程；拟一维 Euler 换成质量、动量、能量三个方程，并引入压力和变截面积；本页增加第二个空间方向，面有了法向，网格坐标也成为设计变量；二维层流 N–S 再加入黏性应力和热传导，面上要用到梯度；二维 RANS–SA 再加一个湍流模型方程。每点的未知量从 1 个增加到 5 个，但每一页的残差仍由面循环组装，前向与伴随也沿同一个循环求导。另有两页与本页相邻：二维标量对流页用的是同一张网格（见它的第 2 节），一维层流 N–S 页在一维的 Euler 方程上加入黏性和热传导。

Figure 1.3 places this page on the route through the site; every box links to Section 1 of its page. Each step to the right mainly adds one kind of ingredient: 1-D Burgers has a single nonlinear convection equation; quasi-1-D Euler replaces it by three equations (mass, momentum and energy) and brings in the pressure and a varying cross-section; this page adds a second space direction, so faces have normals and the mesh coordinates become design variables; 2-D laminar N–S adds viscous stresses and heat conduction, which need gradients at the faces; 2-D RANS–SA adds one turbulence-model equation. The unknowns per point grow from 1 to 5, yet every page still assembles its residual with face loops, and the tangent and the adjoint differentiate the same loops. Two more pages sit next to this one: the 2-D scalar advection page uses the same mesh (see its Section 2), and the 1-D laminar N–S page adds viscosity and heat conduction to the Euler equations in one dimension.

本页要回答的问题

The question this page answers

本页假设读者熟悉微积分、矩阵运算、Newton 法和 Krylov 迭代（GMRES），最好读过二维标量对流页（同一张网格与面循环）和拟一维 Euler 页（同样的热力学与特征边界）；本页用到的内容都会重新给出。

This page assumes calculus, matrix algebra, Newton’s method and a Krylov solver (GMRES), and it helps to have read the 2-D scalar advection page (the same mesh and face loop) and the quasi-1-D Euler page (the same thermodynamics and characteristic boundaries); the parts this page uses are restated here.

本页只回答一个问题：设计变量 $D$ 改变一点，压力合力 $J$ 改变多少？$D$ 有 50 个分量。前向（tangent）路线每个分量解一次线性方程，共 50 次；伴随路线只解一次转置方程，就得到全部分量（第 11 节）。页首的全页路线图把后面各节串成一条线。按页首默认设置实际计数，50 次前向求解共用 2400 次 GMRES 迭代，1 次伴随求解只用 52 次。

This page answers one question: how much does the pressure force $J$ change when the design $D$ changes a little? $D$ has 50 components. The tangent (forward) route solves one linear system per component, 50 in all; the adjoint route solves a single transposed system and obtains every component (Section 11). The page roadmap at the top strings the later sections together. Counted with the default settings at the top of the page, the 50 tangent solves take 2400 GMRES iterations in total and the single adjoint solve only 52.

要在计算机上回答这个问题，先要把连续的方程组变成有限个未知量的代数方程。第 2 节沿用二维标量的网格，决定四个守恒量存放在哪里；第 3 节再把定常的积分形式写成各个面上的通量之和。

To answer it on a computer, the continuous equations must first become algebraic equations in finitely many unknowns. Section 2 reuses the mesh of the 2-D scalar pages and decides where the four conserved variables live; Section 3 then writes the steady integral form as a sum of fluxes over faces.

网格：沿用二维标量的混合 O 型网格

Mesh: the scalar example’s hybrid O-grid

第 1 节的方程组对计算域中的每一点都成立，计算机却只能存放有限个数。本节沿用二维标量页的网格，决定四个守恒量存放在哪里，以及每个面的面积向量怎样计算。

The equations of Section 1 hold at every point of the domain, but a computer stores finitely many numbers. This section reuses the mesh of the 2-D scalar pages and decides where the four conserved variables live and how the area vector of every face is computed.

基础网格的 24 个节点坐标、8 个四边形和 16 个三角形，逐点、逐单元与二维标量一致。节点分布在三圈上，每圈取 8 个角度 $t_j=2\pi j/8$：内圈取在下面的物面曲线上，外圈取在半径 1.9 的圆上，中圈按 0.45 的比例把两者线性混合。数组从 0 开始编号；边界的 $R=-1$ 表示没有右控制体，绝不能拿它访问数组。

The 24 base nodes, 8 quadrilaterals and 16 triangles match the scalar example point for point and cell for cell. The nodes lie on three rings at the 8 angles $t_j=2\pi j/8$: the inner ring on the body curve below, the outer ring on the circle of radius 1.9, and the middle ring blends the two with weight 0.45. Arrays are zero-based. $R=-1$ marks an absent exterior control volume and must never be used as an array index.

图 2.1 把上式与三圈节点画在一起：虚线是曲线本身，实线是计算真正使用的多边形。

Figure 2.1 draws this curve together with the three rings of nodes: the dashed lines are the curves themselves, the solid lines the polygons the computation actually uses.

相邻节点用直线相连，所以计算中的物面是内接于这条曲线的八边形，远场是半径 1.9 的正八边形；本页所有导数都针对这两条多边形边界。第 15 节的加密网格保持这两个八边形不变。

Neighbouring nodes are joined by straight edges, so the computational body is an octagon inscribed in this curve and the far field is a regular octagon of radius 1.9. Every derivative on this page refers to these polygonal boundaries. The refined meshes of Section 15 keep both octagons unchanged.

每个原始单元保存四个守恒分量，单元平均值显示在顶点平均中心。每条边就是一个计算面；逆时针绕单元走时，面积向量指向流出该单元的方向。这里共有 24 个控制体、48 个面、96 个状态自由度。

Each original cell stores four conserved averages, displayed at its vertex-average centre. An edge is a computational face. Walking counterclockwise around a cell gives an outward area vector. There are 24 volumes, 48 faces and 96 state degrees of freedom.

图 2.2 左边给出这 24 个控制体的编号，右边给出它们在状态向量里的位置：控制体 $i$ 的四个分量占第 $4i$ 到 $4i+3$ 个位置。残差 $R$ 和伴随 $\psi$ 按同样的顺序存放，所以三个向量的第 $4i+k$ 个分量指同一个控制体、同一条守恒方程（$k=0,1,2,3$ 依次是质量、$x$ 动量、$y$ 动量、能量）。

The left of Figure 2.2 numbers the 24 control volumes and the right shows where they sit in the state vector: the four components of volume $i$ occupy positions $4i$ to $4i+3$. The residual $R$ and the adjoint $\psi$ are stored in the same order, so component $4i+k$ of all three vectors refers to the same volume and the same conservation equation ($k=0,1,2,3$: mass, $x$ momentum, $y$ momentum, energy).

这里 $(x_s,y_s)$ 与 $(x_e,y_e)$ 是面的起点和终点（格点页的对偶面按各段相加）。面积向量 $\boldsymbol n_f$ 已经包含面长；单位法向 $\hat{\boldsymbol n}_f$ 只用于第 4 节的特征分解，两者不能混用。

Here $(x_s,y_s)$ and $(x_e,y_e)$ are the start and end points of the face (a dual face on the node page sums its segments). The area vector $\boldsymbol n_f$ already includes the face length; the unit normal $\hat{\boldsymbol n}_f$ is used only for the characteristic split of Section 4, and the two must not be mixed.

图 2.3 以一个控制体为例画出上面的约定。

Figure 2.3 illustrates this convention for one control volume.

逐步推导：每个控制体的面积向量之和为零

Step by step: the area vectors of every control volume sum to zero

第一步，单元是一个多边形。沿逆时针方向依次取边界折线的顶点 $\boldsymbol x_1,\boldsymbol x_2,\dots$，绕一圈后回到 $\boldsymbol x_1$；第 $k$ 段从 $\boldsymbol x_k$ 到 $\boldsymbol x_{k+1}$，它的外法向面积向量是把边向量顺时针转 90°：

Step 1, a cell is a polygon. Take the vertices of its boundary counterclockwise, $\boldsymbol x_1,\boldsymbol x_2,\dots$, returning to $\boldsymbol x_1$ after one loop; segment $k$ runs from $\boldsymbol x_k$ to $\boldsymbol x_{k+1}$, and its outward area vector is the edge vector turned clockwise by 90°:

第二步，$\perp$ 是线性运算；绕一圈的边向量首尾相消，和为零：

Step 2, $\perp$ is linear, and the edge vectors of one closed loop telescope to zero:

第三步，代码按面存储 $\boldsymbol n_f$（从 $L$ 指向 $R$），控制体 $i$ 的外法向是 $\boldsymbol n_{if}=\mathcal B_{if}\boldsymbol n_f$，符号与第 5 节的 $\mathcal B_{if}$ 相同。于是 $\sum_f\mathcal B_{if}\boldsymbol n_f=0$。物理通量对面积向量是线性的，所以均匀状态 $U$ 在每个不含边界面的控制体上给出 $\sum_f\mathcal B_{if}F_{\boldsymbol n_f}(U)=F_{\sum_f\mathcal B_{if}\boldsymbol n_f}(U)=0$，耗散项也因 $U_R-U_L=0$ 而为零：均匀状态在所有内部控制体上的残差都是零。第 6 节的对角块和第 7 节的压力合力都用到同一个事实。

Step 3, the code stores $\boldsymbol n_f$ per face (pointing from $L$ to $R$); the outward vector of volume $i$ is $\boldsymbol n_{if}=\mathcal B_{if}\boldsymbol n_f$, with the sign $\mathcal B_{if}$ of Section 5. Hence $\sum_f\mathcal B_{if}\boldsymbol n_f=0$. The physical flux is linear in the area vector, so in every volume without boundary faces a uniform state $U$ gives $\sum_f\mathcal B_{if}F_{\boldsymbol n_f}(U)=F_{\sum_f\mathcal B_{if}\boldsymbol n_f}(U)=0$, and the dissipation vanishes because $U_R-U_L=0$: a uniform state has zero residual in every interior volume. The diagonal blocks of Section 6 and the pressure force of Section 7 use the same fact.

核对：全部 24 个控制体上 $|\sum_f\mathcal B_{if}\boldsymbol n_f|$ 的最大值为 2×10⁻¹⁶，而最长的面长为 1.454（本页代码在构建时计算）。

Check: the largest $|\sum_f\mathcal B_{if}\boldsymbol n_f|$ over all 24 volumes is 2×10⁻¹⁶, while the longest face has length 1.454 (computed by this page’s code at build time).

点击下面面板中的控制体或面，可查看编号、法向和守恒量；上方可切换着色的物理量，色条给出数值范围，图例说明物面（琥珀色）、远场（蓝色）与内部面的画法。

Click a volume or a face in the panel below to inspect its index, normal and state. The selector changes the plotted quantity, the colour bar gives its range, and the legend explains how body faces (amber), far-field faces (blue) and interior faces are drawn.

控制体、面和面积向量都已确定，状态也按控制体排好了顺序。第 3 节把第 1 节的法向通量写到每个面上，得到四分量的数值通量。

The control volumes, the faces and their area vectors are now fixed, and the state is stored volume by volume. Section 3 writes the normal flux of Section 1 on every face and obtains the four-component numerical flux.

面上的四分量通量

The four-component flux on a face

第 1 节的积分形式要求每个控制体边界上的净通量为零。控制体的边界由面组成；第 2 节中面 $f$ 的面积向量 $\boldsymbol n_f=\ell_f\hat{\boldsymbol n}_f$ 已经包含面长，所以面上的物理通量就是第 1 节的 $F_{\hat n}$ 乘以面长：

The integral form of Section 1 asks for zero net flux through the boundary of every control volume. That boundary consists of faces, and the Section 2 area vector $\boldsymbol n_f=\ell_f\hat{\boldsymbol n}_f$ of face $f$ already includes its length, so the physical flux through a face is the $F_{\hat n}$ of Section 1 times the face length:

首版使用一阶分片常值重构与固定耗散系数 $\alpha=2$ 的 Rusanov 形式通量。固定 $\alpha$ 是本例明确选定的数值方法，因此它对状态的导数为零；不是对一个本应依赖状态的最大波速偷偷冻结导数。页面检查 $|u_n|+c<\alpha$。代价是耗散很强：$M_\infty=0.2$ 时 $\alpha$ 是来流速度的 8.5 倍。第 7 节和第 15 节量化这对压力和载荷意味着什么。

This first version uses piecewise-constant reconstruction and a Rusanov-form flux with the fixed dissipation coefficient $\alpha=2$. A constant $\alpha$ is part of the specified method, so its state derivative is zero; it is not an omitted derivative of a state-dependent maximum wave speed. The page checks $|u_n|+c<\alpha$. The price is strong dissipation: at $M_\infty=0.2$, $\alpha$ is 8.5 times the free-stream speed. Sections 7 and 15 quantify what this does to pressures and loads.

图 3.1 把上式拆成代码里的几步：两侧各自从守恒量算出原始变量和物理通量，再与只依赖 $U_R-U_L$ 的耗散项组合。图 3.2 在一个具体的面上给出两部分的数值。

Figure 3.1 splits this formula into the steps of the code: each side computes its primitive variables and physical flux from the conserved variables, which are then combined with a dissipation term that depends only on $U_R-U_L$. Figure 3.2 gives the two parts on one particular face.

逐步推导：压力、声速与它们对 $U$ 的梯度

Step by step: pressure, sound speed and their gradients with respect to $U$

第一步，把速度写成守恒量之比 $u=(\rho u)/\rho$、$v=(\rho v)/\rho$，动能就是 $((\rho u)^2+(\rho v)^2)/(2\rho)$。代入第 1 节的状态方程，压力是四个守恒量的显式函数（代码 prim() 算出的 p）：

Step 1, write the velocity as ratios of conserved variables, $u=(\rho u)/\rho$ and $v=(\rho v)/\rho$; the kinetic energy is $((\rho u)^2+(\rho v)^2)/(2\rho)$. Substituted into the equation of state of Section 1, the pressure is an explicit function of the four conserved variables (the p computed by prim() in the code):

第二步，对每个分量求偏导，其余三个固定。只有动能项依赖 $\rho$、$\rho u$、$\rho v$：

Step 2, differentiate with respect to each component, holding the other three fixed. Only the kinetic energy depends on $\rho$, $\rho u$ and $\rho v$:

这正是第 8 节 $\dot p$ 公式中各项的系数。第三步，声速 $c=\sqrt{\gamma p/\rho}$。对 $c^2=\gamma p/\rho$ 两边求前向导数，$2c\,\dot c=\gamma(\dot p/\rho-p\,\dot\rho/\rho^2)$，再用 $\gamma p/\rho=c^2$：

These are exactly the coefficients of the $\dot p$ formula of Section 8. Step 3, the sound speed is $c=\sqrt{\gamma p/\rho}$. The forward derivative of $c^2=\gamma p/\rho$ is $2c\,\dot c=\gamma(\dot p/\rho-p\,\dot\rho/\rho^2)$, and with $\gamma p/\rho=c^2$:

第 4 节的壁面通量、第 7 节的目标函数、第 8 节的通量 Jacobian 和远场的 Riemann 不变量 $r_\pm=u_n\pm2c/(\gamma-1)$ 都要用这两个梯度。

The wall flux of Section 4, the objective of Section 7, the flux Jacobian of Section 8 and the far-field Riemann invariants $r_\pm=u_n\pm2c/(\gamma-1)$ all use these two gradients.

核对：第 4 节 riemann() 的计算图与残差共用同一个 prim()；由它给出的 $\partial(r_+-r_-)/\partial U$ 与 $\partial S/\partial U$ 按 $c=\tfrac{\gamma-1}4(r_+-r_-)$、$p=\rho^\gamma S$ 换算出 $\partial c/\partial U$ 与 $\partial p/\partial U$，与上面两式在全部 24 个控制体的收敛状态上最大相差 4×10⁻¹⁶（本页代码在构建时计算）。

Check: the computation graph of riemann() in Section 4 shares prim() with the residual; its $\partial(r_+-r_-)/\partial U$ and $\partial S/\partial U$, converted with $c=\tfrac{\gamma-1}4(r_+-r_-)$ and $p=\rho^\gamma S$, give $\partial c/\partial U$ and $\partial p/\partial U$, which differ from the two formulas above by at most 4×10⁻¹⁶ over the converged states of all 24 volumes (computed by this page’s code at build time).

# U = [rho, rho*u, rho*v, rho*E];  n = area vector (length l)
un_l = (U[1]*nx + U[2]*ny)/U[0]              # l*u_n
p = (gamma-1)*(U[3]-(U[1]**2+U[2]**2)/(2*U[0]))
Fn = [U[0]*un_l, U[1]*un_l+p*nx, U[2]*un_l+p*ny, (U[3]+p)*un_l]
H = (Fn(UL,n)+Fn(UR,n))/2 - alpha*length(n)*(UR-UL)/2

内部面两侧都有控制体，$U_L$ 与 $U_R$ 都是未知量。边界面只有一侧：物面外面是物体，远场外面没有网格。第 4 节给这两类面补上通量：物面只传递压力，远场按第 1 节的特征计数，从来流和域内各取一部分信息，重建外侧状态 $U_b$。

An interior face has a control volume on each side, and both $U_L$ and $U_R$ are unknowns. A boundary face has only one: beyond the wall lies the body, and beyond the far field there is no mesh. Section 4 completes the flux on both kinds of face: the wall passes only the pressure, and the far field rebuilds an exterior state $U_b$ from the free stream and the interior, taking from each as many pieces of information as the characteristic count of Section 1 allows.

边界：滑移壁面与亚声速特征远场

Boundaries: slip walls and subsonic characteristic far field

第 1 节按特征速度的符号数过条件：物面只要一个条件 $u_n=0$，远场入流面要三个，出流面要一个。本节把它们写成边界面上的通量。

Section 1 counted the conditions by the signs of the characteristic speeds: the wall needs a single condition, $u_n=0$, a far-field inflow face three and an outflow face one. This section turns them into fluxes on the boundary faces.

这里与二维标量最重要的区别在物面。无穿透只令质量和能量通量为零；压力仍传递两个方向的动量。壁面通量采用内侧压力 $p_L$，不需要把四个守恒量强行赋值。格心和格点都在边界面上弱施加这条通量规则；格点位于物面并不意味着可以把四行全换成 Dirichlet 约束。

The body boundary is the key difference from the scalar example. No penetration removes the mass and energy fluxes, but pressure still transfers momentum. The wall flux uses the interior pressure $p_L$ without assigning all four conserved values. Both schemes impose this flux weakly; placing a node on the wall does not justify replacing all four equations by Dirichlet constraints.

上面的点表示前向（tangent）导数：$\dot x$ 是程序变量 $x$ 沿给定输入方向的导数，这里 $\dot p$ 来自状态，$\dot n_x$、$\dot n_y$ 来自几何。反向（伴随）量写作 $\bar x$。第 8 节完整说明这套自动微分记号。

The dot marks a forward (tangent) derivative: $\dot x$ is the derivative of the program variable $x$ along a given input direction; here $\dot p$ comes from the state and $\dot n_x$, $\dot n_y$ from the geometry. Reverse (adjoint) quantities are written $\bar x$. Section 8 sets out this algorithmic-differentiation notation in full.

图 4.1 在一个物面面上对比这条规则与内部面的公式。

Figure 4.1 compares this rule with the interior-face formula on one wall face.

逐步推导：壁面通量的 Jacobian 与它的转置

Step by step: the Jacobian of the wall flux and its transpose

第一步，壁面通量只通过压力依赖内侧状态，另外依赖面积向量（代码 flux() 中 bc 为 wall 的分支）：

Step 1, the wall flux depends on the interior state only through the pressure, and otherwise on the area vector (the wall branch of flux() in the code):

第二步，用第 3 节的 $\partial p/\partial U$ 做链式法则。第 0 行与第 3 行恒为零，第 1、2 行是同一个行向量乘以 $n_x$、$n_y$，所以这个 $4\times4$ 块的秩是 1：

Step 2, apply the chain rule with $\partial p/\partial U$ from Section 3. Rows 0 and 3 are identically zero and rows 1 and 2 are the same row vector times $n_x$ and $n_y$, so this $4\times4$ block has rank one:

第三步，转置。反向扫描把面通量的伴随 $\bar H$（四个分量）乘以这个块的转置；两行非零的块转置后只剩一个组合：

Step 3, transpose. The reverse sweep multiplies the adjoint $\bar H$ of the face flux (four components) by the transposed block; with two nonzero rows, only one combination survives:

质量和能量两个分量 $\bar H_\rho$、$\bar H_{\rho E}$ 根本不经过物面传回；经过物面的只有动量分量的伴随在法向上的投影 $\bar p=n_x\bar H_{\rho u}+n_y\bar H_{\rho v}$，也就是壁面压力的伴随，它再按 $\partial p/\partial U$ 分给四个守恒量。第 16 节的连续物面条件约束的正是同一个组合 $\lambda_{\rho u}\hat n_x+\lambda_{\rho v}\hat n_y$。第四步，几何导数：状态固定时 $\partial H_{\mathrm{wall}}/\partial\boldsymbol n=p\,(\boldsymbol e_{\rho u},\boldsymbol e_{\rho v})$，即上面的 $\dot H_{\mathrm{wall}}$ 公式中的 $p\,\dot n_x$、$p\,\dot n_y$ 两项。

The mass and energy components $\bar H_\rho$ and $\bar H_{\rho E}$ do not return through the wall at all; only the normal projection of the momentum adjoints, $\bar p=n_x\bar H_{\rho u}+n_y\bar H_{\rho v}$, passes: it is the adjoint of the wall pressure, and $\partial p/\partial U$ distributes it over the four conserved variables. The continuous wall condition of Section 16 constrains exactly the same combination $\lambda_{\rho u}\hat n_x+\lambda_{\rho v}\hat n_y$. Step 4, the geometry derivative: at fixed state $\partial H_{\mathrm{wall}}/\partial\boldsymbol n=p\,(\boldsymbol e_{\rho u},\boldsymbol e_{\rho v})$, which are the terms $p\,\dot n_x$ and $p\,\dot n_y$ of the $\dot H_{\mathrm{wall}}$ formula above.

核对：在全部 8 个物面面上，计算图对单位种子的前向结果与上面的矩阵最大相差 7×10⁻¹⁸；以伴随解 $\psi_L$ 作为 $\bar H$ 时，计算图的逆序结果与第三步的公式最大相差 4×10⁻¹⁹（本页代码在构建时计算）。

Check: on all 8 wall faces the forward results of the computation graph for unit seeds differ from the matrix above by at most 7×10⁻¹⁸; with the adjoint $\psi_L$ as $\bar H$, the reversed graph differs from the formula of Step 3 by at most 4×10⁻¹⁹ (computed by this page’s code at build time).

远场：先旋转到法向与切向

Far field: rotate to normal and tangential components

设 $\hat{\boldsymbol n}$ 为指向域外的单位法向，$\hat{\boldsymbol t}=(-\hat n_y,\hat n_x)$ 为切向。速度分解为 $u_n=\boldsymbol u\cdot\hat{\boldsymbol n}$ 与 $u_t=\boldsymbol u\cdot\hat{\boldsymbol t}$，反过来 $\boldsymbol u=u_n\hat{\boldsymbol n}+u_t\hat{\boldsymbol t}$。

Let $\hat{\boldsymbol n}$ be the unit normal pointing out of the domain and $\hat{\boldsymbol t}=(-\hat n_y,\hat n_x)$ the tangent. Split the velocity into $u_n=\boldsymbol u\cdot\hat{\boldsymbol n}$ and $u_t=\boldsymbol u\cdot\hat{\boldsymbol t}$; conversely $\boldsymbol u=u_n\hat{\boldsymbol n}+u_t\hat{\boldsymbol t}$.

Riemann 不变量从哪里来

Where the Riemann invariants come from

在面附近只保留沿 $\hat{\boldsymbol n}$ 的变化（坐标记为 $\xi$），得到一维 Euler 方程。对等熵流，$dp=c^2d\rho$，而 $c^2=\gamma p/\rho$ 给出 $dc=\tfrac{\gamma-1}2\,c\,d\rho/\rho$，所以 $dp/(\rho c)=2\,dc/(\gamma-1)$。把连续方程乘 $c/\rho$ 后与动量方程相加或相减，得到两个沿特征输运的量；熵函数 $S$ 与切向速度则以速度 $u_n$ 随流输运。

Near the face keep only variations along $\hat{\boldsymbol n}$ (coordinate $\xi$); this gives the one-dimensional Euler equations. For isentropic flow $dp=c^2d\rho$, and $c^2=\gamma p/\rho$ gives $dc=\tfrac{\gamma-1}2\,c\,d\rho/\rho$, so $dp/(\rho c)=2\,dc/(\gamma-1)$. Adding or subtracting $c/\rho$ times the continuity equation to the momentum equation gives two quantities transported along characteristics; the entropy function $S$ and the tangential velocity move with the flow at speed $u_n$.

亚声速时 $u_n+c>0$，$r_+$ 总从域内来；$u_n-c<0$，$r_-$ 总从域外来。$S$ 与 $u_t$ 的速度是 $u_n$：入流面（$u_n<0$）从来流取，出流面（$u_n>0$）从域内取。所以入流面从外部取三项信息、出流面只取一项，不能像标量那样把整个状态一起复制。这个闭合只在局部一维、等熵的意义下是特征关系，不宣称是精确无反射边界。

For subsonic flow $u_n+c>0$, so $r_+$ always comes from inside; $u_n-c<0$, so $r_-$ always comes from outside. $S$ and $u_t$ move at $u_n$: they come from the free stream at an inflow face ($u_n<0$) and from inside at an outflow face ($u_n>0$). An inflow face therefore takes three pieces of outside information and an outflow face only one; copying the whole state as for the scalar is inappropriate. The closure is characteristic only in the local one-dimensional, isentropic sense; it is not claimed to be exactly non-reflecting.

重建边界状态

Rebuilding the boundary state

$\rho_b$ 的公式来自 $c^2=\gamma p/\rho=\gamma S\rho^{\gamma-1}$。远场面的通量仍是第 3 节的 $H(U_L,U_b)$。下表用本页代码逐步算出一个入流面和一个出流面（$M_\infty=0.2$、$\theta=10^\circ$，与页首选择无关；$\gamma=1.4$ 时 $2/(\gamma-1)=5$、$(\gamma-1)/4=1/10$、$1/(\gamma-1)=2.5$）。交互面板可以对任意远场面重复这些步骤。

The $\rho_b$ formula follows from $c^2=\gamma p/\rho=\gamma S\rho^{\gamma-1}$. A far-field face still uses the flux $H(U_L,U_b)$ of Section 3. The table applies these steps with this page’s code to one inflow and one outflow face ($M_\infty=0.2$, $\theta=10^\circ$, independent of the header; with $\gamma=1.4$, $2/(\gamma-1)=5$, $(\gamma-1)/4=1/10$ and $1/(\gamma-1)=2.5$). The interactive panel repeats the steps for any far-field face.