关于multiphaseEulerFoam中固相摩擦黏度的计算问题

Zhy2022

各位前辈，请问在固相摩擦黏度求解过程中，如下处理的意义是什么？

在Wachem2000的论文里，摩擦黏度的定义式也不是这样写的：

此外这行代码对应的公式应该是怎样的呢？

sqrt((1.0/3.0)*sqr(tr(D[celli])) - invariantII(D[celli]))

李东岳

http://dyfluid.com/docs/tensor.html 公式可以在这里翻译

Wachem2000那个是JohnsonJackson固相粘度，你这个是Schaeffer模型，

Zhy2022

@李东岳多谢东岳老师！学生还有一事不太清楚：
invariantII有如下两个版本：
-8
-2006
像8这样定义的话就没算3D的吗？但是tr前面仍是取的1/3的哇——(1.0/3.0)*sqr(tr(D[celli])

Zhy2022

@Zhy2022 不好意思，8中也有非2D的定义-SymmTensorI.H

Zhy2022

请教各位前辈老师：

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

Zhy2022

请教各位前辈一个问题：

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

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

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

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

李东岳

@zhy2022

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

参考 http://dyfluid.cn/无痛苦ns方程笔记.pdf 6.17节

后面那个用phi

另外，你那个0.5找到理论依据了么

Zhy2022

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

李东岳

@zhy2022 我琢磨下帮你处理下这个问题

Zhy2022

@李东岳

Zhy2022

@李东岳东岳老师，从下面几个式子看，固相黏度应该都是半经验表达式，所以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左右（摩擦黏度产生处）报错：

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

CFD中文网

关于multiphaseEulerFoam中固相摩擦黏度的计算问题