参考文献仍然是An accurate adaptive solver for surface-tension-driven interfacial flows

标准的2D规则网格曲率计算：

1. 以欲计算曲率的格子$(i,j,k)$为中心，搞一个$3\times 7$或者$7\times 3$的stencil。具体是哪种取决于之前重建界面的时候MYC算出来的法向朝哪个方向。
2. 在这个stencil里面把每一列的VOF值相加，得到高度函数$y=h(x)$或者$x=h(y)$.
3. 使用有限差分计算曲率$\kappa=\frac{h”}{(1+h’^2)^{3/2}}$.

这样的问题是$3\tim[……]

继续阅读

参考文献：Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry

在流体模拟中，我们可以用VOF（volume of fluid，体积分数）方法追踪两相之间的界面。而VOF场的演化需要对流，做这种对流的最佳方案之一是用所谓的piecewise linear方法，即在每个格子内部估计一个平面作为局部界面，在对流时根据这个界面计算通过格子之间界面的物质通量。

而在这种方法当中，想要估计平面，需要先计算出这个平面的法向，然后再确定其截距。[……]

继续阅读

参考文献：An accurate adaptive solver for surface-tension-driven interfacial flows

表面重建

计算截距

做advection的时候，认为每个格子里的界面都是一个平面（我们默认三维，二维就是直线）。也就是说有一个法向矢量$\bm{m}$和一个常数$\alpha$，那么这个格子的界面就是

$$
\bm{m}\cdot\bm{x}=\alpha.
$$

可以发现，如果这个格子的体积分数$c$（也就是VOF）和$\bm{m}$都已确定，那么$\alpha$就是唯一确定的。确定的方法类似marching cubes，即分类[……]

继续阅读

Conjugate Gradient

对于对称正定矩阵$A$和向量$b$，构造优化目标函数$\phi(x)=\frac{1}{2}x^TAx-b^Tx$.

可以用张量记号写作：$\phi=\frac{1}{2}x_ia_{ij}x_j$.

对$x_i$求导得：

$\frac{\partial\phi}{\partial x_i}=\frac{1}{2}a_{ij}x_j+\frac{1}{2}x_ja_{ji}$.

由于$A$对称，故这两项相等，即$\frac{\partial\phi}{\partial x_i}=a_{ij}x_j=Ax$.

Preconditioner

ht[……]

继续阅读

ConjugateGradient类

Eigen使用ConjugateGradient类进行共轭梯度法计算。一个代码示例为（https://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html ）：

int n = 10000;
VectorXd x(n), b(n);
SparseMatrix<double> A(n,n);
// fill A and b
ConjugateGradient<SparseMatrix<double>, Lower|Upper> cg;
cg.compu[......]
继续阅读

本文提到的方法见于：Tartakovsky, Alexandre M., and Paul Meakin. “A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional Rayleigh–Taylor instability.” Journal of Computational Physics 207.2 (2005): 610-624.

普通SPH形式

密度： $\rho_i=\sum_{j}m_jW_{ij}.$

其中$[……]

继续阅读

值的表示

对于一维情况，现假设我们有一个定义在非负实数域上的标量权值函数$\phi_1(x)$，截断至1，即满足：$\int{\phi_1(x)}\mathrm{d}x=1$，且$\phi_1(x)=0$对$\forall x>1$成立。

我们希望基于它设计一个关于$\vec{r}$的对称核函数$W_1(\vec{r})$，截断至$h$，那么这个核函数的形式应为

$$W_1(\vec{r})=\alpha_1 \phi_1\left(\frac{|\vec{r}|}{h}\right).$$

它的积分应为1，因此

$$\int W_1(x)\mathrm{d}x=h\alpha_1\i[……]

继续阅读

3.10收敛性测试和数值验证

首先是每算例都做三组：

• 固定[latex]\epsilon[/latex]，改变网格尺寸h，检查[latex]h\to 0[/latex]的收敛情况。
• 固定比例关系[latex]\epsilon=\alpha h[/latex]，其中[latex]\alpha[/latex]为一常数，检查[latex]h\to 0[/latex]的收敛情况。
• 交换极限，先算一个[latex]\epsilon\to 0^+[/latex]的内极限，然后检查[latex]h\to 0[/latex]的收敛情况。

所有测试都是在规则正方形/立方体网格上完成的。

3.[……]

继续阅读

文章链接：https://doi.org/10.1016/j.cpc.2021.107849

算例部分在第四章。

4.1基础测试

4.1.1泊肃叶流动

设置一个8*8*2网格，模拟二维正方形水槽内的泊肃叶流动。参数[latex]||\nabla p||=\mu=\rho=1[/latex].在水槽左端右端设置压强边界条件。模拟至流场以1e-3精度稳定，约在1000步后，时刻1达成。使用显式viscous diffusion.上下边界处使用粘性边界条件u=0.

4.1.2二维斯托克斯流

和泊肃叶流动相同，但计算区域中心有一直径为0.5的圆盘。不进行advect[……]

继续阅读

现假设我们有一流体格i，和与之相邻的格子j。我们研究投影矩阵A中，代表(i,j)格子对的元素（下简记为A[i,j]）和i格子所对应的对角元（下简记为A[i,i]）的数值情况。注意，i和j都是格子，本身坐标是一个向量，不要与矩阵坐标相混淆。

现列举j的几种情况。

情况一：j亦是流体格。此时，j格压强也是未知数，A[i,j]=-1，而A[i,i]中有一项1，如下图。（或在porous flow中，面积分数为a，那就是-a和a）

情况二：j亦是流体格，但i和j中间有一Neumann边界条件（虽然这很奇怪，我怀疑它会不会出现）。此时格子i的方程中，i-j相关联的项直接被取消，[……]

继续阅读