[KTG]

Hi everyone,



I am trying to implement a steady state advection diffusion equation. I see that the transient scalarTransport function object can estimate the turbulent diffusivity from the turbulence model, but I would like to know how. I did look in the source code, and could not figure it out. Does anyone know or have a reference for that? I am just interested in the isothermal case at this point.



Thanks

[Tobermory]

Let me give you a few pointers. Start with line 40 of the header file (https://cpp.openfoam.org/v8/scalarTr...8H_source.html) which shows how the diffusivity is built from the laminar and turbulent diffusivities, based on alpha and alphaDt (unless a constant D is chosen):

           D = alphaD*nu + alphaDt*nut

Now move to line 144 of the .C file (https://cpp.openfoam.org/v8/scalarTr...8C_source.html) where you'll see a function read() that either reads in the values of alphaD and alphaDt from the scalarTransport dictionary (ie file), or uses the default values. If there's a D value in the dictionary, then it assumes you want to use a constant D, and it sets the flag constantD_ and reads the value into D_:

 bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
 {
    ...
     constantD_ = dict.readIfPresent("D", D_);
     alphaD_ = dict.lookupOrDefault("alphaD", 1.0);
     alphaDt_ = dict.lookupOrDefault("alphaDt", 1.0);
    ...
}

Now move to line 174 of the same file, and in the execute() function there is:

     // Calculate the diffusivity
     volScalarField D("D" + s_.name(), this->D(phi));

This defines and initialises a field for D using the output from this->D(phi). For the latter, return to line 56 of the same file and you'll see the definition:

 Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
 (
     const surfaceScalarField& phi
 ) const
 {
     ... 
     if (constantD_)
     {
         return volScalarField::New
         (
             Dname,
             mesh_,
             dimensionedScalar(Dname, phi.dimensions()/dimLength, D_)
         );
     }
     else if (mesh_.foundObject<icoModel>(momentumTransportModel::typeName))
     {
         const icoModel& model = mesh_.lookupObject<icoModel>
         (
             momentumTransportModel::typeName
         );
 
         return alphaD_*model.nu() + alphaDt_*model.nut();
     }
     else if (mesh_.foundObject<cmpModel>(momentumTransportModel::typeName))
     {
         const cmpModel& model = mesh_.lookupObject<cmpModel>
         (
             momentumTransportModel::typeName
         );
 
         return alphaD_*model.mu() + alphaDt_*model.mut();
     }
     else
     {
         return volScalarField::New
         (
             Dname,
             mesh_,
             dimensionedScalar(Dname, phi.dimensions()/dimLength, 0)
         );
     }
 }

which checks whether there's a turbulence model in action, and if so returns either the dynamic or kinematic total viscosity, depending on the solver type, using the mut()/nut() etc. functions from the RTS selected turbulence model. Or it returns a field based on the constant D value, if the constantD_ flag is set.

Finally, the D field is used in the laplacian, in the fvMatrix, back in execute().

That's the bones of it. To understand the turbulence model references, you'll need to dig into the RTS system - make a LARGE cup of coffee before doing so, though! Good luck.

[KTG]

Thanks for that explanation of the code. I think the reason I didn't find what I was looking for was because it wasn't there! I was able to get from the documentation that D is constructed from a linear combination of nut and nutilda, but I am still confused as to why the effective diffusivity is equal to nut (in the case of RANS simulations). Doesn't this assume that the schmidt number is 1? I have found many expressions in the literature for estimating the effective diffusion for things like flow in channels (lots of analytical solutions), but none that give an expression of how to get this from a turbulence quantity in the case of doing CFD. I am probably searching for the wrong keywords or something.

[KTG]

Well, now that I posted that, I found the answer on page 68 of Versteeg. "we expect that the value of the turbulent diffusivity is fairly close to that of the turbulent viscosity". The keyword I was missing is "Reynolds analogy".



That said, it makes me wonder why we don't use the Reynolds stresses to make a matrix of eddy diffusivity coefficients from the Reynolds stresses and avoid the assumption that everything is isotropic? Maybe its been done and I just haven't found it...

[Tobermory]

Yes - the alphaD and alphaDt from my earlier post are the reciprocal of the Prandtl and turbulent Prandtl numbers if the scalar is temperature, or the Schmidt and turbulent Schmidt numbers if it's mass concentration. And the Reynolds analogy does indeed suggest Pr_t ~ Sh_t ~ 1. Typically CFD codes will use a value of around 0.9, based on exptal data. The (laminar) Prandtl number for most gases is around 0.7. The laminar Schmidt number is rather more variable, I seem to recall.

Why the isotropic assumption? That's typical for most RANS models - they all are based on an isotropic eddy viscosity model. If you want a better RANS model with non isotropic diffusion, eg for dense gas dispersion, then you need to run some kind of algebraic strain model, or a Reynolds stress model that includes transport equations for the scalar fluxes as well as the Re stresses ... that's a very rare beast, and is not available in OpenFOAM. So, the coding in the scalarTransport just matches the availability of the general RANS models. Hope that helps.

[KTG]

Sorry for my delayed response, that does help. But it raises another question for me. When I run the turbulenceFields function object, I am assuming it outputs the Reynolds stress tensor R by calculating:


\overbar{u_i^prime u_j^\prime = \frac{mu_t}{-\rho} ( \frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i} )




Its not as direct a way to get the R tensor as the other options you mentioned, but wouldn't it be better than just assuming that its isotropic for the purposes of advecting a scalar? I feel like I am making doing some circular logic here and thats why no one has gone in this direction, but I don't know. Maybe this way of calculating R is just not that accurate in practice and I haven't seen the validation cases...