Adding a source term in the momentum equation

Questions and remarks about code_saturne usage
Forum rules
Please read the forum usage recommendations before posting.
andrea28
Posts: 41
Joined: Fri Apr 15, 2016 7:25 pm

Adding a source term in the momentum equation

Post by 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
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Adding a source term in the momentum equation

Post by 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
Posts: 41
Joined: Fri Apr 15, 2016 7:25 pm

Re: Adding a source term in the momentum equation

Post by andrea28 »

Thanks, I found the option in the GUI.
andrea28
Posts: 41
Joined: Fri Apr 15, 2016 7:25 pm

Re: Adding a source term in the momentum equation

Post by 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
Posts: 41
Joined: Fri Apr 15, 2016 7:25 pm

Re: Adding a source term in the momentum equation

Post by andrea28 »

Yvan Fournier wrote:
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.
Attachments
source term.txt
(1.61 KiB) Downloaded 238 times
Luciano Garelli
Posts: 280
Joined: Fri Dec 04, 2015 1:42 pm

Re: Adding a source term in the momentum equation

Post by 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
Results.jpg
Regards.
Attachments
Heat_cavity.tar.gz
(45.82 KiB) Downloaded 228 times
andrea28
Posts: 41
Joined: Fri Apr 15, 2016 7:25 pm

Re: Adding a source term in the momentum equation

Post by andrea28 »

Luciano Garelli wrote:

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
Posts: 4070
Joined: Mon Feb 20, 2012 3:25 pm

Re: Adding a source term in the momentum equation

Post by 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
Posts: 280
Joined: Fri Dec 04, 2015 1:42 pm

Re: Adding a source term in the momentum equation

Post by 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.
Predef_symbols.jpg
Example.jpg
Luciano Garelli
Posts: 280
Joined: Fri Dec 04, 2015 1:42 pm

Re: Adding a source term in the momentum equation

Post by 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"
Post Reply