[Quentin]

Hi Foamers,

I' m setting up (open FOAM v 2.3.0) a case for wave generation and drag calculation using LTSInterFoam. I conduct some basics verification and validation of my set up with a hull, for wich I had the drag values ; and I was able to get a 5% max difference with those value, with a visually coherent wave pattern.

Now I'm trying to apply my set-up wih a new hull - wich had quite the same bow as the mini 6.50 747 magnum. The case brake up as the alpha.water get a very wrong behaviour under the hull : outside the boundary layers, alpha seems quite correct but inside the boundary layer it goes to zero near the hull.

I've attached some picture with a "cylinder-hull" illustrating this, wich i plotted the interface, colored by z ; and the value of alpha.water on the "cylinder-hull". blue color is full water and red is full air.
1 is the initial situation.
1.jpg
2, 3 and 4 are ploted from the iteration 1600, we can see that under the "cylinder-hull" there is a water-air mix, when it should be full water.

3.jpg


2.jpg

the 4 picture show that the matter only occur on the boundary layer - when i put more boundary layers, only the first 4-5 boundary layer on the wall get those wrong alpha.water values.

4.jpg.


I also have upload my test case set up here.

vSolution :

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      0;
        nAlphaSubCycles 1;
        cAlpha          0;
        icAlpha         0;

        MULESCorr       yes;
        nLimiterIter    10;
        alphaApplyPrevCorr  yes;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
        minIter         1;
    }

    pcorr
    {
        solver          PCG;

        preconditioner
        {
            preconditioner  GAMG;

            smoother        GaussSeidel;
            agglomerator    faceAreaPair;
            mergeLevels     1;
            nCellsInCoarsestLevel 100;
            cacheAgglomeration true;

            tolerance       1e-6;
            relTol          0;
        };

        tolerance       1e-6;
        relTol          0;
    };

    p_rgh
    {
        solver          GAMG;

        smoother        DIC;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        nCellsInCoarsestLevel 100;
        cacheAgglomeration true;

        tolerance       1e-8;
        relTol          0.0;

    };

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "(U|k|omega).*"
    {
        solver          smoothSolver;

        smoother        symGaussSeidel;
        nSweeps         1;

        tolerance       1e-7;
        relTol          0;
        minIter         1;

    };
}

PIMPLE
{
    momentumPredictor   yes;

    nOuterCorrectors    1;
    nCorrectors         2;
    nNonOrthogonalCorrectors 0;

    maxCo               5;
    maxAlphaCo          5;

    rDeltaTSmoothingCoeff 0.05;
    rDeltaTDampingCoeff 0.5;
    nAlphaSpreadIter    0;
    nAlphaSweepIter     0;
    maxDeltaT           1;
}

relaxationFactors
{
    fields
    {
    p_rgh     0.3;
    

    }
    equations
    {
        "(U|k|omega).*" 0.7;
    "alpha.water.*" 0.7;
    }
}

fvScheme

ddtSchemes
{
    default         localEuler rDeltaT;
}

gradSchemes
{
    default         cellMDLimited Gauss linear 0.5;//Gauss linear;//
    limitedGrad        cellLimited Gauss linear 1;
}

divSchemes
{
    div(rhoPhi,U)   Gauss linearUpwind grad(U);
    div(phi,alpha)  Gauss limitedLinear 1.0 phi ;//Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression ;
    div(phi,k)      Gauss linearUpwind limitedGrad;
    div(phi,omega)  Gauss linearUpwind limitedGrad;
    div((muEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited 0.777;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 0.777;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha.water;
}

i have already try to change settings there and there, but as I'm unable to diagnose the matter, i've made only hazardous changes.

thank you in advance for your help,

Quentin Coispeau

ps : in my cylinder-test case, the mesh in not high enough in z+ to get the full wave at the bow, but i had the very same alpha problem with my "magnum-like" hull, for wich it was not the case.

[Quentin]

I found out that :

calpha can't be set to 0 : the compression speed used to compressed the interface is calculate as |U_c| = min( calpha * |U|, max |U| ). so if calpha is set to 0, you remove the convective compression flux and alpha equation becomes too diffusive. calpha should be equal to 1 or greater.

I also had to had a alpha corrector and set the alpha div scheme to vanLeer again. And reduce maxCo and maxAlphaCo

[laurent98]

hi Quentin,
did it solve all your pb? i just run the beginning of you case. i call this phenomena 'air under hull'. i got this behevior on a catamaran quiet flat hull. i could solve it in Finemarine with a threshold in pressure. i plan to make the same in OF, not done yet...
http://www.cfd-online.com/Forums/ope...-pressure.html
Laurent

[Quentin]

Hi Laurent,

Try to run the same case with these dictionnaries instead of the previous ones :

fvSolution

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      1;
        nAlphaSubCycles 1;
        cAlpha          1;
        icAlpha         0;

        MULESCorr       yes;
        nLimiterIter    10;
        alphaApplyPrevCorr  yes;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        minIter         1;
    }

    pcorr
    {
        solver          PCG;

        preconditioner
        {
            preconditioner  GAMG;

            smoother        GaussSeidel;
            agglomerator    faceAreaPair;
            mergeLevels     2;
            nCellsInCoarsestLevel 10;
            cacheAgglomeration false;
        nPreSweeps 0;
        nPostSweeps 2;
        nBottomSweeps 2;
            tolerance       1e-5;
            relTol          0;
        };

        tolerance       1e-5;
        relTol          0.01;
    maxIter 100;
    };

    p_rgh
    {
        solver          GAMG;

        smoother        GaussSeidel;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        nCellsInCoarsestLevel 10;
        cacheAgglomeration true;
    nPreSweeps 0;
    nPostSweeps 2;
    nBottomSweeps 2;
        tolerance       1e-5;
        relTol          0.05;

    };

    p_rghFinal
    {
        solver PCG;
    preconditioner
    {
        preconditioner GAMG;
        tolerance 1e-07;
        relTol 0;
        nVcycles 2;
        smoother GaussSeidel;
        nPreSweeps 0;
        nPostSweeps 2;
        nFinestSweeps 2;
        cacheAgglomeration on;
        nCellsInCoarsestLevel 10;
        agglomerator faceAreaPair;
        mergeLevels 1;
    }
    tolerance 1e-07;
    relTol 0;
    maxIter 20;
    }

    "(U|k|omega).*"
    {
        solver          smoothSolver;

        smoother        symGaussSeidel;
        nSweeps         1;

        tolerance       1e-7;
        relTol          0.1;
        minIter         1;

    };
}

PIMPLE
{
    momentumPredictor   yes;

    nOuterCorrectors    1;
    nCorrectors         2;
    nNonOrthogonalCorrectors 0;

    maxCo               0.2499;
    maxAlphaCo          0.2499;

    rDeltaTSmoothingCoeff 0.05;
    rDeltaTDampingCoeff 0.5;
    nAlphaSpreadIter    0;
    nAlphaSweepIter     1;
    maxDeltaT           1;
}

relaxationFactors
{
   /* fields
    {
    p_rgh     0.3;
    

    }
    equations
    {
        "(U|k|omega).*" 0.7;
    "alpha.water.*" 0.7;*/
    }
}

cache
{
    grad(U);
}

fvSchemes

ddtSchemes
{
    default         localEuler rDeltaT;
}

gradSchemes
{
    default         Gauss linear;//cellMDLimited Gauss linear 0.5;//
    limitedGrad        cellLimited Gauss linear 1;
}

divSchemes
{
    div(rhoPhi,U)   Gauss linearUpwind grad(U);
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression ;
    div(phi,k)      Gauss upwind limitedGrad;
    div(phi,omega)  Gauss upwind limitedGrad;
    div((muEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited 0.777;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 0.777;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha.water;
}

Using calpha = 1 permit to get a much sharper interface and quite solve the problem (behind the bow I still have the first layer of the boundary layer which got wrong values as in the pictures below after 6000 times step
6.jpg

7.jpg

But an other matter "the wiggle wave" appears as related in this thread. But as we need the boundary layer to resolve the boundary layer, we can't get a better aspect ratio or am I wrong ?)

If you look in fvSolution, i change maxC0 and maxAlphaCo from 5 to 0.2499, so I need 20 times more timeStep, which is quite longer ; i tried to add more alpha subcycles, but the interface becomes quite crazy after not so much timeStep.

you're idea of pressure threshold seems quite interesting to me, and should probably be implemented as a new boundary condition. I may give a try in it soon if I can't get satisfactory results without. But i fear it is a bit "cheating" with numerical behaviour, which may false the results - at least it needs validation.

I set icalpha to 0.25 because I've seen this value somewhere (don't remember where) but I don't know what's its use and haven't look through the code for it yet.

[Quentin]

Here an illustration of you're pressure treshold idea : the pressure and the interface on the cylinder with the blue/green jump at p = 250 Pa

8.jpg

so it seems possible to implement your pressure condition on alpha quite easely.

[laurent98]

hi,
so, let's say,
if p superior at 250 then put alpha=1
but i'am a very beginner in C++... is anyone have a idea of how to do that???
1- add a variable to declare in transportProperties P0=250
2-add a line in /home/laurent/OpenFOAM/OpenFOAM-2.3.0/applications/solvers/multiphase/interFoam/alphaEqn.H like
if (p > P0)
{
alpha.water=1
}


3- compile
could it be SO simple.... i'm dreaming....
any advice will be very welcome!
thank to all Laurent

[Quentin]

I think it would be much better to implement a new boundary condition which would be something like

if (p < P0)
     zeroGradient;
else 
     fixedValue 1;
end

I think I will do it today because I still have a one layer thick wrong alpha on my real hull simulation

[daveatstyacht]

Dear Quentin,
The behavior you are seeing is called "numerical ventilation". The forcing of the mass fraction alpha to be 1 when the pressure is greater than a certain value like can be done in FINE/Marine is something that should be used with caution. With low speed vessels where real ventilation is unlikely, this can be a way to improve your answer. For high speed craft you should confirm ventilation is numerical instead of real ventilation before applying such a correction. Also be careful about the value of pressure you pick being greater than the stagnation pressure of the air or you will create water in places you don't want it.

I should point out that you can correct the friction forces in post treatment by correcting the rho used. This is only appropriate in the case that wall functions were used. Both have their downsides.

The best is to avoid the issue by refining the mesh near the where the problem starts and minimizing sources of diffusion (schemes, courant number, etc). Easier said than done since they tend to increase the cost.

[Quentin]

Hi Dave,

thank you very much for your answer. knowing the name of the issue will help me find documentation about it.

I had already tried to have less diffusive numerical scheme but it wasn't enough. I'll follow your advice about refining the mesh where the problem appears.

I agree with you saying that the pressure threshold trick is more a bandage than a solution.

I should point out that you can correct the friction forces in post treatment by correcting the rho used. This is only appropriate in the case that wall functions were used. Both have their downsides.

I understand what you mean but i don't understand in which way i could correct rho and to correct what error ? hyper-ventilation ?

EDIT : at the bow, the prismatic boundar layers are oriented ~ 30-40° of the flow, which don't help with the diffusivity. Refining the mesh seems to limit the problem but can't avoid this angle difference between flow and layers... by the way is it correct that the central differencing scheme, i.e, Gauss linear in OF, is the less diffusive scheme (in fact with no numerical diffusion) ?

[laurent98]

hi Quentin,
did you manager to change the OF's code to do the trick, with the pressure? i also did a lots of tryings on better refinement on flat multihull but this trick was working perfectly well, once you know the pressure threshold (mostly depending of the speed ...) could you please help me to make the trick on OF? thanks LL

[Quentin]

Hi Laurent,

yes I manage to do the trick, using the inletOutlet boundary condition as a start. I'll give you my code but i'm not very sure of what I've done and it's definitively not very clean, but it seems to work ; so be careful.

It didn't wok in my case because I have both hyper ventilation of air under the hull and water at the bow ; so i'm working again with Dave's advices.

You have to recompile the inletOutlet chaging all "inletOutlet" for "alphaHull", and modifying the content of
1) alphaHullfvPatchField.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "alphaHullFvPatchField.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class Type>
Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF
)
:
    mixedFvPatchField<Type>(p, iF),
    phiName_("p"),
    pthreshold_(800.0)
{
    this->refValue() = pTraits<Type>::zero;
    this->refGrad() = pTraits<Type>::zero;
    this->valueFraction() = 0.0;
}


template<class Type>
Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField
(
    const alphaHullFvPatchField<Type>& ptf,
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    mixedFvPatchField<Type>(ptf, p, iF, mapper),
    phiName_(ptf.phiName_)
{}


template<class Type>
Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const dictionary& dict
)
:
    mixedFvPatchField<Type>(p, iF),
    phiName_(dict.lookupOrDefault<word>("phi", "p")),
    pthreshold_(dict.lookupOrDefault<scalar>("pthreshold", 800.0))
{
    this->refValue() = Field<Type>("inletValue", dict, p.size());

    fvPatchField<Type>::operator=(this->patchInternalField());

    this->refGrad() = pTraits<Type>::zero;
    
    this->valueFraction() = 1.0;
}


template<class Type>
Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField
(
    const alphaHullFvPatchField<Type>& ptf
)
:
    mixedFvPatchField<Type>(ptf),
    phiName_(ptf.phiName_)
{}


template<class Type>
Foam::alphaHullFvPatchField<Type>::alphaHullFvPatchField
(
    const alphaHullFvPatchField<Type>& ptf,
    const DimensionedField<Type, volMesh>& iF
)
:
    mixedFvPatchField<Type>(ptf, iF),
    phiName_(ptf.phiName_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type>
void Foam::alphaHullFvPatchField<Type>::updateCoeffs()
{
 
    if (this->updated())
    {
        return;
    }

    const Field<scalar>& phip =
    this->patch().template lookupPatchField<volScalarField, scalar>
    (
        "p"
    );
    
    this->valueFraction() = pos(phip - pthreshold_);

    mixedFvPatchField<Type>::updateCoeffs();
}


template<class Type>
void Foam::alphaHullFvPatchField<Type>::write(Ostream& os) const
{
    fvPatchField<Type>::write(os);
    if (phiName_ != "phi")
    {
        os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
    }
    this->refValue().writeEntry("inletValue", os);
    this->writeEntry("value", os);
}


// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
template<class Type>
void Foam::alphaHullFvPatchField<Type>::operator=
(
    const fvPatchField<Type>& ptf
)
{
    fvPatchField<Type>::operator=
    (
        this->valueFraction()*this->refValue()
        + (1 - this->valueFraction())*ptf
    );
}


// ************************************************************************* //

and
2)alphaHullFvPatchField.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::alphaHullFvPatchField

Group
    grpOutletBoundaryConditions

Description
    This boundary condition provides a generic outflow condition, with
    specified inflow for the case of return flow.

    \heading Patch usage

    \table
        Property     | Description             | Required    | Default value
        phi          | flux field name         | no          | phi
        alphaValue   | inlet value for reverse flow | yes    |
    \endtable

    Example of the boundary condition specification:
    \verbatim
    myPatch
    {
        type            alphaHull;
        phi             phi;
        alphaValue      uniform 0;
        value           uniform 0;
    }
    \endverbatim

    The mode of operation is determined by the sign of the flux across the
    patch faces.

Note
    Sign conventions:
    - positive flux (out of domain): apply zero-gradient condition
    - negative flux (into of domain): apply the user-specified fixed value

SeeAlso
    Foam::mixedFvPatchField
    Foam::zeroGradientFvPatchField

SourceFiles
    alphaHullFvPatchField.C

\*---------------------------------------------------------------------------*/

#ifndef alphaHullFvPatchField_H
#define alphaHullFvPatchField_H

#include "mixedFvPatchField.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                   Class alphaHullFvPatchField Declaration
\*---------------------------------------------------------------------------*/

template<class Type>
class alphaHullFvPatchField
:
    public mixedFvPatchField<Type>
{

protected:

    // Protected data

        //- Name of flux field
        word phiName_;

    //- Value of the pressure threshold
    scalar pthreshold_;


public:

    //- Runtime type information
    TypeName("alphaHull");


    // Constructors

        //- Construct from patch and internal field
        alphaHullFvPatchField
        (
            const fvPatch&,
            const DimensionedField<Type, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        alphaHullFvPatchField
        (
            const fvPatch&,
            const DimensionedField<Type, volMesh>&,
            const dictionary&
        );

        //- Construct by mappin"phi"g given alphaHullFvPatchField onto a new patch
        alphaHullFvPatchField
        (
            const alphaHullFvPatchField<Type>&,
            const fvPatch&,
            const DimensionedField<Type, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct as copy
        alphaHullFvPatchField
        (
            const alphaHullFvPatchField<Type>&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchField<Type> > clone() const
        {
            return tmp<fvPatchField<Type> >
            (
                new alphaHullFvPatchField<Type>(*this)
            );
        }

        //- Construct as copy setting internal field reference
        alphaHullFvPatchField
        (
            const alphaHullFvPatchField<Type>&,
            const DimensionedField<Type, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchField<Type> > clone
        (
            const DimensionedField<Type, volMesh>& iF
        ) const
        {
            return tmp<fvPatchField<Type> >
            (
                new alphaHullFvPatchField<Type>(*this, iF)
            );
        }


    // Member functions

        //- Update the coefficients associated with the patch field
        virtual void updateCoeffs();

        //- Write
        virtual void write(Ostream&) const;

    // Member operators

        virtual void operator=(const fvPatchField<Type>& pvf);

};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "alphaHullFvPatchField.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

then recompile and add your library to controlDict dictionnary and the call of the BC would be :

hull
    {
        type            alphaHull;
        inletValue    uniform 1.0;
    }

the important fonction is updatecoeff() in alphaHullFvPatchField.C. I hope it will help you but again my code is not clean at all, It's at test phase.

[Quentin]

Hi foamers, Dave,

With the knowledge of the name of the problem, i found his very intesting publication about resistance calculation for AC33 hulls. In paragraph 4.4 they are speaking about numerical ventilation, explaining how you can correct your skin resistance calculation with rho, but also telling that you could solve the problem adding a source term in the transport equations of ad \alpha_{air} wich would be something like with being an arbitrary distance from the wall and an arbitrary number of iterations:

• in the transport equation of the air phase : if && , then , else
• in the transport equation of the water phase : if && , then , else

As the interFoam solver only solve one transport equation for alpha, i.e :
with ;

I think I could rewrite the source terme as :

if && , then ,
else if && , then ,
else .

but i'm no quite sure of this as if i write the transport equation with a source term for :


this give me if i transform into the transport equation for :

so a gradient of U appears.

What do you think about it ?

EDIT : it's not a gradient, it's a divergence ...

[daveatstyacht]

Dear Quentin,
Yes, this paper was what I had in mind though again it is only good if you have wall functions.

Regarding source terms, I would suggest looking at the cavitatingFoam and/or interPhaseChangeFoam to see how they have handled sources in the alpha equation.

[Quentin]

Dear Dave, Laurent & foamers,

I managed to imlement this source term in LTSinterFoam solver.

1) the source I finally implemented is :
if &&
if &&
else

2) I implemented implicitely he negative part and explicitely the positive one, so the disctretised equation has the form (with Euler explicit time discretisation) :

if &&
if &&
usual form else.

I choose usually n = 200.

3) so Here is the new alphaEqn.H ; please note that i'm not completly sure of what i did in MULS::LTScorrect and MULES::LTSexplicitsolve :

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    // Standard face-flux compression coefficient
    surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));

    // Add the optional isotropic compression contribution
    if (icAlpha > 0)
    {
        phic *= (1.0 - icAlpha);
        phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
    }

    // Do not compress interface at non-coupled boundary faces
    // (inlets, outlets etc.)
    forAll(phic.boundaryField(), patchi)
    {
        fvsPatchScalarField& phicp = phic.boundaryField()[patchi];

        if (!phicp.coupled())
        {
            phicp == 0;
        }
    }

    //** Creating Source Terms
    volScalarField::DimensionedInternalField Su
    (
        IOobject
        (
            "Su",
            runTime.timeName(),
            mesh
        ),
    mesh,
    dimensionedScalar("Su", dimensionSet(0,0, -1,0,0,0,0), 0.0)
    
    );

    volScalarField::DimensionedField Sp 
    (
        IOobject
        (
            "Sp",
            runTime.timeName(),
            mesh
        ),
        mesh,
    dimensionedScalar("Sp", dimensionSet(0,0, -1,0,0,0,0), 0.0)
    );

    //** Evaluating Source Terms
    forAll (alpha1, celli)
    {
    if(y[celli] <= d && alpha1[celli] <= 0.5 && alpha1[celli] >= 0.0)
    {

         Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
    }
    else if(y[celli] <= d && alpha1[celli] >= 0.5 && alpha1[celli] <= 1.0)
    {
         Su[celli] += 1.0/n * rsubDeltaT.field()[celli] ;

         Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
    }
    }
    
    

    tmp<surfaceScalarField> tphiAlpha;

    if (MULESCorr)
    {
    //** Modified Alpha Eqn
        fvScalarMatrix alpha1Eqn
        (
            #ifdef LTSSOLVE
            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
            #else
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            #endif
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phi,
                upwind<scalar>(mesh, phi)
            ).fvmDiv(phi, alpha1)

        );

        solve(alpha1Eqn == Su + fvm::Sp(Sp,alpha1));

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

        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
        tphiAlpha = tmp<surfaceScalarField>
        (
            new surfaceScalarField(tphiAlphaUD())
        );

        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
        {
        //** evaluation of alpha_old for LTScorrect
        volScalarField alpha10(alpha1);
            Info<< "Applying the previous iteration compression flux" << endl;
            #ifdef LTSSOLVE
        //** modified LTScorrect - /!\ unsure
            MULES::LTScorrect
        (
        geometricOneField(), 
        alpha1, 
        tphiAlpha(),
        tphiAlphaCorr0(), 
        Sp, 
        (-Sp*alpha10)(), 
        1, 
        0
        );
            #else
            MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
            #endif

            tphiAlpha() += tphiAlphaCorr0();
        }

        // Cache the upwind-flux
        tphiAlphaCorr0 = tphiAlphaUD;

        alpha2 = 1.0 - alpha1;

        interface.correct();
    }

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        surfaceScalarField phir(phic*interface.nHatf());

        tmp<surfaceScalarField> tphiAlphaUn
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
               -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );
    
        if (MULESCorr)
        {
            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
            volScalarField alpha10(alpha1);

            #ifdef LTSSOLVE
        //** modified LTScorrect - /!\ unsure
            MULES::LTScorrect(
        geometricOneField(), 
        alpha1,
        tphiAlphaUn(), 
        tphiAlphaCorr(), 
        Sp,
        (-Sp*alpha10)(), 
        1, 
        0);
            #else
            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
            #endif

            // Under-relax the correction for all but the 1st corrector
            if (aCorr == 0)
            {
                tphiAlpha() += tphiAlphaCorr();
            }
            else
            {
                alpha1 = 0.5*alpha1 + 0.5*alpha10;
                tphiAlpha() += 0.5*tphiAlphaCorr();
            }
        }
        else
        {
            tphiAlpha = tphiAlphaUn;

            #ifdef LTSSOLVE
        //** modified explicitLTSSolve - /!\ unsure
            MULES::explicitLTSSolve(geometricOneField(), alpha1, phi, tphiAlpha(), Sp, Su, 1, 0);
            #else
            MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0);
            #endif
        }

        alpha2 = 1.0 - alpha1;

        interface.correct();
    }

    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;

    if (alphaApplyPrevCorr && MULESCorr)
    {
        tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0;
    }

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

I also added in createField.H

    //** Wall distance 
    Info<< "Calculating wall distance field" <<endl;
    volScalarField y = wallDist(mesh).y();

and in alphaControls.H, i added

//** Added Source Terms controls

scalar d 
(
    alphaControls.lookupOrDefault<scalar>("d", 0.0)
);

scalar n
(
    alphaControls.lookupOrDefault<scalar>("n", 1.0)
);

so the Source terms n and d can be set in fvSolution in alpha.water solver's parameters.

I havn't tested it with nalphasubcycles > 1, but otherwise, it seems to work for me.

[jameswilson620]

Originally Posted by Quentin

Dear Dave, Laurent & foamers,

I managed to imlement this source term in LTSinterFoam solver.

1) the source I finally implemented is :
if &&
if &&
else

2) I implemented implicitely he negative part and explicitely the positive one, so the disctretised equation has the form (with Euler explicit time discretisation) :

if &&
if &&
usual form else.

I choose usually n = 200.

3) so Here is the new alphaEqn.H ; please note that i'm not completly sure of what i did in MULS::LTScorrect and MULES::LTSexplicitsolve :

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    // Standard face-flux compression coefficient
    surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));

    // Add the optional isotropic compression contribution
    if (icAlpha > 0)
    {
        phic *= (1.0 - icAlpha);
        phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
    }

    // Do not compress interface at non-coupled boundary faces
    // (inlets, outlets etc.)
    forAll(phic.boundaryField(), patchi)
    {
        fvsPatchScalarField& phicp = phic.boundaryField()[patchi];

        if (!phicp.coupled())
        {
            phicp == 0;
        }
    }

    //** Creating Source Terms
    volScalarField::DimensionedInternalField Su
    (
        IOobject
        (
            "Su",
            runTime.timeName(),
            mesh
        ),
    mesh,
    dimensionedScalar("Su", dimensionSet(0,0, -1,0,0,0,0), 0.0)
    
    );

    volScalarField::DimensionedField Sp 
    (
        IOobject
        (
            "Sp",
            runTime.timeName(),
            mesh
        ),
        mesh,
    dimensionedScalar("Sp", dimensionSet(0,0, -1,0,0,0,0), 0.0)
    );

    //** Evaluating Source Terms
    forAll (alpha1, celli)
    {
    if(y[celli] <= d && alpha1[celli] <= 0.5 && alpha1[celli] >= 0.0)
    {

         Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
    }
    else if(y[celli] <= d && alpha1[celli] >= 0.5 && alpha1[celli] <= 1.0)
    {
         Su[celli] += 1.0/n * rsubDeltaT.field()[celli] ;

         Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
    }
    }
    
    

    tmp<surfaceScalarField> tphiAlpha;

    if (MULESCorr)
    {
    //** Modified Alpha Eqn
        fvScalarMatrix alpha1Eqn
        (
            #ifdef LTSSOLVE
            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
            #else
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            #endif
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phi,
                upwind<scalar>(mesh, phi)
            ).fvmDiv(phi, alpha1)

        );

        solve(alpha1Eqn == Su + fvm::Sp(Sp,alpha1));

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

        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
        tphiAlpha = tmp<surfaceScalarField>
        (
            new surfaceScalarField(tphiAlphaUD())
        );

        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
        {
        //** evaluation of alpha_old for LTScorrect
        volScalarField alpha10(alpha1);
            Info<< "Applying the previous iteration compression flux" << endl;
            #ifdef LTSSOLVE
        //** modified LTScorrect - /!\ unsure
            MULES::LTScorrect
        (
        geometricOneField(), 
        alpha1, 
        tphiAlpha(),
        tphiAlphaCorr0(), 
        Sp, 
        (-Sp*alpha10)(), 
        1, 
        0
        );
            #else
            MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
            #endif

            tphiAlpha() += tphiAlphaCorr0();
        }

        // Cache the upwind-flux
        tphiAlphaCorr0 = tphiAlphaUD;

        alpha2 = 1.0 - alpha1;

        interface.correct();
    }

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        surfaceScalarField phir(phic*interface.nHatf());

        tmp<surfaceScalarField> tphiAlphaUn
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
               -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );
    
        if (MULESCorr)
        {
            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
            volScalarField alpha10(alpha1);

            #ifdef LTSSOLVE
        //** modified LTScorrect - /!\ unsure
            MULES::LTScorrect(
        geometricOneField(), 
        alpha1,
        tphiAlphaUn(), 
        tphiAlphaCorr(), 
        Sp,
        (-Sp*alpha10)(), 
        1, 
        0);
            #else
            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
            #endif

            // Under-relax the correction for all but the 1st corrector
            if (aCorr == 0)
            {
                tphiAlpha() += tphiAlphaCorr();
            }
            else
            {
                alpha1 = 0.5*alpha1 + 0.5*alpha10;
                tphiAlpha() += 0.5*tphiAlphaCorr();
            }
        }
        else
        {
            tphiAlpha = tphiAlphaUn;

            #ifdef LTSSOLVE
        //** modified explicitLTSSolve - /!\ unsure
            MULES::explicitLTSSolve(geometricOneField(), alpha1, phi, tphiAlpha(), Sp, Su, 1, 0);
            #else
            MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0);
            #endif
        }

        alpha2 = 1.0 - alpha1;

        interface.correct();
    }

    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;

    if (alphaApplyPrevCorr && MULESCorr)
    {
        tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0;
    }

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

I also added in createField.H

    //** Wall distance 
    Info<< "Calculating wall distance field" <<endl;
    volScalarField y = wallDist(mesh).y();

and in alphaControls.H, i added

//** Added Source Terms controls

scalar d 
(
    alphaControls.lookupOrDefault<scalar>("d", 0.0)
);

scalar n
(
    alphaControls.lookupOrDefault<scalar>("n", 1.0)
);

so the Source terms n and d can be set in fvSolution in alpha.water solver's parameters.

I havn't tested it with nalphasubcycles > 1, but otherwise, it seems to work for me.

Quentin,

Thanks for the info. With little knowledge of OpenFOAM and the discretization of these equations, Implementing source terms is quite the challenge. I have a few questions.

Using these conditions:

"
1) the source I finally implemented is :
if &&
if &&
else
"

I see this code:

 //** Evaluating Source Terms
    forAll (alpha1, celli)
    {
    if(y[celli] <= d && alpha1[celli] <= 0.5 && alpha1[celli] >= 0.0)
    {

         Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
    }
    else if(y[celli] <= d && alpha1[celli] >= 0.5 && alpha1[celli] <= 1.0)
    {
         Su[celli] += 1.0/n * rsubDeltaT.field()[celli] ;

         Sp[celli] -=  1.0/n * rsubDeltaT.field()[celli] ;
    }
    }

What is rsubDeltaT.field()? can you explain the input and output of this field, how its defined or what values and units it takes?

Could you explain for both conditions above (alpha <.5 and alpha > .5) why Sp is defined in both conditions and Su only defined in the latter?

Thank you so much for your help, James

[Quentin]

Hi James,

I havn't work on this for a few month so I don't remember exactly but :

Su is the part of the source term not linked to [MATH]\alpha[MATH\] and Sp is the one linked to.
It allows you to implement the source term explicitely or implicitely (Su => explicit ; Sp => implicit). there are some conditions regarding the stability of the solution whether the implicit part is positiv or negativ ( I don't remember exactly).

I would add that using this source term seems very dangerous to me because it brakes the physic and it was much more correct to activate the compression term already implemented in OpenFOAM through the parameter calpha.

in my case I finally used for the alpha solver :

calpha = 1;
nalphacorr = 0 for the first iterations (before the pression equation was stabilized (~ 50-100iterations) ) then 1
nalphasubcycle = 1

and i wasn't using a PIMPLE loop (noutercorrector = 1)

hopfully it will help you,

Quentin

[Quentin]

rsubdeltaT is the field containing the inverse of local time step (you have rdeltaT which is the inverse of main local time step and r subdeltaT is rDeltaT times bnalphasubcycle).

Regarding Su and Sp, in the "> 0.5" condition, the first part is implemented explicitely (Su) and the second implicitely (Sp). For the "<0.5" one, there is only an implicit part Sp and Su = 0.