c@a c@versb C----------------------------------------------------------------------- C CVERS Code_Saturne version 1.3 C ------------------------ C C This file is part of the Code_Saturne Kernel, element of the C Code_Saturne CFD tool. C C Copyright (C) 1998-2007 EDF S.A., France C C contact: saturne-support@edf.fr C C The Code_Saturne Kernel is free software; you can redistribute it C and/or modify it under the terms of the GNU General Public License C as published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C The Code_Saturne Kernel is distributed in the hope that it will be C useful, but WITHOUT ANY WARRANTY; without even the implied warranty C of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public License C along with the Code_Saturne Kernel; if not, write to the C Free Software Foundation, Inc., C 51 Franklin St, Fifth Floor, C Boston, MA 02110-1301 USA C C----------------------------------------------------------------------- c@verse SUBROUTINE USELCL C ***************** C -------------------------------------------------------------- & ( IDBIA0 , IDBRA0 , & NDIM , NCELET , NCEL , NFAC , NFABOR , NFML , NPRFML , & NNOD , LNDFAC , LNDFBR , NCELBR , & NVAR , NSCAL , NPHAS , & NIDEVE , NRDEVE , NITUSE , NRTUSE , & IFACEL , IFABOR , IFMFBR , IFMCEL , IPRFML , & IPNFAC , NODFAC , IPNFBR , NODFBR , & ICODCL , ITRIFB , ITYPFB , IZFPPP , & IDEVEL , ITUSER , IA , & XYZCEN , SURFAC , SURFBO , CDGFAC , CDGFBO , XYZNOD , VOLUME , & DT , RTP , RTPA , PROPCE , PROPFA , PROPFB , & COEFA , COEFB , RCODCL , & W1 , W2 , W3 , W4 , W5 , W6 , COEFU , & RDEVEL , RTUSER , RA ) C -------------------------------------------------------------- C*********************************************************************** C FONCTION : C -------- c@foncb CFONC CFONC ROUTINE UTILISATEUR POUR PHYSIQUE PARTICULIERE CFONC MODULE ELECTRIQUE CFONC (Effet Joule, Arc Electrique, Conduction ionique) CFONC REMPLISSAGE DU TABLEAU DE CONDITIONS AUX LIMITES CFONC (ICODCL,RCODCL) POUR LES VARIABLES INCONNUES CFONC PENDANT DE USCLIM.F CFONC CFONC CFONC CFONC CE SOUS PROGRAMME UTILISATEUR EST OBLIGATOIRE CFONC ============================================= CFONC CFONC CFONC CFONC INTRODUCTION CFONC ============ CFONC CFONC On donne ici les conditions aux limites par face de bord. CFONC CFONC On balaye les faces de bord et selon un critere on affecte tel CFONC ou tel type de conditions aux limites. Dans l'exemple donne CFONC ci dessous, c'est la couleur (propriete 1 de la famille) CFONC qui permet de distinguer les differents types de bord. On CFONC aurait pu aussi travailler avec les coordonnees du centre CFONC des faces, mais "c'eut ete moins pratique". CFONC CFONC CFONC TYPE DE CONDITIONS AUX LIMITES CFONC ============================== CFONC CFONC On peut affecter les conditions aux limites de deux manieres. CFONC CFONC CFONC Pour les conditions "standard" : CFONC -------------------------------- CFONC CFONC (entree, sortie libre, paroi, symetrie), on donne un code CFONC stocke dans le tableau ITYPFB (dimensionne au nombre de CFONC faces de bord,nombre de phases). Ce code sera ensuite CFONC utilise par un sous-programme non utilisateur pour affecter CFONC les conditions correspondantes (les scalaires, en CFONC particulier, recevront alors les conditions de la phase a CFONC laquelle ils sont associes). Ainsi : CFONC CFONC Type de bord Entree Sortie libre Symetrie Paroi CFONC Code IENTRE ISOLIB ISYMET IPAROI CFONC CFONC Les entiers IENTRE, ISOLIB, ISYMET, IPAROI, CFONC sont affectes par ailleurs (include param.h). Leur valeur CFONC est superieure ou egale a 1 et CFONC inferieure ou egale a NTYPMX (valeur fixee dans paramx.h) CFONC CFONC CFONC CFONC En outre, il faut donner certaines valeurs : CFONC CFONC CFONC - Entree (plus precisement entree/sortie a debit impose, car le debit CFONC peut etre impose sortant) : CFONC CFONC -> Conditions de Dirichlet sur les variables CFONC autres que la pression obligatoire si le flux est entrant, CFONC optionnelle si le flux est sortant (le code affecte flux nul CFONC si aucun Dirichlet n'est specifie) ; ainsi CFONC en face IFAC, sur la variable IVAR : RCODCL(IFAC,IVAR,1) CFONC CFONC CFONC - Paroi : (= solide impermeable, avec frottement) CFONC CFONC -> Valeur de la vitesse de paroi defilante s'il y a lieu CFONC en face IFAC, RCODCL(IFAC,IU,1) CFONC RCODCL(IFAC,IV,1) CFONC RCODCL(IFAC,IW,1) CFONC -> Code specifique et valeur de la temperature imposee CFONC en paroi s'il y a lieu : CFONC en face IFAC, ICODCL(IFAC,IVAR) = 5 CFONC RCODCL(IFAC,IVAR,1) = Temperature imposee CFONC -> Code specifique et valeur du flux imposee CFONC en paroi s'il y a lieu : CFONC en face IFAC, ICODCL(IFAC,IVAR) = 3 CFONC RCODCL(IFAC,IVAR,3) = Flux impose CFONC = CFONC Noter que la condition par defaut pour les scalaires CFONC (hors k et epsilon) est un Neumann homogene CFONC CFONC CFONC CFONC - Symetrie (= paroi impermeable avec glissement) : CFONC CFONC -> Rien a preciser CFONC CFONC CFONC - Sortie libre (plus precisement entree/sortie libre a pression imposee) CFONC CFONC -> Rien a preciser pour la pression et la vitesse CFONC Pour les scalaires et grandeurs turbulentes, une valeur de Dirichlet CFONC peut etre specifiee de maniere optionnelle. Le comportement est le CFONC suivant : CFONC * la pression est toujours traitee en Dirichlet CFONC * si flux de masse entrant : CFONC on retient la vitesse a l'infini CFONC condition de Dirichlet pour les scalaires et grandeurs CFONC turbulentes (ou flux nul si l'utilisateur n'a pas CFONC specifie de Dirichlet) CFONC si flux de masse sortant : CFONC on impose un flux nul sur la vitesse, les scalaires et CFONC les grandeurs turbulentes CFONC CFONC Noter que la pression sera recalee a P0 CFONC sur la premiere face de sortie libre trouvee CFONC CFONC CFONC Pour les conditions "non standard" : CFONC ------------------------------------ CFONC CFONC Autres que (entree, sortie libre, paroi, symetrie), on donne CFONC - d'une part, pour chaque face : CFONC -> une valeur de ITYPFB admissible CFONC ie superieure ou egale a 1 et inferieure ou egale a CFONC NTYPMX (voir sa valeur dans paramx.h). CFONC Les valeurs predefinies dans paramx.h CFONC IENTRE, ISOLIB, ISYMET, IPAROI, sont dans cet CFONC intervalle et il est preferable de ne pas affecter CFONC inconsidrement et par hasard a ITYPFB la valeur CFONC d'un de ces entiers. Pour eviter cela, on peut CFONC utiliser IINDEF si l'on souhaite eviter de verifier CFONC les valeurs dans paramx.h. IINDEF est une valeur CFONC admissible a laquelle n'est attachee aucune condition CFONC limite predefinie. CFONC Noter que le tableau ITYPFB est CFONC reinitialise a chaque pas de temps a la valeur CFONC non admissible 0. Si on oublie de modifier ITYPFB pour CFONC une face, le code s'arretera. CFONC CFONC - et d'autre part pour chaque face et chaque variable : CFONC -> un code ICODCL(IFAC,IVAR) CFONC -> trois reels RCODCL(IFAC,IVAR,1) CFONC RCODCL(IFAC,IVAR,2) CFONC RCODCL(IFAC,IVAR,3) CFONC La valeur de ICODCL est prise parmi les suivantes : CFONC 1 : Dirichlet (utilisable pour toute variable) CFONC 3 : Neumann (utilisable pour toute variable) CFONC 4 : Symetrie (utilisable uniquement pour la vitesse et CFONC les composantes du tenseur Rij) CFONC 5 : Paroi (utilisable pour toute variable sauf la CFONC pression) CFONC 9 : Sortie libre (utilisable uniquement pour la vitesse) CFONC Les valeurs des 3 reels RCODCL sont les suivantes CFONC RCODCL(IFAC,IVAR,1) : CFONC Dirichlet sur la variable si ICODCL(IFAC,IVAR)= 1 CFONC valeur en paroi (defilmnt, temp) si ICODCL(IFAC,IVAR)= 5 CFONC La dimension de RCODCL(IFAC,IVAR,1) est celle de la CFONC variable resolue : ex U (vitesse en m/s), CFONC T (temperature en degres) CFONC H (enthalpie en J/kg) CFONC F (scalaire passif en -) CFONC RCODCL(IFAC,IVAR,2) : CFONC coefficient d'echange "exterieur" (entre la valeur CFONC imposee et la valeur au bord du domaine) CFONC RINFIN = infini par defaut CFONC Pour les vitesses U, en kg/(m2 s) : CFONC RCODCL(IFAC,IVAR,2) = (VISCL+VISCT) / D CFONC Pour la pression P, en s/m : CFONC RCODCL(IFAC,IVAR,2) = DT / D CFONC Pour les temperatures T, en Watt/(m2 degres) : CFONC RCODCL(IFAC,IVAR,2) = CP*(VISCLS+VISCT/SIGMAS) / D CFONC Pour les enthalpies H, en kg /(m2 s) : CFONC RCODCL(IFAC,IVAR,2) = (VISCLS+VISCT/SIGMAS) / D CFONC Pour les autres scalaires F en : CFONC RCODCL(IFAC,IVAR,2) = (VISCLS+VISCT/SIGMAS) / D CFONC (D a la dimension d'une distance en m) CFONC CFONC RCODCL(IFAC,IVAR,3) : CFONC Densite de flux (< 0 si gain, n normale orientee vers l'exterieur) CFONC si ICODCL(IFAC,IVAR)= 3 CFONC Pour les vitesses U, en kg/(m s2) = J : CFONC RCODCL(IFAC,IVAR,3) = -(VISCL+VISCT) * (GRAD U).n CFONC Pour la pression P, en kg/(m2 s) : CFONC RCODCL(IFAC,IVAR,3) = -DT * (GRAD P).n CFONC Pour les temperatures T, en Watt/m2 : CFONC RCODCL(IFAC,IVAR,3) =-CP*(VISCLS+VISCT/SIGMAS) * (GRAD T).n CFONC Pour les enthalpies H, en Watt/m2 : CFONC RCODCL(IFAC,IVAR,3) = -(VISCLS+VISCT/SIGMAS) * (GRAD H).n CFONC Pour les autres scalaires F en : CFONC RCODCL(IFAC,IVAR,3) = -(VISCLS+VISCT/SIGMAS) * (GRAD F).n CFONC CFONC CFONC CFONC Noter bien que si l'utilisateur affecte une valeur a ITYPFB CFONC parmi IENTRE, ISOLIB, ISYMET, IPAROI, CFONC et qu'il ne modifie pas ICODCL (valeur nulle par defaut), CFONC c'est ITYPFB qui imposera la condition limite. CFONC CFONC Par contre, si l'utilisateur impose CFONC ICODCL(IFAC,IVAR) (non nul), CFONC ce sont alors les valeurs de RCODCL qu'il aura fournies CFONC qui sont retenues pour la face et la variable consideree CFONC (s'il ne precise pas RCODCL, ce sont les valeurs CFONC par defaut qui sont retenues pour la face et CFONC la variable consideree soit : CFONC RCODCL(IFAC,IVAR,1) = 0.D0 CFONC RCODCL(IFAC,IVAR,2) = RINFIN CFONC RCODCL(IFAC,IVAR,3) = 0.D0) CFONC En particulier, on peut par exemple CFONC -> donner ITYPFB(IFAC,IPHAS) = IPAROI CFONC qui impose les conditions de paroi par defaut pour toutes CFONC les variables sur la face IFAC, CFONC -> et preciser EN OUTRE pour la variable IVAR sur cette CFONC face IFAC des conditions paarticulieres en imposant CFONC ICODCL(IFAC,IVAR) et les 3 RCODCL. CFONC CFONC CFONC L'utilisateur peut egalement affecter a ITYPFB une valeur CFONC non egale a IENTRE, ISOLIB, ISYMET, IPAROI, IINDEF CFONC mais superieure ou egale a 1 et inferieure ou egale a CFONC NTYPMX (voir les valeurs dans param.h) pour reperer CFONC des groupes de couleurs dans d'autres sous programmes CFONC qui lui seraient propres et ou ITYPFB serait disponible. CFONC Dans ce cas cependant, il faudra CFONC imposer les conditions limites en donnant des valeurs a CFONC ICODCL et aux trois RCODCL (puisque la valeur de ITYPFB CFONC ne sera pas predefinie dans le code). CFONC CFONC CFONC REGLES DE COHERENCE CFONC =================== CFONC CFONC Quelques regles de coherence entre les codes ICODCL CFONC des variables pour les conditions non standard : CFONC CFONC Les codes des vitesses doivent etre identiques CFONC Les codes des Rij doivent etre identiques CFONC Si code (vitesse ou Rij) = 4 CFONC il faut code (vitesse et Rij) = 4 CFONC Si code (vitesse ou turbulence) = 5 CFONC il faut code (vitesse et turbulence) = 5 CFONC Si code scalaire (hormis pression ou fluctuations) = 5 CFONC il faut code vitesse = 5 CFONC CFONC CFONC REMARQUES CFONC ========== CFONC CFONC Attention : pour imposer un flux (non nul) sur les Rij, CFONC la viscosite a prendre en compte est VISCL CFONC meme si VISCT existe (VISCT=RHO CMU K2/EPSILON) CFONC CFONC CFONC On dispose du tableau de tri des faces de bord au pas CFONC de temps precedent (sauf au premier pas de temps, ou CFONC ITRIFB n'a pas ete renseigne). CFONC Le tableau du type des faces de bord ITYPFB a ete CFONC reinitialise avant d'entrer dans le sous programme. CFONC CFONC CFONC CFONC Noter comment acceder a certaines variables : CFONC CFONC Valeurs aux cellules CFONC Soit IEL = IFABOR(IFAC) CFONC CFONC * Masse vol phase IPHAS, cellule IEL : CFONC PROPCE(IEL ,IPPROC(IROM (IPHAS))) CFONC * Viscosite moleculaire dynamique phase IPHAS, cellule IEL : CFONC PROPCE(IEL ,IPPROC(IVISCL(IPHAS))) CFONC * Viscosite turbulente dynamique phase IPHAS, cellule IEL : CFONC PROPCE(IEL ,IPPROC(IVISCT(IPHAS))) CFONC * Chaleur specifique phase IPHAS, cellule IEL : CFONC PROPCE(IEL ,IPPROC(ICP (IPHAS))) CFONC * Diffusivite lambda scalaire ISCAL, cellule IEL : CFONC PROPCE(IEL ,IPPROC(IVISLS(ISCAL))) CFONC CFONC Valeurs aux faces de bord CFONC CFONC * Masse vol phase IPHAS, face de bord IFAC : CFONC PROPFB(IFAC,IPPROB(IROM (IPHAS))) CFONC * Flux de masse relatif a la variable IVAR , face de bord IFAC : CFONC (i.e. le flux de masse servant a la convection de IVAR) CFONC PROPFB(IFAC,IPPROB(IFLUMA(IVAR ))) CFONC * Pour les autres grandeurs a la face de bord IFAC : CFONC prendre pour approximation la valeur dans la cellule IEL CFONC adjacente i.e. comme plus haut avec IEL = IFABOR(IFAC). CFONC c@fonce C----------------------------------------------------------------------- C ARGUMENTS c@argub CARGU .______________.____._____.______________________________________. CARGU ! NOM !TYPE!MODE ! ROLE ! CARGU !______________!____!_____!______________________________________! CARGU ! IDBIA0 ! E ! -> ! NUMERO DE LA 1ERE CASE LIBRE DANS IA ! CARGU ! IDBRA0 ! E ! -> ! NUMERO DE LA 1ERE CASE LIBRE DANS RA ! CARGU ! NDIM ! E ! -> ! DIMENSION DE L'ESPACE ! CARGU ! NCELET ! E ! -> ! NOMBRE D'ELEMENTS HALO COMPRIS ! CARGU ! NCEL ! E ! -> ! NOMBRE D'ELEMENTS ACTIFS ! CARGU ! NFAC ! E ! -> ! NOMBRE DE FACES INTERNES ! CARGU ! NFABOR ! E ! -> ! NOMBRE DE FACES DE BORD ! CARGU ! NFML ! E ! -> ! NOMBRE DE FAMILLES D ENTITES ! CARGU ! NPRFML ! E ! -> ! NOMBRE DE PROPRIETESE DES FAMILLES ! CARGU ! NNOD ! E ! -> ! NOMBRE DE SOMMETS ! CARGU ! LNDFAC ! E ! -> ! LONGUEUR DU TABLEAU NODFAC (OPTIONNEL! CARGU ! LNDFBR ! E ! -> ! LONGUEUR DU TABLEAU NODFBR (OPTIONNEL! CARGU ! NCELBR ! E ! -> ! NOMBRE D'ELEMENTS AYANT AU MOINS UNE ! CARGU ! ! ! ! FACE DE BORD ! CARGU ! NVAR ! E ! -> ! NOMBRE TOTAL DE VARIABLES ! CARGU ! NSCAL ! E ! -> ! NOMBRE TOTAL DE SCALAIRES ! CARGU ! NPHAS ! E ! -> ! NOMBRE DE PHASES ! CARGU ! NIDEVE NRDEVE! E ! -> ! LONGUEUR DE IDEVEL RDEVEL ! CARGU ! NITUSE NRTUSE! E ! -> ! LONGUEUR DE ITUSER RTUSER ! CARGU ! IFACEL ! TE ! -> ! ELEMENTS VOISINS D'UNE FACE INTERNE ! CARGU ! (2, NFAC) ! ! ! ! CARGU ! IFABOR ! TE ! -> ! ELEMENT VOISIN D'UNE FACE DE BORD ! CARGU ! (NFABOR) ! ! ! ! CARGU ! IFMFBR ! TE ! -> ! NUMERO DE FAMILLE D'UNE FACE DE BORD ! CARGU ! (NFABOR) ! ! ! ! CARGU ! IFMCEL ! TE ! -> ! NUMERO DE FAMILLE D'UNE CELLULE ! CARGU ! (NCELET) ! ! ! ! CARGU ! IPRFML ! TE ! -> ! PROPRIETES D'UNE FAMILLE ! CARGU ! NFML ,NPRFML! ! ! ! CARGU ! IPNFAC ! TE ! -> ! POSITION DU PREMIER NOEUD DE CHAQUE ! CARGU ! (LNDFAC) ! ! ! FACE INTERNE DANS NODFAC (OPTIONNEL)! CARGU ! NODFAC ! TE ! -> ! CONNECTIVITE FACES INTERNES/NOEUDS ! CARGU ! (NFAC+1) ! ! ! (OPTIONNEL) ! CARGU ! IPNFBR ! TE ! -> ! POSITION DU PREMIER NOEUD DE CHAQUE ! CARGU ! (LNDFBR) ! ! ! FACE DE BORD DANS NODFBR (OPTIONNEL)! CARGU ! NODFBR ! TE ! -> ! CONNECTIVITE FACES DE BORD/NOEUDS ! CARGU ! (NFABOR+1) ! ! ! (OPTIONNEL) ! CARGU ! ICODCL ! TE ! <- ! CODE DE CONDITION LIMITES AUX FACES ! CARGU ! (NFABOR,NVAR! ! ! DE BORD ! CARGU ! ! ! ! = 1 -> DIRICHLET ! CARGU ! ! ! ! = 3 -> DENSITE DE FLUX ! CARGU ! ! ! ! = 4 -> GLISSEMT ET U.n=0 (VITESSE) ! CARGU ! ! ! ! = 5 -> FROTTEMT ET U.n=0 (VITESSE) ! CARGU ! ! ! ! = 9 -> ENTREE/SORTIE LIBRE (VITESSE! CARGU ! ! ! ! ENTRANTE EVENTUELLE BLOQUEE ! CARGU ! ITRIFB(NFABOR! TE ! -> ! INDIRECTION POUR TRI DES FACES DE BRD! CARGU ! NPHAS )! ! ! ! CARGU ! ITYPFB(NFABOR! TE ! <- ! TYPE DES FACES DE BORD ! CARGU ! NPHAS )! ! ! ! CARGU ! IZFPPP ! TE ! <- ! NUMERO DE ZONE DE LA FACE DE BORD ! CARGU ! (NFABOR) ! ! ! POUR LE MODULE PHYS. PART. ! CARGU ! IDEVEL(NIDEVE! TE ! <-> ! TAB ENTIER COMPLEMENTAIRE DEVELOPEMT ! CARGU ! ITUSER(NITUSE! TE ! <-> ! TAB ENTIER COMPLEMENTAIRE UTILISATEUR! CARGU ! IA(*) ! TR ! - ! MACRO TABLEAU ENTIER ! CARGU ! XYZCEN ! TR ! -> ! POINT ASSOCIES AUX VOLUMES DE CONTROL! CARGU ! (NDIM,NCELET ! ! ! ! CARGU ! SURFAC ! TR ! -> ! VECTEUR SURFACE DES FACES INTERNES ! CARGU ! (NDIM,NFAC) ! ! ! ! CARGU ! SURFBO ! TR ! -> ! VECTEUR SURFACE DES FACES DE BORD ! CARGU ! (NDIM,NFABOR)! ! ! ! CARGU ! CDGFAC ! TR ! -> ! CENTRE DE GRAVITE DES FACES INTERNES ! CARGU ! (NDIM,NFAC) ! ! ! ! CARGU ! CDGFBO ! TR ! -> ! CENTRE DE GRAVITE DES FACES DE BORD ! CARGU ! (NDIM,NFABOR)! ! ! ! CARGU ! XYZNOD ! TR ! -> ! COORDONNES DES NOEUDS (OPTIONNEL) ! CARGU ! (NDIM,NNOD) ! ! ! ! CARGU ! VOLUME ! TR ! -> ! VOLUME D'UN DES NCELET ELEMENTS ! CARGU ! (NCELET ! ! ! ! CARGU ! DT(NCELET) ! TR ! -> ! PAS DE TEMPS ! CARGU ! RTP, RTPA ! TR ! -> ! VARIABLES DE CALCUL AU CENTRE DES ! CARGU ! (NCELET,*) ! ! ! CELLULES (INSTANT COURANT OU PREC)! CARGU ! PROPCE ! TR ! -> ! PROPRIETES PHYSIQUES AU CENTRE DES ! CARGU ! (NCELET,*) ! ! ! CELLULES ! CARGU ! PROPFA ! TR ! -> ! PROPRIETES PHYSIQUES AU CENTRE DES ! CARGU ! (NFAC,*) ! ! ! FACES INTERNES ! CARGU ! PROPFB ! TR ! -> ! PROPRIETES PHYSIQUES AU CENTRE DES ! CARGU ! (NFABOR,*) ! ! ! FACES DE BORD ! CARGU ! COEFA, COEFB ! TR ! <-> ! CONDITIONS AUX LIMITES AUX ! CARGU ! (NFABOR,*) ! ! ! FACES DE BORD ! CARGU ! RCODCL ! TR ! <- ! VALEUR DES CONDITIONS AUX LIMITES ! CARGU ! (NFABOR,NVAR! ! ! AUX FACES DE BORD ! CARGU ! ! ! ! RCODCL(1) = VALEUR DU DIRICHLET ! CARGU ! ! ! ! RCODCL(2) = VALEUR DU COEF. D'ECHANGE! CARGU ! ! ! ! EXT. (INFINIE SI PAS D'ECHANGE) ! CARGU ! ! ! ! RCODCL(3) = VALEUR DE LA DENSITE DE ! CARGU ! ! ! ! FLUX (NEGATIF SI GAIN) W/m2 ! CARGU ! ! ! ! POUR LES VITESSES (VISTL+VISCT)*GRADU! CARGU ! ! ! ! POUR LA PRESSION DT*GRADP! CARGU ! ! ! ! POUR LES SCALAIRES ! CARGU ! ! ! ! CP*(VISCLS+VISCT/SIGMAS)*GRADT! CARGU ! W1,2,3,4,5,6 ! TR ! - ! TABLEAUX DE TRAVAIL ! CARGU ! (NCELET ! ! ! (CALCUL DU GRADIENT DE PRESSION) ! CARGU ! COEFU ! TR ! - ! TAB DE TRAV ! CARGU ! (NFABOR,3) ! ! ! (CALCUL DU GRADIENT DE PRESSION) ! CARGU ! RDEVEL(NRDEVE! TR ! <-> ! TAB REEL COMPLEMENTAIRE DEVELOPEMT ! CARGU ! RTUSER(NRTUSE! TR ! <-> ! TAB REEL COMPLEMENTAIRE UTILISATEUR ! CARGU ! RA(*) ! TR ! - ! MACRO TABLEAU REEL ! CARGU !______________!____!_____!______________________________________! c@argue C c@commb CCOMM COMMONS CCOMM .______________.____._____.______________________________________. CCOMM ! NOM !TYPE!MODE ! ROLE ! CCOMM !______________!____!_____!______________________________________! CCOMM !______________!____!_____!______________________________________! c@comme C C TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU) C L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL) C MODE : -> DONNEE, <- RESULTAT, <-> DONNEE MODIFIEE, C - TABLEAU DE TRAVAIL C*********************************************************************** C IMPLICIT NONE C C*********************************************************************** C DONNEES EN COMMON C*********************************************************************** C INCLUDE "paramx.h" INCLUDE "numvar.h" INCLUDE "optcal.h" INCLUDE "cstphy.h" INCLUDE "cstnum.h" INCLUDE "entsor.h" INCLUDE "ppppar.h" INCLUDE "ppthch.h" INCLUDE "ppincl.h" INCLUDE "elincl.h" C C*********************************************************************** C C ARGUMENTS C INTEGER IDBIA0 , IDBRA0 INTEGER NDIM , NCELET , NCEL , NFAC , NFABOR INTEGER NFML , NPRFML INTEGER NNOD , LNDFAC , LNDFBR , NCELBR INTEGER NVAR , NSCAL , NPHAS INTEGER NIDEVE , NRDEVE , NITUSE , NRTUSE C INTEGER IFACEL(2,NFAC) , IFABOR(NFABOR) INTEGER IFMFBR(NFABOR) , IFMCEL(NCELET) INTEGER IPRFML(NFML,NPRFML) INTEGER IPNFAC(NFAC+1), NODFAC(LNDFAC) INTEGER IPNFBR(NFABOR+1), NODFBR(LNDFBR) INTEGER ICODCL(NFABOR,NVAR) INTEGER ITRIFB(NFABOR,NPHAS), ITYPFB(NFABOR,NPHAS) INTEGER IZFPPP(NFABOR) INTEGER IDEVEL(NIDEVE), ITUSER(NITUSE), IA(*) C 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), RTP(NCELET,*), RTPA(NCELET,*) DOUBLE PRECISION PROPCE(NCELET,*) DOUBLE PRECISION PROPFA(NFAC,*), PROPFB(NFABOR,*) DOUBLE PRECISION COEFA(NFABOR,*), COEFB(NFABOR,*) DOUBLE PRECISION RCODCL(NFABOR,NVAR,3) DOUBLE PRECISION W1(NCELET),W2(NCELET),W3(NCELET) DOUBLE PRECISION W4(NCELET),W5(NCELET),W6(NCELET) DOUBLE PRECISION COEFU(NFABOR,NDIM) DOUBLE PRECISION RDEVEL(NRDEVE), RTUSER(NRTUSE), RA(*) C C VARIABLES LOCALES C INTEGER IDEBIA, IDEBRA INTEGER IFAC, II, ICOUL, IPHAS , IEL INTEGER IFML, IDIM INTEGER IZONE,IESP DOUBLE PRECISION UREF2, D2S3 DOUBLE PRECISION RHOMOY, DHY, USTAR2 DOUBLE PRECISION XKENT, XEENT DOUBLE PRECISION Z1 , Z2 C C*********************************************************************** C C TEST_A_ENLEVER_POUR_UTILISER_LE_SOUS_PROGRAMME_DEBUT C======================================================================= C 0. CE TEST PERMET A L'UTILISATEUR D'ETRE CERTAIN QUE C'EST C SA VERSION DU SOUS PROGRAMME QUI EST UTILISEE C ET NON CELLE DE LA BIBLIOTHEQUE C======================================================================= C IF(1.EQ.1) THEN WRITE(NFECRA,9001) CALL CSEXIT (1) ENDIF C 9001 FORMAT( &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/, &'@ @@ ATTENTION : ARRET LORS DE L''ENTREE DES COND. LIM. ',/, &'@ ********* ',/, &'@ MODULE ELECTRIQUE ',/, &'@ ',/, &'@ LE SOUS-PROGRAMME UTILISATEUR uselcl DOIT ETRE COMPLETE',/, &'@ ',/, &'@ Ce sous-programme utilisateur permet de definir les ',/, &'@ conditions aux limites. Il est indispensable. ',/, &'@ ',/, &'@ Le calcul ne sera pas execute. ',/, &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/) C C C TEST_A_ENLEVER_POUR_UTILISER_LE_SOUS_PROGRAMME_FIN C C======================================================================= C 1. INITIALISATIONS C C======================================================================= C IDEBIA = IDBIA0 IDEBRA = IDBRA0 C D2S3 = 2.D0/3.D0 C C======================================================================= C 2. REMPLISSAGE DU TABLEAU DES CONDITIONS LIMITES C ON BOUCLE SUR LES FACES DE BORD C ON DETERMINE LA FAMILLE ET SES PROPRIETES C ON IMPOSE LA CONDITION LIMITE C C IMPOSER ICI LES CONDITIONS LIMITES SUR LES FACES DE BORD C C======================================================================= C C BOUCLE SUR LES FACES DE BORD C DO IFAC = 1, NFABOR C C FAMILLE ET PREMIERE PROPRIETE C C En general, la couleur des faces est la premiere propriete : C ICOUL = IPRFML(IFML,1). C C Dans le cas plus specifique de maillages contenant des groupes C (ensembles dont l'intersection peut ne pas etre vide), C il est possible que l'on obtienne plusieurs proprietes par C famille. Ainsi, une face appartenant aux groupes 3 et 5 et C portant la couleur 8 sera dans une famille dont les C proprietes seront -3, -5 et 8 (ordre arbitraire). C Il faut donc, si l'on veut savoir si la face IFAC porte C la couleur 8, balayer toutes les proprietes de la C famille : C IFML = IFMFBR(IFAC ) C DO IPROP = 1, NPRFML C IF(IPRFML(IFML,IPROP).EQ.8) THEN C La face porte la couleur 8 ... C ENDIF C ENDDO C C IFML = IFMFBR(IFAC ) ICOUL = IPRFML(IFML,1) IPHAS = 1 C C C ELEMENT ADJACENT A LA FACE DE BORD C C C --- On impose en couleur 1 une entree ; exemple de Cathode C ====================================================== C IF( ICOUL.EQ.1) THEN C ITYPFB(IFAC,IPHAS) = IENTRE C C - Numero de zone (choix du numero de couleur par exemple) IZONE = ICOUL C C - Reperage de la zone a laquelle appartient la face IZFPPP(IFAC) = IZONE C RCODCL(IFAC,IU(IPHAS),1) = 0.D0 RCODCL(IFAC,IV(IPHAS),1) = 0.D0 RCODCL(IFAC,IW(IPHAS),1) = 0.D0 C C Turbulence C C (ITYTUR est un indicateur qui vaut ITURB/10) IF (ITYTUR(IPHAS).EQ.2 .OR. ITYTUR(IPHAS).EQ.3 & .OR. ITURB(IPHAS).EQ.50 .OR. ITURB(IPHAS).EQ.60) THEN C UREF2 = RCODCL(IFAC,IU(IPHAS),1)**2 & +RCODCL(IFAC,IV(IPHAS),1)**2 & +RCODCL(IFAC,IW(IPHAS),1)**2 UREF2 = MAX(UREF2,1.D-12) C C C Exemple de turbulence calculee a partir C de formules valables pour une conduite C C On veillera a specifier le diametre hydraulique C adapte a l'entree courante. C C On s'attachera egalement a utiliser si besoin une formule C plus precise pour la viscosite dynamique utilisee dans le C calcul du nombre de Reynolds (en particulier, lorsqu'elle C est variable, il peut etre utile de reprendre ici la loi C imposee dans USPHYV. On utilise ici par defaut la valeur C VISCL0 donnee dans USINI1 C En ce qui concerne la masse volumique, on dispose directement C de sa valeur aux faces de bord (ROMB) et c'est celle que C utilise donc ici (elle est en particulier coherente avec C le traitement implante dans USPHYV, en cas de masse C volumique variable) C C Diametre hydraulique DHY = 0.075D0 C C Calcul de la vitesse de frottement au carre (USTAR2) C et de k et epsilon en entree (XKENT et XEENT) a partir C de lois standards en conduite circulaire C (leur initialisation est inutile mais plus propre) RHOMOY = PROPFB(IFAC,IPPROB(IROM(IPHAS))) USTAR2 = 0.D0 XKENT = EPZERO XEENT = EPZERO C CALL KEENDB C =========== & ( UREF2, DHY, RHOMOY, VISCL0(IPHAS), CMU, XKAPPA, & USTAR2, XKENT, XEENT ) C IF (ITYTUR(IPHAS).EQ.2) THEN C RCODCL(IFAC,IK(IPHAS),1) = XKENT RCODCL(IFAC,IEP(IPHAS),1) = XEENT C ELSEIF(ITYTUR(IPHAS).EQ.3) THEN C RCODCL(IFAC,IR11(IPHAS),1) = D2S3*XKENT RCODCL(IFAC,IR22(IPHAS),1) = D2S3*XKENT RCODCL(IFAC,IR33(IPHAS),1) = D2S3*XKENT RCODCL(IFAC,IR12(IPHAS),1) = 0.D0 RCODCL(IFAC,IR13(IPHAS),1) = 0.D0 RCODCL(IFAC,IR23(IPHAS),1) = 0.D0 RCODCL(IFAC,IEP(IPHAS),1) = XEENT C ELSEIF (ITURB(IPHAS).EQ.50) THEN C RCODCL(IFAC,IK(IPHAS),1) = XKENT RCODCL(IFAC,IEP(IPHAS),1) = XEENT RCODCL(IFAC,IPHI(IPHAS),1) = D2S3 RCODCL(IFAC,IFB(IPHAS),1) = 0.D0 C ELSEIF (ITURB(IPHAS).EQ.60) THEN C RCODCL(IFAC,IK(IPHAS),1) = XKENT RCODCL(IFAC,IOMG(IPHAS),1) = XEENT/CMU/XKENT C ENDIF C ENDIF C C --- On traite les scalaires C C Enthalpie en J/kg C II = IHM ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = 1.D6 C C Potentiel electrique reel impose a 0. volts (exemple de Cathode en arc) C II = IPOTR ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = 0.D0 C C Fraction massique des (N-1) constituants C IF ( NGAZG .GT. 1 ) THEN DO IESP=1,NGAZG-1 II = IYCOEL(IESP) ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = 0.D0 ENDDO ENDIF C C Specifique Version Effet Joule : C C Potentiel Imaginaire impose a 0 C IF ( IPPMOD(IELJOU).GE. 2 ) THEN II = IPOTI ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = 0.D0 ENDIF C C Specifique Version Arc Electrique : C C Potentiel vecteur : Flux nul C IF ( IPPMOD(IELARC).GE.2 ) THEN DO IDIM= 1,NDIMVE II = IPOTVA(IDIM) ICODCL(IFAC,ISCA(II)) = 3 RCODCL(IFAC,ISCA(II),3) = 0.D0 ENDDO ENDIF C C --- On impose en couleur 5 une entree/sortie ; C ====================================== exemple d'Electrode en Joule C ============================ C ELSE IF( ICOUL.EQ.5) THEN C C SORTIE : FLUX NUL VITESSE ET TEMPERATURE, PRESSION IMPOSEE C Noter que la pression sera recalee a P0 C sur la premiere face de sortie libre (ISOLIB) C ITYPFB(IFAC,IPHAS) = ISOLIB C C - Numero de zone (choix du numero de couleur par exemple) IZONE = ICOUL C C - Reperage de la zone a laquelle appartient la face IZFPPP(IFAC) = IZONE C C --- On traite les scalaires rattaches a la phase courante C C Enthalpie en J/kg (Par defaut flux nul avec ISOLIB) C Rien a faire C C Fraction massique des (N-1) constituants (Par defaut flux nul avec ISOLIB) C Rien a faire C C Specifique Version Effet Joule : C C En effet Joule, C si l'on souhaite faire un calcul en recalant les conditions C aux limites (utiliser IELCOR=1 dans useli1) C pour atteindre la valeur de la puissance PUISIM C (a imposer dans useli1 en Ampere.Volt) C on multiplie la condition limite initiale sur le potentiel C reel (et sur le potentiel imaginaire s'il est pris en C compte) par le coefficient COEJOU. C COEJOU est determine automatiquement pour que la puissance C dissipee par effet Joule (partie reelle et partie C imaginaire si besoin) soit PUISIM C au debut du calcul, COEJOU vaut 1 ; COEJOU est transmis dans C les fichiers suites. C C Si on ne souhaite pas faire un calcul avec recalage, on impose C directement une valeur adaptee. C IF ( IPPMOD(IELJOU).GE. 1 ) THEN II = IPOTR ICODCL(IFAC,ISCA(II)) = 1 IF(IELCOR.EQ.1) THEN RCODCL(IFAC,ISCA(II),1) = 500.D0*COEJOU ELSE RCODCL(IFAC,ISCA(II),1) = 500.D0 ENDIF ENDIF C IF ( IPPMOD(IELJOU).GE. 2 ) THEN II = IPOTI ICODCL(IFAC,ISCA(II)) = 1 IF(IELCOR.EQ.1) THEN RCODCL(IFAC,ISCA(II),1) = SQRT(3.D0)*500.D0*COEJOU ELSE RCODCL(IFAC,ISCA(II),1) = SQRT(3.D0)*500.D0 ENDIF ENDIF C C C --- On impose en couleur 2 une entree/sortie ; C ============================== exemple d'Anode en arc electrique C ================================= C ELSE IF( ICOUL.EQ.2) THEN C C SORTIE : FLUX NUL VITESSE ET TEMPERATURE, PRESSION IMPOSEE C Noter que la pression sera recalee a P0 C sur la premiere face de sortie libre (ISOLIB) C ITYPFB(IFAC,IPHAS) = ISOLIB C C - Numero de zone (choix du numero de couleur par exemple) IZONE = ICOUL C C - Reperage de la zone a laquelle appartient la face IZFPPP(IFAC) = IZONE C C --- On traite les scalaires rattaches a la phase courante C C Enthalpie en J/kg (Par defaut flux nul avec ISOLIB) C Rien a faire C C Potentiel electrique reel C C En arc electrique, C si l'on souhaite faire un calcul en recalant le potentiel C de l'anode (utiliser IELCOR=1 dans useli1) C pour atteindre la valeur du courant COUIMP C (a imposer dans useli1 en Amperes) C on utilise alors la valeur DPOT comme condition limite C DPOT est en effet automatiquement adaptee par le calcul C pour que (j.E Volume/DPOT) = COUIMP C (initialiser DPOT dans useli1 en Volts avec une valeur C representative de la difference de potentiel imposee) C C Si on ne souhaite pas faire un calcul avec recalage, on impose C directement une valeur adaptee au cas C (par exemple, ici 1000 Volts ). C II = IPOTR ICODCL(IFAC,ISCA(II)) = 1 C IF ( IPPMOD(IELARC).GE.1 .AND. IELCOR .EQ.1) THEN RCODCL(IFAC,ISCA(II),1) = DPOT ELSE RCODCL(IFAC,ISCA(II),1) = 1000.D0 ENDIF C C C Fraction massique des (N-1) constituants (Par defaut flux nul avec ISOLIB) C C Specifique Version Arc Electrique : C Potentiel vecteur : flux nul (par defaut) C C C C --- On impose en couleur 3 une paroi C ================================ C ELSE IF( ICOUL.EQ.3) THEN C C PAROI : DEBIT NUL (FLUX NUL POUR LA PRESSION) C FROTTEMENT POUR LES VITESSES (+GRANDEURS TURB) C FLUX NUL SUR LES SCALAIRES C C Pour un calcul arc electrique 3D, on cale le potentiel C vecteur avec une condition de Dirichlet issue des valeurs C du potentiel vecteur au pas de temps precedent C dans une zone de paroi choisie C Par defaut, ailleurs, un flux nul s'applique (paroi isolee). C ITYPFB(IFAC,IPHAS) = IPAROI C C - Numero de zone (choix du numero de couleur par exemple) IZONE = ICOUL C C - Reperage de la zone a laquelle appartient la face IZFPPP(IFAC) = IZONE C IF ( IPPMOD(IELARC).GE.2 ) THEN IF ( CDGFBO(1,IFAC) .LE. 2.249D-2 .OR. & CDGFBO(1,IFAC) .GE. 2.249D-2 .OR. & CDGFBO(3,IFAC) .LE. -2.249D-2 .OR. & CDGFBO(3,IFAC) .GE. 2.249D-2 ) THEN IEL = IFABOR(IFAC) DO IDIM = 1, NDIMVE II = IPOTVA(IDIM) ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = RTPA(IEL,ISCA(II)) ENDDO ENDIF ENDIF C C --- On impose en couleur 50 : anode avec claquage C ============================================= ELSE IF( ICOUL.EQ.51) THEN C ITYPFB(IFAC,IPHAS) = IPAROI C C - Numero de zone (choix du numero de couleur par exemple) IZONE = ICOUL C C - Reperage de la zone a laquelle appartient la face IZFPPP(IFAC) = IZONE C C ---- Enthalpie (J/kg ) : coef echange impose C II=IHM ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = 2.D4 RCODCL(IFAC,ISCA(II),2) = 1.D5 C C Potentiel electrique reel C II = IPOTR ICODCL(IFAC,ISCA(II)) = 1 C IF ( IPPMOD(IELARC).GE.1 .AND. IELCOR .EQ.1) THEN RCODCL(IFAC,ISCA(II),1) = DPOT ELSE RCODCL(IFAC,ISCA(II),1) = 100.D0 ENDIF C C Si CLAQUAGE : a adapter en fonction du cas et du C sous-programme USELRC C IF ( IPPMOD(IELARC).GE.1 .AND. IELCOR .EQ.1) THEN IF(ICLAQ.EQ.1 .AND. NTCABS.LE.NTDCLA+30) THEN C Z1 = ZCLAQ - 2.D-4 IF(Z1.LE.0.D0) Z1 = 0.D0 Z2 = ZCLAQ + 2.D-4 IF(Z2.GE.2.D-2) Z2 = 2.D-2 C IF( CDGFBO(3,IFAC).GE.Z1 .AND. & CDGFBO(3,IFAC).LE.Z2 ) THEN ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = DPOT ELSE ICODCL(IFAC,ISCA(II)) = 3 RCODCL(IFAC,ISCA(II),3) = 0.D0 ENDIF ENDIF ENDIF C C Potentiel vecteur : Flux nul C IF ( IPPMOD(IELARC).GE.2 ) THEN DO IDIM= 1,NDIMVE II = IPOTVA(IDIM) ICODCL(IFAC,ISCA(II)) = 3 RCODCL(IFAC,ISCA(II),3) = 0.D0 ENDDO ENDIF C C --- On impose en couleur 4 une symetrie C =================================== C ELSE IF( ICOUL.EQ.4) THEN C C SYMETRIES C ITYPFB(IFAC,IPHAS) = ISYMET C C - Numero de zone (choix du numero de couleur par exemple) IZONE = ICOUL C C - Reperage de la zone a laquelle appartient la face IZFPPP(IFAC) = IZONE C C Par defaut tous les scalaires (potentiels en particulier) C recoivent une condition de flux nul. C En effet Joule, on peut souhaiter imposer une condition C d'antisymetrie sur le potentiel imaginaire selon la C configuration des electrodes : IF ( IPPMOD(IELJOU).GE. 2 ) THEN II = IPOTI ICODCL(IFAC,ISCA(II)) = 1 RCODCL(IFAC,ISCA(II),1) = 0.D0 ENDIF ENDIF C ENDDO C C---- C FORMATS C---- C C---- C FIN C---- C RETURN END c@z