界面温度超过饱和温度，相变异常

[hongjiewang]

各位老师好，我想要实现的是只在界面处发生相变。下壁面温度398.15K充当热源，先以导热的传热方式提高相界面温度，之后从界面处发生相变。我将模型核边界条件的设置附在下方。

目前存在的问题是：在界面超过饱和温度后会发生相变，刚开始的相变是正常的，但是当整个区域的液相含量达到0.38左右时，相变几乎停止。通过查看温度分布，界面温度是超过373.15K的。


分析下来，最大出问题的地方是我的压强边界的设置。内部场和上壁面的压强都是假设大气压，所以气体可能逸不出去？导致上部气体积聚，增大界面处的界面热阻？所以麻烦各位老师看看我的算例边界条件是否正确。谢谢~

[李东岳]

自己写的求解器？

[hongjiewang]

@李东岳 嗯嗯是的，老师。我是基于interCondensatingEvaporatingFoam写的。
经过我导出一些数据，我发现了存在的问题。之前我的相界面是根据老师提示的用梯度进行判断。
1.在createFields里定义变量

Info<< "a\n" << endl;
volScalarField a
(
    IOobject
    (
      "a",
      runTime.timeName(),
      mesh,
      IOobject::READ_IF_PRESENT,
      IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("a", dimensionSet(0,0,0,0,0,0,0), 1)
); 

surfaceScalarField test=fvc::snGrad(alpha1);

Info<< "test1\n" << endl;
volScalarField test1
(
    IOobject
    (
      "test1",
      runTime.timeName(),
      mesh,
      IOobject::NO_READ,
      IOobject::AUTO_WRITE
    ),
    fvc::surfaceIntegrate(test)
);//识别出相界面
forAll(a,celli)
{
    if (test1[celli] == 0)
    {
        a[celli]= 0;
    }
    else
    {
	a[celli]= 1;
    }
}//判断网格是否处于相界面

2.在.C文件中alphaEqnSubCycle.H后添加updateVarities.H文件，即在相方程计算结束进行test1和a的更新

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

请问老师这段程序哪里可能出现了问题，导致a和test1不会更新，一直处于createFields里的数值。

[hongjiewang]

@李东岳 老师，这个问题已经解决了~不过还是麻烦各位老师看一下我的边界条件是否合适，不是特别清楚对不对~

[hongjiewang]

还有一个问题需要请教~在我用Clausius-Clapeyron 方程
加入饱和压力与饱和温度的影响，编译通过，但是求解算例时候会（在没有考虑饱和温度随压强变化下的算例是可以正常运算的）。

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libpthread.so.0
#3  ? in /lib/x86_64-linux-gnu/libm.so.6
#4  Foam::log(Foam::Field<double>&, Foam::UList<double> const&) at ??:?
#5  void Foam::log<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#6  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::log<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#7  Foam::temperaturePhaseChangeTwoPhaseMixture::TSatLocal() const at ??:?
#8  Foam::temperaturePhaseChangeTwoPhaseMixtures::Lee::mDotAlphal() const at ??:?
#9  Foam::temperaturePhaseChangeTwoPhaseMixture::vDotAlphal() const at ??:?
#10  ? in ~/OpenFOAM/dyfluid-v2006/platforms/linux64GccDPInt32Opt/bin/twoPhaseChangePorousFoam
#11  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#12  ? in ~/OpenFOAM/dyfluid-v2006/platforms/linux64GccDPInt32Opt/bin/twoPhaseChangePorousFoam
Floating point exception (core dumped)

[李东岳]

自己编译求解器很难调试，因为求解的方程都是耦合在一起的，一般发散会互相影响，不好判断问题在哪里

[hongjiewang]

@李东岳 实在是有点困难~请问老师如何让两个场（如图片中的a和vDotAlpha）对应网格的值相乘，这不是内积也不是外积，用forAll循环的话，会报assignment of read-only location错误。

forAll(a,celli)
    {
        vDotcAlphal[celli] = a[celli]*vDotcAlphal0[celli];
        vDotvAlphal[celli] = a[celli]*vDotvAlphal0[celli];
    }