Inlet velocity profile with code_saturne

Hi all,

I'm quite new using code_saturne, and I want to use it to determine the airflow around buildings (this is the general idea). My problem is that at the inlet of the domain, the air is coming with a power law (or a log law) profile, and I don't really know how to specify such a boundary condition in code_saturne.

I've been looking in the tutorials given on the edf-website, and I've found that it should be set in the subroutine usclim.F . So I had a look at this subroutine, and found some infos about it, but it isn't clear in my head.

Basically, I want to set a velocity inlet such as U (ux, uy and uz components of the velocity) function of the height z with a known power law. Could somebody give me a hint about how to do it ?? Or an example of what I should change in the subroutine ? I can't find any of that on the Internet...

Thanks a lot in advance for your help

Alicia

Hi Alicia,

Do you have checked on the Twiki of Code saturne? There are some tutorial and there is another forum.

cfd.mace.manchester.ac.uk/twiki/bin/view/Saturne/WebHome

When you'll find your answer, please post it here. You have a very good question.

Good luck,

Michael Page
IDÉA Innovation PME
www.ideapme.ca

Hi all !

So a bit late but I come back here with an answer which I'm sure will be very helpful for some of you.

To set a inlet with a velocity profile non-constant with saturne :
1-set a constant inlet velocity when defining the boundary condition (which will not be taken into account, but you have to give a value anyway)
2-go to the directory STUDY/CASE/SRC/REFERENCE/BASE, copy the routine usclim.F and paste it in the directory STUDY/CASE/SRC.
3-in the routine, just keep the definition of the boundary which has a complex inlet condition (or outlet, or else). Keep the initialization part too.
Get your boundary with the function GETFBR :

=> call GETFBR('boundary color',NLELT,LSTELT)

where 'boundary color' is the id of your face : it can be a color, a group name, a geometric criterion...
keep the loop on the faces and on the phases :

=> DO IELT=1,NLELT
=> LSTELT(ILELT)
=> IEL=IFABOR(IFAC)
=> DO IPHAS=1, NPHAS

assign the type of boundary to your face :
=> ITYPB(IFAC,IPHAS)=IENTRE
(for an inlet)
then set the expression of your velocity profile for each component. If you have one (or more) velocity components depending on one (or more) coordinates, you need to get the coordinates of the center of the faces of the cells which are on your boundary (am I clear here ??), do that with the function : CDGFBO(i,IFAC) where i=1 if you want the first coordinate (x), 2 if you want y and 3 if you want z.
To assign the value or expression to the velocity components, use the function RCODCL(IFAC,IU(IPHAS),1) IU is for the first component of the velocity, 1 is to assign a Dirichlet condition.
I give you an example of my routine for an inlet velocity profile in which the horizontal velocity u is function of the vertical coordinate z with a power law : u(z)=2,91*z^(0.25) and the other velocity components were 0.

here it is. Before that just keep the part uncommented (name of the subroutine,common data,arguments)

C=======================================================================
C 1. INITIALISATIONS
C
C=======================================================================
C
IDEBIA = IDBIA0
IDEBRA = IDBRA0
C
D2S3 = 2.D0/3.D0
C
C=======================================================================
C 2. REMPLISSAGE DU TABLEAU DES CONDITIONS LIMITES
C ON BOUCLE SUR LES FACES DE BORD
C ON DETERMINE LA FAMILLE ET SES PROPRIETES
C ON IMPOSE LA CONDITION LIMITE
C
C IMPOSER ICI LES CONDITIONS LIMITES SUR LES FACES DE BORD
C
C=======================================================================
C
C
C --- Assign an inlet velocity profile for u with a power law in z C --- to the boundary 'Inlet'
C
CALL GETFBR('Inlet',NLELT,LSTELT)
C ===========
C
DO ILELT = 1, NLELT
C
IFAC = LSTELT(ILELT)
IEL = IFABOR(IFAC)
C
DO IPHAS = 1, NPHAS
C
ITYPFB(IFAC,IPHAS) = IENTRE
C
RCODCL(IFAC,IU(IPHAS),1) = 2.91D0*(CDGFBO(3,IFAC)**0.25D0)
RCODCL(IFAC,IV(IPHAS),1) = 0.0D0
RCODCL(IFAC,IW(IPHAS),1) = 0.0D0

ENDDO
C
ENDDO
C

C
C----
C FORMATS
C----
C
C----
C FIN
C----
C
RETURN
END
c@z


I hope it helps. Of course you can do the same for some variables of the turbulence, but I haven't gone so far yet... !

Enjoy it, and have a good evening everyone !

Alicia

Hi, and thank you very much for this explanation.

Nice answer Alicia! This is indeed the right way to do it, but let me give you a different solution.

Provided you are using the latest version of Code_Saturne (e.g. version 2.0-rc1) with MEI support (MEI stands for Mathematical Expression Interpretor), you can do the same thing within the graphical interface.

In the "inlet" boundary definition, choose the "norm (user law)" item for the velocity setup and then, click on the button on the right. Finally, you can write such formula:[code:1]u_norm = 2.91*z^0.25;[/code:1]

You'll find some explanation and example in the different tabs of this popup.

Hope this will help someone
David

You are right, I made an error in mapping the nodes.


I tried to correct the file (check the attached file). Unfortunately I have no been able to work with Calculix lately, so the error persisted.
Try to substitute the converter with the attached one.
If you vant to play with it, the array inp_to_med_xxx maps the nodes from aster in the order expected from calculix.