@李东岳 东岳老师，从下面几个式子看，固相黏度应该都是半经验表达式，所以0.5可能是个实验获得的参数：
d7cd3790-cd15-45c3-bad9-6dcba069f6f6-image.png （JJ）
bcc611d5-bc11-49c6-beb6-73d08b8496ad-image.png （Schaeffer）
a5a9446d-2284-45ea-b94d-84d8125963b2-image.png （Princeton）
此外，princeton中摩擦应力计算如下：
b334a3b1-6081-4a12-bb32-5f3c04cca839-image.png
3e038a48-fde7-41ec-9125-df2a1283aecb-image.png
对通量的处理：
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左右（摩擦黏度产生处）报错：
149e14c3-5f93-4325-980c-2e15732c8f04-image.png
36f75fde-3108-44a1-88ee-ba7eee7d1aa7-image.png
相分数出现突变，而且无法增大（继续堆积），增大就会报错，不知道是哪个地方的处理有问题，还请东岳老师和其他前辈能解答一二！
