Commit 6f792b6a authored by Alice MAISON's avatar Alice MAISON

Upload src

parent b2d0a4dc
!===============================================================================
! User source terms definition.
!
! 1) Momentum equation (coupled solver)
! 2) Species transport
! 3) Turbulence (k-epsilon, k-omega, Rij-epsilon, v2-f, Spalart-Allmaras)
!===============================================================================
!-------------------------------------------------------------------------------
! Code_Saturne version 6.0-beta
! --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2019 EDF S.A.
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 2 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!-------------------------------------------------------------------------------
!===============================================================================
!> \file cs_user_source_terms.f90
!>
!> \brief User subroutines for additional right-hand side source terms
!>
!> See \subpage cs_user_source_terms and
!> \subpage cs_user_source_terms-scalar_in_a_channel for examples.
!>
!===============================================================================
!> \brief Additional right-hand side source terms for velocity components equation
!> (Navier-Stokes)
!>
!> \section ustsnv_use Usage
!>
!> The additional source term is decomposed into an explicit part (\c crvexp) and
!> an implicit part (\c crvimp) that must be provided here.
!> The resulting equation solved by the code for a velocity is:
!> \f[
!> \rho \norm{\vol{\celli}} \DP{\vect{u}} + ....
!> = \tens{crvimp} \cdot \vect{u} + \vect{crvexp}
!> \f]
!>
!> Note that \c crvexp and \c crvimp are defined after the Finite Volume integration
!> over the cells, so they include the "volume" term. More precisely:
!> - crvexp is expressed in kg.m/s2
!> - crvimp is expressed in kg/s
!>
!> The \c crvexp and \c crvimp arrays are already initialized to 0
!> before entering the
!> the routine. It is not needed to do it in the routine (waste of CPU time).
!>
!> \remark The additional force on \f$ x_i \f$ direction is given by
!> \c crvexp(i, iel) + vel(j, iel)* crvimp(j, i, iel).
!>
!> For stability reasons, Code_Saturne will not add -crvimp directly to the
!> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is
!> treated implicitely only if it strengthens the diagonal of the matrix.
!> However, when using the second-order in time scheme, this limitation cannot
!> be done anymore and -crvimp is added directly. The user should therefore test
!> the negativity of crvimp by himself.
!>
!> When using the second-order in time scheme, one should supply:
!> - crvexp at time n
!> - crvimp at time n+1/2
!>
!> The selection of cells where to apply the source terms is based on a
!> \ref getcel command. For more info on the syntax of the \ref getcel command,
!> refer to the user manual or to the comments on the similar command
!> \ref getfbr in the routine \ref cs_user_boundary_conditions.
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
! mode name role !
!______________________________________________________________________________!
!> \param[in] nvar total number of variables
!> \param[in] nscal total number of scalars
!> \param[in] ncepdp number of cells with head loss terms
!> \param[in] ncesmp number of cells with mass source terms
!> \param[in] ivar index number of the current variable
!> \param[in] icepdc index number of cells with head loss terms
!> \param[in] icetsm index number of cells with mass source terms
!> \param[in] itypsm type of mass source term for each variable
!> (see \ref cs_user_mass_source_terms)
!> \param[in] dt time step (per cell)
!> \param[in] ckupdc head loss coefficient
!> \param[in] smacel value associated to each variable in the mass
!> source terms or mass rate (see
!> \ref cs_user_mass_source_terms)
!> \param[out] crvexp explicit part of the source term
!> \param[out] crvimp implicit part of the source term
!_______________________________________________________________________________
subroutine ustsnv &
( nvar , nscal , ncepdp , ncesmp , &
ivar , &
icepdc , icetsm , itypsm , &
dt , &
ckupdc , smacel , &
crvexp , crvimp )
!===============================================================================
!===============================================================================
! Module files
!===============================================================================
use cstnum
use paramx
use numvar
use entsor
use optcal
use cstphy
use parall
use period
use mesh
use field
use cs_c_bindings
use user_module, only: zone_arbre, lz, alpha
!===============================================================================
implicit none
! Arguments
integer nvar , nscal
integer ncepdp , ncesmp
integer ivar
integer icepdc(ncepdp)
integer icetsm(ncesmp), itypsm(ncesmp,nvar)
double precision dt(ncelet)
double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
double precision crvexp(3,ncelet), crvimp(3,3,ncelet)
! Local variables
logical :: voiture, arbre
double precision meca_voiture, meca_arbre
integer ii, iel, nelt, cellid
double precision CD, normu, tmp_factor
double precision, dimension(:), pointer :: rom
double precision, dimension(:,:), pointer :: vel
integer, allocatable, dimension(:) :: lstelt
!===============================================================================
!===============================================================================
! 1. Initialization
!===============================================================================
! Allocate a temporary array for cells selection
allocate(lstelt(ncel))
! Get density and velocity
call field_get_val_prev_v(ivarfl(iu), vel)
call field_get_val_s(icrom, rom)
meca_arbre = notebook_parameter_value_by_name("meca_arbre")
if (abs(meca_arbre - 0) .lt. epzero) then
arbre = .false.
else
arbre = .true.
endif
! Get CD for trees
CD = notebook_parameter_value_by_name("CD") ! no unit ?
! For trees
if (arbre) then
! Get cells
call getcel(trim(zone_arbre), nelt, lstelt)
! Apply source term
do iel = 1, nelt
cellid = lstelt(iel)
normu = sqrt(vel(1,cellid)**2+vel(2,cellid)**2+vel(3,cellid)**2)
tmp_factor = -rom(cellid)*cell_f_vol(cellid)*alpha*CD
if (normu > epzero) then
crvimp(1,1,cellid) = crvimp(1,1,cellid) &
+ tmp_factor * normu
crvimp(2,2,cellid) = crvimp(2,2,cellid) &
+ tmp_factor * normu
crvimp(3,3,cellid) = crvimp(3,3,cellid) &
+ tmp_factor * normu
endif
enddo
endif
! Deallocate the temporary array
deallocate(lstelt)
return
end subroutine ustsnv
!===============================================================================
!===============================================================================
!> \brief Additional right-hand side source terms for turbulence models
!>
!> \section cs_user_turbulence_source_terms_use Usage
!>
!> The additional source term is decomposed into an explicit part (crvexp) and
!> an implicit part (crvimp) that must be provided here.
!> The resulting equations solved by the code are:
!> \f[
!> \rho \norm{\vol{\celli}} \DP{\varia} + ....
!> = \tens{crvimp} \varia + \vect{crvexp}
!> \f]
!> where \f$ \varia \f$ is the turbulence field of index \c f_id
!>
!> Note that crvexp, crvimp are defined after the Finite Volume
!> integration over the cells, so they include the "volume" term. More precisely:
!> - crvexp is expressed in kg.m2/s2
!> - crvimp is expressed in kg/s
!>
!> The crvexp, crvimp arrays are already initialized to 0 before
!> entering the routine. It is not needed to do it in the routine (waste of CPU time).
!>
!> For stability reasons, Code_Saturne will not add -crvimp directly to the
!> diagonal of the matrix, but Max(-crvimp,0). This way, the crvimp term is
!> treated implicitely only if it strengthens the diagonal of the matrix.
!> However, when using the second-order in time scheme, this limitation cannot
!> be done anymore and -crvimp is added directly. The user should therefore test
!> the negativity of crvimp by himself.
!>
!> When using the second-order in time scheme, one should supply:
!> - crvexp at time n
!> - crvimp at time n+1/2
!>
!> The selection of cells where to apply the source terms is based on a getcel
!> command. For more info on the syntax of the \ref getcel command, refer to the
!> user manual or to the comments on the similar command \ref getfbr in the routine
!> \ref cs_user_boundary_conditions.
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
! mode name role !
!______________________________________________________________________________!
!> \param[in] nvar total number of variables
!> \param[in] nscal total number of scalars
!> \param[in] ncepdp number of cells with head loss terms
!> \param[in] ncesmp number of cells with mass source terms
!> \param[in] f_id field index of the current turbulent variable
!> \param[in] icepdc index number of cells with head loss terms
!> \param[in] icetsm index number of cells with mass source terms
!> \param[in] itypsm type of mass source term for each variable
!> (see \ref cs_user_mass_source_terms)
!> \param[in] ckupdc head loss coefficient
!> \param[in] smacel value associated to each variable in the mass
!> source terms or mass rate (see
!> \ref cs_user_mass_source_terms)
!> \param[out] crvexp explicit part of the source term
!> \param[out] crvimp implicit part of the source term
!_______________________________________________________________________________
subroutine cs_user_turbulence_source_terms &
( nvar , nscal , ncepdp , ncesmp , &
f_id , &
icepdc , icetsm , itypsm , &
ckupdc , smacel , &
crvexp , crvimp )
!===============================================================================
!===============================================================================
! Module files
!===============================================================================
use cstnum
use paramx
use numvar
use entsor
use optcal
use cstphy
use parall
use period
use mesh
use field
use cs_f_interfaces
use cs_c_bindings
use user_module, only: zone_arbre, lz, alpha
!===============================================================================
implicit none
! Arguments
integer nvar , nscal
integer ncepdp , ncesmp
integer f_id
integer icepdc(ncepdp)
integer icetsm(ncesmp), itypsm(ncesmp,nvar)
double precision dt(ncelet)
double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
double precision crvexp(ncelet), crvimp(ncelet)
! Local variables
logical :: arbre
double precision meca_arbre
integer iel, nelt, cellid
character(len=80) :: fname
double precision, dimension(:), pointer :: cpro_rom, cvar_var, cvar_k, cvar_ep
double precision, dimension(:,:), pointer :: vel
double precision normu, CD, tmp_factor
integer, allocatable, dimension(:) :: lstelt
!===============================================================================
!===============================================================================
! 1. Initialization
!===============================================================================
! Allocate a temporary array for cells selection
allocate(lstelt(ncel))
! get field variable, field name and density and velocity and turbulent kinetic energy
call field_get_val_prev_v(ivarfl(iu), vel)
call field_get_val_s(icrom, cpro_rom)
call field_get_val_prev_s(f_id, cvar_var)
call field_get_name(f_id, fname)
call field_get_val_prev_s(ivarfl(ik), cvar_k)
call field_get_val_prev_s(ivarfl(iep), cvar_ep)
!Activate trees
meca_arbre = notebook_parameter_value_by_name("meca_arbre")
if (abs(meca_arbre - 0) .lt. epzero) then
arbre = .false.
else
arbre = .true.
endif
! Get CD for trees
CD = notebook_parameter_value_by_name("CD") ! no unit ?
if (arbre) then
! Get cells
call getcel(trim(zone_arbre), nelt, lstelt)
! Apply k-epsilon model
if (trim(fname).eq.'k') then
do iel = 1, nelt
cellid = lstelt(iel)
normu = sqrt(vel(1,cellid)**2+vel(2,cellid)**2+vel(3,cellid)**2)
tmp_factor = cpro_rom(cellid) * cell_f_vol(cellid) * alpha * CD
crvexp(cellid) = crvexp(cellid) &
+ tmp_factor * 1.d0 * (normu**3) ! Bp=1
crvimp(cellid) = crvimp(cellid) &
- tmp_factor * 5.03d0 * normu ! Bd=5.03
enddo
elseif (trim(fname).eq.'epsilon') then
do iel = 1, nelt
cellid = lstelt(iel)
normu = sqrt(vel(1,cellid)**2+vel(2,cellid)**2+vel(3,cellid)**2)
tmp_factor = cpro_rom(cellid) * cell_f_vol(cellid) * alpha * CD
crvimp(cellid) = crvimp(cellid) &
+ tmp_factor * 0.9d0 * 1.d0 / cvar_k(cellid) * (normu)**3 & ! C4ep=0.9 Bp=1
- tmp_factor * 0.9d0 * 5.03d0 * normu ! C5ep=0.9 Bd=5.03
enddo
endif
endif
! Deallocate the temporary array
deallocate(lstelt)
return
end subroutine cs_user_turbulence_source_terms
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment