!=============================================================================== ! 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 5.0.8 ! -------------------------- ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2018 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 !> !> 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). !> !> 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 character*80 chaine1, chaine2 integer iel, nlfac, ifac, ielt, ii, iflmas integer nlelt, cptfac double precision flowint, flowref, surf, value double precision flown, flown1, dp, rho_ref, mu_ref, hydraulic_d, Re_d, f, tau_w, temp integer, allocatable, dimension(:) :: lstfac, lstelt double precision, dimension(:), pointer :: cpro_rom, cpro_visl, imasfl double precision, dimension(:), pointer :: Re, cf, den double precision, dimension(:,:),pointer :: vel, A_vel save flown, flown1, dp allocate(lstfac(nfac)) allocate(lstelt(nfabor)) chaine1 = 'Inlet' chaine2 = 'RodWall or Wall' !================ ! Variables and flow initialisation !================ !! Target velocity uref = 15.49d0 !! Use a reference density value, since in this case not a variable value rho_ref = 36.00621d0 !! If mometum equation is not solved, then return. Note, to specify the source !! term for the momentum equation, this routine is only called for ivar=iu. if (ivar.ne.iu) return !ipcrom = ipproc(irom) !gives density at the current time step if variable !iflmas = ipprof(ifluma(iu)) !no idea what this does !!=========================== !!Source term calculation !!=========================== !! Impose periodicity using source term method flowref = uref*rho_ref !Mass flux flowint = flown !Mass flux of the previous time step flown = 0.0d0 surf = 0.0d0 cptfac = 0 call getfac(chaine1, nlfac, lstfac) call getfbr(chaine2, nlelt, lstelt) !! This gets the interger for a field associated with given key call field_get_key_int(ivarfl(iw),kimasf,iflmas) !! Access mass flux call field_get_val_s(iflmas, imasfl) call field_get_val_v(ivarfl(iu),vel) !! First loop to determine surface and current flow do ielt = 1, nlfac ifac = lstfac(ielt) !! This sums up surface area surf = surf + surfan(ifac) !! This sums up the total mass flux ! flown = flown + imasfl(ifac) !! Mass flux can be positive or negative, should use the absolute value here flown = flown + abs(imasfl(ifac)) cptfac = cptfac + 1 enddo !! Checking codes !if(ntcabs==30)then ! open(10,file='test.dat') ! do ielt=1, nlfac ! ifac = lstfac(ielt) ! write(10,*)abs(imasfl(ifac)) ! temp = temp+abs(imasfl(ifac)) ! end do ! write(10,*)'temp = ', temp ! write(10,*)'temp/surf = ', temp/surf ! close(10) !end if !! Incase of parralesim sum all the totals from the different processors if (irangp.ge.0) then call parsom(surf) call parsom(flown) call parcpt(cptfac) endif !flown = abs(flown) / surf flown = flown/surf if (ntcabs .eq. 1) then dp = 0.0d0 flown = flowref flown1 = flown else flown1 = flowint endif !dp = dp + (2.0d0*(flown-flowref) + (flown - flown1))/(2.0d0*dt(1)) !! Check masss flow balances dp = dp + (2.0d0*(flown-flowref)-(flown1-flowref))/(2.0d0*dt(1)) !!momentum source on all cells do iel = 1, ncel value = - volume(iel) * dp !write(*,*) dp crvexp(3,iel) = value enddo !!momentum source on one layer cells !do ielt = 1, nlfac ! ifac = lstfac(ielt) ! do ii = 1, 2 ! iel = ifacel(ii, ifac) ! value = - volume(iel) * dp ! !write(*,*) dp ! crvexp(3,iel) = value ! enddo !enddo !do iel = 1, ncel ! if(ntcabs==100)write(*,*)iel, crvexp(3,iel) !end do !! Impose the momentum source of wall frictional force do ielt = 1, nlelt ifac = lstelt(ielt) iel = ifabor(ifac) !! General test ! crvexp(3,iel) = -1.d-5*vel(3,iel)-volume(iel)*dp !! Use global bulk velocity ! tau_w = 0.25d0*f*0.5d0*rho_ref*(uref**2) !! Use local velocity ! tau_w = 0.25d0*f*0.5d0*rho_ref*(vel(3,iel)**2) !! Use local bulk velocity (subchannel velocity) ! tau_w = 0.25d0*f*0.5d0*rho_ref*(A_vel(3,iel)**2) !! Calculate the total source term together with periodic sources ! crvexp(3,iel) = -tau_w*surfbn(ifac)-volume(iel)*dp !! Use for checking results ! if(ntcabs==1)then ! open(10,file='user_monitor.dat',status='replace') ! close(10) ! end if ! open(10,file='user_monitor.dat',status='old', position='append') ! if(iel==780)write(10,*) ttcabs, tau_w, volume(iel)*dp, vel(3,iel) ! close(10) ! write(*,*)xlomlg end do if(mod(ntcabs,100)==0)then write(nfecra,*) '==FLOW RATE CONTROL==' write(nfecra,*) 'current time step number=',ntcabs write(nfecra,*) 'number of inlet faces = ', cptfac write(nfecra,*) 'inlet surface = ', surf write(nfecra,*) 'flow at time n at inlet = ', flown , 'flow forcing = ', value write(nfecra,*) 'pressure gradient = ',-dp write(nfecra,*) 'ref. flow', flowref write(nfecra,*) 'flow at time n-1 at inlet = ', flown1 end if !-------- ! Formats !-------- !< [format_1] !1000 format(' User source termes for variable ',A8,/) !< [format_1] !---- ! End !---- !< [deallocate_1] !! Deallocate the temporary array deallocate(lstfac,lstelt) !< [deallocate_1] return end subroutine ustsnv !=============================================================================== !=============================================================================== ! Purpose: ! ------- !> User subroutine. !> \brief Additional right-hand side source terms for scalar equations (user !> scalars and specific physics scalars). !> !> !> Usage !> ----- !> The routine is called for each scalar, user or specific physisc. It is !> therefore necessary to test the value of the scalar number iscal to separate !> the treatments of the different scalars (if (iscal.eq.p) then ....). !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equation solved by the code for a scalar f is: !> !> \f[ \rho*volume*\frac{df}{dt} + .... = crvimp*f + crvexp \f] !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.[scal]/s, where [scal] is the unit of the scalar !> - crvimp is expressed in kg/s !> !> !> The crvexp and 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 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 \ref getfbr in the routine !> \ref cs_user_boundary_conditions. !> WARNING: If scalar is the temperature, the resulting equation !> solved by the code is: !> !> rho*Cp*volume*dT/dt + .... = crvimp*T + crvexp !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in W !> - crvimp is expressed in W/K !> !> !> STEP SOURCE TERMS !>=================== !> In case of a complex, non-linear source term, say F(f), for scalar f, the !> easiest method is to implement the source term explicitely. !> !> df/dt = .... + F(f(n)) !> where f(n) is the value of f at time tn, the beginning of the time step. !> !> This yields : !> crvexp = volume*F(f(n)) !> crvimp = 0 !> !> However, if the source term is potentially steep, this fully explicit !> method will probably generate instabilities. It is therefore wiser to !> partially implicit the term by writing: !> !> df/dt = .... + dF/df*f(n+1) - dF/df*f(n) + F(f(n)) !> !> This yields: !> crvexp = volume*( F(f(n)) - dF/df*f(n) ) !> crvimp = volume*dF/df !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! 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] iscal index number of the current scalar !> \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 !> \param[in] (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 !> \param[in] 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 ustssc & ( nvar , nscal , ncepdp , ncesmp , & iscal , & 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 iscal integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) double precision dt(ncelet) double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(ncelet), crvimp(ncelet) ! Local variables integer iel double precision, dimension(:), pointer :: imasfl double precision, dimension(:,:), pointer :: vel, A_vel double precision, dimension(:), pointer :: surf_rod, cvol integer, allocatable, dimension(:) :: lstelt integer nlelt, ielt double precision ubulk, surf, xx, sumvol ,rho_ref integer iflmas, ifac, nn character*80 :: chaine double precision :: heat_flux=342000. !< [loc_var_scal] !=============================================================================== !! Impose heat sink for periodicity if (iscal.ne.iscalt) return uref = 15.49d0 rho_ref = 36.00621d0 call field_get_val_v(ivarfl(iu),vel) allocate(lstelt(nfabor)) chaine = 'RodWall' call getfbr(chaine, nlelt, lstelt) surf = 0.d0 do ielt = 1, nlelt ifac = lstelt(ielt) surf = surf + surfbn(ifac) end do sumvol = 0.d0 ubulk = 1.d-15 do iel = 1, ncel sumvol = sumvol+volume(iel) ubulk = ubulk+ vel(3,iel)*volume(iel) end do !! Incase of parralesim sum all the totals from the different processors if (irangp.ge.0) then call parsom(surf) call parsom(sumvol) call parsom(ubulk) end if !! Check area of the boundary faces if(ntcabs==1)then ! surf = 0.d0 ! do ielt = 1, nlelt ! ifac = lstelt(ielt) ! surf = surf + surfbn(ifac) !add total surface area ! write(*,*) surfbn(ifac), ifac ! end do ! write(*,*)surf,nlelt ! write(*,*)sumvol end if !! Calculate heat sink do iel = 1, ncel crvimp(iel) = 0.d0 !! From developers at Code_Saturne user's meeting (05-06/June/2018) !! Check details in '/Working/The University of Sheffield/Other_meeting/ !! Code_Saturne_meeting_05_06_2018/periodicity.doc' (Dell working laptop) crvexp(iel) = -heat_flux*surf*volume(iel)*vel(3,iel)/ubulk !! Heat sink is calculated based on the subchannel rather than the whole domain !! (This is not correct) ! crvexp(iel) = -heat_flux*surf_rod(iel)*volume(iel)/cvol(iel)& ! * vel(3,iel)/(A_vel(3,iel)+1.d-15) end do !-------- ! Formats !-------- !---- ! End !---- return end subroutine ustssc !=============================================================================== !=============================================================================== ! Purpose: ! ------- !> User subroutine. !> \brief Additional right-hand side source terms for vectorial equations !> (user vectors and specific physics vectors). !> !> !> Usage !> ----- !> The routine is called for each vector, user or specific physisc. It is !> therefore necessary to test the value of the vector number iscal to separate !> the treatments of the different vectors (if (iscal.eq.p) then ....). !> !> The additional source term is decomposed into an explicit part (crvexp) and !> an implicit part (crvimp) that must be provided here. !> The resulting equation solved by the code for a vector f is: !> !> \f[ \rho*volume*\frac{d\vect{f}}{dt} + .... = \tens{crvimp}*\vect{f} + !> \vect{crvexp} \f] !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in kg.[scal]/s, where [scal] is vector unit !> - crvimp is expressed in kg/s !> !> !> The crvexp and 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 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 \ref getfbr in the routine !> \ref cs_user_boundary_conditions. !> WARNING: If scalar is the temperature, the resulting equation !> solved by the code is: !> !> rho*Cp*volume*dT/dt + .... = crvimp*T + crvexp !> !> !> Note that crvexp and crvimp are defined after the Finite Volume integration !> over the cells, so they include the "volume" term. More precisely: !> - crvexp is expressed in W !> - crvimp is expressed in W/K !> !> !> STEP SOURCE TERMS !>=================== !> In case of a complex, non-linear source term, say F(f), for scalar f, the !> easiest method is to implement the source term explicitely. !> !> df/dt = .... + F(f(n)) !> where f(n) is the value of f at time tn, the beginning of the time step. !> !> This yields : !> crvexp = volume*F(f(n)) !> crvimp = 0 !> !> However, if the source term is potentially steep, this fully explicit !> method will probably generate instabilities. It is therefore wiser to !> partially implicit the term by writing: !> !> df/dt = .... + dF/df*f(n+1) - dF/df*f(n) + F(f(n)) !> !> This yields: !> crvexp = volume*( F(f(n)) - dF/df*f(n) ) !> crvimp = volume*dF/df !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! 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] iscal index number of the current scalar !> \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 !> \param[in] (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 !> \param[in] 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 ustsvv & ( nvar , nscal , ncepdp , ncesmp , & iscal , & 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 iscal 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, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for variable ',A8,/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine ustsvv !=============================================================================== !=============================================================================== ! Purpose: ! ------- !> \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 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 !=============================================================================== 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(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(ncelet), crvimp(ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for turbulence model',/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine cs_user_turbulence_source_terms !=============================================================================== !=============================================================================== ! Purpose: ! ------- !> \brief Additional right-hand side source terms for turbulence models and !>irijco =1 !> !> \section cs_user_rij_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_terms2 & ( nvar , nscal , ncepdp , ncesmp , & f_id , & icepdc , icetsm , itypsm , & ckupdc , smacel , & crvexp , crvimp ) !=============================================================================== !=============================================================================== ! Module files !=============================================================================== 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 !=============================================================================== 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(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(6,ncelet), crvimp(6,6,ncelet) ! Local variables integer, allocatable, dimension(:) :: lstelt !=============================================================================== !=============================================================================== ! 1. Initialization !=============================================================================== ! Allocate a temporary array for cells selection allocate(lstelt(ncel)) !-------- ! Formats !-------- 1000 format(' User source terms for turbulence model',/) !---- ! End !---- ! Deallocate the temporary array deallocate(lstelt) return end subroutine cs_user_turbulence_source_terms2