CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   InterFoam: add a source term in alpha eq. (https://www.cfd-online.com/Forums/openfoam-programming-development/94845-interfoam-add-source-term-alpha-eq.html)

Alucard November 28, 2011 05:26

InterFoam: add a source term in alpha eq.
 
Goodmorning,

I'm try to simulate the growth of a cristal using,as a basis, interfoam.
As i'm at the beginning, i consider as a first case, a simple 2D circle that grows at a constant velocity "Uoo". Alpha=1 in the circle and 0 outside.

So what i want to "see" it's just a simple circle that grows in time in my square domain.
If the basic alpha transport equation is something like that:

d(alpha)/dt+U*grad(alpha)=0

in my case I add a "second" velocity that is the envelope growth:

d(alpha)/dt+U*grad(alpha)+Ug*grad(alpha)=0

So if I initialise my U=0 (I 've not motion of my circle) and I impose a constant Ug, how it's possible to integrate correctly this term "Ug*grad(alpha)" in the equation?

I tried to replace directly U=Ug BUT the calculation of the phi fluxes is quite complicated and realted to Ueqn so I cannot do that easily.
I know perhaps is a dummy question but it's 10 days I'm trying to figure out how to solve it and I don't see a clever solution to that.


Thank you in advance.

sabin.ceuca November 28, 2011 07:37

Hi,
if i were you i would take a look at the alphaEqn.H of interPhaseChangeFoam and just tailor it to suite your needs.
Although i don't understand your reasons behind it, because the alpha transport equation reads somehow different from what you have posted here.

Alucard November 28, 2011 08:24

Quote:

Originally Posted by sabin.ceuca (Post 333795)
Hi,
if i were you i would take a look at the alphaEqn.H of interPhaseChangeFoam and just tailor it to suite your needs.
Although i don't understand your reasons behind it, because the alpha transport equation reads somehow different from what you have posted here.

Hi,
thank you I'm looking at it right now I guess I've to work on the "Su" term (I don't know exactly how yet!
Why I'm using the interfoam approach?
Cause I have two phases in my study (liquid and solid) and So I need a "marker" for the phase.
So as I need a method that allows me to transport phases without numerical diffusion I thought to interFoam where all the "tricks" in order to have a sharp interface during time are already developped.
So i tried rigid translations and roatation of a circle with success and now I was trying to add the term related to surface growth.
In the complete problem it will be function of Temperature and other thermodynamical parameters but for the moment I impose it.
Why so do you think my approach is not the good one?
Thanks :)

Alucard November 28, 2011 08:39

So a good way to do that could be to write the "Su" term like that?

volScalarField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
fvc::grad(alpha1) & Ug

);


where Ug is the "growth velocity" vector field (I already initalised in my code)?

Thank you again.

olivierG November 28, 2011 08:41

Hello,

I think interFoam is not the best method in cristal grow, because usually the interface transition is very sharp, and with VOF method, you will get a "sharp but not enough" interface, + difficulties to set thermal source / density.

The best method may be Level set, but i don"t know how to deal with in OpenFoam (though i've seen a VOF/ level set coupling related to openFoam somewhere in a workshop pdf.)

Anyway, i would try to setup an hybrid buoyant+InterFoam, i.e thermal with variable properties and VOF, so if you go that way, we can share some part of the solver.

regards,

olivier

sabin.ceuca November 28, 2011 09:20

Hi Alucard,
yes your Su looks fine, give it a try. From the point of view of OF it looks good :), if it makes sense from a point of view of the physics you'll see.
Greets

Alucard November 28, 2011 12:02

Quote:

Originally Posted by sabin.ceuca (Post 333806)
Hi Alucard,
yes your Su looks fine, give it a try. From the point of view of OF it looks good :), if it makes sense from a point of view of the physics you'll see.
Greets

Thank you ,
it worked in OF (compilation+running) but not under a phisical point of view...I don't know why, i also tried to solve explicitly the equation:

" alpha1=alpha1-runTime.deltaT()*(Ug_mod*mag(fvc::grad(alpha1))); "
with Ug_mod, module of the velocity vector
but I'm having strong oscillations close to the circle-interface and alpha that becomes higher than 1.
Do you think something is wrong?

(before I tried to put Su as supposed and i got the same kinda of results)

sabin.ceuca November 29, 2011 04:10

Hmm,
don't worry it never works at the first try:) Did you also implement source and sink terms to the pEqn and UEqn? Because it might be that you create some 1st phase but do not transport it away.
Bon courage!
Sabin

akidess November 29, 2011 08:50

olivierG, the solid-liquid interface is sharp for pure materials, and diffuse for melting of alloys.

Nonetheless I agree interFoam is not the best choice for this problem. Instead, if not too late, I recommend having a look at phase-change solvers which other users have posted:
http://www.cfd-online.com/Forums/ope...g-problem.html

fabian_roesler January 13, 2012 10:16

melting solver
 
Hi all

nice that my melting solver

http://www.cfd-online.com/Forums/ope...g-problem.html

made it into another thread of this forum. :)
I although tried using interFoam for melting and solidification problems. But as you have a solid phase, alpha cannot be transported as done in VOF. My solver uses the enthalpy-porosity-method proposed by Voller. The solver will be improved in the next months.

Regards

Fabian

vak96 April 22, 2017 00:27

Urgent: Addition of source term to alpha equation in interFoam
 
Hi all,

I want to add a source term to the alpha equation in the interFoam solver. I am using OpenFOAM-v1612+ version. I tried to do it in the following way

I added following to alphaEqn.H
Code:

//Creating Source Terms for alpha equation
volScalarField::DimensionedField Su
    (
        IOobject
        (
            "Su",
            runTime.timeName(),
            mesh
        ),
        mesh,
        dimensionedScalar("Su", dimensionSet(0,0,-1,0,0,0,0), scalar(0.0))
    );
   
volScalarField::DimensionedField Sp
    (
        IOobject
        (
            "Sp",
            runTime.timeName(),
            mesh
        ),
        mesh,
    dimensionedScalar("Sp", dimensionSet(0,0, -1,0,0,0,0), 0.0)
    );

Then I added the source term to equation as follows
Code:

  fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
          ==
            Sp
        );

I gave values to this source term in loop like this in interFoam.C, Su is 0 everywhere
Code:

if( (sourceZ-0.05<=centroidZ) & (centroidZ < sourceZ))
                {
                    if(alpha1[cellI] < 0.5)
                    {
                        Sp[cellI] = 1;
                    }
                    else
                    {
                        Sp[cellI] = 0;
                    }
                }
            if( (sourceZ<centroidZ) & (centroidZ<=sourceZ+0.05))
                {
                    if(alpha1[cellI] < 0.5)
                    {
                        Sp[cellI] = -1;
                    }
                    else
                    {
                        Sp[cellI] = 0;
                    }
                }   
            else
                {
                    Sp[cellI] = 0;
                }

I don't understand how to change the MULES::correct and MULES::explicitSolve to include the source terms.

Can anyone please help with this? Whether the way I declared the source term is right or wrong? And the way I added it to the equation?

It is really urgent so please help me if possible.

Thanking you,
Akshay

Mahdi2010 August 3, 2017 02:23

Quote:

Originally Posted by vak96 (Post 645869)
Hi all,

I want to add a source term to the alpha equation in the interFoam solver. I am using OpenFOAM-v1612+ version. I tried to do it in the following way

I added following to alphaEqn.H
Code:

//Creating Source Terms for alpha equation
volScalarField::DimensionedField Su
    (
        IOobject
        (
            "Su",
            runTime.timeName(),
            mesh
        ),
        mesh,
        dimensionedScalar("Su", dimensionSet(0,0,-1,0,0,0,0), scalar(0.0))
    );
   
volScalarField::DimensionedField Sp
    (
        IOobject
        (
            "Sp",
            runTime.timeName(),
            mesh
        ),
        mesh,
    dimensionedScalar("Sp", dimensionSet(0,0, -1,0,0,0,0), 0.0)
    );

Then I added the source term to equation as follows
Code:

  fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
          ==
            Sp
        );

I gave values to this source term in loop like this in interFoam.C, Su is 0 everywhere
Code:

if( (sourceZ-0.05<=centroidZ) & (centroidZ < sourceZ))
                {
                    if(alpha1[cellI] < 0.5)
                    {
                        Sp[cellI] = 1;
                    }
                    else
                    {
                        Sp[cellI] = 0;
                    }
                }
            if( (sourceZ<centroidZ) & (centroidZ<=sourceZ+0.05))
                {
                    if(alpha1[cellI] < 0.5)
                    {
                        Sp[cellI] = -1;
                    }
                    else
                    {
                        Sp[cellI] = 0;
                    }
                }   
            else
                {
                    Sp[cellI] = 0;
                }

I don't understand how to change the MULES::correct and MULES::explicitSolve to include the source terms.

Can anyone please help with this? Whether the way I declared the source term is right or wrong? And the way I added it to the equation?

It is really urgent so please help me if possible.

Thanking you,
Akshay


I am not sure if you have already found the answer, but for others:
You would need to take a look at the source code of MULES as well as the how it is impelemented in interPhaseChangeFoam. Therefore, it seems you should just add:

MULES::correct(geometricOneField(), alpha1, talphaPhiUn(), talphaPhiCorr.ref(), Sp, Su, 1, 0);

and

MULES::explicitSolve(geometricOneField(), alpha1, phiCN, alphaPhi, Sp, Su, 1, 0);

bbita November 5, 2017 19:16

Dear Foamer,

I add a source term to capture interface mass transfer. It works fine except when I change the sign of the term, the result doesn't make any sense.
When transfer is from the reference phase, I can see how my phase shrinks but I cannot model the case when my reference phase is expanding.

Thanks


All times are GMT -4. The time now is 18:37.