code_saturne User's Forum

Hi,

I know this question has already been discussed in an other topic, but I didn't solve my problem.

I simulate a laminar flow through two rough surfaces, and I would like to impose a pressure drop between the inlet and the outlet. I have made several tests.  

My first calculation is carried out with the following usclim.F file in unsteady state : With this boundary condition, calculations diverge whatever the time step used. I also tried to use a time step variable in time. Listing file is given in attachment.

If I use the following condition at the inlet : calculations converge in few iterations (about 200). However, I have some questions. Indeed, from the documentation, if ICODCL(IFAC,IU(IPHAS)) = 9 and if the mass flow is coming in, which is true here, a zero value is imposed to the velocity but not the mass flow. How the velocity can be zero but not the mass flow ? is there a source ? And is there a solution to have a non zero velocity at the inlet ?

Thanks

Hi Christophe,

first a few general comments on the ISOLIB condition.

In Code_Saturne, the mass flux is a separate variable, calculated while solving the mass conservation equation (simplified form : div(rho.u)=0). The mass flux variable is located at the centre of the faces and it is the variable that respects the mass conservation: the sum of the mass flux on all the faces of one element is equal to zero. In the Navier-Stokes resolution, the mass flux is directly used for the convection term.

The velocity variable is located at the centre of the elements.

More precisely, the mass flux is first estimated at the end of the
prediction (Navier-Stokes) phase, using projection of the velocity on
the faces and velocity boundary conditions. But then in the correction
phase (mass conservation), it is modified to respect mass conservation.
 
With the ISOLIB condition, when the mass flux is INGOING (mij<0 since the surface vectors are oriented towards the outside of the domain), the boundary condition of the velocity is indeed put to zero. So after the Navier-Stokes prediction phase your predicted mass-flux will be 0 at the boundary face, but then you go through the pressure correction phase, where the pressure boundary conditions will modify it, so you'll have a (probably) non-null value of the mass flux that respects the div(rho.u)=0. Practically what it does is :
- there is no modifications for non-velocity equations, especially scalars, that will have a correct convection term since it uses the mass flux variable.
- for the velocity equation, this will act as a sink term to prevent the velocity from unstably increasing inwards.
This technique is used for stabilising reasons. It is based on the experience of former codes (in the former code ESTET, the boundary condition velocity could be set to a fraction alpha of the value at the centre, the Code_Saturne method corresponding to alpha=0!).
 
In situations where the flow rate is imposed, there is no problem and you can deal with ISOLIB faces with ingoing flow (for instance a pipe with at one end IENTRE conditions with outgoing flow and at the other end ISOLIB condition that will provide an ingoing flow).

In a case like yours where pressure conditions are imposed and the flow rate adapts itself, I would advise to put the ISOLIB condition where the flow is outgoing. So what you should have is :

* at inlet
  ITYPFB(IFAC,IPHAS)=IENTRE
  ICODCL(IFAC,IPR(IPHAS)) = 1

  RCODCL(IFAC,IPR(IPHAS),1) = 2.D5
  ICODCL(IFAC,IU(IPHAS))=3
  ICODCL(IFAC,IV(IPHAS))=3
  ICODCL(IFAC,IW(IPHAS))=3
  RCODCL(IFAC,IU(IPHAS),3) = 0  (zero Neumann value, it is the default value)
  RCODCL(IFAC,IV(IPHAS),3) = 0
  RCODCL(IFAC,IW(IPHAS),3) = 0
 
* at outlet
  ITYPFB(IFAC,IPHAS)=ISOLIB
  ICODCL(IFAC,IPR(IPHAS)) = 1

RCODCL(IFAC,IPR(IPHAS),1) = 1.D5
  the value ICODCL=9 is set by default for the velocity
 
In this case, the effect will almost identical as your first trial, with ICODCL=3 for velocity at the outlet, but in the transitory convergence phase, and especially at the beginning of the calculation, if you have outlet faces where the mass flux is locally going inwards, the ICODCL=9 "sink term" will stabilise the flow ... hopefully!
 
And of course, be sure to set IPHAS to 1 (it doesn't appear in what you wrote, but it si probably set higher in your routine)
 
Hope it clarifies things. If it still doesn't work, I'll have to look deeper into it.
 
Best regards
 
  --Marc--

Hi Marc,

first of all, thank you for your explanation about ISOLIB condition with ingoing flow. It is a bit clearer for me.

So, as you told me, I put the ISOLIB condition at outlet and there is no difference with the condition ICODCL=3. Calculation diverge after few iterations.

To my mind, the divergence problem comes from inlet.

Cheers

Christophe

Is there any idea to solve my divergence problem ?
 
Thanks
 
Christophe

Hi Christophe,
can you attach the listing of a diverging calculation (check that you write the info at each time step, NTLIST=1), and the fortran routines + xml file, so that I can try to have a look?
 
  --Marc--

Hi Marc,

You will find enclosed all the asked files.
Thanks a lot.

Christophe