流体线性稳定性基础

什么是稳定性⁰

Lyapunov stability：

首先要定义一个平衡态，之后，系统会围绕这个状态受到扰动。如果系统返回到平衡状态，则认为稳定的；如果系统偏离平衡状态系统（或者更准确地说，这个特定的平衡状态）被认为是不稳定的。在Lyapunov稳定性的定义中，一个无限的时间范围允许返回到平衡。

短时间的Lyapunov stability：（放大因子或增益）

$$G(t)=\max _{\mathbf{q}_{0}} \frac{E(\mathbf{q}(t))}{E\left(\mathbf{q}_{0}\right)}$$

E(q) denotes the kinetic energy of the perturbation q.

$$\frac{d}{d t} \mathbf{q}=\mathbf{L} \mathbf{q}+\mathbf{B} \mathbf{f}$$ $$\frac{d}{d t}\left(\begin{array}{l} q_{1} \\ q_{2} \end{array}\right)=\underbrace{\left(\begin{array}{cc} \frac{1}{100}-\frac{1}{R e} & 0 \\ \mu & -\frac{2}{R e} \end{array}\right)}_{\mathrm{A}}\left(\begin{array}{l} q_{1} \\ q_{2} \end{array}\right)$$ $$\mathbf{q}(t)=\exp (t \mathrm{~L}) \mathbf{q}_{0}$$ $$G(t)=\max _{\mathbf{q}_{0}} \frac{\|\mathbf{q}(t)\|_{E}^{2}}{\left\|\mathbf{q}_{0}\right\|_{e}^{2}}=\max _{\mathbf{q}_{0}} \frac{\left\|\exp (t \mathrm{~L}) \mathbf{q}_{0}\right\|_{E}^{2}}{\left\|\mathbf{q}_{0}\right\|_{F}^{2}}=\|\exp (t \mathrm{~L})\|_{E}^{2}$$

singular value decomposition

$$\underbrace{\left[\begin{array}{llll} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \end{array}\right]}_{M}=\underbrace{\left[\begin{array}{lll} \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot \end{array}\right]}_{U}\underbrace{\left[\begin{array}{cccc} \sigma_{1} & & & \ominus \\ & \sigma_{2} & & \ominus \\ & & \sigma_{3} & \ominus \end{array}\right]}_{\Sigma} \underbrace{\left[\begin{array}{cccc} \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \ominus & \ominus & \ominus & \ominus \end{array}\right]}_{V^{*}}$$ $$M=U \Sigma V^{*}=\sum_{i=1}^{m} \sigma_{i} u_{i} v_{i}^{*}$$ $$v_{i}^{*} v_{j}=\delta_{i j}$$ $$M v_{j}=\sum_{i=1}^{m} \sigma_{i} u_{i} v_{i}^{*} v_{j}=\sigma_{j} u_{j}$$

任意函数可以展开成vj基函数的形式

$$a=\sum_{i=1}^{n} v_{i} c_{i}$$ $$\begin{aligned} M a &=\sum_{i=1}^{m} u_{i} \sigma_{i} v_{i}^{*} a \\ &=\sum_{i=1}^{m} u_{i} \sigma_{i} c_{i} \end{aligned}$$

The maximum amplitude of output for a given input amplitude?

$$\sigma_{1}=\max _{a \neq 0} \frac{\|M a\|}{\|a\|}$$

M的作用如图所示，对圆盘进行拉伸和旋转

SVD截断

$$M \simeq M_{r}=\sum_{i=1}^{r} u_{i} \sigma_{i} v_{i}^{*}$$ $$M^{-1}=V \Sigma^{-1} U^{*}$$ $$M^{+}=V_{r} \Sigma_{r}^{-1} U_{r}^{*}$$ $$x^{\prime}=M^{+} b$$

A least squares solution x0 that minimises $|M x-b|$.

全局稳定性分析方法

我们来考虑N-S方程的统一动力学形式：

$${\partial q\over \partial t} =N(q)$$

N 是 Navier-Stokes operator.

定义base flow

$$N(q_b)=0$$

The temporal evolution of perturbations q’

$$q'=q-q_b$$ $${\partial q'\over \partial t}=L(q_b)q'$$

the linearized Navier-Stokes operator about the base flow, also called the Jacobian of the system

$$L(q_b)={\partial N \over \partial q}|_{q_b}$$

A global linear stability analysis consists of searching modes of the form

$$q'=\hat{q}e^{\lambda t}$$

reduces N-S equation into the eigenproblem:

$$L(q_b)\hat{q}=\lambda\hat{q}$$

The corresponding eigenvalues characterize the stability of $q_b$: if an eigenvalue $\lambda_i$ has a positive real part $\delta_i = real(\lambda_i)$ (positive growth rate), then the base flow is unstable. Imaganry part means the oscillating mode at a frequency $w_i$.

Resolvent analysis for the of pseudo-resonance

另一种全局方法，被称为预解分析，包括研究外部作用力对baseflow的影响流动而不是专注于自我维持的振荡。类似在机械振动领域的谐波激励下的稳态响应。这个提法最初Trefethen引入，1993.

Resolvent analysis最早应用是在稳定性和转捩中的nonmetal analysis中运用。例如在寻找线性方程组的矩阵时间指数的范数，对transient growth的研究，e-pseudospectra, 以及定义线性NS在对外部forcing最大反馈的比例。

基本形式可以写成：

$${\partial q'\over \partial t}=L(q_b)q' + f\\\\ f = \hat{f}e^{iwt} \\\\ q = \hat{q}e^{iwt}$$

代入上式可得：

$$\hat{q}= R\hat{f} \\\\ R=(iwI-L(q_b))^{-1}$$

R 被称为resolvent operator, 或者叫放大因子。

The optimal forcing which maximizes the energy gain G

$$G(\hat{f})={||R\hat{f}||^2\over ||\hat{f}||^2}$$

Mean flow stability vs base flow stability analysis

In some cases (e.g. the cylinder flow), the frequency selection process is strongly driven by nonlinear effects. 因此，线性基流稳定性分析不能正确地预测流体的动力学行为，而这产生不满足RPIF(((real part positive) with a frequency (the imaginary part))属性的情况.

最终的结论是mean flow 比base flow 更有效。

整体程序框架

考虑因素：stability analysis、resolvent analysis、Arbitrary Lagrangian Eulerian methods 可考虑流固耦合的统一框架、implicit newton approach for non-linearity、选用成熟的FVM或者FEM开源C++求解器、小扰动假设，主要关心的问题在于流体或者固体起振之初，比如失速起始、无量纲化

固体域$\Omega_s$：

$$M_{s} \frac{\partial^{2} \hat{\xi}}{\partial t}-\hat{\nabla} \cdot \hat{P}(\hat{\xi})=0$$

其中，$M_s$是质量矩阵，P为Kirchhhoff 应力张量.

另外，对于可压固体，我们还需要对其密度进行修正：

$$\rho^*(t)=\hat{J}(\hat{\xi})\rho^*_s(0)$$

流体域：

暂时考虑不可压的N-S方程：

$${\partial \tilde{u}\over \partial t} + (\tilde{\nabla }\tilde{u})\tilde{u} - \tilde{\nabla}\cdot \tilde{\sigma}(\tilde{u},\tilde{p})=0 \\\\ \tilde{\nabla}\cdot \tilde{u}=0$$

$\sigma$ 为Cauchy 应力张量，跟流体的形变和雷诺数又关系。

$$\tilde{\sigma}(\tilde{u},\tilde{p}) = -\tilde{p}I -{1\over R_e}(\tilde{\nabla}\tilde{u}+\tilde{\nabla}\tilde{u}^T)$$

• 固体域使用Lagrangian框架、流体域使用Eulerian框架，这样的好处是表达更加自然美观。

Arbitrary Lagrangian Eulerian methods

在于考虑相互的交界面：

$$\tilde{x}_t = \begin{cases} \hat{A}_t(\hat{x},t), & 流体 \\ \hat{x}+\hat{\xi}(\hat{x},t), & 固体 \end{cases}$$

因此，得到extension displacement in the reference fluid domain:

$$\hat{A}_t(\hat{x},t)=\hat{x}+\hat{\xi}(\hat{x},t)\\\\ \hat{\xi}_e (\hat{x},t)=\tilde{x}_t-\hat{x}=\hat{A}_t(\hat{x},t)-\hat{x}$$

链式法则，

$${\partial \tilde{u}\over \partial t } |_{x^t}={\partial \tilde{u}\over \partial t}|_{\hat{x}} - （\tilde{\nabla}\tilde{u}）\tilde{w}$$

最终，ALE formulation of N-S equation:

$${\partial \tilde{u}\over \partial t}|_{\hat{x}} + (\tilde{\nabla}\tilde{u})(\tilde{u}-\tilde{w})-\tilde{\nabla}\cdot \tilde{\sigma}(\tilde{u},\tilde{p})=0\\\\ \tilde{\nabla}\cdot \tilde{u}=0$$

对比之前的N-S，我们可以发现

$$\tilde{w} = \begin{cases} 0, & 流体 \\ {\partial \hat{A}_t\over \partial t}|_\hat{x} \cdot \hat{A}_t^{-1} , & 靠近固体的流体域\\ 任意， & 两者之间 \end{cases}$$

弱解形式，应用FEM求解

freefem++ 参考代码：https://github.com/ERC-AEROFLEX/ALE

$$\underbrace{\left(\begin{array}{ccc} \hat{\mathscr{T}} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & -\hat{\mathscr{T}} \mathrm{fe}\left(\hat{\boldsymbol{q}}_{\mathrm{f}}, \hat{\boldsymbol{q}}_{\mathrm{e}}\right) & \hat{\mathscr{T}}_{\mathrm{f}}\left(\hat{\boldsymbol{q}}_{\mathrm{e}}\right) \end{array}\right)}_{\hat{\mathscr{T}}_{\mathrm{fsi}}(\hat{\boldsymbol{q}})} \frac{\partial}{\partial t}\left(\begin{array}{l} \hat{\boldsymbol{q}}_{\mathrm{s}} \\ \hat{\boldsymbol{q}}_{\mathrm{e}} \\ \hat{\boldsymbol{q}}_{\mathrm{f}} \end{array}\right)=\underbrace{\left(\begin{array}{c} \hat{\mathscr{S}}\left(\hat{\boldsymbol{q}}_{\mathrm{s}}\right)+\hat{\mathcal{I}}_{\mathrm{fs}}^{\mathrm{T}} \hat{\boldsymbol{q}}_{\mathrm{f}} \\ -\mathscr{A}_{\mathrm{e}} \hat{\boldsymbol{q}}_{\mathrm{e}}+\hat{\mathcal{I}}_{\mathrm{es}} \hat{\boldsymbol{q}}_{\mathrm{s}} \\ \hat{\mathcal{I}}_{\mathrm{fs}} \hat{\boldsymbol{q}}_{\mathrm{s}}+\mathscr{\mathscr { N }}_{\mathrm{f}}\left(\hat{\boldsymbol{q}}_{\mathrm{f}}, \hat{\boldsymbol{q}}_{\mathrm{e}}\right) \end{array}\right)}_{\hat{\mathscr{N}_{fsi}}(\hat{\boldsymbol{q}})}$$

2.1.2公式转化为稳定性问题：

$$\boldsymbol{q}^{\prime}(\boldsymbol{x}, t)=\boldsymbol{q}^{\circ}(\boldsymbol{x}) \exp (\lambda t)+\mathrm{c.c.}$$ $$\lambda\left[\begin{array}{ccc} \mathbf{B}_{\mathrm{s}} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & \mathbf{B}_{\mathrm{fe}} & \mathbf{B}_{\mathrm{f}} \end{array}\right]\left[\begin{array}{l} \boldsymbol{q}_{\mathrm{S}}^{\circ} \\ \boldsymbol{q}_{\mathrm{e}}^{\circ} \\ \boldsymbol{q}_{\mathrm{f}}^{\circ} \end{array}\right]=\left[\begin{array}{lll} \mathbf{A}_{\mathbf{s}} & 0 & \mathbf{A}_{\mathrm{fs}}^{\mathrm{T}} \\ \mathbf{A}_{\mathrm{es}} & \mathbf{A}_{\mathrm{e}} & 0 \\ \mathbf{A}_{\mathrm{fs}} & \mathbf{A}_{\mathrm{fe}} & \mathbf{A}_{\mathrm{f}} \end{array}\right]\left[\begin{array}{l} \boldsymbol{q}_{\mathrm{S}}^{\circ} \\ \boldsymbol{q}_{\mathrm{e}}^{\circ} \\ \boldsymbol{q}_{\mathrm{f}}^{\circ} \end{array}\right]$$

简化矩阵：

$$\lambda\left[\begin{array}{cc} \mathbf{B}_{\mathrm{f}} & -\mathbf{B}_{\mathrm{fe}} \mathbf{A}_{\mathrm{e}}^{-1} \mathbf{A}_{\mathrm{es}} \\ 0 & \mathbf{B}_{\mathrm{s}} \end{array}\right]\left[\begin{array}{l} \boldsymbol{q}_{\mathrm{f}}^{\circ} \\ \boldsymbol{q}_{\mathrm{s}}^{\circ} \end{array}\right]=\left[\begin{array}{cc} \mathbf{A}_{\mathrm{f}} & \mathbf{A}_{\mathrm{fs}}-\mathbf{A}_{\mathrm{fe}} \mathbf{A}_{\mathrm{e}}^{-1} \mathbf{A}_{\mathrm{es}} \\ \mathbf{A}_{\mathrm{fs}}^{\mathrm{T}} & \mathbf{A}_{\mathrm{s}} \end{array}\right]\left[\begin{array}{l} \boldsymbol{q}_{\mathrm{f}}^{\circ} \\ \boldsymbol{q}_{\mathrm{s}}^{\circ} \end{array}\right]$$ $$\lambda\left[\begin{array}{cc} \mathbf{B}_{\mathrm{f}} & \mathbf{B}_{\mathrm{fs}}^{r} \\ 0 & \mathbf{B}_{\mathrm{s}}^{r} \end{array}\right]\left[\begin{array}{l} \boldsymbol{q}_{\mathrm{f}}^{\circ} \\ \boldsymbol{\alpha}_{\mathrm{s}} \end{array}\right]=\left[\begin{array}{ll} \mathbf{A}_{\mathrm{f}} & \mathbf{A}_{\mathrm{fs}}^{r} \\ \mathbf{A}_{\mathrm{sf}}^{r} & \mathbf{A}_{\mathrm{s}}^{r} \end{array}\right]\left[\begin{array}{l} \boldsymbol{q}_{\mathrm{f}}^{\circ} \\ \boldsymbol{\alpha}_{\mathrm{s}} \end{array}\right]$$

问题？ 如何衔接到n-periodic的模型上面去？

$$\begin{array}{c} \mathcal{M}_{s} \frac{\partial^{2} \xi}{\partial t^{2}}-\nabla \cdot(\boldsymbol{F}(\xi) \boldsymbol{S}(\xi))=\mathbf{0} \quad \text { in } \Omega_{s} \\ -\nabla \cdot\left(h_{x} \nabla \xi_{e}\right)=\mathbf{0} \quad \text { in } \Omega_{e} \\ \boldsymbol{\xi}_{e}-\xi=\mathbf{0} \quad \text { on } \Gamma \end{array}$$

Linear Stability Analysis 更深一层次的理解

这是一个Duffing oscillator动力学系统，用相空间的形式表示。

$$\begin{array}{l} \dot{x}=y \\ \dot{y}=-\frac{1}{2} y+x-x^{3} \end{array}$$

如图，可以看到三个fixed point， 或者称为equilibrium points $\mathbf{X}^{*}$。

找到equilibrium points之后，我们需要来判断这个状态是stable还是unstable的。 Lyapunov stability的解释是，过了很长时间能回到equilibrium points。

首先，假设系统的perturbation满足如下动力学方程：

$$\dot{\mathbf{x}}=\mathcal{F}\left(\mathbf{X}^{*}+\mathbf{x}\right)$$

根据Taylor 公式在平衡点展开，

$$\dot{\mathbf{x}}=\mathcal{A} \mathbf{x}$$ $$\mathcal{A}=\partial \mathcal{F} / \partial \mathbf{X}$$

$\mathcal{A}$ 为Jacobian matrix of F.

也就是说，假设给定初始条件x_0, x可以解析得到：

$$\mathbf{x}(t)=\exp (\mathcal{A} t) \mathbf{x}_{0}$$

我们对A进行spectral decomposition：

$$\mathcal{A}=\mathcal{V} \mathbf{\Lambda} \mathcal{V}^{-1}$$

得到x的另外一种表达形式：

$$\mathbf{x}(t)=\mathcal{V} \exp (\boldsymbol{\Lambda} t) \mathcal{V}^{-1} \mathbf{x}_{0}$$

可以发现，其特征值如果按照实部从大到小排列，

$$\lim _{t \rightarrow+\infty} \exp (\mathcal{A} t) \mathbf{x}_{0}=\lim _{t \rightarrow+\infty} \exp \left(\lambda_{1} t\right) \mathbf{v}_{1}$$

也就是说，$\lambda_1$就决定了系统的稳定性，如果实部大于0，那么系统就会呈指数增长，即线性不稳定；如果实数小于0，那么系统就会指数衰减，即线性稳定。实数等于0的情况较为特殊，无法直接判断，需要再继续判断A的二阶Taylor展开。

这里主要针对长时间的稳定性问题，短时间内需要后面再进行理论的探讨。

Non-modal Stability Analysis 更深一层次的理解

给一个例子，

$$\frac{\mathrm{d}}{\mathrm{d} t}\left[\begin{array}{l} x_{1} \\ x_{2} \end{array}\right]=\underbrace{\left[\begin{array}{cc} \frac{1}{100}-\frac{1}{R e} & 0 \\ 1 & -\frac{2}{R e} \end{array}\right]}_{\mathcal{A}}\left[\begin{array}{l} x_{1} \\ x_{2} \end{array}\right]$$

其特征值为

$$\lambda_{1}=\frac{1}{100}-\frac{1}{R e}$$ $$\lambda_{2}=-\frac{2}{R e}$$

根据特征值，可以判断长时间内的衰减趋势。但是短时间内，两者都大幅度的指数放大。这是在之前Stability Analysis无法解释的。

这是由于矩阵A的non-normality所带来的：

$$\mathcal{A}^{\dagger} \mathcal{A} \neq \mathcal{A} \mathcal{A}^{\dagger}$$

$\mathcal{A}^{\dagger} $为伴随矩阵，它具有非常好的性质。关于伴随矩阵和逆矩阵的关系，和其性质，参考链接。

$$A\cdot\operatorname{adj}(A) =\operatorname{adj}(A) \cdot A=|A| . I$$ $$|\operatorname{adj}(\mathrm{A})|=|\mathrm{A}|^{\mathrm{n}-1}$$ $$\operatorname{adj}(A B)=\operatorname{adj}(B) \cdot \operatorname{adi}(A)$$ $$(\operatorname{adj} A)^{\top}=\operatorname{adj}\left(A^{\top}\right)$$ $$\operatorname{adj}(k A)=k^{n-1} \operatorname{adj}(A)$$

As a result of this non-normality, the eigenvectors of A do not form an orthonormal set of vectors。

长时间的过程，都会趋向于线性稳定。但是，non-modal在短时间内可能会造成瞬时复制增大。本质是应为non-modal的矩阵性质，分解得到非正交的特征，物理上表现在更初始值的选取有关。但一旦逼近平衡点，就会按照平衡点的性质稳定下来。

从图中解释，目前这两个特征值是非正交的，初始的选取位置不一样，造成短期的增长效应。但如果是正交的，那么类似于垂直的网格布置，两个特征值的互不影响，这样也就是说，直接取消于稳定点的趋势。

Optimal Perturbation Analysis 紧接着引出

也就是，即使是stable的点逼近，但是短时间内也可能产生大的放大响应，找到最大的放大点是一个非常有意思的问题。

两个策略：

求优化问题

$$\begin{array}{c} \underset{\mathbf{x}_{0}}{\operatorname{maximize}} \mathcal{J}\left(\mathbf{x}_{0}\right)=\|\mathbf{x}(T)\|_{2}^{2} \\ \text { subject to } \dot{\mathbf{x}}-\mathcal{A} \mathbf{x}=0 \\ \left\|\mathbf{x}_{0}\right\|_{2}^{2}-1=0 \end{array}$$

并转化为：

$$\underset{\mathbf{x}, \mathbf{v}, \mu}{\operatorname{maximize}} \mathcal{L}(\mathbf{x}, \mathbf{v}, \mu)$$ $$\mathcal{L}(\mathbf{x}, \mathbf{v}, \mu)=\mathcal{J}\left(\mathbf{x}_{0}\right)+\int_{0}^{T} \mathbf{v}^{T}(\dot{\mathbf{x}}-\mathcal{A} \mathbf{x}) \mathrm{d} t+\mu\left(\left\|\mathbf{x}_{0}\right\|_{2}^{2}-1\right)$$ $$\delta_{\mathbf{x}} \mathcal{L}=\left[\nabla_{\mathbf{x}} \mathcal{J}+\mathbf{v}(T)\right] \cdot \delta \mathbf{x}(0)+\int_{0}^{T}\left[\dot{\mathbf{v}}-\mathcal{A}^{\dagger} \mathbf{v}\right] \cdot \delta \mathbf{x} \mathrm{dt}+\left[2 \mu \mathbf{x}_{0}-\mathbf{v}(0)\right] \cdot \delta \mathbf{x}(0)$$

SVD

$$\underset{\mathbf{x}_{0}}{\operatorname{maximize}} \frac{\|\mathbf{x}(T)\|_{2}^{2}}{\left\|\mathbf{x}_{0}\right\|_{2}^{2}}$$ $$\begin{aligned} \mathcal{G}(T) &=\max _{\mathbf{x}_{0}} \frac{\left\|\exp (\mathcal{A} t) \mathbf{x}_{0}\right\|_{2}^{2}}{\left\|\mathbf{x}_{0}\right\|_{2}^{2}} \\ &=\|\exp (\mathcal{A} T)\|_{2}^{2} \end{aligned}$$

引入SVD分解：

$$exp (\mathcal{A} T)=\mathcal{U} \Sigma \mathcal{V}^{H}$$

那么，最大的特征值即为最大能量响应：

$$\mathcal{G}(T)=\sigma_{1}^{2}$$

其对应的初始x_0 由对应的右奇异特征值给出。

Resolvent Analysis

$$\dot{\mathbf{x}}=\mathcal{A} \mathbf{x}+\mathbf{f}$$

the response of the system to the forcing f(t) 表现为卷积的形式：

$$\mathbf{x}(t)=\int_{0}^{t} \exp (\mathcal{A}(t-\tau)) \mathbf{f}(\tau) \mathrm{d} \tau$$

Consider linear stable systems，influence of the forcing on the current state decays exponentially according to the least stable eigenvalue。

$$\mathbf{f}(t)=\Re\left(\hat{\mathbf{f}} e^{i \omega t}\right)$$

根据卷积和傅立叶变换的性质，可以很容易的得到响应矩阵，这个矩阵也叫做Resolvent operator：

$$\hat{\mathbf{x}}=(i \omega \mathcal{I}-\mathcal{A})^{-1} \hat{\mathbf{f}}$$

类似，寻找最大的响应放大因子：

$$\begin{aligned} \mathcal{R}(\omega) &=\max _{\hat{\mathbf{f}}} \frac{\left\|(i \omega \mathcal{I}-\mathcal{A})^{-1} \hat{\mathbf{f}}\right\|_{2}^{2}}{\|\hat{\mathbf{f}}\|_{2}^{2}} \\ &=\|\mathcal{R}(\omega)\|_{2}^{2} \end{aligned}$$

同样，可以利用SVD寻找到最大的奇异值：

$$\mathcal{R}(\omega)=\sigma_{1}^{2}$$

线化不可压N-S方程

$$\begin{aligned} \partial_{t} \mathbf{u} &=-\nabla p-\mathbf{u} \cdot \nabla \mathbf{u}+\frac{1}{R e} \nabla^{2} \mathbf{u} \\ \boldsymbol{\nabla} \cdot \mathbf{u} &=0 \end{aligned}$$ $$R e=\frac{U_{b u l k} D}{\nu}$$

将均匀项和脉动量分开

$$\begin{aligned} &\left[\partial_{t} \tilde{\mathbf{u}}(\mathbf{x}, t)\right.\\ &+\tilde{\mathbf{u}}(\mathbf{x}, t) \cdot \nabla \tilde{\mathbf{u}}(\mathbf{x}, t)+\overline{\mathbf{u}}(\mathbf{x}) \cdot \nabla \overline{\mathbf{u}}(\mathbf{x}) \\ &+\overline{\mathbf{u}}(\mathbf{x}) \cdot \nabla \tilde{\mathbf{u}}(\mathbf{x}, t)+\tilde{\mathbf{u}}(\mathbf{x}, t) \cdot \nabla \overline{\mathbf{u}}(\mathbf{x}) \\ &+\nabla(\bar{p}(\mathbf{x})+\tilde{p}(\mathbf{x}, t)) \\ &\left.-\frac{1}{R_{0}} \nabla^{2}(\overline{\mathbf{u}}(\mathbf{x})+\tilde{\mathbf{u}}(\mathbf{x}, t))\right] =0 \end{aligned}$$

fluaction equation:

$$\begin{array}{l} i \omega \hat{\mathbf{u}}(\mathbf{x}, \omega)+\overline{\mathbf{u}}(\mathbf{x}) \cdot \nabla \hat{\mathbf{u}}(\mathbf{x}, \omega)+\hat{\mathbf{u}}(\mathbf{x}, \omega) \cdot \nabla \overline{\mathbf{u}}(\mathbf{x}) \\ =-\nabla \hat{p}(\mathbf{x}, \omega)+\frac{1}{R e} \nabla^{2} \hat{\mathbf{u}}(\mathbf{x}, \omega)-\hat{\mathbf{f}}(\mathbf{x}, \omega) \end{array}$$

RANS equation: (w=0)

$$\overline{\mathbf{u}}(\mathbf{x}) \cdot \nabla \overline{\mathbf{u}}(\mathbf{x})=-\nabla \bar{p}(\mathbf{x})+\frac{1}{R e} \nabla^{2} \overline{\mathbf{u}}(\mathbf{x})-\overline{\mathbf{f}}(\mathbf{x})$$

消除压力项：

pressure Poisson equation

$$-\nabla^{2} p(\mathbf{x}, t)=\nabla \cdot(\mathbf{u}(\mathbf{x}, t) \cdot \nabla \mathbf{u}(\mathbf{x}, t))$$

分解得到：

mean component

$$-\nabla^{2} \bar{p}(\mathbf{x})=\nabla \cdot(\overline{\mathbf{u}}(\mathbf{x}) \cdot \nabla \overline{\mathbf{u}}(\mathbf{x}))+\nabla \cdot \overline{\mathbf{f}}(\mathbf{x})$$

frequency component

$$\nabla^{2} \hat{p}(\mathbf{x}, \omega)=-\nabla \cdot L \hat{\mathbf{u}}(\mathbf{x}, \omega)-\nabla \cdot \hat{\mathbf{f}}(\mathbf{x}, \omega)$$ $$L=\overline{\mathbf{u}}(\mathbf{x}) \cdot \nabla+(\nabla \overline{\mathbf{u}}(\mathbf{x}))^{T}$$

代入到前面的fluaction equation，消去压力项：

$$\left(i \omega I+\Pi L-\frac{1}{R e} \nabla^{2}\right) \hat{\mathbf{u}}(\mathbf{x}, \omega)=-\Pi \hat{\mathbf{f}}(\mathbf{x}, \omega)$$ $$\Pi=\left(I-\nabla\left(\nabla^{2}\right)^{-1} \nabla \cdot\right)$$ $$\hat{\mathbf{f}}(\mathbf{x}, \omega)=\int_{-\infty}^{\infty} e^{-i \omega t} \tilde{\mathbf{u}}(\mathbf{x}, t) \cdot \nabla \tilde{\mathbf{u}}(\mathbf{x}, t) d t$$

前面的推导的具体过程参考⁵

$$R(\omega):=-\left(i \omega I+\Pi L-\frac{1}{R e} \nabla^{2}\right)^{-1} \Pi$$ $$\hat{\mathbf{u}}(\mathbf{x}, \omega)=R(\omega) \hat{\mathbf{f}}(\mathbf{x}, \omega)$$

最后总结框图如下：

$$\mathbf{f}=\tilde{\mathbf{u}} \cdot \nabla \tilde{\mathbf{u}}$$

首先，我们将速度和压力分解，

$$\left[\begin{array}{l} \mathbf{u}(x, y, z, t) \\ p(x, y, z, t) \end{array}\right]=\iiint_{-\infty}^{\infty}\left[\begin{array}{l} \mathbf{u}_{\mathbf{k}}(y) \\ p_{\mathbf{k}}(y) \end{array}\right] e^{i\left(\kappa_{x} x+\kappa_{z} z-\omega t\right)} d \kappa_{x} d \kappa_{z} d \omega$$ $$\begin{array}{c} (-i \omega+i k U) \mathbf{u}_{\mathbf{k}}+v_{\mathbf{k}} U^{\prime} \mathbf{e}_{\mathbf{x}}+\nabla p_{\mathbf{k}}-R e^{-1} \nabla^{2} \mathbf{u}_{\mathbf{k}}=\mathbf{f}_{\mathbf{k}} \\ \nabla \cdot \mathbf{u}_{\mathbf{k}}=0 \end{array}$$

预解算子：

$$\begin{aligned} \left[\begin{array}{c} \mathbf{u}_{\mathbf{k}} \\ p_{\mathbf{k}} \end{array}\right] &=\left(-i \omega\left[\begin{array}{ll} \mathbf{I} & \\ & 0 \end{array}\right]-\left[\begin{array}{cc} \mathcal{L}_{\mathbf{k}} & -\nabla_{\mathbf{k}} \\ \nabla_{\mathbf{k}}^{T} & 0 \end{array}\right]\right)^{-1}\left[\begin{array}{l} \mathbf{I} \\ 0 \end{array}\right] \mathbf{f}_{\mathbf{k}} \\ &=\tilde{\mathcal{H}}_{\mathbf{k}} \mathbf{f}_{\mathbf{k}} \end{aligned}$$ $$\nabla_{\mathbf{k}}=\left[i \kappa_{x}, \partial / \partial y, i \kappa_{z}\right]^{T}$$ $$\mathbf{f}_{\mathbf{k}}=(-\mathbf{u} \cdot \mathbf{\nabla} \mathbf{u})_{\mathbf{k}}$$ $$\mathcal{L}_{\mathbf{k}}=\left[\begin{array}{ccc} -i \kappa_{x} U+R e_{\tau}^{-1} \nabla_{\mathbf{k}}^{2} & -\partial U / \partial y & 0 \\ 0 & -i \kappa_{x} U+R e_{\tau}^{-1} \nabla_{\mathbf{k}}^{2} & 0 \\ 0 & 0 & -i \kappa_{x} U+R e_{\tau}^{-1} \nabla_{\mathbf{k}}^{2} \end{array}\right]$$

对管道坐标系来说：

$$\mathcal{L}_{\mathbf{k}}=\left[\begin{array}{ccc} -i k U+\frac{D+r^{-2}}{R e} & -\frac{\partial U}{\partial r} & 0 \\ 0 & -i k U+\frac{D}{R e} & -\frac{2 i n r^{-2}}{R e} \\ 0 & \frac{2 i n r^{-2}}{R e} & -i k U+\frac{D}{R e} \end{array}\right]$$ $$D=-k^{2}-\left(n^{2}+1\right) r^{-2}+\partial_{r}^{2}+r^{-1} \partial_{r}$$

见代码：https://github.com/mluhar/resolvent

McKeon & Sharma (2010), 做的是在管道里面。

预解算子的边界条件

硬壁面无滑移边界条件

对于软壁面边界条件的思路🤔：

Complex wall admittance linking pressure and velocity。²

Impose boundary conditions for Navier-Stokes Resolvent This accounts for a simple compliant surface model, which essentially leads to a harmonic-motion relationship between pressure and velocity。

探究谐波响应激起的最大压力响应，而非速度 ¹

possion equation:

$$\nabla^{2} p_{\mathbf{k}}=-2 i k v_{\mathbf{k}} U^{\prime}+\nabla \cdot \mathbf{f}_{\mathbf{k}}$$

boundary condition at the wall:

$$\left.\frac{\partial p_{\mathbf{k}}}{\partial r}\right|_{r=1}=\frac{1}{R e}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial v_{\mathbf{k}}}{\partial r}\right)\right]_{r=1}$$ $$\begin{array}{ll} \nabla^{2} p_{\mathbf{k}, f}=-2 i k v_{\mathbf{k}} U^{\prime}, & \left.\frac{\partial p_{\mathbf{k}, f}}{\partial r}\right|_{r=1}=0 \\ \nabla^{2} p_{\mathbf{k}, s}=\nabla \cdot \mathbf{f}_{\mathbf{k}}, & \left.\frac{\partial p_{\mathbf{k}, s}}{\partial r}\right|_{r=1}=0 \\ \nabla^{2} p_{\mathbf{k}, s t}=0, & \left.\frac{\partial p_{\mathbf{k}, s t}}{\partial r}\right|_{r=1}=\frac{1}{R e}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial v_{\mathbf{k}}}{\partial r}\right)\right]_{r=1} \end{array}$$ $$\begin{array}{c} \nabla^{2} G_{\mathbf{k}}=\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial G_{\mathbf{k}}}{\partial r}\right)-\left(k^{2}+\frac{n^{2}}{r^{2}}\right) G_{\mathbf{k}}=\frac{1}{r} \delta\left(r-r^{\prime}\right) \\ \left.\frac{\partial G_{\mathbf{k}}}{\partial r}\right|_{r=1}=0 \end{array}$$

构建格林公式满足：

$$\begin{array}{c} \nabla^{2} G_{\mathbf{k}}=\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial G_{\mathbf{k}}}{\partial r}\right)-\left(k^{2}+\frac{n^{2}}{r^{2}}\right) G_{\mathbf{k}}=\frac{1}{r} \delta\left(r-r^{\prime}\right) \\ \left.\frac{\partial G_{\mathbf{k}}}{\partial r}\right|_{r=1}=0 \end{array}$$ $$\begin{aligned} &G_{\mathbf{k}}\left(r, r^{\prime}\right)=\left\{\begin{array}{ll} -A_{\mathbf{k}} I_{n}\left(k r^{\prime}\right) I_{n}(k r)-I_{n}\left(k r^{\prime}\right) K_{n}(k r), & 1 \geq r \geq r^{\prime} \\ -A_{\mathbf{k}} I_{n}\left(k r^{\prime}\right) I_{n}(k r)-I_{n}(k r) K_{n}\left(k r^{\prime}\right), & 0 \leq r \leq r^{\prime} \end{array}\right.\\ &\text { with }\\ &A_{\mathbf{k}}=\frac{K_{n-1}(k)+K_{n+1}(k)}{I_{n-1}(k)+I_{n+1}(k)} \end{aligned}$$

可以解析得到压力对格林公式的表达式：

$$\begin{array}{c} p_{\mathbf{k}, f}(r)=\int_{0}^{1} G_{\mathbf{k}}\left(r, r^{\prime}\right)\left(-\left.2 i k \hat{v}_{\mathbf{k}} \frac{\partial U}{\partial r}\right|_{r^{\prime}}\right) r^{\prime} d r^{\prime} \\ p_{\mathbf{k}, s}(r)=\int_{0}^{1} G_{\mathbf{k}}\left(r, r^{\prime}\right)\left(\nabla \cdot \mathbf{f}_{\mathbf{k}}\right) r^{\prime} d r^{\prime} \end{array}$$ $$p_{\mathbf{k}, s t}(r)=\left[\frac{I_{n}(k r)}{I_{n-1}(k)+I_{n+1}(k)}\right] \frac{2}{k R e}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial v_{\mathbf{k}}}{\partial r}\right)\right]_{r=1}$$

加强筋边界-channel

Solid obstructions within the fluid domain—riblets in our case—are modeled as a spatially varying permeability function Kx; y; z,

which appears as an additional body force in the momentum equations. ³

【方程1】

$$\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \cdot \nabla \boldsymbol{u}=-\nabla p+\frac{1}{\operatorname{Re}_{\tau}} \nabla^{2} \boldsymbol{u}-K^{-1} \boldsymbol{u}, \quad \nabla \cdot \boldsymbol{u}=0$$

在流体域，K趋向于∞；在固体域，K趋向于0；

【方程2】

$$K^{-1}(y, z)=\sum_{n=-\infty}^{\infty} a_{n}(y) \exp \left\{i n \kappa_{s} z\right\}$$

同样的，meanflow也会呈现类似的周期性特征：

【方程3】

$$U(y, z)=\sum_{n=-\infty}^{\infty} U_{n}(y) \exp \left\{i n \kappa_{s} z\right\}$$

在此之上，脉动小量被提取出来：可以看到，均匀流与x的方向无关

【方程4】

$$\boldsymbol{u}(\boldsymbol{x}, t)=\boldsymbol{U}(y, z)+\boldsymbol{u}^{\prime}(\boldsymbol{x}, t)$$

代入动量方程，均匀量满足关系：

【方程5】

$$-\frac{\partial P}{\partial x}+\frac{1}{R e_{\tau}} \nabla^{2} U-K^{-1} U-\left(\nabla \cdot\left(\overline{u^{\prime} u^{\prime}}\right)\right)_{x}=0$$

脉动量满足：

【方程6】

$$\begin{array}{l} \frac{\partial \boldsymbol{u}^{\prime}}{\partial t}+\boldsymbol{U} \cdot \nabla \boldsymbol{u}^{\prime}+\boldsymbol{u}^{\prime} \cdot \nabla \boldsymbol{U}+\nabla \cdot\left(\boldsymbol{u}^{\prime} \boldsymbol{u}^{\prime}-\overline{\boldsymbol{u}^{\prime} \boldsymbol{u}^{\prime}}\right) \\ \quad=-\nabla p^{\prime}+\frac{1}{\boldsymbol{R} e_{\tau}} \nabla^{2} \boldsymbol{u}^{\prime}-K^{-1} \boldsymbol{u}^{\prime} \end{array}$$ $$\nabla \cdot \boldsymbol{u}^{\prime}=0$$

将总结的均匀流的周期性特性【方程2、3】代入到【方程5】，得：

【方程7】

$$\begin{array}{l} \frac{1}{\operatorname{Re}_{\tau}}\left(\left(\nu_{e} \frac{d^{2}}{d y^{2}}+\frac{d \nu_{e}}{d y} \frac{d}{d y}\right) \delta(n)+\frac{d^{2}}{d y^{2}}-n^{2} \kappa_{s}^{2}\right) U_{n} \\ -\sum_{p+q=n} U_{p} a_{q}=\frac{\partial P}{\partial x} \delta(n) \end{array}$$

脉动量方程6的分解形式：

【方程8】

$$\begin{array}{l} \left(-i \omega+i \kappa_{x} U_{0}+a_{0}-\frac{1}{R e_{\tau}} \tilde{\nabla}^{2}\right) \boldsymbol{u}_{\boldsymbol{k}}+\tilde{\nabla} p_{\boldsymbol{k}}+v_{\boldsymbol{k}} \frac{\partial U_{0}}{\partial y} \boldsymbol{e}_{x} \\ \quad+\sum_{n=-\infty \atop n \neq 0}^{\infty}\left(a_{n}+i \kappa_{x} U_{n}\right) \boldsymbol{u}_{k-n \kappa_{s}} \\ \quad+\sum_{n=-\infty \atop n \neq 0}^{\infty}\left(v_{k-n \kappa_{s}} \frac{\partial}{\partial y}+i n \kappa_{s} w_{k-n \kappa_{s}}\right) U_{n} \boldsymbol{e}_{x}=\boldsymbol{f}_{\boldsymbol{k}} \end{array}$$

方程8可以整理成矩阵的形式：

【方程9】

$$\begin{array}{l} \left(-i \omega\left[\begin{array}{ll} \boldsymbol{I} & \\ & 0 \end{array}\right]-\left[\begin{array}{cc} \boldsymbol{L}_{k}-a_{0} \boldsymbol{I} & -\tilde{\nabla} \\ \tilde{\nabla}^{T} & 0 \end{array}\right]\right)\left[\begin{array}{l} \boldsymbol{u}_{k} \\ p_{\boldsymbol{k}} \end{array}\right]+\sum_{n=-\infty, \atop n \neq 0}^{\infty} \boldsymbol{F}_{n}\left[\begin{array}{l} \boldsymbol{u}_{k-n \kappa_{s}} \\ p_{\boldsymbol{k}-n \kappa_{s}} \end{array}\right] \\ =\left[\begin{array}{l} \boldsymbol{I} \\ 0 \end{array}\right] \boldsymbol{f}_{\boldsymbol{k}} \end{array}$$ $$\boldsymbol{L}_{\boldsymbol{k}}=\left[\begin{array}{ccc} -i \kappa_{x} U_{0}+R e_{\tau}^{-1} \tilde{\nabla}^{2} & -d U_{0} / d y & 0 \\ 0 & -i \kappa_{x} U_{0}+R e_{\tau}^{-1} \tilde{\nabla}^{2} & 0 \\ 0 & 0 & -i k_{x} U_{0}+R e_{\tau}^{-1} \tilde{\nabla}^{2} \end{array}\right]$$ $$\boldsymbol{F}_{n}=\left[\begin{array}{cccc} a_{n}+i \kappa_{x} U_{n} & d U_{n} / d y & i n \kappa_{s} U_{n} & 0 \\ & a_{n}+i k_{x} U_{n} & 0 & 0 \\ & & a_{n}+i k_{x} U_{n} & 0 \\ 0 & & & 0 \end{array}\right]$$

可以发现，k和k_s空间波数耦合。k_s在压气机中类似叶片扫过的均匀流。nk_s 表示其多阶谐波。

将方程9整理成输入和输出的形式：

$$\tilde{\boldsymbol{I}}=\left[\begin{array}{ll} \boldsymbol{I} & \\ & 0 \end{array}\right] \quad \text { and } \tilde{\boldsymbol{L}}_{\boldsymbol{k}}=-\left[\begin{array}{cc} \boldsymbol{L}_{\boldsymbol{k}}-a_{0} \boldsymbol{I} & -\tilde{\nabla} \\ \tilde{\nabla}^{T} & 0 \end{array}\right]$$ $$\boldsymbol{u}_{k}^{s}=\left[\boldsymbol{u}_{\boldsymbol{k}}, p_{\boldsymbol{k}}\right]^{T}$$ $$\boldsymbol{f}_{k}^{s}=\left[f_{k}, 0\right]^{T}$$

最终，写成简洁的形式：

$$\boldsymbol{u}_{\boldsymbol{k}+}^{s}=\boldsymbol{H}_{\boldsymbol{k}+} \boldsymbol{f}_{\boldsymbol{k}+}^{s}$$ $$\begin{array}{r} \boldsymbol{u}_{\boldsymbol{k}+}^{s}=\left[\boldsymbol{u}_{k+n_{h} \kappa_{s}}^{s}, \quad \cdots \quad \boldsymbol{u}_{\boldsymbol{k}}^{s}, \quad \cdots \quad \boldsymbol{u}_{k-n_{h} \kappa_{s}}^{s}\right]^{T} \\ \text { and } \boldsymbol{f}_{\boldsymbol{k}+}^{s}=\left[\boldsymbol{f}_{\boldsymbol{k}+n_{h} \kappa_{s}}^{s}, \quad \cdots \quad \boldsymbol{f}_{\boldsymbol{k}}^{s}, \quad \cdots \quad \boldsymbol{f}_{\boldsymbol{k}-n_{h} \kappa_{s}}^{s}\right]^{T} \end{array}$$

SVD分解

对预解算子重新整理：

$$\left[\begin{array}{ll} W_{\mathbf{u}} & 0 \end{array}\right]\left[\begin{array}{c} \mathbf{u}_{\mathbf{k}} \\ p_{\mathbf{k}} \end{array}\right]=\left(\left[\begin{array}{ll} W_{\mathbf{u}} & 0 \end{array}\right] \tilde{\mathcal{H}}_{\mathbf{k}} W_{\mathbf{f}}^{-1}\right) W_{\mathbf{f}} \mathbf{f}_{\mathbf{k}}$$

f_k仅仅于速度有关系，因此也可以写成：

$$W_{\mathbf{u}} \mathbf{u}_{\mathbf{k}}=\tilde{\mathcal{H}}_{\mathbf{k}}^{S} W_{\mathbf{f}} \mathbf{f}_{\mathbf{k}}$$

svd分解：

$$\tilde{\mathcal{H}}_{\mathbf{k}}^{S}=\sum_{m} \psi_{\mathbf{k}, m} \sigma_{\mathbf{k}, m} \phi_{\mathbf{k}, m}^{*}$$

对应的force mode：

$$\mathbf{f}_{\mathbf{k}, m}=W_{\mathbf{f}}^{-1} \phi_{\mathbf{k}, m}$$

对应的response mode：

$$\mathbf{u}_{\mathbf{k}, m}=W_{\mathbf{u}}^{-1} \psi_{\mathbf{k}, m}$$

基函数满足正交性质：

$$\int \mathbf{f}_{\mathbf{k}, l}^{*} \mathbf{f}_{\mathbf{k}, m} d y=\delta_{l m}, \int \mathbf{u}_{\mathbf{k}, l}^{*} \mathbf{u}_{\mathbf{k}, m} d y=\delta_{l m}$$

响应的放大程度表现在奇异值的大小。

整体目标：把握信号的特点和机理的关系，构建模型

验证模型1:

2d，周期性变换，均匀流，环形

优点：定向分析，目前论文分析流肋条大小对减阻对效果

缺点：未考虑湍流的影响

验证模型2:

循环矩阵，应用到上述的模型1中，DNS（湍流模型）

优点：模型1的局部可以得到涡流的细节，利用peter schmid 循环矩阵的思路，再进行周期化，并做预解分析

验证模型3:

https://su2code.github.io/tutorials/Unsteady_NACA0012/

https://su2code.github.io/tutorials/Inc_Turbulent_NACA0012/

https://su2code.github.io/tutorials/Unsteady_Shape_Opt_NACA0012/

http://www.cemeai.icmc.usp.br/projetos/item/965-simulacoes-numericas-de-turbulencia-em-asas-de-aeronaves

可压的resolvent anlaysis N-S model

参考教材⁶

可压LES model （Mach > 0.2，Mach < 6）

基本方程推导

$$\frac{\partial \rho}{\partial t}+\frac{\partial \rho u_{j}}{\partial x_{j}}=0$$ $$\frac{\partial \rho u_{i}}{\partial t}+\frac{\partial \rho u_{i} u_{j}}{\partial x_{j}}+\frac{\partial p}{\partial x_{i}}=\frac{\partial \sigma_{i j}}{\partial x_{j}}$$ $$\frac{\partial \rho E}{\partial t}+\frac{\partial(\rho E+p) u_{j}}{\partial x_{j}}=\frac{\partial \sigma_{i j} u_{i}}{\partial x_{j}}-\frac{\partial}{\partial x_{j}} q_{j}$$ $$\rho E=\frac{p}{\gamma-1}+\frac{1}{2} \rho u_{i} u_{i}$$ $$p=\rho R T$$ $$R=C_{p}-C_{v}$$ $$\sigma_{i j}=2 \mu(T) S_{i j}-\frac{2}{3} \mu(T) \delta_{i j} S_{k k}$$ $$S_{i j}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)$$ $$\frac{\mu(T)}{\mu\left(T_{0}\right)}=\left(\frac{T}{T_{0}}\right)^{0.76}$$ $$q_{j}=-\kappa \frac{\partial T}{\partial x_{j}}$$ $$\kappa=\mu C_{p} / \operatorname{Pr}$$ $$\nu=\mu / \rho$$

推导

能量方程可以写成焓的表达形式：

$$\frac{\partial \rho h}{\partial t}+\frac{\partial \rho h u_{j}}{\partial x_{j}}=\frac{\partial p}{\partial t}+u_{j} \frac{\partial p}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial}{\partial x_{j}} T\right)+\Phi$$ $$h=C_{p} T$$

viscous dissipation

$$\Phi=\sigma_{i j} \frac{\partial u_{i}}{\partial x_{j}}$$

能量方程可以写成internal energy的表达形式：

$$\frac{\partial}{\partial t}\left(\rho C_{v} T\right)+\frac{\partial}{\partial x_{j}}\left(\rho u_{j} C_{v} T\right)=-p \frac{\partial u_{j}}{\partial x_{j}}+\Phi+\frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)$$

能量方程可以写成pressure的表达形式：

$$\frac{\partial p}{\partial t}+u_{j} \frac{\partial p}{\partial x_{j}}+\gamma p \frac{\partial u_{j}}{\partial x_{j}}=(\gamma-1) \Phi+(\gamma-1) \frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)$$ $$\gamma=C_{p} / C_{v}$$

能量方程可以写成entropy的表达形式：

$$\frac{\partial \rho s}{\partial t}+\frac{\partial \rho s u_{j}}{\partial x_{j}}=\frac{1}{T}\left[\Phi+\frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)\right]$$ $$s=C_{v} \ln \left(p \rho^{-\gamma}\right)$$

LES的滤波概念：尺度分离

假设湍流各向同性、均匀性。滤波器特性与坐标系无关。非各向同性可以利用非中心滤波器。滤波本质是卷积的思想。无非类似空间的格林函数。

$$\bar{\phi}(\mathbf{x}, t)=\frac{1}{\Delta} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G\left(\frac{\mathbf{x}-\boldsymbol{\xi}}{\Delta}, t-t^{\prime}\right) \phi\left(\xi, t^{\prime}\right) d t^{\prime} d^{3} \boldsymbol{\xi}$$

其中，cut-off scale in space and time，$\Delta$ ,$\tau_{c}$.

该方程可写成时域的简化形式：

$$\bar{\phi}=G \star \phi$$

转化为频域，

$$\overline{\hat{\phi}}(\mathbf{k}, \omega)=\hat{G}(\mathbf{k}, \omega) \hat{\phi}(\mathbf{k}, \omega)$$

或者：

$$\overline{\hat{\phi}}=\hat{G} \hat{\phi}$$

non-resolved part：

$$\phi^{\prime}(\mathbf{x}, t)=\phi(\mathbf{x}, t)-\bar{\phi}(\mathbf{x}, t)$$ $$\phi^{\prime}=(1-G) \star \phi$$

常用的三种滤波器类型：

$$\begin{array}{lll} \hline \text { Filter } & \text { Kernel in physical space } G & \text { Kernel in spectral space } \hat{G} \\ \hline \text { Box filter } & G(x-\xi)=\left\{\begin{array}{cc} \frac{1}{\Delta} & \text { if }|x-\xi| \leq \frac{\Delta}{2} \\ 0 & \text { otherwise } \end{array}\right. & \hat{G}(k)=\frac{\sin (k \Delta / 2)}{k \Delta / 2} \\ \text { Gaussian } & G(x-\xi)=\left(\frac{\varsigma}{\pi \Delta^{2}}\right)^{1 / 2} \exp \left(\frac{-\varsigma(x-\xi)^{2}}{\Delta^{2}}\right) & \hat{G}(k)=e^{-\left(\Delta^{2} k^{2}\right) / 4 \varsigma} \\ \text { Sharp cutoff } & G(x-\xi)=\frac{\sin \left(k_{c}(x-\xi)\right)}{k_{c}(x-\xi)} \text { with } k_{c}=\frac{\pi}{\Delta} & \hat{G}(k)=\left\{\begin{array}{cl} 1 & \text { if }|k|<k_{c} \\ 0 & \text { otherwise } \end{array}\right. \\ \hline \end{array}$$

利用taylor公式：

$$\phi(x-\Delta z, t)=\phi(x, t)+\sum_{l=1}^{\infty} \frac{(-1)^{l}(\Delta z)^{l}}{l !} \frac{\partial^{l} \phi(x, t)}{\partial x^{l}}$$ $$\begin{aligned} \bar{\phi}(x, t) &=\frac{1}{\Delta} \int_{x-\Delta \beta}^{x-\Delta \alpha} G\left(\frac{x-\xi}{\Delta}\right) \phi(\xi, t) d \xi \\ &=\int_{\alpha}^{\beta} G(z) \phi(x-\Delta z, t) d z \end{aligned}$$

代入公式后，得

$$\begin{aligned} \bar{\phi}(x, t) &=\int_{\alpha}^{\beta} G(z) \phi(x, t) d z+\int_{\alpha}^{\beta} \sum_{l=1}^{\infty} \frac{(-1)^{l}}{l !} \Delta^{l} z^{l} G(z) \frac{\partial^{l} \phi(x, t)}{\partial x^{l}} d z \\ &=\phi(x, t)+\sum_{l=1}^{\infty} \frac{(-1)^{l}}{l !} \Delta^{l} M_{l} \frac{\partial^{l} \phi(x, t)}{\partial x^{l}} \end{aligned}$$ $$M_{l}=\int_{\alpha}^{\beta} z^{l} G(z) d z$$

我们的目标当然是整理成矩阵的离散形式：

$$\bar{\phi}_{i}=G \phi_{i}=\sum_{n=-M}^{N} a_{n} \phi_{i+n}$$ $$\hat{G}(k)=\sum_{n=-M}^{N} a_{n} e^{j k n \Delta x}, \quad \text { with } j^{2}=-1$$

Commutation Error

$$B(a, b)=a b$$ $$\overline{\phi \psi}=\bar{\phi} \bar{\psi}+[G \star, B](\phi, \psi)$$

这里是建模的核心。

LES 方程

利用滤波器Favre Filtering

$$\overline{\rho \phi}=\bar{\rho} \tilde{\phi}$$ $$\phi=\tilde{\phi}+\phi^{\prime \prime}$$

分别分解为低频项和高频项

主控方程组：

【方程1】

$$\frac{\partial \rho}{\partial t}+\frac{\partial \rho u_{j}}{\partial x_{j}}=0$$

【方程1滤波形式】：方程整体滤波，卷积核函数操作，不可避免引入Commutation Error，通常密度最稳定，分离误差较小。

$$\frac{\partial \bar{\rho}}{\partial t}+\frac{\partial \bar{\rho} \tilde{u}_{j}}{\partial x_{j}}=0$$

【方程2-焓的形式】

$$\frac{\partial \rho h}{\partial t}+\frac{\partial \rho h u_{j}}{\partial x_{j}}=\frac{\partial p}{\partial t}+u_{j} \frac{\partial p}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial}{\partial x_{j}} T\right)+\Phi$$

viscous dissipation

$$\Phi=\sigma_{i j} \frac{\partial u_{i}}{\partial x_{j}}$$

【方程2-焓的滤波形式】：

$$\begin{array}{l} \frac{\partial \bar{\rho} \tilde{h}}{\partial t}+\frac{\partial \bar{\rho} \tilde{h} \tilde{u}_{j}}{\partial x_{j}}-\frac{\partial \bar{p}}{\partial t}-\tilde{u}_{j} \frac{\partial \bar{p}}{\partial x_{j}}+\frac{\partial \breve{q}_{j}}{\partial x_{j}}-\check{\Phi} \\ =-\left[\frac{\partial \bar{\rho} C_{p}\left(\widetilde{T u_{j}}-\tilde{T} \tilde{u}_{j}\right)}{\partial x_{j}}-\left(\overline{u_{j} \frac{\partial p}{\partial x_{j}}}-\tilde{u}_{j} \frac{\partial \bar{p}}{\partial x_{j}}\right)-(\bar{\Phi}-\check{\Phi})+\frac{\partial\left(\overline{q_{j}}-\check{q}_{j}\right)}{\partial x_{j}}\right] \end{array}$$ $$\check{\Phi}=\check{\sigma}_{i j} \frac{\partial \tilde{u}_{i}}{\partial x_{j}}$$ $$\check{\sigma}_{i j}=\mu(\tilde{T})\left(2 \tilde{S}_{i j}-\frac{2}{3} \delta_{i j} \tilde{S}_{k k}\right)$$ $$\tilde{S}_{i j}=\frac{1}{2}\left(\frac{\partial \tilde{u}_{i}}{\partial x_{j}}+\frac{\partial \tilde{u}_{j}}{\partial x_{i}}\right)$$ $$\check{q}_{j}=-\kappa(\tilde{T}) \frac{\partial \tilde{T}}{\partial x_{j}}$$

最终写成（此处漏掉了一些推导过程）：

$$\begin{array}{c} \frac{\partial \bar{\rho} \tilde{h}}{\partial t}+\frac{\partial \bar{\rho} \tilde{h} \tilde{u}_{j}}{\partial x_{j}}-\frac{\partial \bar{p}}{\partial t}-\tilde{u}_{j} \frac{\partial \bar{p}}{\partial x_{j}}+\frac{\partial \check{q}_{j}}{\partial x_{j}}-\check{\Phi} \\ =-\left[\frac{\partial C_{v} Q_{j}}{\partial x_{j}}+\Pi_{d i l}-\epsilon_{v}+\frac{\partial\left(\overline{q_{j}}-\check{q}_{j}\right)}{\partial x_{j}}\right] \end{array}$$ $$Q_{j}=\bar{\rho}\left(\widetilde{u_{j} T}-\tilde{u}_{j} \tilde{T}\right)$$ $$\Pi_{d i l}=\overline{p \frac{\partial u_{j}}{\partial x_{j}}}-\bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}$$ $$\epsilon_{v}=\bar{\Phi}-\check{\Phi}$$

【方程2-温度的形式】

$$\frac{\partial}{\partial t}\left(\rho C_{v} T\right)+\frac{\partial}{\partial x_{j}}\left(\rho u_{j} C_{v} T\right)=-p \frac{\partial u_{j}}{\partial x_{j}}+\Phi+\frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)$$

【方程2-温度的滤波形式】

$$\begin{array}{l} \frac{\partial \bar{\rho} C_{v} \tilde{T}}{\partial t}+\frac{\partial \bar{\rho} \tilde{u}_{j} C_{v} \tilde{T}}{\partial x_{j}}+\bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}-\check{\Phi}+\frac{\partial \check{q}_{j}}{\partial x_{j}} \\ =-\left[\frac{\partial \overline{\rho u_{j} C_{v} T}}{\partial x_{j}}-\frac{\partial \bar{\rho} \tilde{u}_{j} C_{v} \tilde{T}}{\partial x_{j}}+\overline{p \frac{\partial u_{j}}{\partial x_{j}}}-\bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}-(\bar{\Phi}-\check{\Phi})+\frac{\partial\left(\overline{q_{j}}-\check{q}_{j}\right)}{\partial x_{j}}\right] \end{array}$$

简化成：

$$\begin{array}{c} \frac{\partial \bar{\rho} C_{v} \tilde{T}}{\partial t}+\frac{\partial \bar{\rho} \tilde{u}_{j} C_{v} \tilde{T}}{\partial x_{j}}+\bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}-\check{\Phi}+\frac{\partial \check{q}_{j}}{\partial x_{j}} \\ =-\left[\frac{\partial C_{v} Q_{j}}{\partial x_{j}}+\Pi_{d i l}-\epsilon_{v}+\frac{\partial\left(\overline{q_{j}}-\check{q}_{j}\right)}{\partial x_{j}}\right] \end{array}$$

【方程2-压力的形式】

$$\frac{\partial p}{\partial t}+u_{j} \frac{\partial p}{\partial x_{j}}+\gamma p \frac{\partial u_{j}}{\partial x_{j}}=(\gamma-1) \Phi+(\gamma-1) \frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)\frac{\partial p}{\partial t}+u_{j} \frac{\partial p}{\partial x_{j}}+\gamma p \frac{\partial u_{j}}{\partial x_{j}}=(\gamma-1) \Phi+(\gamma-1) \frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)$$

【方程2-压力的滤波形式】

$$\begin{aligned} \frac{\partial \bar{p}}{\partial t}+& \tilde{u}_{j} \frac{\partial \bar{p}}{\partial x_{j}}+\gamma \bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}-(\gamma-1) \check{\Phi}+(\gamma-1) \frac{\partial \check{q}_{j}}{\partial x_{j}} \\ =&-\left[\overline{u_{j} \frac{\partial p}{\partial x_{j}}}-\tilde{u}_{j} \frac{\partial \bar{p}}{\partial x_{j}}+\gamma \overline{p \frac{\partial u_{j}}{\partial x_{j}}}-\gamma \bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}-(\gamma-1)(\bar{\Phi}-\check{\Phi})\right.\\ &\left.+(\gamma-1) \frac{\partial\left(\overline{q_{j}}-\check{q}_{j}\right)}{\partial x_{j}}\right] \end{aligned}$$

简化：

$$\begin{array}{c} \frac{\partial \bar{p}}{\partial t}+\tilde{u}_{j} \frac{\partial \bar{p}}{\partial x_{j}}+\gamma \bar{p} \frac{\partial \tilde{u}_{j}}{\partial x_{j}}-(\gamma-1) \check{\Phi}+(\gamma-1) \frac{\partial \check{q}_{j}}{\partial x_{j}} \\ =-R \frac{\partial Q_{j}}{\partial x_{j}}-(\gamma-1)\left[\Pi_{d i l}-\epsilon_{v}+\frac{\partial\left(\overline{q_{j}}-\check{q}_{j}\right)}{\partial x_{j}}\right] \end{array}$$

【方程2-熵的形式】

$$\frac{\partial \rho s}{\partial t}+\frac{\partial \rho s u_{j}}{\partial x_{j}}=\frac{1}{T}\left[\Phi+\frac{\partial}{\partial x_{j}}\left(\kappa \frac{\partial T}{\partial x_{j}}\right)\right]$$ $$s=C_{v} \ln \left(p \rho^{-\gamma}\right)$$

【方程2-熵的滤波形式】

$$\begin{aligned} \frac{\partial \bar{\rho} \tilde{s}}{\partial t} &+\frac{\partial \bar{\rho} \tilde{s} \tilde{u}_{j}}{\partial x_{j}}-\frac{1}{\tilde{T}}\left[\check{\Phi}+\frac{\partial}{\partial x_{j}}\left(\kappa(\tilde{T}) \frac{\partial \tilde{T}}{\partial x_{j}}\right)\right] \\ =&-\frac{\partial \bar{\rho}\left(\overline{s u_{j}}-\tilde{s} \tilde{u}_{j}\right)}{\partial x_{j}}+\left(\frac{\bar{\Phi}}{T}-\frac{\check{\Phi}}{\tilde{T}}\right)+\frac{1}{T} \frac{\partial}{\partial x_{j}}\left(\kappa(T) \frac{\partial T}{\partial x_{j}}\right) \\ &-\frac{1}{\tilde{T}} \frac{\partial}{\partial x_{j}}\left(\kappa(\tilde{T}) \frac{\partial \tilde{T}}{\partial x_{j}}\right) \end{aligned}$$ $$\tilde{s}=\frac{1}{\bar{\rho}} \overline{\rho C_{v} \ln \left(p \rho^{-\gamma}\right)}$$ $$\check{s}=C_{v} \ln \left(\bar{p} \bar{\rho}^{-\gamma}\right)$$

【方程3-状态方程】

$$\rho E=\frac{p}{\gamma-1}+\frac{1}{2} \rho u_{i} u_{i}$$

【方程3-状态方程的滤波形式】

$$\bar{\rho} \tilde{E}=\frac{\bar{p}}{\gamma-1}+\frac{1}{2} \bar{\rho} \widetilde{u_{i} u_{i}}$$

第二项没有办法直接计算，因此，需要引入模型。

【模型】：还有一些变参数模型，不具体介绍

$$\bar{\rho} \tilde{E}=\frac{\bar{p}}{\gamma-1}+\frac{1}{2} \bar{\rho} \tilde{u}_{i} \tilde{u}_{i}+\frac{\tau_{i i}}{2}$$

等价：

$$\bar{p}=(\gamma-1)\left[\bar{\rho} \tilde{E}-\frac{1}{2} \bar{\rho} \tilde{u}_{i} \tilde{u}_{i}-\frac{\tau_{i i}}{2}\right]$$

和

$$\tilde{T}=\frac{(\gamma-1)}{R}\left[\tilde{E}-\frac{1}{2} \tilde{u}_{i} \tilde{u}_{i}-\frac{\tau_{i i}}{2 \bar{\rho}}\right]$$

【方程4-Total Energy】

$$\frac{\partial \rho E}{\partial t}+\frac{\partial(\rho E+p) u_{j}}{\partial x_{j}}=\frac{\partial \sigma_{i j} u_{i}}{\partial x_{j}}-\frac{\partial}{\partial x_{j}} q_{j}$$

【方程4-Total Energy滤波形式】

$$\begin{array}{l} \frac{\partial \bar{\rho} \tilde{E}}{\partial t}+\frac{\partial(\bar{\rho} \tilde{E}+\bar{p}) \tilde{u}_{j}}{\partial x_{j}}-\frac{\partial \check{\sigma}_{i j} \tilde{u}_{i}}{\partial x_{j}}+\frac{\partial \breve{q}_{j}}{\partial x_{j}} \\ =-\frac{\partial}{\partial r}\left[\left(\overline{\rho u_{j} E}-\bar{\rho} \tilde{u}_{j} \tilde{E}\right)+\left(\overline{u_{j} p}-\tilde{u}_{j} \bar{p}\right)-\left(\overline{\sigma_{i j} u_{j}}-\check{\sigma}_{i j} \tilde{u}_{j}\right)-\left(\overline{q_{j}}-\check{q}_{j}\right)\right] \end{array}$$

简化：

$$\begin{array}{c} \frac{\partial \bar{\rho} \tilde{E}}{\partial t}+\frac{\partial(\bar{\rho} \tilde{E}+\bar{p}) \tilde{u}_{j}}{\partial x_{j}}-\frac{\partial \check{\sigma}_{i j} \tilde{u}_{i}}{\partial x_{j}}+\frac{\partial \check{q}_{j}}{\partial x_{j}} \\ =-\frac{\partial}{\partial x_{j}}\left[C_{p} Q_{j}+\mathcal{J}_{j}-\mathcal{D}_{j}-\left(\overline{q_{j}}-\check{q}_{j}\right)\right] \end{array}$$ $$\mathcal{D}_{j}=\bar{\sigma}_{i j u_{j}}-\check{\sigma}_{i j} \tilde{u}_{j}$$ $$\mathcal{J}_{j}=\frac{1}{2}\left(\bar{\rho} u_{j \bar{u}_{i}} u_{i}-\bar{\rho} \tilde{u}_{j} \widetilde{u_{i} u_{i}}\right)=\frac{1}{2}\left(\bar{\rho} u_{j \overline{u_{i}} u_{i}}-\bar{\rho} \tilde{u}_{j} \tilde{u}_{i} \tilde{u}_{i}-\tau_{i i}\right)$$

还有一些类似的模型，具体见📖。

【方程5-动量方程】

$$\frac{\partial \rho u_{i}}{\partial t}+\frac{\partial \rho u_{i} u_{j}}{\partial x_{j}}+\frac{\partial p}{\partial x_{i}}=\frac{\partial \sigma_{i j}}{\partial x_{j}}$$

【方程5-动量方程滤波形式】

$$\frac{\partial \bar{\rho} \tilde{u}_{i}}{\partial t}+\frac{\partial \bar{\rho} \tilde{u}_{i} \tilde{u}_{j}}{\partial x_{j}}+\frac{\partial \bar{p}}{\partial x_{i}}-\frac{\partial \check{\sigma}_{i j}}{\partial x_{j}}=-\frac{\partial \tau_{i j}}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left(\overline{\sigma_{i j}}-\check{\sigma}_{i j}\right)$$ $$\tau_{i j}=[G \star, H]\left(\rho, \rho u_{i}, \rho u_{j}\right)=\bar{\rho}\left(\widetilde{u_{i} u_{j}}-\tilde{u}_{i} \tilde{u}_{j}\right)$$

这里，模型还略去了压力的脉动量，也就是说，声学小扰动项不在les作为重点考虑。

Resolvent analysis of compressible N-S equation

可压N-S方程：

$$\left.\begin{array}{rl} \partial_{t} \rho+\mathbf{u} \cdot \nabla \rho+\rho \nabla \cdot \mathbf{u} & =0 \\ \rho \partial_{t} \mathbf{u}+\rho \mathbf{u} \cdot \nabla \mathbf{u}+\frac{1}{\gamma M a^{2}} \nabla p-\frac{1}{R e} \nabla \cdot \tau(\mathbf{u}) & =\mathbf{0} \\ \rho \partial_{t} T+\rho \mathbf{u} \cdot \nabla T+(\gamma-1) p \nabla \cdot \mathbf{u} & = \\ =\gamma(\gamma-1) \frac{M a^{2}}{R e} \tau(\mathbf{u}): \mathbf{d}(\mathbf{u})+\frac{\gamma}{\operatorname{Pr} R e} \Delta_{2} T & \\ \rho T-1-\gamma M^{2} p=0 \end{array}\right\}$$

strain and stress tensor：

$$\mathbf{d}(\mathbf{u})=\frac{1}{2}\left(\nabla \mathbf{u}+\nabla \mathbf{u}^{T}\right), \quad \tau(\mathbf{u})=\left[2 \mathbf{d}(\mathbf{u})-\frac{2}{3}(\nabla \cdot \mathbf{u}) \mathbf{I}\right]$$ $$\begin{aligned} \boldsymbol{u}(\boldsymbol{x}, t) &=\boldsymbol{U}_{b}(\boldsymbol{x})+\epsilon \boldsymbol{u}^{\prime}(\boldsymbol{x}, t) \\ p(\boldsymbol{x}, t) &=p_{b}(\boldsymbol{x})+\epsilon p^{\prime}(\boldsymbol{x}, t) \\ T(\boldsymbol{x}, t) &=T_{b}(\boldsymbol{x})+\epsilon T^{\prime}(\boldsymbol{x}, t) \\ \rho(\boldsymbol{x}, t) &=\rho_{b}(\boldsymbol{x})+\epsilon \rho^{\prime}(\boldsymbol{x}, t), \end{aligned}$$ $$\left[\boldsymbol{u}^{\prime}, p^{\prime}, T^{\prime}, \rho^{\prime}\right](\boldsymbol{x}, t)=\sum\left[\hat{\boldsymbol{u}}_{n}, \hat{p}_{n}, \hat{T}_{n}, \hat{\rho}_{n}\right](\boldsymbol{x}) \exp \left\{\lambda_{n} t\right\}+c . c .$$

steady version of compressible Navier- Stokes equations：

$$\begin{array}{c} \mathbf{U}_{b} \cdot \nabla \rho_{b}+\rho_{b} \nabla \cdot \mathbf{U}_{b}=0 \\ \rho_{b} \mathbf{U}_{b} \cdot \nabla \mathbf{U}_{b}+\nabla P_{b}-\frac{1}{R e} \nabla \cdot \tau\left(\mathbf{U}_{b}\right)=\mathbf{0} \\ \rho_{b} \mathbf{U}_{b} \cdot \nabla T_{b}+(\gamma-1) \rho_{b} T_{b} \nabla \cdot \mathbf{U}_{b}-\gamma(\gamma-1) \frac{M^{2}}{R e} \tau\left(\mathbf{U}_{b}\right) ：\\ \mathbf{d}\left(\mathbf{U}_{b}\right)-\frac{\gamma}{\operatorname{Pr} R e} \nabla^{2} T_{b}=0 \\ \rho_{b} T_{b}-1-\gamma M^{2} P_{b}=0 \end{array}$$

脉动量：

$$\begin{array}{c} \lambda_{n} \hat{\rho}_{n}+\mathbf{U}_{b} \cdot \nabla \hat{\rho}_{n}+\hat{\mathbf{u}}_{n} \cdot \nabla \rho_{b}+\rho_{b} \nabla \cdot \hat{\mathbf{u}}_{n}+\hat{\rho}_{n} \nabla \cdot \mathbf{U}_{b}=0, \\ \lambda_{n} \rho_{b} \hat{\mathbf{u}}_{n}+\hat{\rho}_{n} \mathbf{U}_{b} \cdot \nabla \mathbf{U}_{b}+\rho_{b} \hat{\mathbf{u}}_{n} \cdot \nabla \mathbf{U}_{b}+\rho_{b} \mathbf{U}_{b} \cdot \nabla \hat{\mathbf{u}}_{n} \\ +\nabla \hat{p}_{n}-\frac{1}{R e} \nabla \cdot \tau\left(\hat{\mathbf{u}}_{n}\right)=\mathbf{0} \\ \rho_{b} \hat{T}_{n} \lambda_{n}+\hat{\rho}_{n} \mathbf{U}_{b} \cdot \nabla T_{b}+\rho_{b} \hat{\mathbf{u}}_{n} \cdot \nabla T_{b}+\rho_{b} \mathbf{U}_{b} \cdot \nabla \hat{T}_{n} \\ +(\gamma-1)\left(\hat{\rho}_{n} T_{b} \nabla \cdot \mathbf{U}_{b}+\rho_{b} \hat{T}_{n} \nabla \cdot \mathbf{U}_{b}+\rho_{b} T_{b} \nabla \cdot \hat{\mathbf{u}}_{n}\right) \\ -\gamma(\gamma-1) \frac{M^{2}}{R e}\left[\tau\left(\hat{\mathbf{u}}_{n}\right): \mathbf{d}\left(\mathbf{U}_{b}\right)+\tau\left(\mathbf{U}_{b}\right): \mathbf{d}\left(\hat{\mathbf{u}}_{n}\right)\right] \\ -\frac{\gamma}{\operatorname{Pr} R e} \nabla^{2} \hat{T}_{n}=0 \\ \hat{\rho}_{n} T_{b}+\rho_{b} \hat{T}_{n}-\gamma M^{2} \hat{P}_{n}=0 \end{array}$$

上述方程是基于baseflow的。一些论文结果表明，在线性稳定性分析过程中，其结果不理想。

接下来，用真正的meanflow再分析一波：周期性满足下，对t的导数为0；

$$\bar{\alpha}=\frac{1}{\widetilde{T}} \int_{0}^{\tilde{T}} \alpha(t)$$ $$\begin{array}{c} \overline{\mathbf{u}} \cdot \nabla \bar{\rho}+\bar{\rho} \nabla \cdot \overline{\mathbf{u}}=-\overline{\mathbf{u}^{\prime} \cdot \nabla \rho^{\prime}}-\overline{\rho^{\prime} \nabla \cdot \mathbf{u}^{\prime}} \\ \overline{\rho \mathbf{u}} \cdot \nabla \overline{\mathbf{u}}+\nabla \bar{p}-\frac{1}{R e} \nabla \cdot \tau(\overline{\mathbf{u}})=-\overline{\rho^{\prime} \frac{\partial \mathbf{u}^{\prime}}{\partial t}}-\overline{\rho^{\prime} \overline{\mathbf{u}} \cdot \nabla \mathbf{u}^{\prime}}-\overline{\bar{\rho} \mathbf{u}^{\prime} \cdot \nabla \mathbf{u}^{\prime}}-\overline{\rho^{\prime} \mathbf{u}^{\prime} \cdot \nabla \overline{\mathbf{u}}}-\overline{\rho^{\prime} \mathbf{u}^{\prime} \cdot \nabla \mathbf{u}^{\prime}} \\ \overline{\rho \mathbf{u}} \cdot \nabla \bar{T}+(\gamma-1) \bar{\rho} \bar{T} \nabla \cdot \overline{\mathbf{u}}-\gamma(\gamma-1) \frac{M^{2}}{R e} \tau(\overline{\mathbf{u}}): \mathbf{d}(\overline{\mathbf{u}})-\frac{\gamma}{\operatorname{Pr} R e} \nabla^{2} \bar{T} \\ \bar{\rho} \bar{T}-1-\gamma M^{2} \bar{P}=-\overline{\rho^{\prime} T^{\prime}} \end{array}$$

代入

$$\mathbf{q}^{\prime}(\mathbf{x}, t)=\tilde{\mathbf{q}} \exp [\eta t]+\tilde{\mathbf{q}}^{*} \exp \left[\eta^{*} t\right]$$

可以写成：

$$\begin{array}{c} \overline{\mathbf{u}} \cdot \nabla \bar{\rho}+\bar{\rho} \nabla \cdot \overline{\mathbf{u}}=-2 A^{2}\left[\mathfrak{R}\left\{\tilde{\tilde{\mathbf{u}}}^{*} \cdot \nabla \tilde{\rho}\right\}+\mathfrak{R}\left\{\rho^{*} \nabla \cdot \tilde{\mathbf{u}}\right\}\right] \\ \overline{\rho \mathbf{u}} \cdot \nabla \overline{\mathbf{u}}+\nabla \bar{p}-\frac{1}{\operatorname{Re}} \nabla \cdot \tau(\overline{\mathbf{u}})=-2 A^{2}\left[\mathfrak{R}\left\{\tilde{\rho}^{*} i \eta_{i} \tilde{\mathbf{u}}\right\}+\mathfrak{R}\left\{\tilde{\rho}^{*} \tilde{\mathbf{u}} \cdot \nabla \overline{\mathbf{u}}\right\}+\mathfrak{R}\left\{\tilde{\rho}^{*} \overline{\mathbf{u}} \cdot \nabla \tilde{\mathbf{u}}\right\}+\mathfrak{R}\left\{\left.\bar{\rho} \tilde{\mathbf{u}}\right|^{*} \cdot \nabla \tilde{\mathbf{u}}\right\}\right] \\ \overline{\rho \mathbf{u}} \cdot \nabla \bar{T}+(\gamma-1) \bar{\rho} \bar{T} \nabla \cdot \overline{\mathbf{u}}-\gamma(\gamma-1) \frac{M^{2}}{R e} \tau(\overline{\mathbf{u}}): \mathbf{d}(\overline{\mathbf{u}})-\frac{\gamma}{\operatorname{Pr} R e} \nabla^{2} \bar{T} \\ =-2 A^{2}\left[\mathfrak{R}\{\bar{\rho} \tilde{\mathbf{u}} \cdot \nabla \tilde{T}\}+\mathfrak{R}\left\{\tilde{\rho}^{*} \overline{\mathbf{u}} \cdot \nabla \tilde{T}\right\}+\mathfrak{R}\left\{\tilde{\rho}^{*} \tilde{\mathbf{u}} \cdot \nabla \bar{T}\right\}\right]+-2 A^{2}(\gamma-1)\left[\mathfrak{R}\left\{\bar{\rho} \tilde{T}^{*} \nabla \cdot \tilde{\mathbf{u}}\right\}+\mathfrak{R}\left\{\tilde{\rho}^{*} \tilde{T} \nabla \cdot \overline{\mathbf{u}}\right\}\right. \\ \left.+\mathfrak{R}\left\{\tilde{\rho}^{*} \bar{T} \nabla \cdot \tilde{\mathbf{u}}\right\}\right]+-2 A^{2}\left[\mathfrak{R}\left\{\tilde{\rho}^{*} i \eta_{i} \tilde{T}\right\}-\gamma(\gamma-1) \frac{M^{2}}{R e} \mathfrak{R}\left\{\tau(\tilde{\mathbf{u}})^{*}: \mathbf{d}(\tilde{\mathbf{u}})\right\}\right] \\ \bar{\rho} \bar{T}-1-\gamma M^{2} \bar{P}=-2 A^{2} \mathfrak{K}\left\{\tilde{\rho}^{*} \tilde{T}\right\} \end{array}$$ $$A=\int_{\Omega}\left\langle\bar{\rho} u^{\prime 2}+\bar{\rho} T^{\prime 2}+\rho^{\prime 2}\right\rangle d \Omega$$

Weak form of incompressible N-S equation

$$\begin{array}{c} \partial_{t}[\mathbf{u}, p]=\mathcal{N} S([\mathbf{u}, p]) \equiv-\mathbf{u} \cdot \nabla \mathbf{u}-\nabla p+\frac{2}{\operatorname{Re}} \nabla \cdot \mathbf{D}(\mathbf{u}) \\ \nabla \cdot \mathbf{u}=0 \end{array}$$ $$\mathrm{D}(\mathbf{u})=1 / 2\left(\nabla \mathbf{u}+\nabla \mathbf{u}^{T}\right)$$

其弱解形式：

$$\forall[\mathbf{v}, q], \quad \partial_{t}\langle\mathbf{v}, \mathbf{u}\rangle=\langle\mathbf{v}, \mathcal{N} S(\mathbf{u}, p)\rangle+\langle q, \nabla \cdot \mathbf{u}\rangle$$

整个过程需要考虑边界条件：

$$\begin{array}{l} \forall[\mathbf{v}, q], \quad \partial_{t}\langle\mathbf{v}, \mathbf{u}\rangle=\langle\mathbf{v}, \mathcal{N} S(\mathbf{u}, p)\rangle+\langle q, \nabla \cdot \mathbf{u}\rangle \\ +\frac{1}{\varepsilon}\left(\langle\mathbf{v}, \mathbf{u}\rangle_{\Gamma_{\mathrm{cyl}}}+\left\langle\mathbf{v},\left(\mathbf{u}-\mathbf{e}_{x}\right)\right\rangle_{\Gamma_{\mathrm{in}}}+\left\langle v_{y}, u_{y}\right\rangle_{\Gamma_{\mathrm{axis}}}\right) \\ +\langle\mathbf{v}, \sigma \cdot \mathbf{n}\rangle_{\Gamma_{\mathrm{out}}} \cup \Gamma_{\mathrm{lat}} \cup \Gamma_{\mathrm{axis}} \end{array}$$

最终，整理出易于编写代码的形式：（引入分步积分）

$$\begin{aligned} \forall[\mathbf{v}, q], \partial_{t}\langle\mathbf{v}, \mathbf{u}\rangle=&-\langle\mathbf{v}, \mathcal{C}(\mathbf{u}, \mathbf{u}) / 2\rangle-2 \operatorname{Re}^{-1}\langle\mathbf{D}(\mathbf{v}): \mathbf{D}(\mathbf{u})\rangle \\ &+\langle\nabla \cdot \mathbf{v}, p\rangle+\langle q, \nabla \cdot \mathbf{u}\rangle \\ &+\frac{1}{\varepsilon}\left(\langle\mathbf{v}, \mathbf{u}\rangle_{\Gamma_{\mathrm{cyl}}}+\left\langle\mathbf{v},\left(\mathbf{u}-\mathbf{e}_{x}\right)\right\rangle_{\Gamma_{\mathrm{in}}}+\left\langle v_{y}, u_{y}\right\rangle_{\Gamma_{\mathrm{axis}}}\right) \end{aligned}$$

通过Newton Iteration求解base-flow

$$\mathcal{N} S\left(\mathbf{u}_{b}, p_{b}\right)=0$$

首先猜测初始值，然后不断的求梯度逼近，

$$\left[\mathbf{u}_{b}, p_{b}\right]=\left[\mathbf{u}_{b}^{g}, p_{b}^{g}\right]+\left[\delta \mathbf{u}_{b}, \delta p_{b}\right]$$

代入上面的式子，类似taylor分解操作：

$$\mathcal{N} S\left(\mathbf{u}_{b}^{g}, p_{b}^{g}\right)+\mathcal{L} N S_{\mathbf{u}_{s}^{g}}\left(\delta \mathbf{u}_{b}, \delta p_{b}\right)=0$$

同样整理成弱解形式：

$$\begin{array}{l} \left\langle\mathbf{v}, \mathcal{N} S\left(\mathbf{u}_{b}^{g}, p_{b}^{g}\right)\right\rangle+\left\langle q, \nabla \cdot \mathbf{u}_{b}^{g}\right\rangle+\left\langle\mathbf{v}, \mathcal{L} N S_{\mathbf{u}_{b}^{g}}\left(\delta \mathbf{u}_{b}, \delta p_{b}\right)\right\rangle \\ \quad+\left\langle q, \nabla \cdot \delta \mathbf{u}_{b}\right\rangle=0 \end{array}$$

其中，

$$\mathcal{L} N S_{\mathbf{U}}(\mathbf{u}, p)=-\mathcal{C}(\mathbf{U}, \mathbf{u})-\nabla p+\frac{2}{\operatorname{Re}} \nabla \cdot \mathbf{D}(\mathbf{u})$$ $$\mathcal{C}(\mathbf{U}, \mathbf{u})=(\mathbf{U} \cdot \nabla) \mathbf{u}+(\mathbf{u} \cdot \nabla) \mathbf{U}$$

最后整理出如下形式：

$$A \cdot \delta X=Y$$ $$\begin{aligned} \forall[\mathbf{v}, q], 0=&-\langle\mathbf{v}, \mathcal{C}(\mathbf{u}, \mathbf{u}) / 2\rangle-2 \operatorname{Re}^{-1}\langle\mathbf{D}(\mathbf{v}): \mathbf{D}(\mathbf{u})\rangle \\ &+\langle\nabla \cdot \mathbf{v}, p\rangle+\langle q, \nabla \cdot \mathbf{u}\rangle \\ &+\frac{1}{\varepsilon}\left(\langle\mathbf{v}, \mathbf{u}\rangle_{\Gamma_{\mathrm{cyl}}}+\left\langle\mathbf{v},\left(\mathbf{u}-\mathbf{e}_{x}\right)\right\rangle_{\Gamma_{\mathrm{in}}}+\left\langle v_{y}, u_{y}\right\rangle_{\Gamma_{\mathrm{axis}}}\right) \end{aligned}$$

Linear Stability

将参数分解为baseflow和脉动量的表达形式，导入N-S方程

$$\mathbf{u}(\mathbf{x}, t)=\mathbf{u}_{b}(\mathbf{x})+\varepsilon \hat{\mathbf{u}}(\mathbf{x}) e^{\lambda t}, \quad p(\mathbf{x}, t)=p_{b}(\mathbf{x})+\varepsilon \hat{p}(\mathbf{x}) e^{\lambda t}$$ $$\lambda \hat{\mathbf{u}}=\mathcal{L} \mathcal{N} \mathcal{S}_{\mathbf{u}_{b}}(\hat{\mathbf{u}}, \hat{p})$$

其弱解形式也可以写出：

$$\lambda\langle\mathbf{v}, \hat{\mathbf{u}}\rangle=\left\langle\mathbf{v}, \mathcal{L} \mathcal{N} \mathcal{S}_{\mathbf{u}_{b}}(\hat{\mathbf{u}}, \hat{p})\right\rangle+\langle q, \nabla \cdot \hat{\mathbf{u}}\rangle$$

最终整理出离散的矩阵形式：

$$\lambda \mathrm{B} \hat{X}=\mathrm{A} \hat{X}$$

Adjoint Eigenvalue Problem

$$\begin{array}{l} \mathbf{A} \mathbf{x}=\sigma \mathbf{B} \mathbf{x} \\ \mathbf{v} \mathbf{A}=\sigma \mathbf{v} \mathbf{B} \end{array}$$

Connection between adjoint and resolvent

Adjoint Eigenvalue Problem （resolvent anlaysis的再次不同理解）

前面已经提到线性N-S方程其弱解形式也可以写出：

$$\lambda\langle\mathbf{v}, \hat{\mathbf{u}}\rangle=\left\langle\mathbf{v}, \mathcal{L} \mathcal{N} \mathcal{S}_{\mathbf{u}_{b}}(\hat{\mathbf{u}}, \hat{p})\right\rangle+\langle q, \nabla \cdot \hat{\mathbf{u}}\rangle$$

Adjoint 形式：

$$\begin{array}{l} \forall(\mathbf{u}, p ; \mathbf{v}, q), \quad\left\langle\mathcal{L} \mathcal{N} \mathcal{S}_{\mathbf{U}}^{\dagger}(\mathbf{v}, q), \mathbf{u}\right\rangle+\langle\nabla \cdot \mathbf{v}, p\rangle \\ =\left\langle\mathbf{v}, \mathcal{L} \mathcal{N} \mathcal{S}_{\mathbf{U}}(\mathbf{u}, p)\right\rangle+\langle q, \nabla \cdot \mathbf{u}\rangle \end{array}$$

代入脉动量，最后得到频域表达式：

$$\forall(\mathbf{u}, p), \quad \lambda^{\dagger}\langle\hat{\mathbf{v}}, \mathbf{u}\rangle=\left\langle\mathcal{L} \mathcal{N} \mathcal{S}_{\mathbf{U}}^{\dagger}(\hat{\mathbf{v}}, \hat{q}), \mathbf{u}\right\rangle+\langle\nabla \cdot \hat{\mathbf{v}}, p\rangle$$

伴随特征值和特征值互为复共轭关系。

$$\bar{\lambda}^{\dagger} B \hat{X}^{\dagger}=A^{T} \hat{X}^{\dagger}$$

若

$$A=\left[\begin{array}{cc} 3+i & 5 \\ 2-2 i & i \end{array}\right]$$

则

$$A^{\dagger}=\left[\begin{array}{cc} 3-i & 2+2 i \\ 5 & -i \end{array}\right]$$ $$\mathbf{S}(\mathbf{x})=\frac{\hat{\mathbf{u}}^{\dagger} \otimes \hat{\mathbf{u}}}{\left\langle\hat{\mathbf{u}}^{\dagger}, \hat{\mathbf{u}}\right\rangle}$$ $$S_{w}(\mathbf{x})=\|\mathbf{S}(\mathbf{x})\|_{\infty} \equiv \frac{\left\|\hat{\mathbf{u}}^{\dagger}(\mathbf{x})\right\|\|\hat{\mathbf{u}}(\mathbf{x})\|}{\left\langle\hat{\mathbf{u}}^{\dagger}, \hat{\mathbf{u}}\right\rangle}$$

Iterative Methods for Eigenvalue Computations

$$X^{n}=\left(A-\lambda_{\text {shift }} B\right)^{-1} B X^{n-1}$$

Compressible N-S equation

$$\partial_{t} \mathcal{B}[\mathbf{u}, p, \rho, T]=\mathcal{N} \mathcal{S}[\mathbf{u}, p, \rho, T]$$

baseflow solution：

$$\mathcal{N} \mathcal{S}\left[\mathbf{u}_{b}, p_{b}, \rho_{b}, T_{b}\right]=0$$

eigenvalue：

$$\lambda \mathcal{B}[\hat{\mathbf{u}}, \hat{p}, \hat{\rho}, \hat{T}]^{T}=\mathcal{L} \mathcal{N} \mathcal{S}[\hat{\mathbf{u}}, \hat{p}, \hat{\rho}, \hat{T}]^{T}$$

Meanflow

在非线性域，baseflow引入的线性小扰动假设不再成立。因此，有必要重新定义新的”baseflow“，称为:meanflow

$$\mathbf{u}_{m}(\mathbf{x})=\frac{1}{T} \int_{0}^{T} \mathbf{u}(\mathbf{x}, t) d t$$ $$\mathbf{u}^{\prime}(\mathbf{x}, t)=\mathbf{u}(\mathbf{x}, t)-\mathbf{u}_{m}(\mathbf{x})$$

Ajoint ⁹

$$\dot{\mathbf{q}}=\mathbf{f}(\mathbf{q}, \mathbf{u}) \quad \text { PDE-constrained optimization }$$ $$\mathcal{J}=h(\mathbf{q}, \mathbf{u}) \rightarrow \text { opt/ cost fucntion }$$ $$\mathcal{L}=h(\mathbf{q}, \mathbf{u})-\left\langle\mathbf{q}^{+}, \dot{\mathbf{q}}-\mathbf{f}(\mathbf{q}, \mathbf{u})\right\rangle \rightarrow \text { opt }$$

q plus叫做ajoint或者Langrange multiplier.

优化问题 转化为求极值问题：

$$\frac{\delta \mathcal{L}}{\delta \mathbf{q}^{+}}=0 \quad \frac{\delta \mathcal{L}}{\delta \mathbf{q}}=0 \quad \frac{\delta \mathcal{L}}{\delta \mathbf{u}}=0$$

这三项分别可以写成：

$$\begin{aligned} &\frac{\delta \mathcal{L}}{\delta \mathbf{q}^{+}}=0 \rightarrow \dot{\mathbf{q}}-\mathbf{f}(\mathbf{q}, \mathbf{u})=0 \quad gov. equ.\\ &\frac{\delta \mathcal{L}}{\delta \mathbf{q}}=0 \rightarrow \dot{\mathbf{q}}^{+}+\left(\frac{\partial \mathbf{f}}{\partial \mathbf{q}}\right)^{H} \mathbf{q}^{+}+\left(\frac{\partial h}{\partial \mathbf{q}}\right)^{H}=0 \quad adj. equ. \\ &\frac{\delta \mathcal{L}}{\delta \mathbf{u}}=0 \rightarrow\left(\frac{\partial h}{\partial \mathbf{u}}\right)^{H}=-\left(\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\right)^{H} \mathbf{q}^{+} \quad opt. cond. (KKT) \end{aligned}$$

Ajoint优化问题的求解思路（梯度下降）：

第1个方程： u_guess, 然后求解；将q代入到第2方程分别求解q_plus；然后再代到第3个方程里面，即可得到gradient，或者叫做sentisitivity。作为优化问题，我们明确了迭代的方向，通过优化u，再重新循环迭代。直到，梯度为0。表示达到了极值条件。

在深度学习思路中，同样有类似的cost fucntion, 最终殊途同归，利用梯度下降法，求极值问题。

根据此思路，我们可以了解系统的stability, receptivity, sensitivity.

The role of ajoint

• 有自己独立的系统方程

how does a perturbution of the governing equations influence the output measure? （引入小扰动）

外部扰动（加扰动）

$$\mathcal{L}+\delta \mathcal{L}=\mathcal{J}+\delta \mathcal{J}-\left\langle\mathbf{q}^{+}, \dot{\mathbf{q}}-\mathbf{f}(\mathbf{q}, \mathbf{u})+\delta \xi\right\rangle$$

内部扰动（比如，将雷诺数上调一些，或其他参数或形状稍微变化一些）

$$\mathcal{L}+\delta \mathcal{L}=\mathcal{J}+\delta \mathcal{J}-\left\langle\mathbf{q}^{+}, \dot{\mathbf{q}}-\mathbf{f}(\mathbf{q}, \mathbf{u})-\delta \mathbf{f} \mathbf{q}\right\rangle$$

将其代入到上面的变分方程中，得到的类似的线性小扰动方程（线性化过程）：

外部扰动

$$\delta \mathcal{J} \sim\left\langle\mathbf{q}^{+}, \delta \xi\right\rangle$$

内部扰动

$$\delta \mathcal{J} \sim\left\langle\mathbf{q}^{+}, \delta \mathbf{f} \mathbf{q}\right\rangle$$

总结：表明包含 ajoint系统的sensitivity信息，可以用来判断：流场哪些区域，或者哪些参数对系统比较敏感，哪些对系统不敏感。

• 流动中最敏感区域是什么？
• 当我们改变平衡点、几何结构或参数时，稳定性/接受性等是如何变化的？
• 不稳定的“根源”在哪里？
• 我们应该在哪里放置控制元素来有效地操纵流？

比如临近失速的压气机系统，哪些位置的扰动会是系统立刻恶化，该位置用于布置传感器进行失速的检测。也可以对系统进行控制或者扩稳。

adjoint执行

$$\begin{aligned} &\frac{\delta \mathcal{L}}{\delta \mathbf{q}^{+}}=0 \rightarrow \dot{\mathbf{q}}-\mathbf{f}(\mathbf{q}, \mathbf{u})=0 \quad gov. equ.\\ &\frac{\delta \mathcal{L}}{\delta \mathbf{q}}=0 \rightarrow \dot{\mathbf{q}}^{+}+\left(\frac{\partial \mathbf{f}}{\partial \mathbf{q}}\right)^{H} \mathbf{q}^{+}+\left(\frac{\partial h}{\partial \mathbf{q}}\right)^{H}=0 \quad adj. equ. \\ &\frac{\delta \mathcal{L}}{\delta \mathbf{u}}=0 \rightarrow\left(\frac{\partial h}{\partial \mathbf{u}}\right)^{H}=-\left(\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\right)^{H} \mathbf{q}^{+} \quad opt. cond. (KKT) \end{aligned}$$

在执行中，我们必须要获得Jabobian matrix

$$\left(\frac{\partial \mathbf{f}}{\partial \mathbf{q}}\right)^{H}$$

对于incompressible，已经有很多这方面的工作。但对于compressible，大家都还在努力。

ajoint可以采用automatic differentiation, 但目前效率可能不高。

adjoint的好处：¹⁰

正向：100个参数需要求解100次10000维 整个流场（n-s equ）；逆向：2个参数需要求解2次10000维 整个流场（ajoint方程）；两者得到同样的结果。

一个论文里面经常出现的案例：圆柱扰流

SU2中的adjoint¹¹

可压RANS equation:

$$\left\{\begin{array}{ll} \mathcal{R}(U)=\frac{\partial U}{\partial t}+\nabla \cdot \vec{F}_{a l e}^{c}-\nabla \cdot\left(\mu_{t o t}^{1} \vec{F}^{v 1}+\mu_{t o t}^{2} \vec{F}^{v 2}\right)-\mathcal{Q}=0, & \text { in } \Omega, \quad t>0 \\ \vec{v}=\vec{u}_{\Omega}, & \text { on } S, \\ \partial_{n} T=0, & \text { on } S, \\ (W)_{+}=W_{\infty}, & \text { on } \Gamma_{\infty} \end{array}\right.$$ $$U=\left\{\begin{array}{c} \rho \\ \rho \vec{v} \\ \rho E \end{array}\right\}, \vec{F}_{a l e}^{c}=\left\{\begin{array}{c} \rho\left(\vec{v}-\vec{u}_{\Omega}\right) \\ \rho \vec{v} \otimes\left(\vec{v}-\vec{u}_{\Omega}\right)+\bar{I} p \\ \rho E\left(\vec{v}-\vec{u}_{\Omega}\right)+p \vec{v} \end{array}\right\}, \vec{F}^{v 1}=\left\{\begin{array}{c} \cdot \\ \overline{\bar{\tau}} \\ \overline{\bar{\tau}} \cdot \vec{v} \end{array}\right\}, \vec{F}^{v 2}=\left\{\begin{array}{c} q_{\rho} \\ \cdot \\ c_{p} \nabla T \end{array}\right\}, \mathcal{Q}=\left\{\begin{array}{l} \vec{q}_{\rho \vec{v}} \\ q_{\rho E} \end{array}\right\}$$ $$\begin{aligned} &\text { Minimize } \quad J(S)=\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} j(\vec{f}, \vec{n}) d s d t\\ &\text { such that } \quad \mathcal{R}(U)=0 \end{aligned}$$

写出Langrange multiplier

$$\mathcal{J}=\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} j(\vec{f}, \vec{n}) d s d t+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\Omega} \Psi^{\top} \mathcal{R}(U) d \Omega d t$$

变分，小扰动

$$\delta \mathcal{J}=\delta J+\frac{1}{\mathbb{T}} \int_{t_{0}}^{t_{f}} \int_{\Omega} \Psi^{\top} \delta \mathcal{R}(U) d \Omega d t$$ $$\delta J=\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\delta S} j(\vec{f}, \vec{n}) d s d t+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} \delta j(\vec{f}, \vec{n}) d s d t$$

第一项为原始状态下几何结构和标量函数值的变化，而第二项则依赖于原始几何和变形引起的标量函数的变化。

$$\begin{aligned} \int_{\delta S} j(\vec{f}, \vec{n}) d s &=\int_{S}\left(\frac{\partial j}{\partial \vec{f}} \partial_{n} \vec{f}-2 H_{m} j\right) \delta S d s \\ \int_{S} \delta j(\vec{f}, \vec{n}) d s &=\int_{S} \frac{\partial j}{\partial \vec{f}} \delta \vec{f}-\frac{\partial j}{\partial \vec{n}} \cdot \nabla_{S}(\delta S) d s \\ &=\int_{S} \frac{\partial j}{\partial \vec{f}} \cdot(\delta p \vec{n}-\delta \overline{\bar{\sigma}} \cdot \vec{n})-\left(\frac{\partial j}{\partial \vec{n}}+\frac{\partial j}{\partial \vec{f}} p-\frac{\partial j}{\partial \vec{f}} \cdot \overline{\bar{\sigma}}\right) \cdot \nabla_{S}(\delta S) d s \end{aligned}$$

阻力或升力优化问题：

$$j(\vec{f})=\vec{f} \cdot \vec{d}$$

最终，变分可以写出：

$$\delta \mathcal{J}=\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} \vec{d} \cdot(\delta p \vec{n}-\delta \overline{\bar{\sigma}} \cdot \vec{n}) d s d t+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\Omega} \Psi^{\top} \delta \mathcal{R}(U) d \Omega d t$$

第二项线性化：

$$\begin{aligned} \delta \mathcal{R}(U, \nabla U) &=\delta\left[\frac{\partial U}{\partial t}+\nabla \cdot \vec{F}_{a l e}^{c}-\nabla \cdot \mu_{t o t}^{k} \vec{F}^{v k}-\mathcal{Q}\right] \\ &=\delta\left[\frac{\partial U}{\partial t}+\nabla \cdot \vec{F}^{c}-\nabla \cdot\left(U \otimes \vec{u}_{\Omega}\right)-\nabla \cdot \mu_{t o t}^{k} \vec{F}^{v k}-\mathcal{Q}\right] \\ &=\frac{\partial}{\partial t}(\delta U)+\nabla \cdot\left(\frac{\partial \vec{F}^{c}}{\partial U} \delta U\right)-\nabla \cdot\left[\frac{\partial\left(U \otimes \vec{u}_{\Omega}\right)}{\partial U} \delta U\right]-\nabla \cdot \mu_{t o t}^{k}\left[\frac{\partial \vec{F}^{v k}}{\partial U} \delta U+\frac{\partial \vec{F}^{v k}}{\partial(\nabla U)} \delta(\nabla U)\right]-\frac{\partial \mathcal{Q}}{\partial U} \delta U \\ &=\frac{\partial}{\partial t}(\delta U)+\nabla \cdot\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}-\mu_{t o t}^{k} \vec{A}^{v k}\right) \delta U-\nabla \cdot \mu_{t o t}^{k} \bar{D}^{v k} \delta(\nabla U)-\frac{\partial \mathcal{Q}}{\partial U} \delta U \end{aligned}$$ $$\left.\begin{array}{ll} \vec{A}^{c}=\left(A_{x}^{c}, A_{y}^{c}, A_{z}^{c}\right), & A_{i}^{c}=\left.\frac{\partial \vec{F}_{i}^{c}}{\partial U}\right|_{U(x, y, z)} \\ \vec{A}^{v k}=\left(A_{x}^{v k}, A_{y}^{v k}, A_{z}^{v k}\right), & A_{i}^{v k}=\left.\frac{\partial \vec{F}_{i}^{v k}}{\partial U}\right|_{U(x, y, z)} \\ \overline{\bar{D}}^{v k}=\left(\begin{array}{ccc} D_{x x}^{v k} & D_{x y}^{v k} & D_{x z}^{v k} \\ D_{y x}^{v k} & D_{y y}^{v k} & D_{y z}^{v k} \\ D_{z x}^{v k} & D_{z y}^{v k} & D_{z z}^{v k} \end{array}\right), \quad D_{i j}^{v k}=\left.\frac{\partial \vec{F}_{i}^{v k}}{\partial\left(\partial_{j} U\right)}\right|_{U(x, y, z)} \end{array}\right\} \quad i, j=1 \ldots 3, \quad k=1,2,$$

linearized Navier-Stokes equations,

$$\left\{\begin{array}{ll} \delta \mathcal{R}(U)=\frac{\partial}{\partial t}(\delta U)+\nabla \cdot\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}-\mu_{t o t}^{k} \vec{A}^{v k}\right) \delta U-\nabla \cdot \mu_{t o t}^{k} \overline{\bar{D}}^{v k} \delta(\nabla U)-\frac{\partial \mathcal{Q}}{\partial U} \delta U=0 & \text { in } \Omega, \quad t>0 \\ \delta \vec{v}=-\partial_{n}\left(\vec{v}-\vec{u}_{\Omega}\right) \delta S & \text { on } S, \\ \partial_{n}(\delta T)=\nabla T \cdot \nabla_{S}(\delta S)-\partial_{n}^{2}(T) \delta S & \text { on } S, \\ (\delta W)_{+}=0 & \text { on } \Gamma_{\infty} \end{array}\right.$$

代回到原来的式子

$$\begin{aligned} \delta \mathcal{J}=\delta J+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\Omega} \Psi^{\top} \frac{\partial}{\partial t}(\delta U) d \Omega d t+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\Omega} \Psi^{\top} \nabla \cdot\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}-\mu_{t o t}^{k} \vec{A}^{v k}\right) \delta U d \Omega d t \\ -\frac{1}{\mathbb{T}} \int_{+}^{t_{f}} \int_{\Omega} \Psi^{\top} \nabla \cdot \mu_{t o t}^{k} \overline{\bar{D}}^{v k} \delta(\nabla U) d \Omega d t-\frac{1}{\mathbb{T}} \int_{+}^{t_{f}} \int_{\Omega} \Psi^{\top} \frac{\partial \mathcal{Q}}{\partial U} \delta U d \Omega d t \end{aligned}$$

利用分步积分，整理

$$\begin{aligned} \delta \mathcal{J}=\delta J &+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\partial \Omega}\left(B_{1}-B_{2}+B_{3}\right) d s d t \\ &+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\Omega}\left[-\frac{\partial \Psi^{\top}}{\partial t}-\nabla \Psi^{\top} \cdot\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}-\mu_{t o t}^{k} \vec{A}^{v k}\right)-\nabla \cdot\left(\nabla \Psi^{\top} \cdot \mu_{t o t}^{k} \overline{\bar{D}}^{v k}\right)-\Psi^{\top} \frac{\partial \mathcal{Q}}{\partial U}\right] \delta U d \Omega d t \end{aligned}$$ $$\begin{array}{l} B_{1}=\Psi^{\top}\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}\right) \delta U \cdot \vec{n} \\ B_{2}=\Psi^{\top} \mu_{t o t}^{k} \vec{A}^{v k} \delta U \cdot \vec{n}+\Psi^{\top} \mu_{t o t}^{k} \overline{\bar{D}}^{v k} \cdot \nabla(\delta U) \cdot \vec{n} \\ B_{3}=\nabla \Psi^{\top} \cdot \mu_{t o t}^{k} \overline{\bar{D}}^{v k} \delta U \cdot \vec{n} \end{array}$$

如果不考虑边界条件的影响

最终，整理得到

$$\begin{aligned} \delta \mathcal{J} &=\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} \vec{d} \cdot(\delta p \vec{n}-\delta \overline{\bar{\sigma}} \cdot \vec{n}) d s d t-\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S}\left(\vec{\varphi}+\psi_{\rho E} \vec{v}\right) \cdot(\delta p \vec{n}-\delta \overline{\bar{\sigma}} \cdot \vec{n}) d s d t \\ &-\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} \mu_{t o t}^{2} c_{p} \partial_{n}\left(\psi_{\rho E}\right) \delta T d s d t-\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{S} \frac{\partial \mathcal{J}}{\partial S} \delta S d s d t \\ &+\frac{1}{\mathbb{T}} \int_{t_{o}}^{t_{f}} \int_{\Omega}\left[-\frac{\partial \Psi^{\top}}{\partial t}-\nabla \Psi^{\top} \cdot\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}-\mu_{t o t}^{k} \vec{A}^{v k}\right)-\nabla \cdot\left(\nabla \Psi^{\top} \cdot \mu_{t o t}^{k} \overline{\bar{D}}^{v k}\right)-\Psi^{\top} \frac{\partial \mathcal{Q}}{\partial U}\right] \delta U d \Omega d t \end{aligned}$$

至此，adjoint N-S equation可以得到：

$$\left\{\begin{array}{ll} -\frac{\partial \Psi^{\top}}{\partial t}-\nabla \Psi^{\top} \cdot\left(\vec{A}^{c}-\overline{\bar{I}} \vec{u}_{\Omega}-\mu_{t o t}^{k} \vec{A}^{v k}\right)-\nabla \cdot\left(\nabla \Psi^{\top} \cdot \mu_{t o t}^{k} \overline{\bar{D}}^{v k}\right)-\Psi^{\top} \frac{\partial \mathcal{Q}}{\partial U}=0, & \text { in } \Omega, \quad t>0 \\ \vec{\varphi}=\vec{d}-\psi_{\rho E} \vec{v}, & \text { on } S, \\ \partial_{n}\left(\psi_{\rho E}\right)=0, & \text { on } S . \end{array}\right.$$

⁰. Analysis of fluid systems: stability, receptivity, sensitivity ↩

¹. On the structure and origin of pressureuctuations in wall turbulence: predictionsbased on the resolvent analysis ↩

². A framework for studying the eect of compliant surfaces on wall turbulence ↩

³. Resolvent Analysis for Turbulent Channel Flow with Riblets ↩

⁴. On the existence of steady-state resonant waves in experiments ↩

⁵. Elements of resolvent methods in uid mechanics: notes for an introductory short course v0.3 ↩

⁶. Large Eddy Simulation for Compressible Flows ↩

⁷. phd thesis: Fluid dynamic instabilities in complexflow systems ↩

⁸. Computation of the blu -body sound generation by a self-consistent mean flow formulation ↩

⁹. https://www.youtube.com/watch?v=JFPTQ-ALDcM ↩

¹⁰. http://www2.eng.cam.ac.uk/~mpj1001/MJ_ABCD.html ↩

¹¹. Viscous Continuous Adjoint Approach for the Design of Rotating Engineering Applications ↩