[finn_amann]

Hello friends,

I found a momentum source code based on meanVelocityForce in Lee et al (2021) (Code written by Y. C. Xu, an if-condition added by me).

In the scenario of a two-phase open-channel flow with water and air, you could imagine that you might want to limit a momentum source to the water phase. With the basic meanVelocityForce, this proves to be difficult, as one can only select a static set of cells, which is insufficient for a freely deforming water surface.

The new code, vofMeanVelocityForce, simply multiplies certain variables with the alpha coefficient that represents the two phases (alpha = 1 -> water --- alpha = 0 -> air).
The three important functions (magUBarAve(...), correct(...), addSup(...)) in vofMeanVelocityForce.C look like this:

Foam::scalar Foam::fv::vofMeanVelocityForce::magUbarAve
(
    const volVectorField& U
) const
{
    // Read IB cell information

     // Here we get a hold of the alpha field.
    volScalarField& alpha = const_cast<volScalarField&>
        (mesh_.lookupObject<volScalarField>("alpha.water"));

    scalar Vtotal = 0.0;
    scalar magUbarAve = 0.0;

    const scalarField& cv = mesh_.V();

    forAll(cells_, i)
    {
        label cellI = cells_[i];
        scalar volCell = cv[cellI];
        // Only account for cells that are filled with water.
        Vtotal += volCell*alpha[cellI]; 
        magUbarAve += (flowDir_ & U[cellI])*volCell*alpha[cellI]; 
    }

    reduce(magUbarAve, sumOp<scalar>());
    reduce(Vtotal, sumOp<scalar>());

    magUbarAve /= Vtotal;

    return magUbarAve;
}

In interFoam, this function is called as fvOptions.correct(U) in the final iteration of the non-orthogonal correction in pEqn.H. It is also called in the momentumPredictor condition in UEqn.H, but I work without the momentum predictor.

 void Foam::fv::vofMeanVelocityForce::correct(volVectorField& U)
{
    // Read IB cell information

    volScalarField& alpha = const_cast<volScalarField&>
        (mesh_.lookupObject<volScalarField>("alpha.water"));
    scalar Vtotal = 0.0;

    const scalarField& rAU = rAPtr_().internalField();

    // Integrate flow variables over cell set
    scalar rAUave = 0.0;
    const scalarField& cv = mesh_.V();
    forAll(cells_, i)
    {
        label cellI = cells_[i];
        scalar volCell = cv[cellI];
        // Again, only account for cells that are filled with water.
        Vtotal += volCell*alpha[cellI];
        rAUave += rAU[cellI]*volCell*alpha[cellI];
    }

    // Collect across all processors
    reduce(rAUave, sumOp<scalar>());
    reduce(Vtotal, sumOp<scalar>());

    // Volume averages
    rAUave /= Vtotal;

    scalar magUbarAve = this->magUbarAve(U);

    // Calculate the pressure gradient increment needed to adjust the average
    // flow-rate to the desired value
    dGradP_ = relaxation_*(mag(Ubar_) - magUbarAve)/rAUave;

    // Apply correction to velocity field
    forAll(cells_, i)
    {
        label cellI = cells_[i];
        // Only apply the change to cells that are filled with water (the if-condition was added by me):
        if (alpha[cellI] > 0.9) {
            U[cellI] += flowDir_*rAU[cellI]*dGradP_;
        }
    }
    scalar gradP = gradP0_ + dGradP_;

    Info<< "Pressure gradient source: uncorrected Ubar = " << magUbarAve
        << ", pressure gradient = " << gradP << endl;

    writeProps(gradP);
}

In interFoam, as far as I understand, the function fvOptions(rho, U), which is called on the right-hand side of UEqn in UEqn.H, is actually calling the following addSup(...) and applies a pressure gradient to the cells.

 void Foam::fv::vofMeanVelocityForce::addSup
(
    fvMatrix<vector>& eqn,
    const label fieldI
)
{
    volVectorField::Internal Su
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "Sup",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedVector(eqn.dimensions()/dimVolume, Zero)
    );

    volScalarField& alpha = const_cast<volScalarField&>
        (mesh_.lookupObject<volScalarField>("alpha.water"));

    scalar gradP = gradP0_ + dGradP_;
    scalarField alpha_d (alpha, cells_);
    forAll(alpha_d, cellI)
     { 
        // set alpha to 1 in cells with water and 0 in cells with air.
        if(alpha[cellI]>0.5)
        {
            alpha_d[cellI] = 1.0;
        }
        else
        {
            alpha_d[cellI] = 0;
        }
    }

    // Only add the pressure gradient to cells with water, otherwise add zero.
    UIndirectList<vector>(Su, cells_) = flowDir_*gradP*alpha_d;

    eqn += Su;
}

Now on to the problem I am encountering with this momentum source.

I created a 3D domain with cyclic sides. interFoam is used to simulate flow in x-direction.

If I apply the source to a cellSet, let's say a box containing the water surface, i.e. cells with water and cells with air, there will be some unexplained velocities in the air phase of that cellSet. I expected that those cells would not be affected by the source term since the it is always multiplied by zero if it contains air.

If I modify the cellSet to only contain the water cells (by shortening the height of the box to coincide with the water surface for example), the source seems to work fine and there are no weird velocities in the air phase.
I added two images of slices that cut the domain and the cellSet to illustrate this a little bit. Images were taken after 0.1 s of runtime. Initial velocity everywhere is 0.1 m/s and the source term should keep the velocity inside the cellset to 0.1 m/s.

In the image of the velocity, the cells in the air phase of the cellSet develop curious velocities, while the rest behave the way I expect. As I said, if I limit the cellSet to just the water phase, the cells do not exhibit these velocities.

So what am I missing here? How can I make sure, that the source term is only affecting the cells that are filled with water? Why does it not suffice to add "nothing" to the UEqn and U for cells in the air phase?

Have a nice weekend!


Finn