请教高手关于模拟堆积床内焦炭颗粒燃烧，加上源项后出现错误“pressure correction”

yang111457

大家好，我想用fluent模拟堆积床内焦炭颗粒（直径0.005m）燃烧过程传热传质，采用气固双能量方程，组分方程。
遇到的问题：不加气固能量源项和气体组分方程，给定固体初始温度1000K，气体初始温度300K，模拟气固换热是没有问题的；但是加上气固能量源项后，点开始计算，fluent就会提示错误“pressure correction”，是不是能量源项太大？需要对源项进行线性化处理，具体怎么处理呢？还请fluent高手给予指点，小弟不胜感激。

模拟焦炭燃烧过程传热传质

焦炭燃烧反应：C+O2→CO2
焦炭燃点：923K
反应速率：Rc=3.14*（dc）^2*kc *CO2*Nc
       式中，dc是焦炭颗粒直径，kc是综合反应速率常数，CO2是氧气摩尔浓度，Nc是单位体积焦炭颗粒数；
初始条件和边界条件：
Tg，0=1473K，Ts=340K，Co2,0= Co2，Co2,0= Cco2，Cn2,0= Cn2
点火温度1473K(Tg=1473K)，时间60s；保温温度1043K(Tg=1043K)，时间90s；接下来抽风温度300K(Tg=300K)（压力出口，-10KPa）

以下是我的部分udf代码：
real kfc;                               //传质系数,m/s；
    real Sh,Re,Sc;                           //舍伍德数，颗粒雷诺数，施密特准数；
    real mu_lam=C_MU_L(c,t);                //气体黏度；
    real den_gas=C_R(c,t);                    //气体密度，kg/m3；
    real Nc=9.92*pow(10.0,5.0);               //单位体积焦炭颗粒数，1/m3;
    real mole_total;                         //气体总的摩尔数
    real P_total;                            //气体总压,pa
    real P_O2;                             //氧气分压,pa
    real DO2;                              //氧气扩散系数,m2/s;
    real Co;                                //氧气摩尔浓度，mol/m3;
    real O2_mole_frac;                       //氧气摩尔分数；
    real kc1;                               //综合反应速率常数，m/s;
    real k1;                                //反应速率常数,m/s;
    real R1;                                //焦炭完全燃烧的反应速率，mol/m3s;

    P_total=ABS_P(C_P(c,t),op_pres);           //气体总压?
    mole_total=C_YI(c,t,0)/32.0+C_YI(c,t,1)/44.0+C_YI(c,t,2)/28.0; //
    O2_mole_frac=C_YI(c,t,0)/(32.0*mole_total);
    P_O2=P_total*O2_mole_frac;              //氧气分压
    Co=P_O2/R*C_T(c,t);

    DO2=2.00*pow(10.0,-4.0)*pow(0.001,1.78)/P_total;
    Sc=mu_lam/(den_gas*DO2);              //施密特准数舍伍德准数
    Re=ND_MAG(C_U(c,t),C_V(c,t),C_W(c,t))*DIAMETER*den_gas/mu_lam;   //雷诺数
    Sh=2.0/POROUSITY+0.75*pow(Re,0.5)*pow(Sc,0.333)/POROUSITY;      //舍伍德数
    kfc=Sh*DO2/DIAMETER;                                     //传质系数
    k1=6.53*pow(10.0,7.0)*pow(C_UDSI(c,t,0),0.5)*exp(-44000.0/(R*C_UDSI(c,t,0))); //反应速率常数
    kc1=(kfc+k1)/kfc*k1;

    if (C_UDSI(c,t,0)<923.0)

        R1=0;

    else

        R1=3.14*pow(DIAMETER,2.0)*kc1*Co*Nc;

return R1;       //
DEFINE_SOURCE(source_energy_gas,c,t,dS,eqn)
{
    real source;
    real delta_H=-2.79*pow(10.0,5.0);
    real a=0.5;  //反应热分配系数
    source=htc(c,t)*(C_UDSI(c,t,0)-C_T(c,t))-(1-a)*(delta_H*RC(c,t));
    dS[eqn]=0;
    return source;
}

DEFINE_SOURCE(source_energy_solid,c,t,dS,eqn)
{
    real source;
    real delta_H=-2.79*pow(10.0,5.0);
    real a=0.5;  //反应热分配系数
    source=htc(c,t)*(C_T(c,t)-C_UDSI(c,t,0))-a*(delta_H*RC(c,t));
    dS[eqn]=0;
    return source;

aaa-1234

减小时间步长试试