这个是当我把水的物性从常数改为多项式后遇到的问题
Courant Number mean: 0.134868655 max: 0.237479859
deltaT = 3.4820232e-06
Time = 0.00378807

PIMPLE: iteration 1
MULES: Solving for alpha.gas
MULES: Solving for alpha.liquid
gas fraction, min, max = 0.00155200072 1e-15 0.0672465558
liquid fraction, min, max = 0.998447999 0.932753444 1
Phase-sum volume fraction, min, max = 1 0.999994574 1.00026614
MULES: Solving for alpha.gas
MULES: Solving for alpha.liquid
gas fraction, min, max = 0.00155951826 1e-15 0.0674835707
liquid fraction, min, max = 0.998440482 0.932516429 1
Phase-sum volume fraction, min, max = 1 0.999994794 1.00025528
MULES: Solving for alpha.gas
MULES: Solving for alpha.liquid
gas fraction, min, max = 0.00156703579 1e-15 0.0677203557
liquid fraction, min, max = 0.998432964 0.932279644 1
Phase-sum volume fraction, min, max = 1 0.999995014 1.00025493
Constructing momentum equations
Tf.gasAndLiquid: min = 350.131001, mean = 456.035361, max = 549.778465
iDmdt.gasAndLiquid: min = -0.959612341, mean = 11.8287784, max = 447.419943, integral = 0.000182163188
wDmdt.gasAndLiquid: min = 0, mean = 138.094307, max = 5727.35791, integral = 0.00212665234
smoothSolver: Solving for e.gas, Initial residual = 0.104953394, Final residual = 1.40390231e-09, No Iterations 2
smoothSolver: Solving for e.liquid, Initial residual = 0.282448508, Final residual = 6.30884149e-18, No Iterations 2
[stack trace]

#1 Foam::sigFpe::sigHandler(int) in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#2 ? in /lib/x86_64-linux-gnu/libpthread.so.0
#3 Foam::sqrt(Foam::Field<double>&, Foam::UList<double> const&) in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#4 void Foam::sqrt<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so
#5 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::sqrt<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so
#6 Foam::heatTransferModels::RanzMarshall::K(double) const in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libreactingMultiphaseSystem.so
#7 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::BlendedInterfacialModelFoam::heatTransferModel::evaluate<double, Foam::fvPatchField, Foam::volMesh, double>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > (Foam::heatTransferModel::*)(double) const, Foam::word const&, Foam::dimensionSet const&, bool, double) const in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libreactingMultiphaseSystem.so
#8 Foam::BlendedInterfacialModelFoam::heatTransferModel::K(double) const in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libreactingMultiphaseSystem.so
#9 Foam::ThermalPhaseChangePhaseSystem<Foam::PhaseTransferPhaseSystem<Foam::TwoResistanceHeatTransferPhaseSystem<Foam::MomentumTransferPhaseSystemFoam::multiphaseSystem > > >::correctInterfaceThermo() in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libreactingMultiphaseSystem.so
#10 ? in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/bin/reactingMultiphaseEulerFoam
#11 __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#12 ? in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/bin/reactingMultiphaseEulerFoam

最近我在使用reactingMultiphaseEulerFoam计算模拟垂直管道内水蒸气壁面沸腾，长度为 2m，直径为 15.4 mm，管壁热通量为 570000w.m-2，入口处的平均速度为 1m.s-1，工作压力为 45 bar。两个阶段的湍流模型都是 k-epsilon。我在选用系统自带饱和模型的时候，计算一段时间就会崩溃，只有我将饱和模型选为常数的时候才可以模拟成功，但是模拟得到的结果偏差较大，这种情况怎么解决呢？以下是我选用的三种饱和模型：
saturationModel
{
type constant;
pSat 100000;
Tsat 373;
};

/*
{
type polynomial;
C<8>
(
308.0422
0.0015096
-1.61589e-8
1.114106e-13
-4.52216e-19
1.05192e-24
-1.2953e-30
6.5365e-37
);
};
/
/
{
type polynomial;
C<8>
(
440.9049
2.5741e-5
-1.28758e-12
-1.2014e-28
2.45488e-35
-2.94338e-42
1.92021e-49
-5.26545e-57
);
};*/

我用系统自带的网格（ 1 75 1）可以算出来，但是我改成（30 2000 1）就算崩溃了，打算算高2m宽15.4mm的圆管截面

我去尝试修改一下粒径和曳力模型试一试，还有一个问题就是我在计算同一个模型网格量很少的时候就可以计算成功，但是我增大网格密度，就不会算收敛，这是为什么呢

我使用原始算例不修改任何参数，计算出来的内部场分布均匀，但是当我修改网格或者修改一些工况之后就立马变成这个情况，那reactingTwophaseEulerFoam这个求解器会好一些吗？

壁面热流为0.57MW

我是根据reactingMultiphaseEulerFoam/RAS/wallboiling算例修改的，修改了工况为4.5MPA,进口470k，重新拟合了压力与饱和温度关系式，计算出来的结果和文献中想复现的结果完全不同，

第一张图是我自己计算出来的alpha.gas，第二张图是文献中的计算云图