[trl0007]

I've got a simple 2-D test case where one mole bubble of H2 is allowed to diffuse into a surrounding mole of ambient Air. Obviously, because I care about mass diffusion exclusively, I have set gravity to zero. I was doing this simulation to compare to Fick's Law of Diffusion in polar coordinates, just to see how close an FDM model of Fick's Second Law agrees with OpenFoam. I'm using rhoReactingBuoyantFoam because in the future I plan on adding in the effects of gravity and ultimately I'm going to look at a combustion event. But for now, I'd like to see that mass diffusion can be handled correctly.

However, I am noticing that the temperature profile is not well behaved in my simulation. Everything starts at 298 K, but instantly after starting the simulation the temperature rises at the mixing boundary. After simulating for 20 seconds, the maximum temperature can be as high at ~330 K and the inner portions of the H2 bubble cool down to about ~230 K. I've attached pictures of these results. I'm using second order accurate schemes spatially as well as temporally and have checked for convergence with very fine meshes and the issue persists throughout.

I expected some small temperature fluctuations if the gas compresses at all, but nothing of this order of magnitude.

Sorry, I am very new to OpenFOAM, so any advice given would be much appreciated.

The script that I'm using that I anticipated would turn off the effects of combustion is below:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      combustionProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

combustionModel  noCombustion<rhoThermoCombustion>;

active  true;

noCombustionCoeffs
{
}


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

My fvSchemes are below:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         CrankNicolson 0.9;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,h)      Gauss linearUpwind grad(h);
    div(phi,e)      Gauss linearUpwind grad(e);
    div(phi,k)      Gauss linearUpwind grad(k);
    div(phi,epsilon) Gauss linearUpwind grad(epsilon);
    div(phi,R)      Gauss linearUpwind grad(R);
    div(phi,K)      Gauss linearUpwind grad(K);
    div(phi,Ekp)    Gauss linearUpwind grad(Ekp);
    div(phi,Yi_h)   Gauss limitedLinear 0.5;
    div(R)          Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.5;
}

fluxRequired
{
    default         no;
    p_rgh;
}

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

temperature.jpg

concentration.jpg

[Swagga5aur]

Hello Tyler,
This issue seems quite similar to that reported in this thread: having trouble using reactingFoam with reactions turned off
I would like to take a look at your case if you would attach the files.

[trl0007]

I had to reduce the mesh size substantially to meet the file size limits on attachments, but the rest of the files I've left the same as my last run.

Any help on this matter would be much appreciated. The temperature increase always occurs on the outer edge of the mixing boundary, and cooling occurs directly inside of the mixing boundary.

I've tried compiling my own solver and getting rid of the reaction term inside of EEqn.H and it didn't change this at all.

[Swagga5aur]

Thak you I'll look into it in the near future, quick question I tried to solve it just now and none of the mass fractions where specified so I am unable to run the case with rhoReactingBuoyantFoam as no species are specified in the domain.

What is the internalField consisting of when considering mass fractions?

[trl0007]

I appreciate the help! That's curious that it won't run for you, I just took those folders and ran setFields and then rhoReactingBuoyantFoam ran on my end. I create the bubble of H2 inside of the setFieldsDict file such that it's 100% H2 inside of the bubble and 100% air outside of it.

[Swagga5aur]

Ah sorry did not run setfields that explains it.

Sent from my A0001 using CFD Online Forum mobile app

[Swagga5aur]

Hello again Tyler,
I tried enforcing the alphaEff -> muEff, as suggested in the thread I mentioned earlier: having trouble using reactingFoam with reactions turned off

Additionally, noticed some inconsistencies in the boundary conditions which I altered such that the outlet are the atmosphere patches.

This resulted in a constant temperature field of 298K, and I have attached the new solver and case to this post.

Just compile the solver with wmake and run the case with the Allrun script, creating the blockMesh, setFields and solves with the new solver.

Let me know if you have any questions.

[trl0007]

Sure enough, that works perfectly well. I appreciate your help on the matter and have only a parting question:

If I understand this correctly, the solver works under the assumption that the Lewis number and the Schmidt numbers are both one (i.e. thermal diffusivity is the same as momentum diffusivity is the same as mass diffusivity)

This should introduce some manner of error into the solution, yes?

Either way, you've been a huge help. I'm still getting my feet wet with OpenFOAM in general, so this kind of thing aids with the learning curve.

[Swagga5aur]

No problem, and yes this approach is based on those assumptions, however, the original energy equation for the combustion model already assumes a Lewis number of 1 to be derived.
Lewis number of 1 is valid for small scale turbulence -> dominating turbulence transport, which is invalid for such cases as yours, and should be solved with diffusion based libraries.

Schmidt number of 1 corresponds to equal momentum and diffusion mass transfer is in the same order, which is applicable for the majority of gases.

Hope this makes sense else feel free to let me know.

[trl0007]

Yeah, that all makes sense. Many thanks again!

[Attavino]

Yes, the solver is valid only under the assumption the the Lewis number is equal to 1, you should have an error also in a laminar simulation. To avoid that error you should add a term related to the light gas diffusion in the energy equation.