求解过程中出现的错误

史浩

在多相流求解器中，植入温度方程进行求解，求解一段时间后，会出现下面的错误：

Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solve() at ??:?
#8  ? at ??:?
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  ? at ??:?

让算例续算，发现求解其可以越过该问题时刻继续算下去，但是过了某一段时间后又出现相同的错误。
根据log文件发现问题出现在我的植入的温度方程求解的过程中。温度方程求解代码如下：

    fvScalarMatrix TEqn
    (
        fvm::ddt(rhoCp, T)
      + fvm::div(fvc::interpolate(rhoCp)*phi, T)
      - fvm::laplacian(alphat, T)
     ==
        fvc::Su(jouleHeat,T)
    );

    TEqn.relax();

    fvOptions.constrain(TEqn);

    TEqn.solve();

    fvOptions.correct(T);

其中，jouleHeat是电流产生的焦耳热源项，该项恒定。
各位觉得问题可能出在哪？

danao

写成fvm::Sp()?

史浩

@danao 这个是显示源项，是个常数，用Su函数应该没有问题。如果是函数的问题，应该一开始就会报错，应该算一会儿才报错

史浩

有的时候报错报的不完整，完整的错误如下：

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PBiCGStab::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#6  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solve() at ??:?
#8  ? at ??:?
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  ? at ??:?
Floating point exception (core dumped)

李东岳

T不与任何变量耦合？

史浩

@东岳在速度场中有一个浮升力项包含了T，其他的没有了，不知道这算不算耦合

李东岳

把速度方程中这一项去掉试试?

史浩

@东岳多谢东岳。经过一天半的追踪，我发现错误出在了边界条件上。一个边界条件是我重写的边界条件，这个边界条件和mixed很类似。问题应该是在这个函数上：

const scalarField& k_ = 
    patch().lookupPatchField<volScalarField, scalar>("alphat");

因为我的边界条件和另外一个场量alphat有关，在边界条件设置的时候，需要获得这个变量的边界上的场值。我再fvPatch.H下找到了这个函数，写在了程序中。这个函数在非动网格算例中运行很稳定，没有出现错误，但是在拓扑变化的网格中就不太稳定，时不时的会出错，但是出错的地方并不是网格拓扑变化的时候，而是随时都有可能出错。当我把这个函数去掉，程序就可以稳定运行了。
请问东岳老师，还有没有其他的方式获得边界上的场变量？

李东岳

@史浩在求解过程中出现的错误中说：

当我把这个函数去掉，程序就可以稳定运行了。

你之前只是声明，没调用?

史浩

@东岳有调用。我说的去掉是指不用这个函数，随便给这个k_赋值，他就没有问题。但是这样明显不符合物理意义，只是单纯的程序能运行而已

void Foam::robinTemperatureFvPatchScalarField::evaluate(const Pstream::commsTypes)
{
    if (!this->updated())
    {
        updateCoeffs();
    }

    scalarField k_ = //patch().deltaCoeffs();  //随便给个值
        patch().lookupPatchField<volScalarField, scalar>("alphat");

    valueFraction() =
        1.0/
        (
            1.0
          + k_*patch().deltaCoeffs()/(h_ + 1e-10)
        );

    mixedFvPatchScalarField::evaluate();
}

robinTemperatureFvPatchScalarField是mixedFvPatchScalarField的子类

史浩

顺便补充一下，这个函数的执行没什么问题，而是在solve的时候出现的错误。但是把这个去掉，求解就不会有问题。

李东岳

@史浩在求解过程中出现的错误中说：

函数的执行没什么问题，而是在solve的时候出现的错

valueFraction() =
        1.0/
        (
            1.0
          + k_*patch().deltaCoeffs()/(h_ + 1e-10)
        );

会不会是这个的问题呢？