VOF两相流中表面曲率的计算

参考文献仍然是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[……]

继续阅读

Mixed Youngs-Centered格式笔记

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

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

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

继续阅读

Piecewise-Linear VOF Advection笔记

参考文献: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[……]

继续阅读

Eigen库的共轭梯度法(Conjugate Gradient)代码分析

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[......]

继续阅读

数值密度SPH(Number Density SPH)的一些推导

本文提到的方法见于: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}.$

其中$[……]

继续阅读

关于SPH通用核函数的一些推导

值的表示

对于一维情况,现假设我们有一个定义在非负实数域上的标量权值函数$\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[……]

继续阅读

Robert Saye:VIIM论文算例简介

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.[……]

继续阅读

PARIS模拟器算例介绍

文章链接: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[……]

继续阅读

关于MAC网格不可压流体求解中投影矩阵的笔记

现假设我们有一流体格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相关联的项直接被取消,[……]

继续阅读