Problem with basic simulation

[d.charreron]

hello.

I am trying to run a basic simulation of an air flow in a long corridor. This corridor have a inlet, walls and outlet. Until now everything is alright. However, I have trouble when I try to add an second inlet inside the domain. The second inlet, at the center of the corridor is at 5 meters above the ground. I use the usclim.f90 to specify the boundary conditions of the inlet. I want this inlet to have one side as a wall and the other side as an 1m/s inlet. Here is the condition I use:

call getfbr("inlet_hot", nlelt, lstelt)
!==========
do ilelt = 1, nlelt

ifac = lstelt(ilelt)
iel = ifabor(ifac)

itypfb(ifac) = ientre

if (surfbo(1, ifac) > 0) then


ifac = lstelt(ilelt)
itypfb(ifac) = iparoi

else

rcodcl(ifac, iu, 1) = 0.d0
rcodcl(ifac, iv, 1) = 0.d0
rcodcl(ifac, iw, 1) = 1.d0

endif

uref2 = rcodcl(ifac,iu,1)**2 &
+rcodcl(ifac,iv,1)**2 &
+rcodcl(ifac,iw,1)**2
uref2 = max(uref2,1.d-12)

! Hydraulic diameter

dh = 3.3d0
! Turbulence intensity
xintur = 0.02d0

xkent = epzero
xeent = epzero

call keenin &
!==========
( uref2, xintur, dh, cmu, xkappa, xkent, xeent )

! itytur is a flag equal to iturb/10
if (itytur.eq.2) then

rcodcl(ifac,ik,1) = xkent
rcodcl(ifac,iep,1) = xeent

elseif(itytur.eq.3) then

rcodcl(ifac,ir11,1) = d2s3*xkent
rcodcl(ifac,ir22,1) = d2s3*xkent
rcodcl(ifac,ir33,1) = d2s3*xkent
rcodcl(ifac,ir12,1) = 0.d0
rcodcl(ifac,ir13,1) = 0.d0
rcodcl(ifac,ir23,1) = 0.d0
rcodcl(ifac,iep,1) = xeent

elseif(iturb.eq.50) then

rcodcl(ifac,ik,1) = xkent
rcodcl(ifac,iep,1) = xeent
rcodcl(ifac,iphi,1) = d2s3
rcodcl(ifac,ifb,1) = 0.d0

elseif(iturb.eq.60) then

rcodcl(ifac,ik,1) = xkent
rcodcl(ifac,iomg,1) = xeent/cmu/xkent

elseif(iturb.eq.70) then

rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent

endif

if(nscal.gt.0) then
do ii = 1, nscal

rcodcl(ifac,isca(ii),1) = 30.d0

enddo
endif

enddo

!----
! Formats
!----

!----
! End
!----

This doesn't work properly, on Paraview I can see that the flow seems to go through the inlet. Also, the mass flow obtained is 0.

The second problem is about showing a volume with head losses on Paraview. I can see it when I open the "check mech" file but not in my "result.case".

Thanks,


damien

[Yvan Fournier]

Hello,

What version of the code are you using ?

Are you sure "inlet_hot" is a boundary (thin wall) and not a group of interior faces (the mesh check will allow you to check this under ParaView, but even the "listing" file from versions 2.1 or above will give you a count of interior and boundary faces of each group)

Also, in your test, you should only define itypfb for the wall side. The rest (icodcl, rcodcl) should be defined only for the inlet side (move the "endif" much further).

Finally, you do not need to use a user subroutine for your boundary conditions. Using the gui, you can simply boundary condition regions using:

"inlet_hot and normal[1, 0, 0, 0.1]"

and:

"inlet_hot and normal[-1, 0, 0, 0.1]"

selection criteria for the wall and inlet sides respectively.

Best regards,

Yvan

[d.charreron]

Hello,

Indeed, "inlet_hot" appears as an interior face in the listing file. However, I am using the subroutine cs_user_mesh_thinwall.c which permitted to change it into a boundary (thin wall). According to our earliest conversation I understood that using the subroutines are the only solution.

Code saturne is seeing the "inlet_hot" group as an interior face, only after when I start the calculation, the subroutine is changing it into a thin wall. Also, I think It's impossible to define the boundary condition of "inlet_hot" in the code. No?

Regards,

Damien

[d.charreron]

Also, if the problem comes from the contain of the subroutine, I don't really understand the meaning of the following code:
if (surfbo(1, ifac) > 0) then
ifac = lstelt(ilelt)
itypfb(ifac) = iparoi
else
rcodcl(ifac, iu, 1) = 0.d0
rcodcl(ifac, iv, 1) = 0.d0
rcodcl(ifac, iw, 1) = 1.d0
endif

(surfbo(1, ifac) > 0 refers to the side of the surface?

[Yvan Fournier]

Hello,

surfbo(1, ifac) is the X-component of the face's normal (not normalized, proportional to the face surace). Surfbo(2, ifac) and surfbo(3, ifac) are for the Y and Z-components.

So the test is for the orientation of the face normal, to discern between the "inside" and "outside" face.

Best regards,

Yvan

[d.charreron]

Well, thank you for this information. However, it doesn't solve my problem.
Programming is not my strenght, could it be possible to get the subroutine done?
Just with one inlet on one side and wall on other one. It is essencial for every simulation and I am really not able to do it by myself.

Thanks in advance.


Damien

[Yvan Fournier]

Hello,

If programming is not your strength, you (almost) do not need to program. In one of the previous posts, I gave you a solution using only the GUI for boundary conditions. You still need a user subroutine (cs_user_mesh_thinwall.c) to define thin walls, but all the rest can be done with the GUI, so you do not need to program usclim.f90.

Note: with Code_Saturne version 2.2, cs_user_mesh_thinwall becomes part of cs_user_mesh.c (and usclim.f90 is renamed to cs_user_boundary_conditions.f90, but you do not need it).

Best regards,

Yvan

[d.charreron]

Thank you for your answers; finally the both solutions are working fine.
My problem came effectively from the surfbo variable.
I have now another little problem when running the calculation, as it is possible to see on the picture attached, there is a great recirculation of air from the outlet, which leads to stop the calculations after more iterations. I have tried putting a initial velocity put it doesn't help much. What can solve this?

Best regards,

Damien CHARRERON

[Jacques Fontaine]

Hello,

This is a know problem - it is very difficult to simulate an inlet-outlet correctly - .

To reduce/correct this problem:
-you must have a layer of hexahedron (or prisms) on your output/input (maybe it's already the case)
-you can try to move your outlet and inlet away from your obstacle (enlarge your calculation domain)
-you can try to change the boundary condition for the velocity by imposing du/dn = 0 (disable the switch du/dn = 0, if flumas>0 (outlet) and u = 0 else (outlet transformed into an inlet))

Best regards,

JF