CFD Online Discussion Forums - alpha1 wrong behaviour with LTSInterFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - alpha1 wrong behaviour with LTSInterFoam (https://www.cfd-online.com/Forums/openfoam-solving/138324-alpha1-wrong-behaviour-ltsinterfoam.html)

alpha1 wrong behaviour with LTSInterFoam [solved]

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.
Attachment 32069
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.

Attachment 32071

Attachment 32070

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.

Attachment 32072.

I also have upload my test case set up here.

vSolution :

Code:

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

Code:

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.

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

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

Hi Laurent,

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

fvSolution

Code:

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

Code:

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
Attachment 32105

Attachment 32106

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.

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

Attachment 32107

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

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

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

Code:

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

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.

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.

Quote:

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) ?

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

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

Code:

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

  =========                 |

  \\      /  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

Code:

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

  =========                 |

  \\      /  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 :

Code:

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.

Adding source term to avoid numerical ventilation

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 $\alpha_{water}$ 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 $\alpha_{water} > 0.5$ && , then $S = -\frac{\alpha_{air}}{n}$ , else
in the transport equation of the water phase : if $\alpha_{air} > 0.5$ && , then $S = -\frac{\alpha_{water}}{n}$ , else

As the interFoam solver only solve one transport equation for alpha, i.e :
$\partial_t \alpha + \nabla \cdot ( U \alpha) + \nabla \cdot (U_r(1-\alpha)\alpha) = S$ with

;

I think I could rewrite the source terme as :

if $\alpha_{water} > 0.5$ &&

, then $S = - \frac{1 - \alpha_{water}}{n}$ ,
else if $\alpha_{water} < 0.5$ &&

, then $S = -\frac{\alpha_{water}}{n}$ ,
else

.

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

$\partial_t (1 - \alpha ) + \nabla \cdot ( U (1-\alpha)) = S$
this give me if i transform into the transport equation for $\alpha$ :
$\partial_t (\alpha ) + \nabla \cdot ( U \alpha) = - S + \nabla \cdot U$
so a gradient of U appears.

What do you think about it ?

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

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.

Dear Dave, Laurent & foamers,

I managed to imlement this source term in LTSinterFoam solver.

1) the source I finally implemented is :
$S = - \frac{\alpha}{n \, \Delta T}$ if $0 < \alpha < 0.5$ &&

$S =\frac{1 - \alpha }{n \, \Delta T}$ if $0.5 < \alpha < 1$ &&

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) :

$\alpha^n = \frac{ n}{1+n} \Big( \alpha^o + \Delta T \textrm{fvm::div(U,alpha)} \Big)$ if $0 < \alpha < 0.5$ &&

$\alpha^n = \frac{ n}{1+n} \Big( \alpha^o + \frac{1}{n} + \Delta T \textrm{fvm::div(U,alpha)} \Big)$ if $0.5 < \alpha < 1.0$ &&

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 :

Code:

{

    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

HTML Code:

    //** Wall distance 

    Info<< "Calculating wall distance field" <<endl;

    volScalarField y = wallDist(mesh).y();

and in alphaControls.H, i added

Code:

//** 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.

Quote:

Originally Posted by Quentin (Post 503703)

Dear Dave, Laurent & foamers,

I managed to imlement this source term in LTSinterFoam solver.

1) the source I finally implemented is :
$S = - \frac{\alpha}{n \, \Delta T}$ if $0 < \alpha < 0.5$ &&

$S =\frac{1 - \alpha }{n \, \Delta T}$ if $0.5 < \alpha < 1$ &&

else

$\alpha^n = \frac{ n}{1+n} \Big( \alpha^o + \frac{1}{n} + \Delta T \textrm{fvm::div(U,alpha)} \Big)$ if $0.5 < \alpha < 1.0$ &&

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 :

Code:

{

    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

HTML Code:

    //** Wall distance 

    Info<< "Calculating wall distance field" <<endl;

    volScalarField y = wallDist(mesh).y();

and in alphaControls.H, i added

Code:

//** 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 :
$S = - \frac{\alpha}{n \, \Delta T}$ if $0 < \alpha < 0.5$ &&

$S =\frac{1 - \alpha }{n \, \Delta T}$ if $0.5 < \alpha < 1$ &&

else

"

I see this code:

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

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

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.