[lihaoyan]

Hi, guys.
I built a threePhasechangefoam based on InterphaseChangefoam. But the calculation results are very strange, the liquid-mixture is flowing to the higher-pressure-direction by the calculated results.
The initial velocity of liquid is 0, in momentum equation, the abs-value is (-gradient P) is much more than source item(rho*Volc*U). But the velocity of liquid is always positive and to the higher-pressure-direction.
Here are the codes in alphaEqn.H:
{

word alpha1Scheme("div(phi,alpha1)");
word alpha2Scheme("div(phi,alpha2)");
word alpha3Scheme("div(phi,alpha3)");



if (MULES1Corr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi1,
upwind<scalar>(mesh, phi1)
).fvmDiv(phi1, alpha1)
==
m1
);

alpha1Eqn.solve();

Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;

talphaPhi1 = alpha1Eqn.flux();
}

if (MULES2Corr)
{
fvScalarMatrix alpha2Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha2)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi2,
upwind<scalar>(mesh, phi2)
).fvmDiv(phi2, alpha2)
==
m2
);

alpha2Eqn.solve();

Info<< "Phase-2 volume fraction = "
<< alpha2.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha2.name() << ") = " << min(alpha2).value()
<< " Max(" << alpha2.name() << ") = " << max(alpha2).value()
<< endl;

talphaPhi2 = alpha2Eqn.flux();
}

if (MULES3Corr)
{
fvScalarMatrix alpha3Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha3)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi3,
upwind<scalar>(mesh, phi3)
).fvmDiv(phi3, alpha3)
==
m3
);

alpha3Eqn.solve();

Info<< "Phase-3 volume fraction = "
<< alpha3.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha3.name() << ") = " << min(alpha3).value()
<< " Max(" << alpha3.name() << ") = " << max(alpha3).value()
<< endl;

talphaPhi3 = alpha3Eqn.flux();

}
//alpha3 == 1- alpha1 - alpha2;
rhoPhi = talphaPhi1*rho1 + talphaPhi2*rho2 + talphaPhi3*rho3;
rho == alpha1*rho1 + alpha2*rho2 + alpha3*rho3;
rho1v == alpha1*rho1;
rho2v == alpha2*rho2;
rho3v == alpha3*rho3;
n1 == rho1v/rho;
n2 == rho2v/rho;
n3 == rho3v/rho;
mu == alpha1*mu1_ + alpha2*mu2_ + alpha3*mu3_;
nu == mu/rho;
w2 == ((rho2-rho)*pow(d1_,2)/(18*mu))*g;
w3 == ((rho3-rho)*pow(d2_,2)/(18*mu))*g;
m1 == (-pow(alpha1*rho1,2)/(rho*rho1)) *(K1_+K2_);
m2 == (pow(alpha1*rho1,2)/(rho*rho2)) *(K1_);
m3 == (pow(alpha1*rho1,2)/(rho*rho3)) *(K2_);
Volc == (m1 + m2 + m3);
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
Info<< "Phase-2 volume fraction = "
<< alpha2.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha2.name() << ") = " << min(alpha2).value()
<< " Max(" << alpha2.name() << ") = " << max(alpha2).value()
<< endl;
Info<< "Phase-3 volume fraction = "
<< alpha3.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha3.name() << ") = " << min(alpha3).value()
<< " Max(" << alpha3.name() << ") = " << max(alpha3).value()
<< endl;

}

And the codes in UEqn.H:
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
//- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
- fvm::laplacian(mu, U) + rho*Volc*U
//
);

UEqn.relax();

if (pimple.momentumPredictor())
{
//Info<< "Solving p\n" << endl;
//Info<< fvc::snGrad(p_rgh)<<endl;
Info<< " Min(bU) = " << min(U).value()<< " Max(bU) = " << max(U).value()<<endl;
solve
(UEqn ==
fvc::reconstruct
(
(
//interface.surfaceTensionForce()
//- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
Info<< " Min(U) = " << min(U).value()
<< " Max(U) = " << max(U).value()<<endl;
}