[Jost K]

Dear All,

I am still fairly new to this forum so please do not get mad if I am doing this all wrong….

I am working on a project including the simulation of a jet (velocity 2[m/s]) of cold depth water (temperature=10[°C], salinity 35.4 [g/kg]) into the ocean surface region (temperature=25[°C], salinity 36.1 [g/kg]).
After some research in this forum, I found a few threads asking questions about salinity induced buoyant flow in OpenFOAM. However, no solutions have been presented so far.
For this reason, I would like to put together this thread in order to concentrate the efforts in this field and, if possible, derive general instructions on how to model salinity induced buoyant flow in OpenFOAM.

First, I will give a short summary on what has been said in this forum so far. Since there is a high probability that will forget important statements on this matter, I kindly ask all of you to help put together this review.

[1] Incompressible isothermal buoyancy is on isothermal salinity induced buoyant flow. There is a Suggestion to develop a solver based on the buoyantBoussinesq solvers, using a passive scalar transport equation for salinity, which is inspired by the temperature equation implemented in those solvers and calculate a density change and corresponding Buoyancy from this.

[2] Introducing a scalar field for buoyancy production in k-epsilon is on buoyancy driven production of turbulence induced by salinity (if I understand correctly).

[3] add concentration states that temperature and salinity can be modeled using the same equations. However, it is said that modeling salinity is computationally more demanding due the small diffusion coefficients. It is further suggested to model the influences of temperature and salinity in two separate simulations.

We can conclude from this that:

1. A solver based on buoyantBoussinesqPimpleFoam might be the best starting point for the development of our solver. (A thorough documentation for buoyantBoussinesqPimpleFoam can be found here: http://openfoamwiki.net/index.php/Bu...sinesqPisoFoam ).
2. Salinity can be modeled as a passive scalar.

A passive scalar can easily be introduced into a solver (e.g. following this tutorial https://www.youtube.com/watch?v=L94iYGvZr9Q). However, we will need a more sophisticated transport equation for salinity, which might be very similar to buoyantBoussinesqPimpleFoams Teqtn.
Modifying buoyantBoussinesqPimpleFoams linear Boussinesq approximation, we can introduce a term for the salinity influence on density.

Seems easy so far...
However, Ferziger and Peric state that, in water Boussinesqs linear approximation for the temperature density coupling does merely hold for temperature differences up to 2[°C], which is by far not enough for many environmental flows (including mine -see above-).
Even though a higher order polynomial for the temperature density coupling can easily be derived (e. g. based on TEOS-10 http://www.teos-10.org) and implemented, it appears to me that buoyantBoussinesqPimpleFoam is not applicable to problems with large temperature differences, since

1. it becomes very unstable
2. it implements further approximations (like a constant Prandl number) which do not hold for this type of flow.

It thus appears to me that the compressible solver buoyantPimpleFoam has to be used for the simulation of typical environmental flows including salinity.
A short description of the theory is given here: https://caefn.com/openfoam/solvers-buoyantpimplefoam .
Using the icoPolynomial equation of state, buoyantPimpleFoam can be used in “incompressible mode” (i.e. density is only depended on temperature not on pressure), which seems appropriate.
But, how do we implement the density salinity coupling now?
Do we:

1. Implement a new, salinity dependent equation of state?
2. Modify the density on the solver level?

Both ways do not seem very generic to me.

Also, not only the density is salinity dependent. Especially the heat capacity Cp, but also thermal conductivity and dynamic viscosity show some (minor) dependence on salinty.

I will thus stop my monologue here and give the ball to you.
Maybe more of a coding question at this stage but I am sure we will face many more challenges on our way towards a solver for buoyant environmental flow including salinity in OpenFOAM.

Best regards from Kiel(Germany),
Jost

[mAlletto]

This paper may help https://www.researchgate.net/publica..._the_k-e_Model

[Jost K]

Dear Mr. Alletto,

thanks a lot for your input!
Having moved forward a bit on the problem of implementing the density salinity coupling, I now see that turbulence modeling is definitely a big topic when dealing with buoyant environmental flow.

Also based on recent progress, I have to revise a statement I made on buoyantBoussinesqPimpleFoam.

Even though a higher order polynomial for the temperature density coupling can easily be derived (e. g. based on TEOS-10 http://www.teos-10.org) and implemented, it appears to me that buoyantBoussinesqPimpleFoam is not applicable to problems with large temperature differences, since

• it becomes very unstable

Adding a higher order polynomial to buoyantBoussinesqPimpleFoam works fine.
I use the following polynomial:

where:
is Temperature.
is absolute salinity, which is modeled as a passive scalar.
is the standard density.
and the are scalar coefficients.

This approach delivers good results. However, the problem that buoyantBoussinesqPimpleFoam does not have a real thermophysical model, to account for viscosity and thermal conductivity changes due to temperature and salinity, does still remain. It thus still seems attractive to use the solver buoyantPimpleFoam for my simulations.
I have implemented a new salinity dependent thermophysical model based on OpenFOAMs polynomal based thermophysical models. The polynomials I implemented are similar to the ones shown above.

The code structure of the new thermophysical model is shown in the attached pdf (salineThermoProperties.pdf), which is derived from this schematic.

With this recent progress, I now have two solvers for buoyant flow including salinity. However my transport equation for salinity is still very basic.

I would really appreciate if anyone could give me some input on a proper passive scalar transport equation for salinity.

Cheers,
Jost

[mAlletto]

Turbulent diffusion in rans models are usually modelled assuming gradient diffusion. The ratio between mass and momentum diffusivity is described by the turbulent Schmidt number which is a function of stratification. See e.g. https://agupubs.onlinelibrary.wiley....9/2008GL036056. So using a similar equation as for the temperature equation will be OK.

Best

Michael

[TomS]

Dear Jost,

I was interested on how you programmed the polynomial density-salinity reltionship in OpenFOAM. Would it be possible to advise how you went about this? Also, where does this density-salinity reltionship come from?

Thanks.

Tom

www.mts-cfd.com

[Aidah]

Greeting, am kindly requesting for your help. I have been modifying the boussinesq equation of state to add in the salinity in openfoam7. But am failing can you please advice

[Jost K]

Hi,

what problems exactly do you run into?

Originally Posted by Jost K

Adding a higher order polynomial to buoyantBoussinesqPimpleFoam works fine.
I use the following polynomial:

where:
is Temperature.
is absolute salinity, which is modeled as a passive scalar.
is the standard density.
and the are scalar coefficients.

I am sure for many applications a simpler approach is sufficient but this is what I am using.
The code looks like this

rhok = ( C_0 + CT_1 * T + CT_2 * T*T + CT_3 * T*T*T + CSA_1 * salinity ) / rho_0;

Since rhok needs to be dimensionless, the coefficients have funny dimensions.
In my case, they look like this

//Coeffeicients for density approximaton
dimensionedScalar C_0("C_0", dimless, laminarTransport);
dimensionedScalar CT_1("CT_1", dimless/dimTemperature, laminarTransport);
dimensionedScalar CT_2("CT_2", dimless/(dimTemperature * dimTemperature), laminarTransport);
dimensionedScalar CT_3("CT_3", dimless/(dimTemperature * dimTemperature * dimTemperature), laminarTransport);
dimensionedScalar CSA_1("CSA_1", dimless, laminarTransport);

//Reference temperature for rho_0 in K
dimensionedScalar T_0("T_0", dimTemperature, laminarTransport);
//Reference salinity for rho_0
dimensionedScalar SA_0("SA_0", dimless, laminarTransport);
//Calculate reference density
dimensionedScalar rho_0
(
   "rho_0",
   dimless, 
   C_0.value() + CT_1.value() * T_0.value() + CT_2.value() * T_0.value() *T_0.value() + CT_3.value() * T_0.value() *T_0.value() *T_0.value() + CSA_1.value() * SA_0.value()
);

This is probably a slightly messy way to do it but it worked for me and might be a good starting point for you.
Maybe someone in this forum knows a nicer way program this...


I hope this helps with your problem.

Cheers,
Jost

[Aidah]

Thanks alot dear,

Am working on problem involving heattransfer within a solarpond and am using buoyantpimplefoam in openfoam7. The openfoam7 comprise of boussinesq equation of state rho = rho0 = 1 - beta*(T - T0) but depends on only temperature. The problem am working on assumes that the density follows boussinesq approximation rho0 = rhoO*( 1 - beta*(T - T0) + betas*( S - So)) where S is the salinity and the So is the reference salinity depending on both temperature and salinity. On trying to modify the boussinesq equation of state within the openfoam7 adding in salinity and on compling am obtaining an error of S not declared within the scope.

How can I declare S within my equation of state.

Will be grateful for your help

[mAlletto]

Did you include a field for the salinity in the file creatFields

[Aidah]

Yes I did within the solver but the error is appearing in the compling of the thermophysicalModel under specie directory.

[mAlletto]

In the solver you are using the thermodynamic package is responsible for updating the density.

[mAlletto]

the function which calculates rho for compressible pressure based solvers is in the equation of state class.

-

$FOAM_SRC/thermophysicalModels/specie/equationOfState/

[Aidah]

Thanks Alletto,

$FOAM_SRC/thermophysicalModels/specie/equationOfState - this is true. In $FOAM_SRC/thermophysicalModels/specie/equationOfState/boussinesq is what I modified and renamed it myBoussinesq and I modified rho0 = 1 - beta*(T - T0) to rho0 = rhoO*( 1 - beta*(T - T0) + betas*( S - So)). And on compling in the $FOAM_SRC/thermophysicalModels/specie/, it's when the error S is not declared in the scope of the myBoussinesq equation of state appears. Is it that I didn't define S with in EOS.

[TomS]

Hi Aidah,

Please find attached my coding for this which I hope will help.

"C" corresponds to salinity.

Regards,

Tom

www.mts-cfd.com

[Aidah]

thanks Tom


but am using buoyantPimpleFoam which involues the energy equation not temperature equation and the rho is coded in the Equation of state files. please see the added file of my solver and the EOS( where i include in salinity S which am denoting with C). but when i compile in $FOAM_SRC/thermophysicalModels/basic/ - the error




In file included from /home/aidah/OpenFOAM/root-7/run/src/thermophysicalModels/specie/lnInclude/myBoussinesq.H:243:0,
from rhoThermo/rhoThermos.C:33:
/home/aidah/OpenFOAM/root-7/run/src/thermophysicalModels/specie/lnInclude/myBoussinesqI.H: In member function ‘Foam::scalar Foam::myBoussinesq<Specie>::rho(Foam::scalar, Foam::scalar) const’:
/home/aidah/OpenFOAM/root-7/run/src/thermophysicalModels/specie/lnInclude/myBoussinesqI.H:96:51: error: ‘C’ was not declared in this scope
return rho0_*(1.0 - beta_*(T - T0_) + betas_*(C - C0_));
^
/opt/openfoam7/wmake/rules/General/transform:25: recipe for target 'Make/linux64GccDPInt32Opt/rhoThermo/rhoThermos.o' failed

[Aidah]

sorry i had not attached the files but for the solver i have excluded the pressure and velocity equation files (UEqn.H and PEqn.H)


myBoussinesq.zip

solver.zip

[mAlletto]

Not declared in this scope means that your field c is not seen where you want to use it

[Aidah]

Am not so good at programming, How can I really declare it in my file

[Aidah]

Dear TomS


kindly need your help.


am running my simulation but am getting an error of negative temperature when i add scalrcodedsource fvOptions to the simulation.
1) when i use the temperature limit it runs smoothly but the results are wrong
2) i have set the boundary condition of the top surface of the solar pond as zeroGradient please help