Adding diffusion term to interFoam transport equation

Gearb0x · August 30, 2017, 06:01

Dear Foamers,

I would like to add a diffusion term to the interFoam transport equation. Indeed, I want to simulate an air-hydrogen mixture and diffusion due to concentration gradient is normally significant for such mixture.

So, here is my thought process of what to do. Since I have never modified an OpenFoam solver before, I want to be sure that thatI am understanding each step and that I am doing things correctly.

So the current transport equation for $\alpha_1$ (one of the two phases/components) equation in interFoam is found in alphaEqn.H and equals to :

$\frac{\partial \alpha_1}{\partial t} + \nabla \cdot (\alpha_1 \bold{U}) + \nabla\cdot (\alpha_1 \alpha_2 \bold{U_r})$

and I understand that in this equation $\rho$ has been simplified because it is a solver for incompressible fluids.

I would like to add the diffusion due to concentration. Based on the generic transport equation :

$\frac{\partial \rho \alpha_1}{\partial t} + \nabla \cdot (\rho \alpha_1 \bold{U}) = \nabla \cdot (\Gamma \nabla \alpha_1) + S_{\alpha_1}$

So I understand that for interFoam solver, I must somehow include the following term, assuming diffusion term $\Gamma$ is constant and remembering to divide by $\rho$ due to its simplification in interFoam equation sinceit is for incompressible fluids

$\frac{1}{\rho}\nabla \cdot (\Gamma_{\alpha1} \nabla \alpha_1) = \frac{\Gamma_{\alpha1}}{\rho} \nabla \cdot (\nabla \alpha_1) = \frac{\Gamma_{\alpha1}}{\rho} \Delta \alpha_1)$

I copy/pasted the interFoam solver directory to a personal directory and I understand that the code for [/math] \alpha_1 [math] transport equation is found in alphaEqn.H and I assume I should do the following modification :

Code:

    
if (MULESCorr)
    {
        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)
           
          - fvm:: laplacian(Gamma1/rho1, alpha1)

        );

based on the code I found in scalarTransportFoam.C.

Now in scalarTransportFoam.C:

it is not fvm:: and not fv:: which is different than the current functions implemented in alphaEqn.H
sclarTransportFoam is a .C file and not a .H file

it is not fv:: but fvm::

Unfortunately, I do not master C language sufficiently nor OpenFoam structure/classes/functions to understand such subtleties and I would like to be sure I do the right thing and that I understand what I am doing.

Also I see that for example, the function gaussConvectionScheme<scalar> requires additional parameters to work

Code:

   fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)

and I don't know if such things are necessary for the function laplacian.

Eventually I see that

Code:

 
LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)

must be some kind of management of variable types. I am not sure I have to manage this thing also for my added code sample ?

Once alphaEqn.H is fine, I have to add the new constant

to CreateFileds.H

Code:

Info<< "Reading transportProperties\n" << endl;
immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);

volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());

const dimensionedScalar& rho1 = mixture.rho1();
const dimensionedScalar& rho2 = mixture.rho2();
const dimensionedScalar& Gamma1= mixture.Gamma1();

Thanks in advance for confirmation and support.

Once those changes have been confirmed to be fine, I will pursue with the modifications of the OpenFoam dictionnaries as instructed in this tutorial

Dipsomaniac · August 30, 2017, 06:37

Dear GearB0x

This seems like a case where interFoam, which uses interface compression and the MULES algorithm to decrease the effects of numerical diffusion at the interface between two immiscible phases, and probably will not be the ideal solver for the mixing of your two gases. Have you considered looking at a two phase Euler-Euler solver like twoPhaseEulerFoam or reactingTwoPhaseEulerFoam?

Regards,
Brian

Gearb0x · August 30, 2017, 08:51

I want to use interFoam as a first step because it is a simpler solver i.e. incompressible, no heat transfer.

Beans8 · February 14, 2023, 05:16

Did you manage to add the diffusion term?

August 30, 2017, 06:01	Adding diffusion term to interFoam transport equation	#1
Gearb0x Senior Member Join Date: Nov 2009 Posts: 111 Rep Power: 16	Dear Foamers, I would like to add a diffusion term to the interFoam transport equation. Indeed, I want to simulate an air-hydrogen mixture and diffusion due to concentration gradient is normally significant for such mixture. So, here is my thought process of what to do. Since I have never modified an OpenFoam solver before, I want to be sure that thatI am understanding each step and that I am doing things correctly. So the current transport equation for $\alpha_1$ (one of the two phases/components) equation in interFoam is found in alphaEqn.H and equals to : $\frac{\partial \alpha_1}{\partial t} + \nabla \cdot (\alpha_1 \bold{U}) + \nabla\cdot (\alpha_1 \alpha_2 \bold{U_r})$ and I understand that in this equation $\rho$ has been simplified because it is a solver for incompressible fluids. I would like to add the diffusion due to concentration. Based on the generic transport equation : $\frac{\partial \rho \alpha_1}{\partial t} + \nabla \cdot (\rho \alpha_1 \bold{U}) = \nabla \cdot (\Gamma \nabla \alpha_1) + S_{\alpha_1}$ So I understand that for interFoam solver, I must somehow include the following term, assuming diffusion term $\Gamma$ is constant and remembering to divide by $\rho$ due to its simplification in interFoam equation sinceit is for incompressible fluids $\frac{1}{\rho}\nabla \cdot (\Gamma_{\alpha1} \nabla \alpha_1) = \frac{\Gamma_{\alpha1}}{\rho} \nabla \cdot (\nabla \alpha_1) = \frac{\Gamma_{\alpha1}}{\rho} \Delta \alpha_1)$ I copy/pasted the interFoam solver directory to a personal directory and I understand that the code for [/math] \alpha_1 [math] transport equation is found in alphaEqn.H and I assume I should do the following modification : Code: if (MULESCorr) { 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) - fvm:: laplacian(Gamma1/rho1, alpha1) ); based on the code I found in scalarTransportFoam.C. Now in scalarTransportFoam.C: it is not fvm:: and not fv:: which is different than the current functions implemented in alphaEqn.H sclarTransportFoam is a .C file and not a .H file it is not fv:: but fvm:: Unfortunately, I do not master C language sufficiently nor OpenFoam structure/classes/functions to understand such subtleties and I would like to be sure I do the right thing and that I understand what I am doing. Also I see that for example, the function gaussConvectionScheme<scalar> requires additional parameters to work Code: fv::gaussConvectionScheme<scalar> ( mesh, phiCN, upwind<scalar>(mesh, phiCN) ).fvmDiv(phiCN, alpha1) and I don't know if such things are necessary for the function laplacian. Eventually I see that Code: LTS ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) must be some kind of management of variable types. I am not sure I have to manage this thing also for my added code sample ? Once alphaEqn.H is fine, I have to add the new constant to CreateFileds.H Code: Info<< "Reading transportProperties\n" << endl; immiscibleIncompressibleTwoPhaseMixture mixture(U, phi); volScalarField& alpha1(mixture.alpha1()); volScalarField& alpha2(mixture.alpha2()); const dimensionedScalar& rho1 = mixture.rho1(); const dimensionedScalar& rho2 = mixture.rho2(); const dimensionedScalar& Gamma1= mixture.Gamma1(); Thanks in advance for confirmation and support. Once those changes have been confirmed to be fine, I will pursue with the modifications of the OpenFoam dictionnaries as instructed in this tutorial tonnykz, the_ichthyologist and erinsam like this. Last edited by Gearb0x; August 30, 2017 at 08:52.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Source Term due to evaporation in energy transport equation	styleworker	OpenFOAM Programming & Development	3	September 7, 2022 04:09
adding a constant volumetric source term to transport equation in a particular region	cfdonline2mohsen	OpenFOAM Programming & Development	15	February 16, 2017 10:55
A simple question about adding Div and Grad of a source term to the momentum equation	adkar	Main CFD Forum	5	May 23, 2016 18:31
add source term to transport equation	m003020002	Fluent UDF and Scheme Programming	0	September 10, 2013 04:09
Adding a source term to concervation of species equation	bernarde	STAR-CCM+	0	June 4, 2009 04:17

August 30, 2017, 06:37		#2
Dipsomaniac Member Brian Willis Join Date: Mar 2011 Location: Cape Town, South Africa Posts: 58 Rep Power: 15	Dear GearB0x This seems like a case where interFoam, which uses interface compression and the MULES algorithm to decrease the effects of numerical diffusion at the interface between two immiscible phases, and probably will not be the ideal solver for the mixing of your two gases. Have you considered looking at a two phase Euler-Euler solver like twoPhaseEulerFoam or reactingTwoPhaseEulerFoam? Regards, Brian

August 30, 2017, 08:51		#3
Gearb0x Senior Member Join Date: Nov 2009 Posts: 111 Rep Power: 16	I want to use interFoam as a first step because it is a simpler solver i.e. incompressible, no heat transfer.

February 14, 2023, 05:16		#4
Beans8 New Member Join Date: Oct 2022 Posts: 24 Rep Power: 4	Did you manage to add the diffusion term?