[ZX-zha]

I know for sure that OpenFOAM can be used in 1D calculation. So I want to simulate the 1D two-phase pipe flow using the multiphaseEulerFoam solver, but some modifications are also needed. So I ignored the lift model, the virtual mass model often used in the 3D simulation and developed a drag model used in 1D simulation, the basic equations should be:


where is the depth of water.

So I thought the only thing I needed to do was write a drag model and add the water depth term. I also found a paper that uses the same idea as mine.

The K() member function of the drag model is:

Foam::tmp<Foam::volScalarField> Foam::dragModels::PipeLayer::K() const
{
    const fvMesh& mesh(interface_.phase1().mesh());
    const volScalarField& alpha2(interface_.phase2());

    const volScalarField& rho1(interface_.phase1().rho());

    tmp<volScalarField> tnu1(interface_.phase1().thermo().nu());

    const volScalarField& nu1(tnu1());

    Info<< "Create the Circle Object"<< endl;
    Circle cir(d_);

    wordList pbt = alpha2.boundaryField().types();
    forAll(pbt,pbi)
    {
        if (pbt[pbi] != "empty")
        {
            pbt[pbi] = "zeroGradient";
        }
    }

    volScalarField relativeHeight_
    (
        IOobject
        (
            "calc",
            mesh.time().timeName(),
            mesh
        ),
        mesh,
        dimensionedScalar("zero",dimLength,0.0),
        pbt
    );
    relativeHeight_.primitiveFieldRef() = cir.HeightRef(alpha2);
    relativeHeight_.correctBoundaryConditions();

    volScalarField hydroDiameter(cir.hydroDiameterG(alpha2,relativeHeight_));

    volScalarField Re1(Foam::max(hydroDiameter*interface_.magUr()/nu1,Foam::SMALL));

    volScalarField f1(Foam::max(16/Re1,0.014));
    volScalarField AG = Foam::max(cir.AG(alpha2)/pow(d_,2),Foam::SMALL)*pow(d_,2);

    return 0.5*f1*rho1*hydroDiameter*interface_.magUr()/AG;
}

And the water depth term(in UEqn.H file):

        UEqns.set
        (
            phase.index(),
            new fvVectorMatrix
            (
                // fvm::ddt(alpha, rho, U)
                // + fvm::div(alphaRhoPhi, U)
                // + fvm::SuSp(-this->continuityError(), U)
                phase.UEqn()
                ==
                *momentumTransfer[phase.name()]
                + fvModels.source(alpha, rho, U)
                - cosTheta*alpha*rho*Foam::mag(g)*fvc::grad(height)
            )
        );

But for me, it didn't work.

I know there is something wrong, but I can't find the mistake(s). It would be very helpful if anyone could take a look.

[ZX-zha]

I also took some tests:
Test1(without source term):
in this paper:Study of flow at the air-Water interface

The result shows:


Test2(without source term):
in this paper:Hybrid flux-splitting schemes for a common two-fluid model

The result shows:

(The left is the paper result and the right is my result)

These tests show that this idea is feasible, so either my drag model or the water depth term is wrong.