在2025年之后，每一次的LCO课程将会从CFD算法文章库中进行筛选并在课堂上现场介绍并植入。

CFD算法文章库

如何保证高时间精度同时加快计算效率？

该相关算法爆了大量JCP[AF20, ATF23, KS23, KSS21, LM91]。或者说该方面研究已经可以单分一块。不可压NS方程求解最耗费时间的是压力方程。常规有限体积法求解充其量获得的是空间二阶、时间一阶的瞬态结果。龙格库塔（RK）时间步推进可以获得关于时间的高精度结果。大名鼎鼎的Moin以及Le看起来最早做了相关的尝试并发表在1991年的JCP上[LM91]，全文非常短，引用文献只有6个。单多步的RK推进需要求解更多的压力方程导致效率低下。一些文章力图实现的就是在多步RK推进下，仅仅推进一次压力方程来增加计算速度[AF20]。在OpenFOAM中植入RK算法并提升计算效率，看起来非常简单并且可以用于各种不同的工况验证（多篇sci就这样水出来了）。

_images/rk.JPG

在练手阶段，植入瞬态投影法、稳态SIMPLE算法、普适性密度基可压缩算法是比较合适的。通过这些经典算法，你会透彻的了解NS方程速度压力耦合求解的思想。侧重算法的理解与编程。此部分内容是编程课的热身内容。3个求解器都是先从算法讲解开始，然后手搓求解器。另外，下图中红圈处会发现震荡。课堂上需要通过算法处理这个问题。

_images/density.jpg

什么样的数值格式是动能守恒的？

考虑\(\frac{\p\frac{\phi^2}{2}}{\p t}\)，有：

\[ \frac{\p\frac{\phi^2}{2}}{\p t}=\phi\frac{\p\rho\phi}{\p t}+\frac{\phi^2}{2}\frac{\p\rho}{\p t} \]

将\(\frac{\p\rho\phi}{\p t}\)写为\(-\mathcal{C}\)，如果\(\frac{\p\rho\phi}{\p t}=-\mathcal{C}=-\nabla\cdot(\rho\phi\bfU)\)是成立的，同时定义中间变量\(\mathcal{M}\)：\(\frac{\p\rho}{\p t}=\mathcal{M}=-\nabla\cdot(\rho\bfU)\)，上式可以写为：

\[ \frac{\p\frac{\phi^2}{2}}{\p t}=- \left(\phi\mathcal{C} - \frac{\phi^2}{2}\mathcal{M}\right) \]

也即

\[ \frac{\p\frac{\phi^2}{2}}{\p t}=-\phi\nabla\cdot(\rho\phi\bfU) - \frac{\phi^2}{2} \nabla\cdot(\rho\bfU)=\nabla\cdot\left(\rho\frac{\phi^2}{2}\bfU\right) \]

这表明了如果\(\mathcal{C}\)和\(\mathcal{M}\)具有满足散度形式的方程，那么对于二次变量\(\phi^2\)也存在一个散度形式的守恒方程。这通常意味着整个系统是能量守恒的（\(\bfU^2\)能量相关）。相关的恒等关系在Coppola等发表在JCP的文章中被认为是二次变量恒等式[CCPdL19]。对于NS方程系统，所有的对流项都可以被认为是满足二次变量恒等（能量守恒）。但粘性项以及压力梯度项会导致能量的耗散。问题是，二次变量恒等在推导过程中，使用了散度的链条法则，这对于连续形式是没有问题的，但是在离散过程中并不适用。因此，对流项离散后往往不能保证二次变量恒等（能量守恒）。同样的，对于\(\mathcal{C}\)同样也具有不同的形式，例如其可以写为\(\nabla\cdot(\rho\phi\bfU)\)，也可以写成\(\phi\nabla\cdot(\rho\bfU)+\rho\bfU\cdot\nabla\phi\)。这里面牵涉的概念是使用什么样的变量形式去进行离散求解。Coppola等的文章中的公式7-11将其写成5种形式[CCPdL19]，并推导出满足能量守恒的系数法则。在后人的工作中被称之为动能守恒类方法（KEP）。Michele等在KEP方法基础之上进行改进，新的方法在2023年发表在JCP上[DMC23]。Ranocha等人在2022年发表文章，标题很有意思：如果你使用基于熵的高阶格式，即使压力不震荡，也不能解决局部不稳定问题[RG21]，也是基于该方法，对该方法进行植入并测试看起来是有意思的。

稀疏气固跨音速算法？

欧拉方程比较有意思，是因为欧拉方程存在间断。多维的欧拉方程求解存在激波，结果看起来也非常酷炫。在单相欧拉方程的基础上继续考虑颗粒，如果假定颗粒是稀疏的，变构成了跨音速的稀疏气固流动系统。在使用这一套方法用于经典激波管案例时候，变量的型线与单相欧拉方程存在区别[SMT03]：

_images/saito2003.JPG

进一步的考虑方程系统，有：

\[ \frac{\p \rho}{\p t}+\nabla\cdot(\rho\bfU_g)=0 \]
\[ \frac{\p \rho\bfU_g}{\p t}+\nabla\cdot(\rho\bfU_g\bfU_g + p)=-\frac{\sigma}{m}D \]
\[ \frac{\p \rho E_g}{\p t}+\nabla\cdot( (\rho E_g + p)\bfU_g)=-\frac{\sigma}{m} (Q+\bfU_d D) \]
\[ \frac{\p \sigma}{\p t}+\nabla\cdot(\sigma\bfU_d)=0 \]
\[ \frac{\p \sigma\bfU_d}{\p t}+\nabla\cdot(\rho\bfU_d\bfU_d)=\frac{\sigma}{m}D \]
\[ \frac{\p \rho E_d}{\p t}+\nabla\cdot( \rho E_d \bfU_g)=\frac{\sigma}{m} (Q+\bfU_d D) \]

由于固相颗粒是稀疏的因此并不存在固相的压力。这个看起来很有意思因为稀疏气固激波管的流动型线看起来很特别，值得植入研究一下。

尝试写一个高精度、相场法的多相流求解器？

多相流的计算方法可以从NS方程进行简单的拓展。最简单的多相流计算方法就是直接模拟，在CFD领域，单相流的直接模拟也简单，因为不需要附加湍流模型。简单的模型的缺点就是比较贵。OpenFOAM目前显存了VOF方法。但植入一个新的方法看起来是有意思的，也比较适合练手。因此通过本算例植入的代码，主要研究1）相场法的思想、2）FVM的高精度格式、3）解的有界性、4）upwind与downwind等。如下图，采用不同的方法来计算zalesak disk的情况下，不同的格式会预测不同的结果。即使使用TVD格式也会产生数值耗散。对于普通的标量数值耗散并没有太大的问题。但某些情况下我们需要非常明晰的界面。这就需要不同的数值方法来处理。算法的难点在于：1）使用高精度算法如何保证有界性？2）如何保证界面的尖锐？3）不会产生一些奇怪的结果（如下图第三个，解的分布有点奇怪） 。

随后，可以尝试句许植入一个两相流耦合求解器，算法主要考虑体积力对速度压力耦合的影响。CFD中的体积力有3种植入方式：1）隐性、2）显性、3）压力方程。可以对比不同的方法预测的结果。下面的结果就是计算出来的多相流不稳定性。经验来看，在入门CFD的过程中，很多人植入的算法都会发生数值震荡。

_images/aza.png

下面可以进一步的进行拓展。很明显上图的结果是比较耗散的。下一步可以尝试把两相流算法和高精度相场方法相结合，制定一个新的求解器。新的求解器求解后的结果如下图。可以很明显的看出求解的结果非常尖锐！在这个基础上，如果附加上表面张力，既可以形成一个全新的VOF求解器。

_images/high-order.JPG

隐藏流体力学？

数据驱动CFD里面独树一帜的是PINN。通过PINN来练手看起来是个很好的选择，并且可以用来测试全连接网络、残差神经网络植入、以及新的激活函数等。使用PINN求解cavity flow的代码已经开源。不得不承认的是，PINN作为一个新兴的计算方法，用来计算流场实在是太慢了。因此使用PINN来计算流场，不像是走正常的套路。因此应该尝试将PINN用于正道，也即使用PINN来处理病态问题，用其处理经典CFD方法不能处理的问题，这个应该才是PINN真正的用处。这种应用也即之前刊发在Science正刊上的隐藏流体力学。否则用PINN来硬钢经典FVM，PINN必输无疑，太慢了，基本用不上。本案例的设计理念，就是把PINN用在正道上。在下图中，求解的ODE是未知的，但是PINN成功的发现了未知的ODE，并且数据与解析解完美吻合。

_images/hdf.png

植入一个全新的湍流模型

从NASA收录的湍流模型来看，NASA认为最近10年能打的只有2个湍流模型，一个是2018年提出的代数应力模型k-kl，一个是2015年的Wray-Agarwal。2018年的代数应力模型引用少的可怜，应该是NASA自家的模型，因此被自家收录。因此，最近10年，能打的应该就是Wray-Agarwal湍流模型。把Wray-Agarwal湍流模型进行植入看起来挺有意思。相关的方程如下：

\[ S=\sqrt{2}|\bfS|, \bfS=\frac{\nabla\bfU+\nabla\bfU^T}{2}, W=\sqrt{2}|{\bf{W}}|, {\bf{W}}=\frac{\nabla\bfU-\nabla\bfU^T}{2} \]
\[ \eta=S \max\left(1, \left|\frac{W}{S}\right|\right) ,k=\frac{\nu_t S}{\sqrt{C_\mu}},\omega=\frac{S}{\sqrt{C_\mu}}, \]
\[ \mathrm{arg_1}=\frac{v+R}{2}\frac{\eta^2}{C_\mu k\omega},f_1=\mathrm{tanh}(\mathrm{arg_1}^4) \]
\[ C_1=f_1(C_{1kw} - C_{1k\varepsilon}) + C_{1k\varepsilon}, \sigma_R=f_1(\sigma_{kw}-\sigma_{k\varepsilon}) + \sigma_{k\varepsilon} \]
\[ C_{2kw}=\frac{C_{1kw}}{\kappa^2}+\sigma_{kw}, C_{2k\varepsilon}=\frac{C_{1k\varepsilon}}{\kappa^2}+\sigma_{k\varepsilon} \]
\[\begin{split} \begin{split} \frac{\p R}{\p t}+\nabla\cdot(\bfU R) &-\nabla\cdot((\sigma_R R+\nu)\nabla R)= \\\\ & C_1 R S + f_1 C_{2kw}\frac{R}{S}(\nabla R)\cdot(\nabla S)-(1-f_1)\min\left(C_{2ke}R^2 \frac{|\nabla S|^2}{S^2}, C_m|\nabla R|^2\right) \end{split} \end{split}\]
\[ \chi=R/v,f_\mu=\frac{\chi^3}{\chi^2+C_w^3},\nu_t= f_\mu R, \]
\[ \sigma_{kw}=0.72,\sigma_{k\varepsilon}=1.0,C_{1kw}=0.0829, C_{1k\varepsilon}=0.1284 \]
\[ \kappa=0.41,C_w=8.54,C_\mu=0.09, C_m=8.0 \]

植入数据驱动湍流模型

数据驱动湍流模型是目前国际的研究顶流。所有的课题组都对这个蠢蠢欲动。然而其实这是个非常简单的研究领域。如果说在OpenFOAM里面植入MULES算法可能需要一个博士，但植入数据驱动LES模型一个硕士生完全没问题。比如咱们进行的这个算例，采用经典神经网络架构，共训练2000个批次的数据，由于训练量比较小，使用CPU就可以进行训练。但是结果相当好。在下图中，数据驱动LES与经典CFD-LES预测的数据云图一致（局部差异基本可忽略）。在观测时间平均变量时候，结果非常完美。对于比较难预测的雷诺应力相关量，数据驱动LES预测的结果与经典CFD-LES吻合度非常高。类似的文章在2023年，2024年都刊发在了JFM。本算例的结果足以支撑一篇sci的数据。该算例部分代码开源在ML: 代码实例-数据驱动LES，该算例是以后LCO课程的教学案例。

_images/ddles.JPG
_images/data_driven_LES.svg

数据驱动前沿热点

经典CFD经历了60多年的发展后，同位网格SIMPLE类算法在不可压缩领域统一江湖，可压缩领域通通都是密度基。数据驱动CFD目前的发展状况类似CFD在1970年代的情况，群魔乱舞百家齐放。目前还没有一个大统一的寡头算法出现。目前的态势是每家都有自己的特色，那么在10年后20年后数据驱动CFD，哪个算法会统一江湖目前还不确定。目前这一块研究的特色是杂！尤其在结合各种机器学习的网络后，更杂！但可以简单的进行归类。在《无痛苦NS方程笔记》里面，机器学习与CFD相应的结合点主要分为

  • 流场重组（卷积网络、Unet、解码器、超分辨率、扩散模型与GAN）

  • 求解加速（全连接网络、卷积网络、残差网络等）

  • 数据驱动湍流模型（全连接网络、卷积网络、极限学习机等）

  • 数据驱动多相流

  • 循环神经网络与瞬态流场（循环神经网络、LSTM）

  • 正向逆向PINN、隐藏流体力学

  • POD、DMD、流场降阶预测方法（不需要神经网络）

  • 强化学习与过程控制

上述内容全部在《无痛苦NS方程笔记》里面有介绍，课堂上也会进一步的进行介绍。

参考文献

[AF20] (1,2)

Abhiram B Aithal and Antonino Ferrante. A fast pressure-correction method for incompressible flows over curved walls. Journal of Computational Physics, 421:109693, 2020.

[ATF23]

Abhiram B Aithal, Mira Tipirneni, and Antonino Ferrante. Temporal accuracy of fastrk3. Journal of Computational Physics, 475:111853, 2023.

[CCPdL19] (1,2)

Gennaro Coppola, Francesco Capuano, Sergio Pirozzoli, and Luigi de Luca. Numerically stable formulations of convective terms for turbulent compressible flows. Journal of Computational Physics, 382:86–104, 2019.

[DMC23]

Carlo De Michele and Gennaro Coppola. Asymptotically entropy-conservative and kinetic-energy preserving numerical fluxes for compressible euler equations. Journal of Computational Physics, 492:112439, 2023.

[KS23]

Mokbel Karam and Tony Saad. On the theory of fast projection methods for high-order navier-stokes solvers. Journal of Computational Physics, 495:112557, 2023.

[KSS21]

Mokbel Karam, James C Sutherland, and Tony Saad. Low-cost runge-kutta integrators for incompressible flow simulations. Journal of Computational Physics, 443:110518, 2021.

[LM91] (1,2)

Hung Le and Parviz Moin. An improvement of fractional step methods for the incompressible navier-stokes equations. Journal of computational physics, 92(2):369–379, 1991.

[RG21]

Hendrik Ranocha and Gregor J Gassner. Preventing pressure oscillations does not fix local linear stability issues of entropy-based split-form high-order schemes. Communications on Applied Mathematics and Computation, pages 1–24, 2021.

[SMT03]

Tsutomu Saito, M Marumoto, and K Takayama. Numerical investigations of shock waves in gas-particle mixtures: evaluation of numerical methods for dusty-gas shock wave phenomena. Shock Waves, 13(4):299–322, 2003.