您好，请问这个问题您是否有好的解决办法了吗？
强制设为0可能不太对，反映出来的可能不是压实效果。。

@李东岳 李老师您客气了！一直记得

感谢李老师搭建的这么好的平台，论坛的其他老师们也都乐意帮助像我这样的新手答疑解惑！师者，传道授业解惑也。李老师虽不在高校，但是还是要说一声：祝李老师节日快乐！也祝论坛的其他老师节日快乐！

@李东岳 东岳老师，从下面几个式子看，固相黏度应该都是半经验表达式，所以0.5可能是个实验获得的参数：
（JJ）
（Schaeffer）
（Princeton）

此外，princeton中摩擦应力计算如下：

对通量的处理：

const volVectorField& U = phase.U();
surfaceScalarField phiU = fvc::interpolate(U) & U.mesh().Sf();
 volScalarField Ud = fvc::div(phiU);

对压力计算：

const volScalarField pc = Fr_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_)
                            /pow(max(alphaMax - alpha, alphaDeltaMin_), p_);
......

//对体
forAll(n,celli) 
    {
        if (Ud[celli] < 0.0)
        {
            n[celli] = 1.03;
        }
        else
        {
            n[celli] = 0.5*sqrt(3.0)*sin(phi_.value());            
        }
        m[celli] = n[celli]-1;
    
        //..No negative base. 对不对？
        pt[celli] =  1.0-Ud[celli]/( n[celli]*sqrt(2.0)*sin(phi_.value())* (sqrt(Us[celli]) + small) );
        if (!(pt[celli] > 0.0))
        {
            pt[celli] = small;
        }
        
        //..No negative exponent.
        if (m[celli] > 0)
        {
            //n-1 power
            pt[celli] = pc[celli]*pow(pt[celli], m[celli]);
        }
        else
        {
            pt[celli] = pc[celli]*pow(1/pt[celli], mag(m[celli]));
        }
    }              
......

//对面
forAll(currPatch,facei)
        {
            if(UdBf[patchi][facei] < 0.0)
            {
               nBf[patchi][facei] = 1.03;
            }
            else
            {
               nBf[patchi][facei] = 0.5*sqrt(3.0)*sin(phi_.value());
            }
            mBf[patchi][facei] = nBf[patchi][facei]-1;

            ptBf[patchi][facei] = 1.0-UdBf[patchi][facei]/( nBf[patchi][facei]*sqrt(2.0)*sin(phi_.value())* (sqrt(UsBf[patchi][facei]) + small) );
            
            if (!(ptBf[patchi][facei] > 0.0))
            {
               ptBf[patchi][facei] = small;
            }
        
            //..No negative exponent.
            if (mBf[patchi][facei] > 0)
            {
                ptBf[patchi][facei] = pcBf[patchi][facei]*pow(ptBf[patchi][facei], mBf[patchi][facei]);
            }
            else
            {
                ptBf[patchi][facei] = pcBf[patchi][facei]*pow(1/ptBf[patchi][facei], mag(mBf[patchi][facei]));
            }
        } 
    }  
    
    // Correct BCs
    pt.correctBoundaryConditions();
    Us.correctBoundaryConditions();
    Ud.correctBoundaryConditions();

黏度计算：

//对体
forAll(Ud,celli) 
    {
        if (Ud[celli] < 0)
        {
            n[celli] = 1.03;
        }
        else
        {
            n[celli] = 0.5*sqrt(3.0)*sin(phi_.value());            
        }
        m[celli] = n[celli]-1;

        if (!(pc[celli] > 0))
        {
            pc[celli] = small;
        }
        
        //..pf is a const.
       
        //calculation of nu
        if (m[celli] > 0)
        {            
            nu[celli] = sqrt(2.0)*pf[celli]*sin(phi_.value())*
                        ( n[celli]-(n[celli]-1.0)*pow(pf[celli]/(pc[celli]), 1/m[celli]) )
                        /(sqrt(Us[celli]) + small);
        }
        else
        {
            nu[celli] = sqrt(2.0)*pf[celli]*sin(phi_.value())*
                        ( n[celli]-(n[celli]-1.0)*pow(pc[celli]/(pf[celli]+small), mag(1/m[celli])) )
                        /(sqrt(Us[celli]) + small);
        }
    }
......


//对面
 forAll(patches, patchi)
    {
     
      if (!patches[patchi].coupled())
      {
        const fvPatch& currPatch = patches[patchi];
        forAll(currPatch,facei)
        {
            if(UdBf[patchi][facei] < 0.0)
            {
               nBf[patchi][facei] = 1.03;
            }
            else
            {
               nBf[patchi][facei] = 0.5*sqrt(3.0)*sin(phi_.value());
            }
            mBf[patchi][facei] = nBf[patchi][facei]-1;

            if (!(pcBf[patchi][facei] > 0.0))
            {
                pcBf[patchi][facei] = small;
            }

                              
            if (mBf[patchi][facei] > 0)
            {
                nuBf[patchi][facei] = sqrt(2.0)*pf.boundaryField()[patchi][facei]*sin(phi_.value())*
                                      ( nBf[patchi][facei]-(nBf[patchi][facei]-1.0)*pow(pfBf[patchi][facei]/pcBf[patchi][facei], 1/mBf[patchi][facei]) )
                                      /(sqrt(UsBf[patchi][facei]) + small);
            }
            else
            {
                nuBf[patchi][facei] = sqrt(2.0)*pf.boundaryField()[patchi][facei]*sin(phi_.value())*
                                      ( nBf[patchi][facei]-(nBf[patchi][facei]-1.0)*pow(pcBf[patchi][facei]/(pfBf[patchi][facei]+small), mag(1/mBf[patchi][facei])) )
                                      /(sqrt(UsBf[patchi][facei]) + small);
            }
        }
      } 
    } 
    
    nu.correctBoundaryConditions();
    Us.correctBoundaryConditions();
    pc.correctBoundaryConditions(); 
    Ud.correctBoundaryConditions(); 

结果在相分数0.5左右（摩擦黏度产生处）报错：

相分数出现突变，而且无法增大（继续堆积），增大就会报错，不知道是哪个地方的处理有问题，还请东岳老师和其他前辈能解答一二！

@李东岳

@李东岳 谢谢李老师！0.5暂时还没有找到足够有说服力的依据，只是看到有些文献里直接用的μ=0.5pfsin(phi)；
但我觉得wachem2000里面公式不太对，因为他引用的JJ的文章是这样定义的：

请教各位前辈一个问题：

phi是一个surfaceScalarField，和volScalarField不矛盾吗？

此外，现在需要对速度散度判断，以确定n的取值：

n是一个volScalarField或者scalar也可以，要计算：

写循环的话，应该用div（U），还是div（phi）？（ps和pc都是volScalarField）

@Yongbo 谢谢

@李东岳 李老师，请问发散后继续计算又不发散了，比如如下报错：

或者报错与solver有关（gauss-seidel），这证明植入的模型还是不容易收敛是吧？

@veen 还在整理研究结果 ，但是也有了一点新认识：
1）沉降效果与曳力模型和颗粒粒径有关；
2）砂堤回缩与边界条件和固相应力模型（猜测，用JJ比JJShaffer要好，且正在植入其他模型）有关；
3）达不到alphaMax我猜想可能也和固相应力模型有关。
争取早日把成果研究整理分享出来