CFD中文网

@李东岳 李老师，直接foamRun还是报一样的错误

@李东岳 李老师，问题解决了，是我把option文件改错了个地方。。新的问题又来了，我在OF11中（仍是该虚拟机）对bed案例进行foamRun或foamRun -solver multiphaseEuler时，为何报错：```
--> FOAM FATAL ERROR:
solvers table is empty

From function static Foam::autoPtr<Foam::solver> Foam::solver::New(const Foam::word&, Foam::fvMesh&)
in file solver/solverNew.C at line 48.

想请问李老师：我使用全系列虚拟机（您已提前配置好），再运行第6步这个简单例子，报错未找到torchFoam这个命令，是啥原因造成的呢 编译过程中有很多warning和note应该不影响吧

@hitsc30 稳定状态感觉没什么影响，但是瞬态的话影响很大。我一般都还是保持在1以下。个人愚见库朗数30+还不如把网格画粗一点。不过结果对得上最重要

@李东岳 东岳老师，请问这个速度能不能从10-2下降到10-4或者-6？

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

有关tfm中spurious velocity的问题，不知各位前辈有没有好的解决方法

@李东岳 东岳老师，从下面几个式子看，固相黏度应该都是半经验表达式，所以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我猜想可能也和固相应力模型有关。
争取早日把成果研究整理分享出来

请教各位前辈老师：


JJ模型中的0.5是哪来的呢？

请教各位前辈：

针对成员函数：

定义：

返回值：

报错：

但直接返回：

计算没问题，不知问题出在哪里

在dargmodel.C这个衍生类中没有mesh这个成员对象，不知道该如何引入

如果改变写法：

Foam::dragModels::wallEffect::wallEffect
(
    const dictionary& dict,
    const phasePair& pair,
    const bool registerObject
   
)
:
    dragModel(dict, pair, registerObject),
    residualRe_("residualRe", dimless, dict),
    y_(nullptr)
       
{}

Foam::dragModels::wallEffect::wallEffect
(
    const dictionary& dict,
    const phasePair& pair,
    const bool registerObject,
    const fvMesh& mesh
   
)
:
    dragModel(dict, pair, registerObject),
    residualRe_("residualRe", dimless, dict),
    y_(&wallDist::New(mesh).y())
       
{}

计算的时候会报错：

@李东岳 也烦请东岳老师帮忙看看这个问题如何实现~ 三月再您的课上还问过您这个问题，但是下来自己搞没搞通