Efficient monolithic immersed boundary projection method for incompressible... 阅读笔记

Efficient monolithic immersed boundary projection method for incompressible flows with heat transfer

https://doi.org/10.1016/j.jcp.2023.111929

时间格式相关的处理

这处理是什么意思完全看不懂

就是说前面一项是都是 n+1 后面一项都是 n

然后现在改成 n+1 和 n 组合，n 和 n+1 组合

这就是所谓的交错时间？

得到他这个矩阵形式的控制方程

NS 方程

$$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}+\frac{1}{2}({\mathbf{u}}^{n+1}\cdot \mathcal{G}{\mathbf{u}}^{n+1}+{\mathbf{u}}^n\cdot \mathcal{G}{\mathbf{u}}^n)=-\mathcal{G}{p}^{n+1/2}+\frac{1}{2\mathrm{Re}}(\mathcal{L}{\mathbf{u}}^{n+1}+\mathcal{L}{\mathbf{u}}^n)+{\mathcal{H}}_m{\mathbf{F}}^n+\mathcal{F}({\theta }^{n+1/2})+{\mathbf{mbc}}^{n+1/2}$$

代入他那个线性化

$$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}+\frac{1}{2}{\mathcal{N}}_m({\mathbf{u}}^{n+1})=-\mathcal{G}{p}^{n+1/2}+\frac{1}{2\mathrm{Re}}(\mathcal{L}{\mathbf{u}}^{n+1}+\mathcal{L}{\mathbf{u}}^n)+{\mathcal{H}}_m{\mathbf{F}}^n+\mathcal{F}({\theta }^{n+1/2})+{\mathbf{mbc}}^{n+1/2}$$

移项，验证一下左端项……

$$-\mathcal{F}({\theta }^{n+1/2})+\frac{1}{\Delta t}{\mathbf{u}}^{n+1}+\frac{1}{2}{\mathcal{N}}_m({\mathbf{u}}^{n+1})-\frac{1}{2\mathrm{Re}}\mathcal{L}{\mathbf{u}}^{n+1}-{\mathcal{H}}_m{\mathbf{F}}^n+\mathcal{G}{p}^{n+1/2}$$

真的跟他是一样的

无敌了

他这个的构造的第一行是能量方程，第二行是动量方程，第三行是能量外力的插值方程，第四行是动能外力的插值方程，都是 IBM 的经典核函数插值

然后他这个 cbc mbc ebc 都是什么啊，我无语了

他只说了边界条件会纳入 cbc mbc ebc

但是 cbc mbc ebc 它本身都是什么啊，无语

我大概理解 bc 是 boundary condition 的缩写，m 是动量的缩写，e 是能量的缩写

然后这里明确了速度的散度是 cbc

分离速度梯度

他这个很神奇啊……

最初构造一个大的线代方程组

第一行是能量方程，第二行是动量方程

第三行是能量外力的插值方程，第四行是动能外力的插值方程，都是 IBM 的经典核函数插值

第五行是速度散度 = cbc 的 n+1 时刻

然后他对这个矩阵 A 分块，分块之后的矩阵再通过一个近似，把左上角的那个块 A1 的逆视为 dt 乘以单位矩阵，使得分块之后的矩阵 A 又可以被拆成两个矩阵 B 乘 C

B 里面就可以出现散度算子*梯度算子，把 C 直接和向量相乘，得到新的向量，新的向量里面还是那个压强变量

于是 B 乘以这个新向量就有散度算子*梯度算子*压强 = cbc = 速度散度 也就是泊松方程

而且如果细想的话，他这个泊松方程里面就是多一个 lambda* 的那个东西

A lambda = R 的话，lambda 本身就是 温度 速度 能量外力 动量外力 的向量

lambda* 通过 A1 lambda* = R1 解出来，A1 是 A 的左上角，R1 是 R 的前四项，那 lambda* 似乎也应该具有类似的物理意义

然后 lambda* 的散度出现在泊松方程，就……奇妙

解耦温度和速度

他这里就是再次用了一次 A2 的逆的近似，使得能够 LU 分解

但是我看不出来最后得到这个分解式之后会有什么效果呢？

他最后就说一句，因为是交错时间离散方案，所以温度和速度自动解耦，32 33 很容易解

有点东西……

就，它两次 LU 分解

第一次分解出了压力泊松方程和应用压力梯度两个公式

第二次分解出了温度的对流扩散，速度的对流扩散加浮力项

就很神奇

求解步骤

他最终得到的步骤就很清晰

确实印证了我之前的想法，就是，如果他这里写成了单独的步骤，那么就是说明是可以单独求解

公式推导

动量方程

NS 方程

$$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}+\frac{1}{2}({\mathbf{u}}^{n+1}\cdot \mathcal{G}{\mathbf{u}}^{n+1}+{\mathbf{u}}^n\cdot \mathcal{G}{\mathbf{u}}^n)=-\mathcal{G}{p}^{n+1/2}+\frac{1}{2\mathrm{Re}}(\mathcal{L}{\mathbf{u}}^{n+1}+\mathcal{L}{\mathbf{u}}^n)+{\mathcal{H}}_m{\mathbf{F}}_m+\mathcal{F}({\theta }^{n+1/2})+{\mathbf{mbc}}^{n+1/2}$$

其中 $\mathcal{G}$ 表示梯度算子，$\mathcal{L}$ 表示拉普拉斯算子，${\mathcal{H}}_m$ 表示将 IB 力插值回流场的算子，${\mathbf{F}}_m$ 表示动能外力，$\mathcal{F}$ 表示浮力项的算子，$\theta$ 是温度，${\mathbf{mbc}}^{n+1/2}$ 表示速度的边界条件

他认为，能量方程中的对流项已经因为交错时间离散方案而自动线性化了，所以剩下唯一要做的就是动量方程中的对流项的线性化

$${\mathbf{u}}^{n+1}\cdot \mathcal{G}{\mathbf{u}}^{n+1}+{\mathbf{u}}^n\cdot \mathcal{G}{\mathbf{u}}^n\approx {\mathbf{u}}^{n+1}\cdot \mathcal{G}{\mathbf{u}}^n+{\mathbf{u}}^n\cdot \mathcal{G}{\mathbf{u}}^{n+1}:={\mathcal{N}}_m({\mathbf{u}}^{n+1})$$

这使得 ${\mathcal{N}}_m$ 是一个线性算子，由上式定义，表示速度的对流项

代入速度的对流项的线性化，得

$$\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\Delta t}+\frac{1}{2}{\mathcal{N}}_m({\mathbf{u}}^{n+1})=-\mathcal{G}{p}^{n+1/2}+\frac{1}{2\mathrm{Re}}(\mathcal{L}{\mathbf{u}}^{n+1}+\mathcal{L}{\mathbf{u}}^n)+{\mathcal{H}}_m{\mathbf{F}}_m+\mathcal{F}({\theta }^{n+1/2})+{\mathbf{mbc}}^{n+1/2}$$

移项

$$-\mathcal{F}({\theta }^{n+1/2})+\frac{1}{\Delta t}{\mathbf{u}}^{n+1}+\frac{1}{2}{\mathcal{N}}_m({\mathbf{u}}^{n+1})-\frac{1}{2\mathrm{Re}}\mathcal{L}{\mathbf{u}}^{n+1}-{\mathcal{H}}_m{\mathbf{F}}_m+\mathcal{G}{p}^{n+1/2}=\frac{1}{\Delta t}{\mathbf{u}}^n+\frac{1}{2\mathrm{Re}}\mathcal{L}{\mathbf{u}}^n+{\mathbf{mbc}}^{n+1/2}$$

为了获得压强的差 $\delta p=p^{n+1/2}-p^{n-1/2}$，变为

$$-\mathcal{F}({\theta }^{n+1/2})+\frac{1}{\Delta t}{\mathbf{u}}^{n+1}+\frac{1}{2}{\mathcal{N}}_m({\mathbf{u}}^{n+1})-\frac{1}{2\mathrm{Re}}\mathcal{L}{\mathbf{u}}^{n+1}-{\mathcal{H}}_m{\mathbf{F}}_m+\mathcal{G}\delta p=\frac{1}{\Delta t}{\mathbf{u}}^n+\frac{1}{2\mathrm{Re}}\mathcal{L}{\mathbf{u}}^n-\mathcal{G}{p}^{n-1/2}+{\mathbf{mbc}}^{n+1/2}$$

令

$$\begin{array}{r}{\mathcal{A}}_m=\frac{1}{\Delta t}\mathcal{I}+{\mathcal{M}}_m,\\ {\mathcal{M}}_m=\frac{1}{2}{\mathcal{N}}_m-\frac{1}{2\mathrm{Re}}\mathcal{L},\end{array}$$

$$\begin{array}{rl}{r}_m& =\frac{1}{\Delta t}{\mathbf{u}}^n+\frac{1}{2\mathrm{Re}}\mathcal{L}{\mathbf{u}}^n-\mathcal{G}{p}^{n-1/2}+{\mathbf{mbc}}^{n+1/2}\\ & =(\frac{1}{\Delta t}\mathcal{I}+\frac{1}{2\mathrm{Re}}\mathcal{L}){\mathbf{u}}^n-\mathcal{G}{p}^{n-1/2}+{\mathbf{mbc}}^{n+1/2},\end{array}$$

得

$$-\mathcal{F}({\theta }^{n+1/2})+{\mathcal{A}}_m{\mathbf{u}}^{n+1}-{\mathcal{H}}_m{\mathbf{F}}_m+\mathcal{G}\delta p=r_m$$

对应了原文式 18 中系数矩阵的第二行

这个 ${\mathcal{M}}_m={\mathcal{N}}_m-\frac{1}{2\mathrm{Re}}\mathcal{L}$ 是原文的记法

不知道是不是错了，看上去应该是 ${\mathcal{M}}_m=\frac{1}{2}{\mathcal{N}}_m-\frac{1}{2\mathrm{Re}}\mathcal{L}$

温度方程

温度方程

$$\frac{{\theta }^{n+1/2}-{\theta }^{n-1/2}}{\Delta t}+\frac{1}{2}{\mathbf{u}}^n\cdot (\mathcal{G}{\theta }^{n+1/2}+\mathcal{G}{\theta }^{n-1/2})=\frac{1}{2\mathrm{Pe}}(\mathcal{L}{\theta }^{n+1/2}+\mathcal{L}{\theta }^{n-1/2})+{\mathcal{H}}_e{F}_e+{\mathrm{ebc}}^n$$

其中 ${\mathcal{H}}_e$ 是将温度外力插值回流场的算子，$F_e$ 是温度外力，${\mathrm{ebc}}^n$ 是温度边界条件

原文前面在处理动量方程的时候只说了温度方程的对流项自动线性化，但是他没有把线性化的那个公式给出来，相关的符号也就没有给出说明

我猜测他对温度方程的对流项的线性算子的记法是

$${\mathcal{N}}_{\theta }(\theta )={\mathbf{u}}^n\mathcal{G}\theta$$

代入温度的对流项的线性化，得

$$\frac{{\theta }^{n+1/2}-{\theta }^{n-1/2}}{\Delta t}+\frac{1}{2}({\mathcal{N}}_{\theta }{\theta }^{n+1/2}+{\mathcal{N}}_{\theta }{\theta }^{n-1/2})=\frac{1}{2\mathrm{Pe}}(\mathcal{L}{\theta }^{n+1/2}+\mathcal{L}{\theta }^{n-1/2})+{\mathcal{H}}_e{F}_e+{\mathrm{ebc}}^n$$

移项，得到

$$\frac{1}{\Delta t}{\theta }^{n+1/2}+\frac{1}{2}{\mathcal{N}}_{\theta }{\theta }^{n+1/2}-\frac{1}{2\mathrm{Pe}}\mathcal{L}{\theta }^{n+1/2}-{\mathcal{H}}_e{F}_e=\frac{1}{\Delta t}{\theta }^{n-1/2}-\frac{1}{2}{\mathcal{N}}_{\theta }{\theta }^{n-1/2}+\frac{1}{2\mathrm{Pe}}\mathcal{L}{\theta }^{n-1/2}+{\mathrm{ebc}}^n$$

令

$$\begin{array}{r}{\mathcal{A}}_e=\frac{1}{\Delta t}\mathcal{I}+{\mathcal{M}}_e,\\ {\mathcal{M}}_e=\frac{1}{2}{\mathcal{N}}_{\theta }-\frac{1}{2\mathrm{Re}}\mathcal{L},\end{array}$$

$$\begin{array}{rl}{r}_e& =\frac{1}{\Delta t}{\theta }^{n-1/2}-\frac{1}{2}{\mathcal{N}}_{\theta }{\theta }^{n-1/2}+\frac{1}{2\mathrm{Pe}}\mathcal{L}{\theta }^{n-1/2}+{\mathrm{ebc}}^n\\ & =(\frac{1}{\Delta t}\mathcal{I}-{\mathcal{M}}_e){\theta }^{n-1/2}+{\mathrm{ebc}}^n,\end{array}$$

得

$${\mathcal{A}}_e{\theta }^{n+1/2}-{\mathcal{H}}_e{F}_e=r_e$$

对应了原文式 18 中系数矩阵的第一行

实验

1.加热圆柱
入口速度 狄利克雷边界条件
出口速度 对流边界条件
下壁面和上壁面 诺伊曼边界 表示无限远
没说温度边界条件，猜测是四个面都是绝热边界
只说了圆柱体温度为 1，流体温度为 0
猜测这个说的是初始状态
猜测是保持圆柱体温度不变

2.加热圆柱与方形外壳之间的热对流
没说速度边界条件，猜测是四个面都是壁面边界
没说温度边界条件，猜测是四个面都是等温边界
没说圆柱体初始温度，外壳初始温度
猜测圆柱体温度恒为 1，外壳恒为 0
与有限元结果相对比
以最细网格作为参考解，验证了时间精度和空间精度
误差是全场所有点各个物理量与参考解的 l2 范数

3.加热圆球与方形外壳之间的热对流
速度边界都是壁面条件
没说温度边界条件，猜测是所有面都是绝热边界
只说了圆球温度为 1，流体温度为 0
猜测圆球温度恒为 1，外壳恒为 0
与别人的努塞尔数结果对比

4.单个粒子沉降，有外壳，没有热对流，三维
没说速度边界条件，“外壳”应该默认都是壁面边界
粒子有向下的速度 w，实验对比速度 w 的数值结果与实验结果

5.两个粒子沉降，有外壳，没有热对流，三维
初始的时候两个粒子错开方式，随着沉降的进行，他们之间会有碰撞，然后会有一个粒子陷入另一个粒子的尾流，表述为“kissing”；两个粒子彼此远离，表述为“tumbling”

6.多粒子沉降
做了三个 case，粒子初始全是热的；初始全是冷的；等温粒子
三个 case 的粒子放在相同的高度，比较这三种情况的粒子运动的不同：沉降速度，散射图案，努塞尔数

其他

他这个文章似乎没有说他是怎么实现的，有限差分还是有限体积？

难道说他是要做一个通用的算法？

算精度的时候，计算误差是和最高的分辨率比较