!------------------------------------------------------------------------------- ! Code_Saturne version 2.0.0-rc1 ! -------------------------- ! This file is part of the Code_Saturne Kernel, element of the ! Code_Saturne CFD tool. ! Copyright (C) 1998-2009 EDF S.A., France ! contact: saturne-support@edf.fr ! The Code_Saturne Kernel 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. ! The Code_Saturne Kernel 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 the Code_Saturne Kernel; if not, write to the ! Free Software Foundation, Inc., ! 51 Franklin St, Fifth Floor, ! Boston, MA 02110-1301 USA !------------------------------------------------------------------------------- subroutine ustsns & !================ ( idbia0 , idbra0 , & ndim , ncelet , ncel , nfac , nfabor , nfml , nprfml , & nnod , lndfac , lndfbr , ncelbr , & nvar , nscal , nphas , ncepdp , ncesmp , & nideve , nrdeve , nituse , nrtuse , & ivar , iphas , & ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , & ipnfac , nodfac , ipnfbr , nodfbr , & icepdc , icetsm , itypsm , & idevel , ituser , ia , & xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , & dt , rtpa , propce , propfa , propfb , & coefa , coefb , ckupdc , smacel , & crvexp , crvimp , & dam , xam , & w1 , w2 , w3 , w4 , w5 , w6 , & rdevel , rtuser , ra ) !=============================================================================== ! Purpose: ! ------- ! User subroutine. ! Additional right-hand side source terms for velocity components equation ! (Navier-Stokes) ! ! Usage ! ----- ! The routine is called for each velocity component. It is therefore necessary ! to test the value of the variable ivar to separate the treatments of the ! components iu(iphas), iv(iphas) or iw(iphas). ! ! 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 velocity component u is: ! ! rho*volume*du/dt + .... = crvimp*u + 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 kg.m/s2 ! - 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 getfbr in the routine ! usclim. !------------------------------------------------------------------------------- ! Arguments !__________________.____._____.________________________________________________. ! name !type!mode ! role ! !__________________!____!_____!________________________________________________! ! idbia0 ! i ! <-- ! number of first free position in ia ! ! idbra0 ! i ! <-- ! number of first free position in ra ! ! ndim ! i ! <-- ! spatial dimension ! ! ncelet ! i ! <-- ! number of extended (real + ghost) cells ! ! ncel ! i ! <-- ! number of cells ! ! nfac ! i ! <-- ! number of interior faces ! ! nfabor ! i ! <-- ! number of boundary faces ! ! nfml ! i ! <-- ! number of families (group classes) ! ! nprfml ! i ! <-- ! number of properties per family (group class) ! ! nnod ! i ! <-- ! number of vertices ! ! lndfac ! i ! <-- ! size of nodfac indexed array ! ! lndfbr ! i ! <-- ! size of nodfbr indexed array ! ! ncelbr ! i ! <-- ! number of cells with faces on boundary ! ! nvar ! i ! <-- ! total number of variables ! ! nscal ! i ! <-- ! total number of scalars ! ! nphas ! i ! <-- ! number of phases ! ! ncepdp ! i ! <-- ! number of cells with head loss terms ! ! ncssmp ! i ! <-- ! number of cells with mass source terms ! ! nideve, nrdeve ! i ! <-- ! sizes of idevel and rdevel arrays ! ! nituse, nrtuse ! i ! <-- ! sizes of ituser and rtuser arrays ! ! ivar ! i ! <-- ! index number of the current variable ! ! iphas ! i ! <-- ! index number of the current phase ! ! ifacel(2, nfac) ! ia ! <-- ! interior faces -> cells connectivity ! ! ifabor(nfabor) ! ia ! <-- ! boundary faces -> cells connectivity ! ! ifmfbr(nfabor) ! ia ! <-- ! boundary face family numbers ! ! ifmcel(ncelet) ! ia ! <-- ! cell family numbers ! ! iprfml ! ia ! <-- ! property numbers per family ! ! (nfml, nprfml) ! ! ! ! ! maxelt ! i ! <-- ! max number of cells and faces (int/boundary) ! ! lstelt(maxelt) ! ia ! --- ! work array ! ! ipnfac(nfac+1) ! ia ! <-- ! interior faces -> vertices index (optional) ! ! nodfac(lndfac) ! ia ! <-- ! interior faces -> vertices list (optional) ! ! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional) ! ! nodfbr(lndfbr) ! ia ! <-- ! boundary faces -> vertices list (optional) ! ! icepdc(ncepdp) ! ia ! <-- ! index number of cells with head loss terms ! ! icetsm(ncesmp) ! ia ! <-- ! index number of cells with mass source terms ! ! itypsm ! ia ! <-- ! type of mass source term for each variable ! ! (ncesmp,nvar) ! ! ! (see ustsma) ! ! idevel(nideve) ! ia ! <-- ! integer work array for temporary developpement ! ! ituser(nituse ! ia ! <-- ! user-reserved integer work array ! ! ia(*) ! ia ! --- ! main integer work array ! ! xyzcen ! ra ! <-- ! cell centers ! ! (ndim, ncelet) ! ! ! ! ! surfac ! ra ! <-- ! interior faces surface vectors ! ! (ndim, nfac) ! ! ! ! ! surfbo ! ra ! <-- ! boundary faces surface vectors ! ! (ndim, nfavor) ! ! ! ! ! cdgfac ! ra ! <-- ! interior faces centers of gravity ! ! (ndim, nfac) ! ! ! ! ! cdgfbo ! ra ! <-- ! boundary faces centers of gravity ! ! (ndim, nfabor) ! ! ! ! ! xyznod ! ra ! <-- ! vertex coordinates (optional) ! ! (ndim, nnod) ! ! ! ! ! volume(ncelet) ! ra ! <-- ! cell volumes ! ! 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 ! ! coefa, coefb ! ra ! <-- ! boundary conditions ! ! (nfabor, *) ! ! ! ! ! ckupdc(ncepdp,6) ! ra ! <-- ! head loss coefficient ! ! smacel ! ra ! <-- ! value associated to each variable in the mass ! ! (ncesmp,nvar) ! ! ! source terms or mass rate (see ustsma) ! ! crvexp ! ra ! --> ! explicit part of the source term ! ! crvimp ! ra ! --> ! implicit part of the source term ! ! dam(ncelet) ! ra ! --- ! work array ! ! xam(nfac,2) ! ra ! --- ! work array ! ! w1,2,3,4,5,6 ! ra ! --- ! work arrays ! ! (ncelet) ! ! ! (computation of pressure gradient) ! ! rdevel(nrdeve) ! ra ! <-> ! real work array for temporary developpement ! ! rtuser(nituse ! ra ! <-- ! user-reserved real work array ! ! ra(*) ! ra ! --- ! main real work array ! !__________________!____!_____!________________________________________________! ! 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 !=============================================================================== implicit none !=============================================================================== ! Common blocks !=============================================================================== include "paramx.h" include "pointe.h" include "numvar.h" include "entsor.h" include "optcal.h" include "cstphy.h" include "parall.h" include "period.h" include "inccha.h" !=============================================================================== ! Arguments integer idbia0 , idbra0 integer ndim , ncelet , ncel , nfac , nfabor integer nfml , nprfml integer nnod , lndfac , lndfbr , ncelbr integer nvar , nscal , nphas integer ncepdp , ncesmp integer nideve , nrdeve , nituse , nrtuse integer ivar , iphas integer ifacel(2,nfac) , ifabor(nfabor) integer ifmfbr(nfabor) , ifmcel(ncelet) integer iprfml(nfml,nprfml), maxelt, lstelt(maxelt) integer ipnfac(nfac+1), nodfac(lndfac) integer ipnfbr(nfabor+1), nodfbr(lndfbr) integer icepdc(ncepdp) integer icetsm(ncesmp), itypsm(ncesmp,nvar) integer idevel(nideve), ituser(nituse), ia(*) double precision xyzcen(ndim,ncelet) double precision surfac(ndim,nfac), surfbo(ndim,nfabor) double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor) double precision xyznod(ndim,nnod), volume(ncelet) double precision dt(ncelet), rtpa(ncelet,*) double precision propce(ncelet,*) double precision propfa(nfac,*), propfb(nfabor,*) double precision coefa(nfabor,*), coefb(nfabor,*) double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar) double precision crvexp(ncelet), crvimp(ncelet) double precision dam(ncelet ),xam(nfac ,2) double precision w1(ncelet),w2(ncelet),w3(ncelet) double precision w4(ncelet),w5(ncelet),w6(ncelet) double precision rdevel(nrdeve), rtuser(nrtuse), ra(*) ! Local variables character*80 chaine integer idebia, idebra integer iel, ipcrom, ipp, iutile double precision ckp, qdm integer ipass integer ifac,iifac,iunit,iuiph,iflmas double precision debitint,surf double precision debitref double precision surfx,surfy,surfz,surfan data ipass /0/ save ipass !=============================================================================== !=============================================================================== ! 1. INITIALISATION !=============================================================================== idebia = idbia0 idebra = idbra0 ipp = ipprtp(ivar) iuiph = iu(iphas) if(ivar.eq.iuiph) then if(iwarni(ivar).ge.1) then chaine = nomvar(ipp) write(nfecra,1000) chaine(1:8) endif ipcrom = ipproc(irom (iphas)) iflmas = ipprof(ifluma(iuiph)) !=============================================================================== ! 2. EXEMPLE DE TERME SOURCE ARBITRAIRE, POUR LA VARIABLE U : ! S = A * U + B ! APPARAISSANT DANS LES EQUATIONS SOUS LA FORME : ! RHO VOLUME D(U)/Dt = VOLUME*S ! CE TERME A UNE PARTIE QU'ON VEUT IMPLICITER : A ! ET UNE PARTIE QU'ON VA TRAITER EN EXPLICITE : B ! ICI PAR EXEMPLE ARBITRAIREMENT : ! A = - RHO CKP ! B = QDM ! AVEC ! CKP = 10.D0 [1/s ] (COEFFICIENT DE RAPPEL SUR U) ! QDM = 100.D0 [kg/(m2 s2)] (PRODUCTION DE QUANTITE DE MOUVEMENT ! PAR UNITE DE VOLUME ET PAR UNITE DE TEMPS) ! ON A ALORS ! CRVIMP(IEL) = VOLUME(IEL)* A = - VOLUME(IEL) (RHO CKP ) ! CRVEXP(IEL) = VOLUME(IEL)* B = VOLUME(IEL) (QDM ) ! ON REMPLIT CI-DESSOUS CRVIMP ET CRVEXP CORRESPONDANTS. !=============================================================================== iunit = impusr(10) debitref = uref(iphas)*propce(1,ipcrom) debitint = debitn debitn = 0.0d0 surf = 0.0d0 iifac = 0 do ifac=1,nfac surfx = surfac(1,ifac) surfy = surfac(2,ifac) surfz = surfac(3,ifac) surfan = sqrt(surfx**2 + surfy**2 + surfz**2) if(abs(cdgfac(3,ifac)-(-11.1d0)).lt.0.000001d0) then iifac = iifac + 1 debitn = debitn + propfa(ifac,iflmas) surf = surf + surfan endif enddo debitn = abs(debitn) / surf if(isuite.eq.0.and.ipass.eq.1) then dp = 0.0d0 debitn = debitref debitn1 = debitn ! ELSEIF(ISUITE.NE.0.AND.IPASS.EQ.1) THEN ! OPEN(FILE='fluxes',UNIT=IUNIT) ! REWIND(IUNIT) ! READ(IUNIT,*) DEBITN1 ! READ(IUNIT,*) DP ! CLOSE(IUNIT) else debitn1 = debitint endif ! OPEN(FILE='fluxes',UNIT=IUNIT) ! REWIND(IUNIT) ! WRITE(IUNIT,*) DEBITN ! WRITE(IUNIT,*) DP ! CLOSE(IUNIT) ! FORCE DE RAPPEL ! VERIFICATIONS ! Attention ceci n'est valable qu'en pas de temps constant dp = dp + (2.0d0*(debitn-debitref)-(debitn1-debitref))/(2.0d0*dt(1)) WRITE(NFECRA,*) 'NOMBRE DE FACES A L ENTREE = ',IIFAC WRITE(NFECRA,*) 'SURFACE DE L ENTREE = ',SURF WRITE(NFECRA,*) 'DEBIT AU TEMPS N A L ENTREE = ',DEBITN WRITE(NFECRA,*) 'DEBIT AU TEMPS N-1 A L ENTREE = ',DEBITN1 WRITE(NFECRA,*) 'FORCAGE DU DEBIT = ',DP do iel = 1, ncel crvexp(iel) = - volume(iel) * dp enddo endif !-------- ! FORMATS !-------- 1000 format(' TERMES SOURCES UTILISATEURS POUR LA VARIABLE ',A8,/) !---- ! FIN !---- return end