!=============================================================================== ! User source terms definition. ! ! 1) Momentum equation (segegrated solver) ! 2) Momentum equation (coupled solver) ! 3) Species transport ! 4) Turbulence (k-epsilon, k-omega, Rij-epsilon, v2-f, Spalart-Allmaras) !=============================================================================== !------------------------------------------------------------------------------- !VERS ! 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. !------------------------------------------------------------------------------- !=============================================================================== ! 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 use Usage !> !> The routine is called if the coulped solving of the velocity components is !> turned on (\c ivelco=1). !> !> 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} \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). !> !> 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 ustsma) !> \param[in] dt time step (per cell) !> \param[in] rtpa calculated variables at cell centers !> (preceding time steps) !> \param[in] propce physical properties at cell centers !> \param[in] propfa physical properties at interior face centers !> \param[in] propfb physical properties at boundary face centers !> \param[in] ckupdc head loss coefficient !> \param[in] smacel value associated to each variable in the mass !> source terms or mass rate (see \ref ustsma) !> \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 , & dt , rtpa , rtp, propce , propfa , propfb , & crvexp ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== use paramx use numvar use entsor use optcal use cstphy use parall use period use mesh !=============================================================================== implicit none ! Arguments integer nvar , nscal integer ncepdp , ncesmp integer ivar double precision dt(ncelet), rtpa(ncelet,*), rtp(ncelet,*) double precision propce(ncelet,*) double precision propfa(nfac,*), propfb(nfabor,*) double precision crvexp(3,ncelet) ! Local variables integer ipcrom, ifac, iifac, ielt, iflmas, ii integer compt, iel double precision surf, dptmp double precision flowref, flown, dp, flown1, flowint save flown, flown1, dp, compt integer, allocatable, dimension(:) :: lstfac !============================================================================== compt = compt + 1 print *, compt allocate(lstfac(nfac)) !=============================================================================== ! 1. Initialization !=============================================================================== uref = 10.0d0 if (ivar.ne.iu) return ! To adapt if need of iv and iw (orientation of the tube) if (iwarni(ivar).ge.1) then write(nfecra,*) 'User source terms for variable ', nomvar(ipprtp(ivar)) endif iflmas = ipprof(ifluma(iu)) ipcrom = ipproc(irom) !=============================================================================== ! 2. Example of arbitrary source term for component u: ! S = A * u + B ! appearing in the equation under the form: ! rho*du/dt = S (+ standard Navier-Stokes terms) !In the following example: ! A = -rho*CKP ! B = XMMT ! !with: ! CKP = 1.D0 [1/s ] (return term on velocity) ! MMT = 100.D0 [kg/m2/s2] (momentum production by volume and time unit) ! !which yields: ! crvimp(iel) = volume(iel)* A = - volume(iel)*(rho*CKP ) ! crvexp(iel) = volume(iel)* B = volume(iel)*(XMMT ) ! ---------------------------------------------- ! 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. ! ---------------------------------------------- ! Initialisation flowref = uref*996.25d0 flowint = flown flown = 0.0d0 surf = 0.0d0 dp = 0.0d0 call getfac('inlet', ifac, lstfac) ! Determine surface and current flow do ielt = 1, ifac iifac = lstfac(ifac) surf = surf + surfan(iifac) flown = flown + propfa(iifac, iflmas) enddo if (irangp.ge.0) then call parsom(surf) call parsom(flown) endif ! Convert flow to surface flow flown = abs(flown) / surf if (ntcabs.eq.1) then dp = 0.0d0 flown = flowref flown1 = flown else flown1 = flowint endif dp = dp + (2.0d0*(flown-flowref)-(flown1-flowref))/2.0d0 do ielt = 1, ifac iifac = lstfac(ielt) do ii = 1, 2 iel = ifacel(ii, iifac) crvexp(1, iel) = volume(iel) * dp / propce(iel, ipcrom) print *, propce(iel, ipcrom) enddo enddo !-------- ! Formats !-------- !---- ! End !---- ! Deallocate the temporary array deallocate(lstfac) return end subroutine ustsnv