!=============================================================================== ! 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 4.0.1 ! -------------------------- ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2015 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. !------------------------------------------------------------------------------- !=============================================================================== ! Purpose: ! ------- !> \file cs_user_source_terms.f90 !> !> \brief Additional right-hand side source terms !> !> \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). !> !> 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 paramx use numvar use entsor use optcal use cstphy use parall use period use mesh use field !=============================================================================== 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(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(3,ncelet), crvimp(3,3,ncelet) ! Local variables integer ipcrom integer ii, ielt, ifac, nlfac, cptfac, iflmas,icel1,icel2 integer, dimension(:), allocatable :: lstfac double precision flowint, flowref, surf double precision, save :: flown, flown1, dp, interr double precision, dimension(:,:), pointer :: vel double precision, dimension(:), pointer :: cpro_rom, imasfl character(30):: inletselection logical :: writeonlisting double precision :: propamplifactor, derivamplifactor, integramplifactor, corr !save flown, flown1, dp !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! TO BE EDITED AT EACH USE flowref = 0.2230308875d3 ! reference mass flux: set your value here inletselection = 'CHANNELBOTTOMBOUNDARY' ! selection string for inlet writeonlisting = .true. ! whether to write on the listing or not propamplifactor = 1 derivamplifactor = 0 !=============================================================================== if (ivar.ne.iu) return ! adapt if you need also iv and iu allocate(lstfac(nfac)) open(799,FILE="pd.txt",status="old") read(799,*) propamplifactor, derivamplifactor close(799) !=============================================================================== ! 2. Body !=============================================================================== ! Is it the first iteration? then set everything to initial vaues if (ntcabs .eq. 1) then dp = 0d0 flown = flowref flown1 = flown interr = 0 else flown1 = flown endif flowint = flown flown = 0.0d0 surf = 0.0d0 cptfac = 0 call getfac(inletselection, nlfac, lstfac) call field_get_key_int(ivarfl(iu),kimasf,iflmas) call field_get_val_s(iflmas, imasfl) do ielt = 1, nlfac ifac = lstfac(ielt) surf = surf + surfan(ifac) flown = flown + imasfl(ifac) cptfac = cptfac + 1 end do if (irangp.ge.0) then call parsom(surf) call parsom(flown) call parcpt(cptfac) end if flown = abs(flown)/surf flowref = abs(flowref) interr = interr + (flown-flowref)*dt(1) corr = (propamplifactor*(flown-flowref)) + & (derivamplifactor*(flown-flown1))/dt(1) !dp = dp + corr dp = dp + (2.0d0*(flown-flowref)-(flown1-flowref))/(2.0d0*dt(1)) do ielt = 1, nlfac ifac = lstfac(ielt) do ii = 1, 2 crvexp(3,ifacel(ii,ifac)) = -volume(ifacel(ii,ifac)) * dp enddo enddo write(nfecra,*) '==FLOW RATE CONTROL==' write(nfecra,*) 'number of inlet faces = ', cptfac write(nfecra,*) 'inlet surface = ', surf write(nfecra,*) 'flow at time n at inlet = ', flown , 'flow forcing = ', dp, 'flowref corr', flowref, corr write(nfecra,*) 'ref. flow', flowref write(nfecra,*) 'flow at time n-1 at inlet = ', flown1 !-------- ! Formats !-------- !---- ! End !---- ! Deallocate the temporary array deallocate(lstfac) return end subroutine ustsnv