界面相变

[李东岳]

应该通过相分数的梯度来判断

[hongjiewang]

@李东岳 嗯嗯好的，我会尝试一下的。
有一位前辈推荐了我一个utility，或许也有帮助，贴在下方
https://github.com/daidezhi/alphaInitializerFoam

[hongjiewang]

@李东岳
老师，请问下面的语句应该怎么进行修改一下~a是一个volScalarField类型的量

forAll(a,celli)
    {
        if (fvc::snGrad(alpha1[celli]) == 0)
	{
	a[celli]= 0;
	}
	else
	{
	a[celli]= 1;
	}
    }

[李东岳]

volScalarField test = fvc::snGrad(alpha);

forAll(test, celli)
    {
        if (test[celli] == 0)
	{

	}
	else
	{

	}
    }

[hongjiewang]

@李东岳 老师，我按照下面的进行了修改

surfaceScalarField test = fvc::snGrad(alpha1);
forAll(a,celli)
    {
        if (test[celli]) == 0)
	{
	a[celli]= 0;
	}
	else
	{
	a[celli]= 1;
	}
    }

并且看a的输出文件，确实有a=1的网格（即相梯度不等于0）。但是查看运行结果，他并不是在相梯度不等于0的地方相变，还是和之前一样温度高了就相变，并没有抑制作用。我的处理方法是在和相变相关的量前乘以a进行判断，请问老师这样处理哪里出现了问题~

[李东岳]

@hongjiewang 在 界面相变 中说：

还是和之前一样

如果完全一样，应该是代码某些部分没起作用。更深的你需要一个专门的1D或者2D算例做debug

[史浩]

@hongjiewang 你的a是一个volScalarField，而test是一个surfaceScalarField，两个变量不匹配。想办法把surfaeScalarField转化成一个volScalarField试一下

[队长别开枪]

@hongjiewang 在 界面相变 中说：

https://github.com/daidezhi/alphaInitializerFoam

通过 https://github.com/daidezhi/alphaInitializerFoam 精确初始化相场后，使用$\alpha ( 1- \alpha )$的值进行判断，在$[0,1]$区间其值位于$[0,0.25]$之间，为了剔除相分数十分接近0或者1(体单元接近empty或者full，有些单元由于数值计算的误差累计会出现该为empty的时候其相分数是个十分小的数值，full也有这种可能)的情况，加入阈值判断。代码可以参考

// Tolerance
const scalar tol(1e-7);

const volScalarField test(alpha1*(scalar(1)-alpha1));

forAll(test, cellI)
{
    if (mag(test[cellI]) > tol)
    {
        // do something ...
    }
    else
    {
        // do something else ...
    }
}

[hongjiewang]

@史浩 确实如老师所说，我对这部分进行修改之后，已经可以算出预期的效果。感谢老师。

[hongjiewang]

@队长别开枪 感谢老师的提醒，之后我会对代码进一步的完善。非常感谢~

[史浩]

@hongjiewang 能帮到你我也很开心，继续加油~

[hongjiewang]

@史浩 老师，还得麻烦您再帮我看看是哪里出了问题~我是在createFields里定义了a、test、test1这三个变量识别相界面的。并在updateVariables.H（位于相含量求解后）中进行更新

test=fvc::snGrad(alpha1);
test1=fvc::surfaceIntegrate(test);
forAll(a,celli)
{
    if (test1[celli] == 0)
    {
        a[celli]=0;
    }
    else
    {
        a[celli]=1;
    }
}

用a去乘以质量传递速率去控制在该网格是否发生相变，比如
a*fvm::Sp(vDotvmcAlphal, alpha1) + a*vDotcAlphal

结果发现，识别相界面的变量是正确的，但是他的相变进程依然是按照温度分布在进行，并没有考虑从相界面开始。

[史浩]

@hongjiewang 在 界面相变 中说：

afvm::Sp(vDotvmcAlphal, alpha1) + avDotcAlphal
从你给的代码里面，没看出什么问题，可能是其他方面可以出问题了
相含量的变化不一定只和源项有关，还可能和宏观流动、求解误差等相关
建议跟踪一下你给的源项这个东西的值

a*fvm::Sp(vDotvmcAlphal, alpha1) + a*vDotcAlphal

如果这个值没有问题，那看一下你写的代码，看看你代码计算的过程是不是和现实中的物理过程有偏差。比如说这个地方相分数为1，但是你的源项可能是一个正值。因为你用他的散度去表示界面处的相变过程，这个表示本身就是不严谨的。
话说回来，你为啥不用相分数，来判断是否为相界面？代码如下：

forAll(a,celli)
{
    if (alpha1[celli] > 0.001 && alpha1[celli] < 0.999)
    {
        a[celli]=1;
    }
    else
    {
        a[celli]=0;
    }
}

[hongjiewang]

非常感谢老师的回复

@史浩 在 界面相变 中说：

相含量的变化不一定只和源项有关，还可能和宏观流动、求解误差等相关

从运算出来的test1看，即使规定它不相变，网格内的相含量也是有变化的，这也导致我采用梯度是否为0判断相界面的不准确性，和老师所说的一致。目前我将test1==0修改为一个范围值，也不知道是否正确。

不用相分数，来判断相界面

整个计算空间是1✖50✖1mm（网格划分为50✖250✖50），初始时利用setFields将1✖50✖0.2mm区域设置为液体区，即0时刻，每个网格的液体含量非0即1，所以可能没办法利用相分数判断0时刻的相界面。

[hongjiewang]

@史浩 按照老师说的用alpha识别相界面，不管我设置的那个横截面是否把网格切分，比如我设置z轴0.0002m一下为液体区域，当z轴方向50个网格时，初始计算他的液相含量是0.2；当我z轴方向51个网格时，初始计算他的液相含量变成了0.196078，导致0文件中alpha非0即1，因此并没有识别出相界面。