Adding a source term in the momentum equation

[andrea28]

Hi,
it's been a couple of week that i have been using CS, it is the first time I'm using an open source code so I'm still figuring out how does it work.
My idea now is to add a source term in the mom equation, specifically the gravitational term for the Boussinesq approximation. I have found a tutorial on this matter, but at the end it creates a subroutine to make the density variable with the temperature so it is not my case.
Reading the reference file for source terms I can't figure out how to put the volume in my source term or how to select the right cells...
Sorry for being a noob, please help!

[Yvan Fournier]

Hello,

To make it easier to help you, please provide the info recommendend in the forum usage guidelines (especially the code version).

User subroutines for source terms are not the simplest, as they require 3 different passes, and I believe for a Boussinesq-type term, the MEI (mathematical expression interpreter) type expressions enabled through the GUI should be enough, so I recommend trying to do this with the GUI. Chec the "examples" tab in the GUI expression definition box itself.

If you are stuck, post again, but tell us exactly how you express you source terms (with explicit constants, variables, etc. as it may be defined in several slightly different ways).

You ma already have used selection criteria fo boundary conditions, but for more info on how to use those (to select a zone), you may search for "selection criteria" in the user documentation pdf.

Regards,

Yvan

[andrea28]

Thanks, I found the option in the GUI.

[andrea28]

Sorry but in the MEI for the momentum source term, I don't know how to specify the temperature, I can't use the TempK identifier.

[andrea28]

If you are stuck, post again, but tell us exactly how you express you source terms (with explicit constants, variables, etc. as it may be defined in several slightly different ways).

I am stuck with the MEI, so I tried to write down the subroutine following the reference file.
My term depends on the temperature and is like:
b * rho * ( T_ref - T ) * g
and should be added only in the y component of the momentum equation.
Now b, rho, T_ref, g are constant T depends on the cell.
In the attachment there is my attempt.

[Luciano Garelli]

Hello,

As Yvan says, you can write your expression (Boussinesq-type term) through the MEI as

Code: Select all

#air expansion coefficient (1/K)
b=0.00343;
#Ref temperature
tref=300;
# Density of air
density = rho0*(1-b*(temperature-tref));

Then in "Gravity" you have to define the acceleration direction.

I attach an example of the heat transfer in a square cavity, it can be helpful

Regards.

[andrea28]

Code: Select all

#air expansion coefficient (1/K)
b=0.00343;
#Ref temperature
tref=300;
# Density of air
density = rho0*(1-b*(temperature-tref));

Thanks for the answer, but that's not what I am looking for. The density in my case remain constant, but there is a source term in the momentum equation:
rho0*Du/Dt = ... + beta * rho0 * (T_ref - T ) * g
Where:
beta=0.000207 1/K
rho0=998.21kg/m3
T_ref=293 K
g=9.81 m/s2 (in the y direction)
Now the question is does the subroutine I posted earlier make any sense?

[Yvan Fournier]

Hello,

I took a quick look at your subroutine, and it seems to make sense (i.e. I did not notice any obvious bug, but I read id fast).

I believed you could also easily adapt the option suggested by Luciano:

Code: Select all

#air expansion coefficient (1/K)
b=0.00343;
#Ref temperature
tref=300;

Su = 0;
Sv = 0;
Sw = b*rho0*(temperature-tref);

dSudu = 0;
dSudv = 0;
dSudw = 0;

dSvdu = 0;
dSvdv = 0;
dSvdw = 0;

dSwdu = 0;
dSwdv = 0;
dSwdw = 0;

But unfortunately, the temperature does not appear amid the predefined symbols for that MEI expression, so this would not currently work (I consider this a GUI/data input bug, and we need to fix it).

So the best solution is probably to test your user subroutine.

Regards,

Yvan

[Luciano Garelli]

Hello,

In the predefined symbols the temperature doesn't appear, but is used in the example tab at least in CS 4.0. In the example that I have posted I only use the GUI and it seems to works.

Additionally, you mention that the density is a constant, but you are using an expansion coefficient, so the density will change with the temperature. Please excuse me if I am not understanding your question.

[Luciano Garelli]

In the subroutine that you post, you use

Code: Select all

crvexp(3,ncelet)
.......
do iel = 1, ncel
     crvexp(iel) =   volume(iel)*beta*rho*(Tref-Temp)*gi
  enddo

And in the documetation crvexp has two argument crvexp(i,iel) for vectorial fields,
http://code-saturne.org/doxygen/src/cs ... 8f90.html
where "The additional force on direction is given by crvexp(i, iel) + vel(j, iel)* crvimp(j, i)."

Also you can check this post
http://code-saturne.org/forum/viewtopic.php?f=2&t=1096
and the file "cs_user_source_terms.f90"