!------------------------------------------------------------------------------- ! Code_Saturne version 3.0.0 ! -------------------------- ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2013 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. !------------------------------------------------------------------------------- subroutine ustsma & !================ ( nvar , nscal , ncepdp , & ncesmp , iappel , & icepdc , icetsm , itypsm , izctsm , & dt , rtpa , propce , propfa , propfb , & ckupdc , smacel ) !=============================================================================== ! Purpose: ! ------- ! User subroutine. ! Mass source term ! The subroutine ustsma is called at three different stages in the code ! (iappel = 1, 2 or 3) ! iappel = 1 ! Calculation of the number of cells where a mass source term is ! imposed: ncesmp ! Called once at the beginning of the calculation ! iappel = 2 ! Identification of the cells where a mass source term is imposed: ! array icesmp(ncesmp) ! Called once at the beginning of the calculation ! iappel = 3 ! Calculation of the values of the mass source term ! Called at each time step ! The equation for mass conservation becomes ! d(rho)/dt + div(rho u) = gamma ! The equation for a variable f becomes ! d(f)/dt = ..... + gamma*(f_i - f) ! discretized as ! rho*(f^(n+1) - f^(n))/dt = ..... ! + gamma*(f_i - f^(n+1)) ! f_i is the value of f associated to the injected mass. ! Two options are available: ! - the mass flux is injected with the local value of variable f ! --> f_i = f^(n+1) ! (the equation for f is therefore not modified) ! ! - the mass flux is injected with a specific value for f ! --> f_i is specified by the user ! Variables to be specified by the user ! ===================================== ! ncesmp: number of cells where a mass source term is imposed ! icetsm(ieltsm): identification of the cells where a mass source ! term is imposed. ! For each cell where a mass source term is imposed ! (ielstm in [1;ncesmp]), icetsm(ieltsm) is the ! global index number of the corresponding cell ! (icestm(ieltsm) in [1;ncel]) ! smacel(ieltsm,ipr): value of the injection mass rate gamma (kg/m3/s) ! in the ieltsm cell with mass source term ! itypsm(ieltsm,ivar): type of treatment for variable ivar in the ! ieltsm cell with mass source term. ! * itypsm = 0 --> injection of ivar at local value ! * itypsm = 1 --> injection of ivar at user ! specified value ! smacel(ieltsm,ivar): specified value for variable ivar associated ! to the injected mass in the ieltsm cell with ! a mass source term ! except for ivar=ipr ! ! Remarks ! ======= ! ! - if itypsm(ieltsm,ivar)=0, smacel(ieltsm,ivar) is not used ! - if smacel(ieltsm,ipr)<0, mass is removed from the system, ! therefore Code_Saturna automatically considers f_i=f^(n+1), ! whatever the values of itypsm or smacel specified by the user ! - if a value ivar is not linked for a mass source ! term is imposed, no source term will be taen into account. ! - if a scalar doesn't evolve following the standard equation ! d(rho f)/dt + d(rho U f)/dx = ... ! (alternate convective field for instance), the source term ! set by this routine will nto be correct (except in case of ! injection at the local value of the variable). The proper source ! term should be added directly in ustssc. ! Identification of cells ! ======================= ! The selection of cells where to apply the source terms is based on a getcel ! command. For more info on the syntax of the getcel command, refer to the ! user manual or to the comments on the similar command getfbr in the routine ! cs_user_boundary_conditions. !------------------------------------------------------------------------------- ! Arguments !__________________.____._____.________________________________________________. ! name !type!mode ! role ! !__________________!____!_____!________________________________________________! ! nvar ! i ! <-- ! total number of variables ! ! nscal ! i ! <-- ! total number of scalars ! ! ncepdp ! i ! <-- ! number of cells with head loss terms ! ! ncssmp ! i ! <-- ! number of cells with mass source terms ! ! iappel ! i ! <-- ! indicates which at which stage the routine is ! ! ! ! ! is called ! ! icepdc(ncepdp) ! ia ! <-- ! index number of cells with head loss terms ! ! ! ! ! (usable only for iappel > 1) ! ! icetsm(ncesmp) ! ia ! <-- ! index number of cells with mass source terms ! ! itypsm ! ia ! <-- ! type of mass source term for each variable ! ! (ncesmp,nvar) ! ! ! (see uttsma.f90) ! ! izctsm(ncelet) ! ia ! <-- ! cells zone for mass source terms definition ! ! dt(ncelet) ! ra ! <-- ! time step (per cell) ! ! rtpa ! ra ! <-- ! calculated variables at cell centers ! ! (ncelet, *) ! ! ! (preceding time steps) ! ! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers ! ! propfa(nfac, *) ! ra ! <-- ! physical properties at interior face centers ! ! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers ! ! ckupdc(ncepdp,6) ! ra ! <-- ! head loss coefficient ! ! smacel ! ra ! <-- ! value associated to each variable in the mass ! ! (ncesmp,nvar) ! ! ! source terms or mass rate ! !__________________!____!_____!________________________________________________! ! Type: i (integer), r (real), s (string), a (array), l (logical), ! and composite types (ex: ra real array) ! mode: <-- input, --> output, <-> modifies data, --- work array !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use cstnum use parall use period use mesh !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer iappel integer icepdc(*) integer icetsm(ncesmp), itypsm(ncesmp,nvar) integer izctsm(ncel) double precision dt(ncelet), rtpa(ncelet,*) double precision propce(ncelet,*) double precision propfa(nfac,*), propfb(nfabor,*) double precision ckupdc(ncepdp,6) double precision smacel(ncesmp,nvar) ! Local variables integer ieltsm integer ifac, ii, nlfac integer ilelt, nlelt integer izone integer ipass integer iact integer iactlim data ipass /0/ save ipass save iact character*10 namezone double precision flucel, flutemp double precision vtot , gamma double precision surf, surfx, surfy, surfz double precision flufac, flumaref, flumadd, flumap integer, allocatable, dimension(:) :: lstelt integer, allocatable, dimension(:) :: lstfac !=============================================================================== !=============================================================================== ! PARAMETERS SETTINGS !=============================================================================== flumaref = -0.5 ! Desired mass flow rate (if negative, computed through uref) namezone = "inlet" ! Zone where the source terms are applied iactlim = 10 ! Nombre d'itérations entre deux recalibrages !=============================================================================== !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) allocate(lstfac(nfac)) if (iappel.eq.1.or.iappel.eq.2) then !=============================================================================== ! 1. One or two calls ! First call: ! ! iappel = 1: ncesmp: calculation of the number of cells with ! mass source term ! Second call (if ncesmp>0): ! iappel = 2: icetsm: index number of cells with mass source terms ! WARNINGS ! ======== ! Do not use smacel in this section (it is set on the third call, iappel=3) ! Do not use icetsm in this section on the first call (iappel=1) ! This section (iappel=1 or 2) is only accessed at the beginning of a ! calculation. Should the localization of the mass source terms evolve ! in time, the user must identify at the beginning all cells that can ! potentially becomea mass source term. !=============================================================================== ! 1.1 To be completed by the user: cell selection ! ----------------------------------------------- ! Example 2 : Mass source term one in the cells that ! have a boundary face of color 3 and the cells ! with a coordinate X between 2.5 and 5. ! ! In this test in two parts, one mut pay attention not to count ! the cells twice (a cell with a boundary face of color 3 can ! also have a coordinate X between 2.5 and 5). ! One should also pay attention that, on the first call, the ! array icetsm doesn't exist yet. It mustn't be used outside ! of tests (iappel.eq.2). ! It is quite frequent to forget to remove this example when it is ! not needed. Therefore the following test is designed to prevent ! any bad surprise. izone = 0 ieltsm = 0 ! Cells within the face group (if cells, use getcel, if boundary faces, use getfbr) call getfac(namezone,nlfac,lstfac) izone = izone + 1 do ilelt = 1, nlfac ifac = lstfac(ilelt) ii = ifacel(2,ifac) izctsm(ii) = izone ieltsm = ieltsm + 1 if (iappel.eq.2) icetsm(ieltsm) = ii enddo ! 1.2 Generic subsection: do not modify ! ------------------------------------- ! --- For iappel = 1, ! Specification of ncesmp. This block is valid for both examples. if (iappel.eq.1) then ncesmp = ieltsm endif !------------------------------------------------------------------------------- elseif (iappel.eq.3) then !=============================================================================== ! 2. For ncesmp > 0 , third call ! iappel = 3 : itypsm : type of mass source term ! smacel : mass source term ! Remark ! ====== ! If itypsm(ieltsm,ivar) is set to 1, smacel(ieltsm,ivar) must be set. !=============================================================================== ! 2.1 To be completed by the user: itypsm and smacel ! -------------------------------------------------- !------------------------------------------------------------------------------- ! Example 2 : simulation of a suction (by a pump for instance) with a ! total rate of 80 000 kg/s. ! The suction rate is supposed to be uniformly distributed ! on all the cells selected above. ! It is quite frequent to forget to remove this example when it is ! not needed. Therefore the following test is designed to prevent ! any bad surprise. ! Calculation of the total volume of the area where the mass source ! term is imposed (the case of parallel computing is taken into ! account with the call to parsom). if (ntcabs.gt.1.and.iact.eq.iactlim) then vtot = 0.d0 do ieltsm = 1, ncesmp vtot = vtot + volume(icetsm(ieltsm)) enddo if (irangp.ge.0) then call parsom (vtot) endif ! Précaution: le volume doit être plus grand que 0 if (vtot.gt.0.d0) then ! Boucle pour calcul du flux massique et de la surface call getfac('inlet',nlfac,lstfac) surf = 0.d0 do ilelt = 1, nlfac ifac = lstfac(ilelt) flufac = flufac + propfa(ifac,ipprof(ifluma(ipr))) surf = surf + surfan(ifac) enddo ! Si pas de donnée pour flumaref (flumaref < 0), calcul via uref if (flumaref.lt.0) then flumaref = -uref*propce(1,ipproc(irom))*surf endif ! Comparaison du flux massique de référence avec le flux massique actuel flumadd = flumaref - flufac ! Correction gamma = flumadd/vtot else write(nfecra,9000) vtot call csexit (1) endif ! Mise à jour du terme source dans toutes les cellules flucel = 0.d0 do ieltsm = 1, ncesmp smacel(ieltsm,ipr) = gamma flucel = flucel+ & volume(icetsm(ieltsm))*smacel(ieltsm,ipr) enddo ! Si calcul en parallèle if (irangp.ge.0) then call parsom (flucel) endif ! Affichage dans le listing write(nfecra,2000) flucel, namezone, (flufac+flucel) iact = 0 endif print *, ' ' iact = iact + 1 endif !-------- ! Formats !-------- 1000 format(/,'Mass rate generated in the domain: ',E14.5,/) 2000 format(/,'Mass flux rate generated in the domain: ',E14.5,' kg/s/m³'/, & ' distributed in the zone: ',A,/, & 'Total mass flux in the zone: ',E14.8,' kg/s'/) 9000 format(/,'Error in ustsma ',/, & ' the volume of the mass suction area is = ',E14.5,/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine ustsma