Discussion Forum » Code_Saturne » ustsns for testcase5 : Vertical Heated Pipe


Hello, I am a new user of code-saturne and I am very interesting by the test case5.

Is there already done the file ustsns.f90 ? And may be with some explanations.

My problem is the following: I have to perform a simulation in a non-plane channel with periodic condition in the flow direction. A constant heat flux is imposed to the walls.

What I wish to do is to have a temperature of 300K at the beginning of the domain and a mass flow rate of 30 m3s-1.

So I have to add a source term for the momentum equation and another for the energy equation.

Any help or advice on this problem (or on the solution of the test case5) would be appreciated. Thanks,


-- SylvainSerra - 2010-09-29

Hi If you have a look at the diffuser test case here you will see that the inlet conditions are calculated from a periodic channel (IN-KEP, IN-PHI etc ) which uses the subroutine ustsns.F to impose a fixed mass flux (rather than a fixed dP/dx). This is done by computing the mass flux at each iteration and comparing it with the reference (UREF). %BR%

For the temperature use the subroutine ustssc.f90

-- JuanUribe - 2010-09-30

Thank you Juan for your answer,

I try to understand what is done in the ustsns file of this test case but I do not understand all...

specialy, this part:

do ifac=1,nfac surfx = surfac(1,ifac) surfy = surfac(2,ifac) surfz = surfac(3,ifac) surfan = sqrt(surfx**2 + surfy**2 + surfz**2) if(abs(cdgfac(1,ifac)-(-11.1d0)).lt.0.000001d0) then iifac = iifac + 1 debitn = debitn + propfa(ifac,iflmas) surf = surf + surfan endif enddo

what does test the if (abs(cdgfac(1,ifac)-(-11.1d0)).lt.0.000001d0) ?

and how can I know the surface used for the calcul ?

Many thanks for your help,


-- SylvainSerra - 2010-10-01

Hi Sylvain, that loop is trying to find the area which you want to use to compute the mass flux.

  • NFAC is the number of internal faces
  • cdgfac() is an array with the coordinate of the centre of gravity of that face
  • sufac() is the normal vector whose length is the area of the face.
  • propfa(ifac,iflmas) has the mass flux through the face

the test tries to select all faces that are close to x=0 which is where the mass flux is computed for that particular case. In your case you should find the corresponding plane(s) where you want to compute it and modify that test.

I hope this helps

-- JuanUribe - 2010-10-07

Hi Juan, YES, it help me. But have met a problem, I only change the following line: if (abs(cdgfac(1,ifac)).lt.0.000001d0) because in my case, the origin is located at (0.0.0)

if I change only the location in the x direction, the result of the simulation give me the good number of face and the good area but increase the velocity more over than 70000 m s-1. In my case, the "u ref" is 1.3 m s-1.

the mass flow rate solved at the time step N et N-1 are around 1E-7 whereas the reference is 1364.

so it did not works and I do not understand why.

the other test I did is to not do the test of "if (abs(cdgfac(1,ifac)-(-11.1d0)).lt.0.000001d0)"

In that case, the number of faces and the area are not good but the two mass flow rate (N and N-1) are in agreement with the reference. In that case, the velocity is already to height but converge in 8.5 m s-1.

May be you already met this problem? Now I think that I understand wath is done by the ustsns file but I do not understant why the results are not good...

for more information, my case is in 3D, the same kind of simulation than the DNS of kim and Moin (0<x<0.09 ; 0<y<0.03 ; 0<z<0.06) with u Ref = & m S-& and the density 1050; the walls are in y=0 and y=0.03

-- SylvainSerra - 2010-10-07

Hi Sylvain, if you are trying to compute the channel flow to compare with DNS, it might be better to look at the channel flow case:

where you will see that the ustsns.f90 subroutine is used to impose a fixed dP which gives you the friction velocity u_* that you want (i.e. instead of fixing the mass flux, fix u_*).

If you still want to fix the mass flux, what you are doing seems to be fine but there might be a problem with the tolerance. If your reference plane is at 0, try using a test (abs(cdgfac(1,ifac)).lt.xeps) where xeps is a small number. Smaller than the distance between two cells in the X direction but bigger than the difference between the X coordinates of the cells in the reference plane. Try changing xpes until the number of faces the subroutine finds is the the same as in the reference plane.

-- JuanUribe - 2010-10-18

Hi juan,

The ustsns file gives the good surface and the good number of faces...

in fact, my final job will be on a complex geometrie in parallel and I think that it could be faster to impose the mass flow than to solve the friction velocity... but may be it is false.

I have a question about the initialisation of the domain. Is there something special to do for the initialisation when we use periodic boundary? In fact, when I look at the mass flow calculated at each iteration, the value is around 1 10-15 so I have the feeling that my initial conditions are not taken into account.

Is there the ustsns file for the testcase5 and solution that the source term is imposed into the implicit part?

thank you very much for all your answers

-- SylvainSerra - 2010-10-22

Hi, I have other questions, the ustsns file developped for the diffuser test case is done for a geometry where there are only on cell in the y direction and one cell in the z direction. In my case, they are many cells in each direction. Is this diffrence could explain the fact that this file does not works for my geometry?

Is there some diffrence between the version that was used for the test case of the diffuseur and the version 2RC1 (that I use)?

what about the command "getfac", is this command can be use to solve the mass flow rate on the periodic boundary at each itération.

May be the anwer of my problem is in these questions...

I think that my problem is how to calclate the mass flow rate at each iteration on the periodic boundary, because the one solve actually is very far from the reference.

If somebody have an idea, I'd appreciate.


-- SylvainSerra - 2010-10-26

Hi, the problem I had come from the mesh. It was created using gmsh and in that case, the command ifluma gives always zero.... I don't know why...

At the end, I create a new mesh using gambit and it's works now.

I have always a question about the file ustsns.f90 of the diffuseur test case. I don't well understand the following ligne:

dp = dp + & (2.0d0*(debitn-debitref)-(debitn1-debitref))/(2.0d0*dt(1))

is that a kind of time scheme?.. or a relaxation factor?

Many thanks for your help,


-- SylvainSerra - 2011-01-11

Hi Sylvain, nice to hear you could solve the problem. it is a relaxation factor. You want to increase (or decrease) the source term based on the difference between the reference bulk velocity (debitref) and the calculated bulk velocity (debitn, at time step n or debitn1 at timestep n-1). You don't need to use this form which in fact is only usefull for unsteady calculations with constant delta t. You could also write something like $\Delta p = \alpha* (U_{ref}-U_b^n) + (1-\alpha)*(U_b^{n-1}-U_b^n) $

-- JuanUribe - 2011-01-14

You can also check the fortrans here: http://cfd.mace.manchester.ac.uk/twiki/pub/CfdTm/TestCase013Res000/ustsns_INLET.F which is simpler to understand.

-- JuanUribe - 2011-01-14


Current Tags:
create new tag
, view all tags
Topic revision: r10 - 2011-01-14 - 17:07:07 - JuanUribe

Computational Fluid Dynamics and Turbulence Mechanics
@ the University of Manchester
Copyright © by the contributing authors. Unless noted otherwise, all material on this web site is the property of the contributing authors.

%TMPL:DEF{"top:toolbarbuttons"}Raw editAttachPrint version