Non-standard boundary conditions : setting a pressure drop

Questions and remarks about code_saturne usage
Forum rules
Please read the forum usage recommendations before posting.
Post Reply
Christophe Vallet

Non-standard boundary conditions : setting a pressure drop

Post by Christophe Vallet »

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 :

Code: Select all

      CALL GETFBR('entree',NLELT,LSTELT) 
C 
      DO ILELT = 1, NLELT 
        IFAC = LSTELT(ILELT) 
        ITYPFB(IFAC,1) = IENTRE 
CC   ON IMPOSE UNE CONDTION DE TYPE DIRICHLET EN ENTREE SUR LA PRESSION 
        ICODCL(IFAC,IPR(IPHAS)) = 1 
        RCODCL(IFAC,IPR(IPHAS),1) = 2.D5 
C CL DE TYPE NEUMANN POUR LES AUTRES VARIABLES 
        ICODCL(IFAC,IU(IPHAS)) = 3 
        ICODCL(IFAC,IV(IPHAS)) = 3 
        ICODCL(IFAC,IW(IPHAS)) = 3 
      ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  C   CL en sortie CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC       CALL GETFBR('sortie',NLELT,LSTELT) 
CC 
      DO ILELT = 1, NLELT 
        IFAC = LSTELT(ILELT) 
        ITYPFB(IFAC,1) = ISOLIB 
        ICODCL(IFAC,IPR(IPHAS)) = 1 
        RCODCL(IFAC,IPR(IPHAS),1) = 1.D5 
        ICODCL(IFAC,IU(IPHAS)) = 3 
        ICODCL(IFAC,IV(IPHAS)) = 3 
        ICODCL(IFAC,IW(IPHAS)) = 3 
      ENDDO 
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 :

Code: Select all

         ICODCL(IFAC,IU(IPHAS)) = 9 
 
        ICODCL(IFAC,IV(IPHAS)) = 9 
 
        ICODCL(IFAC,IW(IPHAS)) = 9 
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
Attachments
listing.txt
(85.74 KiB) Downloaded 171 times
Marc Sakiz

Re: [ask] Non-standard boundary conditions : setting a pressure drop

Post by Marc Sakiz »

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--
Christophe Vallet

Re: Non-standard boundary conditions : setting a pressure drop

Post by Christophe Vallet »

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
Christophe Vallet

Re: Non-standard boundary conditions : setting a pressure drop

Post by Christophe Vallet »

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

Re: Non-standard boundary conditions : setting a pressure drop

Post by Marc Sakiz »

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--
Christophe Vallet

Re: Non-standard boundary conditions : setting a pressure drop

Post by Christophe Vallet »

Hi Marc,

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

Christophe
Attachments
files-tar.bz2
(84.1 KiB) Downloaded 180 times
Post Reply