Specify a rotation velocity on boundary

[Ruonan]

Hello everyone,

I am simulating rotating flows in an enclosed rotor-stator cavity. The boundary conditions are shown below.

I want to specify a rotation velocity (for example 100 rad/s) on the two boundaries: rotor and hub. Can I do this in the latest version of Saturne GUI, or I need to define it in the user subroutine /SRC/REFERENCE/cs_user_boundary_conditions.f90?

In Code_Saturne 6.0.5, I found when I choose "Sliding wall" in "Boundary conditions", the U, V, W are all in Cartesian frame.

As for the setting in user subroutines, because I am not familiar with it, is there any tutorials or examples I can follow please?

Best regards,
Ruonan

[Yvan Fournier]

Hello,

I do not think there are specific examples for this, but you may find general examples in the documentation:
https://www.code-saturne.org/cms/sites/ ... mples.html

What you need here is simply a Dirichlet value for the velocity whose vecto is oriented based on the current position.

If your case is a for a turbomachinery application, also check the turbomachinery features. There is a tutorial for more complex and transient meshes(https://www.code-saturne.org/cms/sites/ ... SALOME.pdf) but in your case, the "frozen rotor" option is probably more similar to what you are trying to do (so you might not need user-defined functions after all).

Best regards,

Yvan

[Ruonan]

Hello Yvan,

Thank you very much for your reply.

Yes, the case is for a turbomachinery application, but I think "forzen rotor" may not be the best option. I tried to use frozen rotor for this case, found I can only specify rotation velocity on cells, not on faces. But I want to let the boundary faces rotate, while leave the whole domain cells in stationary frame.

I agree that I need to specify a Dirichlet value for the velocity on these two boundary faces "Rotor and Hub". But in order to specify the rotation speed on the faces, I need to define velocity components u v w in Cylindrical frame.

I think I need to:

(1) Define a variable "theta" using y and z (here x is the rotation axis).

(2) Define v, w using theta (predefined) and rotation speed (which is a constant).

Here are the codes I added in cs_user_boundary_conditions.f90:

Code: Select all

! INSERT_VARIABLE_DEFINITIONS_HERE

real::theta
theta=atan(z/y)

integer, allocatable, dimension(:) :: lstelt

!===============================================================================

......


! INSERT_MAIN_CODE_HERE

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

call getfbr('ROTOR and HUB', nlelt, lstelt)       !choose two boundary faces named ROTOR and HUB respectively

do ilelt = 1, nlelt
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   
!Dirichlet B.C for velocity:

    icodcl(ifac,iu) = 1               ! 1 means Dirichlet
    rcodcl(ifac,iu,1) = 0                           !ux

    icodcl(ifac,iv) = 1               ! 1 means Dirichlet
    rcodcl(ifac,iv,1) = 766*sin(theta)              !uy

    icodcl(ifac,iw) = 1              ! 1 means Dirichlet
    rcodcl(ifac,iw,1) = -766*cos(theta)          !uz

enddo

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

But when I trying to compile it in SRC folder, many errors appeared.

I think there must be many mistakes in the codes, so any help will be much appreciated. I am new to Code_Saturne and Fortran, trying to make progress.

Best regards,
Ruonan

[Yvan Fournier]

Hello,

To compute theta, you need to define z and y first. And base them on each face's coordinates. So using :

Code: Select all

theta = atan(cdgfbo(3,ifac) / cdgfbo(2,ifac))

inside the loop is better.

Also, if you have compilation errors you cannot solve, you can post them here. But in most cases, reading the error messages carefully should provide hints to what you need to do. If you are not familiar at all with Fortran or C programming, you can find many tutorials or courses on the web. I could suggest a few.

Best regards,

Yvan

[Ruonan]

Hello Yvan,

Thank you very much for your help. Yes I am learning Fortran online with some tutorials recently.

I revised my codes following your suggestions, it improved a lot, but there are still two errors which I didn't manage to solve it. Could you please have a look at it?

The subroutine file are attached below.

Thank you very much in advance!

Best regards,
Ruonan

[Yvan Fournier]

Hello,

According the the error message, you cannot initialize the theta value at declaration time.

You need to simply declare it outside the loop, and assign the value as you have done inside the loop on faces.

Regards,

Yvan

[Ruonan]

Hello Yvan,

Thank you very much for your reply and help. Finally I succeed in specifying a rotation velocity on boundary face. I would like to post my codes here.

Code: Select all

subroutine cs_f_user_boundary_conditions &
 ( nvar   , nscal  ,                                              &
   icodcl , itrifb , itypfb , izfppp ,                            &
   dt     ,                                                       &
   rcodcl )

!===============================================================================

!===============================================================================
! Module files
!===============================================================================

use paramx
use numvar
use optcal
use cstphy
use cstnum
use entsor
use parall
use period
use ihmpre
use ppppar
use ppthch
use coincl
use cpincl
use ppincl
use ppcpfu
use atincl
use atsoil
use ctincl
use cs_fuel_incl
use mesh
use field
use turbomachinery
use iso_c_binding
use cs_c_bindings

!===============================================================================

implicit none

! Arguments

integer          nvar   , nscal

integer          icodcl(nfabor,nvar)
integer          itrifb(nfabor), itypfb(nfabor)
integer          izfppp(nfabor)

double precision dt(ncelet)
double precision rcodcl(nfabor,nvar,3)


! Local variables
integer          ifac, iel, ii
integer          ilelt, nlelt
double precision theta
double precision x , y , z
double precision r

integer, allocatable, dimension(:) :: lstelt

! INSERT_VARIABLE_DEFINITIONS_HERE





!===============================================================================


!===============================================================================
! Initialization
!===============================================================================



allocate(lstelt(nfabor))  ! temporary array for boundary faces selection


! INSERT_ADDITIONAL_INITIALIZATION_CODE_HERE

!===============================================================================
! Assign boundary conditions to boundary faces here

! For each subset:
! - use selection criteria to filter boundary faces of a given subset
! - loop on faces from a subset
!   - set the boundary condition for each face
!===============================================================================

! INSERT_MAIN_CODE_HERE

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

call getfbr('ROTOR or HUB', nlelt, lstelt) !choose two boundary faces named ROTOR and DOWN respectively

do ilelt = 1, nlelt
   ifac = lstelt(ilelt)

x=cdgfbo(1,ifac) !cdgfbo means centers of faces on boundary
y=cdgfbo(2,ifac)
z=cdgfbo(3,ifac)

theta=atan(y/z)    ! define the theta

r=sqrt(y*y+z*z)  ! define the radius r

!   iel = ifabor(ifac)
!   itypfb(ifac) = iparoi !itypfb is used to define the type of faces, iparoi means smooth solid wall face, impermeable and with friction

!Dirichlet B.C for velocity:

    icodcl(ifac,iu) = 1 ! icodcl defines the type of boundary condition for the variable iu on the face ifac, 1 means Dirichlet
    rcodcl(ifac,iu,1) = 0.0d0 !rcodcl gives the numerical values associated with the type of boundary conditions, ux

    icodcl(ifac,iv) = 1 ! 1 means Dirichlet
    rcodcl(ifac,iv,1) = 7.66d0*r*cos(theta) !uy

    icodcl(ifac,iw) = 1 ! 1 means Dirichlet
    rcodcl(ifac,iw,1) = (-1)*7.66d0*r*sin(theta) !uz

enddo

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

deallocate(lstelt)  ! temporary array for boundary faces selection

return
end subroutine cs_f_user_boundary_conditions

Thank you again for your time and kind help!

Best regards,
Ruonan

[Ruonan]

Hello,

Sorry in last post there are some errors in the codes. Here are the correct ones:

Code: Select all

subroutine cs_f_user_boundary_conditions &
 ( nvar   , nscal  ,                                              &
   icodcl , itrifb , itypfb , izfppp ,                            &
   dt     ,                                                       &
   rcodcl )

!===============================================================================

!===============================================================================
! Module files
!===============================================================================

use paramx
use numvar
use optcal
use cstphy
use cstnum
use entsor
use parall
use period
use ihmpre
use ppppar
use ppthch
use coincl
use cpincl
use ppincl
use ppcpfu
use atincl
use atsoil
use ctincl
use cs_fuel_incl
use mesh
use field
use turbomachinery
use iso_c_binding
use cs_c_bindings

!===============================================================================

implicit none

! Arguments

integer          nvar   , nscal

integer          icodcl(nfabor,nvar)
integer          itrifb(nfabor), itypfb(nfabor)
integer          izfppp(nfabor)

double precision dt(ncelet)
double precision rcodcl(nfabor,nvar,3)


! Local variables
integer          ifac, iel, ii
integer          ilelt, nlelt
double precision theta
double precision x , y , z
double precision r

integer, allocatable, dimension(:) :: lstelt

! INSERT_VARIABLE_DEFINITIONS_HERE





!===============================================================================


!===============================================================================
! Initialization
!===============================================================================



allocate(lstelt(nfabor))  ! temporary array for boundary faces selection


! INSERT_ADDITIONAL_INITIALIZATION_CODE_HERE

!===============================================================================
! Assign boundary conditions to boundary faces here

! For each subset:
! - use selection criteria to filter boundary faces of a given subset
! - loop on faces from a subset
!   - set the boundary condition for each face
!===============================================================================

! INSERT_MAIN_CODE_HERE

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

call getfbr('ROTOR or HUB', nlelt, lstelt) !choose two boundary faces named ROTOR and DOWN respectively

do ilelt = 1, nlelt
   ifac = lstelt(ilelt)
   iel = ifabor(ifac)   
   itypfb(ifac) = iparoi !itypfb is used to define the type of faces, iparoi means smooth solid wall face, impermeable and with friction

x=cdgfbo(1,ifac) !cdgfbo means centers of faces on boundary
y=cdgfbo(2,ifac)
z=cdgfbo(3,ifac)

theta=atan(y/z)    ! define the theta

r=sqrt(y*y+z*z)  ! define the radius r

!Dirichlet B.C for velocity:

    icodcl(ifac,iv) = 1 ! 1 means Dirichlet
    rcodcl(ifac,iv,1) = 766.0d0*r*cos(theta) !uy

    icodcl(ifac,iw) = 1 ! 1 means Dirichlet
    rcodcl(ifac,iw,1) = (-1)*766.0d0*r*sin(theta) !uz

enddo

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

deallocate(lstelt)  ! temporary array for boundary faces selection

return
end subroutine cs_f_user_boundary_conditions

Best regards,
Ruonan