!------------------------------------------------------------------------------- ! Code_Saturne version 2.0.0-beta2 ! -------------------------- ! This file is part of the Code_Saturne Kernel, element of the ! Code_Saturne CFD tool. ! Copyright (C) 1998-2008 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 ) !=============================================================================== ! FONCTION : ! ---------- ! ROUTINE UTILSATEUR ON PRECISE LES TERMES SOURCES UTILISATEURS ! POUR LA COMPOSANTE DE VITESSE IVAR ! ON RESOUT RHO*VOLUME*D(VAR)/DT = CRVIMP*VAR + CRVEXP ! ON FOURNIT ICI CRVIMP ET CRVEXP (ILS CONTIENNENT RHO*VOLUME) ! CRVEXP en kg variable/s : ! ex : pour la vitesse kg m/s2 ! CRVIMP en kg /s : ! VEILLER A UTILISER UN CRVIMP NEGATIF ! (ON IMPLICITERA CRVIMP ! IE SUR LA DIAGONALE DE LA MATRICE, LE CODE AJOUTERA : ! MAX(-CRVIMP,0) EN SCHEMA STANDARD EN TEMPS ! -CRVIMP SI LES TERMES SOURCES SONT A L'ORDRE 2 ! CES TABLEAUX SONT INITIALISES A ZERO AVANT APPEL A CE SOUS ! PROGRAMME ET AJOUTES ENSUITE AUX TABLEAUX PRIS EN COMPTE ! POUR LA RESOLUTION ! EN CAS D'ORDRE 2 DEMANDE SUR LES TERMES SOURCES, ON DOIT ! FOURNIR CRVEXP A L'INSTANT N (IL SERA EXTRAPOLE) ET ! CRVIMP A L'INSTANT N+1/2 (IL EST DANS LA MATRICE, ! ON LE SUPPOSE NEGATIF) ! L'identification des cellules concernees peut s'appuyer sur la ! commande GETCEL. ! GETCEL(CHAINE,NLELT,LSTELT) : ! - CHAINE est une chaine de caractere fournie par l'utilisateur ! qui donne les criteres de selection ! - NLTELT est renvoye par la commande. C'est un entier qui ! correspond au nombre de cellules trouveees repondant au ! critere ! - LSTELT est renvoye par la commande. C'est un tableau d'entiers ! de taille NLTELT donnant la liste des cellules trouvees ! repondant au critere. ! CHAINE peut etre constitue de : ! - references de couleurs (ex. : 1, 8, 26, ... ! - references de groupes (ex. : entrees, groupe1, ...) ! - criteres geometriques (ex. X<0.1, Y>=0.25, ...) ! Ces criteres peuvent etre combines par des operateurs logiques ! (AND et OR) et des parentheses ! ex. : '1 AND (groupe2 OR groupe3) AND Y<1' permettra de recuperer ! les cellules de couleur 1, appartenant aux groupes 'groupe2' ! ou 'groupe3' et de coordonnee Y inferieure a 1. !------------------------------------------------------------------------------- !ARGU ARGUMENTS !__________________.____._____.________________________________________________. ! nom !type!mode ! role ! !__________________!____!_____!________________________________________________! ! idbia0 ! e ! <-- ! numero de la 1ere case libre dans ia ! ! idbra0 ! e ! <-- ! numero de la 1ere case libre dans ra ! ! ndim ! e ! <-- ! dimension de l'espace ! ! ncelet ! e ! <-- ! nombre d'elements halo compris ! ! ncel ! e ! <-- ! nombre d'elements actifs ! ! nfac ! e ! <-- ! nombre de faces internes ! ! nfabor ! e ! <-- ! nombre de faces de bord ! ! nfml ! e ! <-- ! nombre de familles d entites ! ! nprfml ! e ! <-- ! nombre de proprietese des familles ! ! nnod ! e ! <-- ! nombre de sommets ! ! lndfac ! e ! <-- ! longueur du tableau nodfac (optionnel ! ! lndfbr ! e ! <-- ! longueur du tableau nodfbr (optionnel ! ! ncelbr ! e ! <-- ! nombre d'elements ayant au moins une ! ! ! ! ! face de bord ! ! nvar ! e ! <-- ! nombre total de variables ! ! nscal ! e ! <-- ! nombre total de scalaires ! ! nphas ! e ! <-- ! nombre de phases ! ! ncepdp ! e ! <-- ! nombre de cellules avec pdc ! ! ncesmp ! e ! <-- ! nombre de cellules a source de masse ! ! nideve nrdeve ! e ! <-- ! longueur de idevel rdevel ! ! nituse nrtuse ! e ! <-- ! longueur de ituser rtuser ! ! ivar ! e ! <-- ! numero de variable ! ! iphas ! e ! <-- ! numero de phase ! ! ifacel ! te ! <-- ! elements voisins d'une face interne ! ! (2, nfac) ! ! ! ! ! ifabor ! te ! <-- ! element voisin d'une face de bord ! ! (nfabor) ! ! ! ! ! ifmfbr ! te ! <-- ! numero de famille d'une face de bord ! ! (nfabor) ! ! ! ! ! ifmcel ! te ! <-- ! numero de famille d'une cellule ! ! (ncelet) ! ! ! ! ! iprfml ! te ! <-- ! proprietes d'une famille ! ! nfml ,nprfml ! ! ! ! ! maxelt ! e ! <-- ! nb max d'elements (cell,fac,fbr) ! ! lstelt(maxelt) te ! --- ! tableau de travail ! ! ipnfac ! te ! <-- ! position du premier noeud de chaque ! ! (lndfac) ! ! ! face interne dans nodfac (optionnel) ! ! nodfac ! te ! <-- ! connectivite faces internes/noeuds ! ! (nfac+1) ! ! ! (optionnel) ! ! ipnfbr ! te ! <-- ! position du premier noeud de chaque ! ! (lndfbr) ! ! ! face de bord dans nodfbr (optionnel) ! ! nodfbr ! te ! <-- ! connectivite faces de bord/noeuds ! ! (nfabor+1) ! ! ! (optionnel) ! ! icepdc(ncelet ! te ! <-- ! numero des ncepdp cellules avec pdc ! ! icetsm(ncesmp ! te ! <-- ! numero des cellules a source de masse ! ! itypsm ! te ! <-- ! type de source de masse pour les ! ! (ncesmp,nvar) ! ! ! variables (cf. ustsma) ! ! idevel(nideve ! te ! <-- ! tab entier complementaire developemt ! ! ituser(nituse ! te ! <-- ! tab entier complementaire utilisateur ! ! ia(*) ! tr ! --- ! macro tableau entier ! ! xyzcen ! tr ! <-- ! point associes aux volumes de control ! ! (ndim,ncelet ! ! ! ! ! surfac ! tr ! <-- ! vecteur surface des faces internes ! ! (ndim,nfac) ! ! ! ! ! surfbo ! tr ! <-- ! vecteur surface des faces de bord ! ! (ndim,nfabor) ! ! ! ! ! cdgfac ! tr ! <-- ! centre de gravite des faces internes ! ! (ndim,nfac) ! ! ! ! ! cdgfbo ! tr ! <-- ! centre de gravite des faces de bord ! ! (ndim,nfabor) ! ! ! ! ! xyznod ! tr ! <-- ! coordonnes des noeuds (optionnel) ! ! (ndim,nnod) ! ! ! ! ! volume ! tr ! <-- ! volume d'un des ncelet elements ! ! (ncelet ! ! ! ! ! dt(ncelet) ! tr ! <-- ! pas de temps ! ! rtpa ! tr ! <-- ! variables de calcul au centre des ! ! (ncelet,*) ! ! ! cellules (instant prec) ! ! propce ! tr ! <-- ! proprietes physiques au centre des ! ! (ncelet,*) ! ! ! cellules ! ! propfa ! tr ! <-- ! proprietes physiques au centre des ! ! (nfac,*) ! ! ! faces internes ! ! propfb ! tr ! <-- ! proprietes physiques au centre des ! ! (nfabor,*) ! ! ! faces de bord ! ! coefa, coefb ! tr ! <-- ! conditions aux limites aux ! ! (nfabor,*) ! ! ! faces de bord ! ! ckupdc ! tr ! <-- ! tableau de travail pour pdc ! ! (ncepdp,6) ! ! ! ! ! smacel ! tr ! <-- ! valeur des variables associee a la ! ! (ncesmp,* ) ! ! ! source de masse ! ! ! ! ! pour ivar=ipr, smacel=flux de masse ! ! crvexp(ncelet ! tr ! --> ! tableau de travail pour part explicit ! ! crvimp(ncelet ! tr ! --> ! tableau de travail pour part implicit ! ! dam(ncelet ! tr ! --- ! tableau de travail pour matrice ! ! xam(nfac,*) ! tr ! --- ! tableau de travail pour matrice ! ! w1...6(ncelet ! tr ! --- ! tableau de travail ! ! rdevel(nrdeve ! tr ! <-- ! tab reel complementaire developemt ! ! rtuser(nrtuse ! tr ! <-- ! tab reel complementaire utilisateur ! ! ra(*) ! tr ! --- ! macro tableau reel ! !__________________!____!_____!________________________________________________! ! TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU) ! L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL) ! MODE : <-- donnee, --> resultat, <-> Donnee modifiee ! --- tableau de travail !------------------------------------------------------------------------------- !=============================================================================== implicit none !=============================================================================== ! DONNEES EN COMMON !=============================================================================== 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(*) ! VARIABLES LOCALES character*80 chaine integer idebia, idebra integer iel, ipcrom, ipp, iutile double precision ckp, qdm integer ipass, nlelt, ilelt 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) WRITE(NFECRA,*) 'DEBIT de ref = ',debitref debitint = debitn debitn = 0.0d0 surf = 0.0d0 iifac = 0 call getfac ('entree', nlelt, lstelt) do ilelt = 1, nlelt ifac=lstelt(ilelt) iifac = iifac + 1 debitn = debitn + propfa(ifac,iflmas) enddo if (irangp.ge.0) then call parcpt (iifac) call parsom (debitn) endif ! ici, 0.000028 est la surface d'entrŽe debitn = abs(debitn) / 0.0012 if(isuite.eq.0.and.ipass.eq.1) then dp = 0.0d0 debitn = debitref debitn1 = debitn else debitn1 = debitint endif ! 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)) if (ntcabs.le.ntpabs+2) then dp = - 6591.0d0 endif WRITE(NFECRA,*) 'NOMBRE DE FACES A L ENTREE (ustsns)= ',IIFAC 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 (Pa/m) = ',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