如何提高方程求解的稳定性

[Do1975]

下面是我需要计算的方程：

fvScalarMatrix TbEqn
(
    - fvm::div(-phi, Tb) 
    + fvm::laplacian(DT, Tb)
);
TbEqn -= Q;
TbEqn.solve();

其中，DT为热传导系数，Q为数值1的标量场，phi通过SIMPLE求解。
设置出口边界条件为：

outlet
{
    type                 groovyBC;
    variables            "Tk=DT;h=U&normal();Tinf=0;f=1/(1+Tk/(h*mag(delta())));";
    valueExpression      "Tinf";
    gradientExpression   "0";
    fractionExpression   "f";
    value                uniform 0;
}

设置的divSchemes为：

div(-phi,Tb)     Gauss linearUpwind grad(Tb);

使用的solvers为：

"(Tb)"
{
     solver                  PBiCGStab;
     preconditioner          DILU;
     tolerance               1e-6;
     relTol                  0;
 }

上述方程是用于拓扑优化求解的，但是优化过程中，经常出现浮点数错误“Signal: Floating point exception (8)”，如何有效解决这个问题呢？
尝试过在h*mag(delta()中添加一个小数，但是也同样出现浮点数报错。