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 USLABO C ***************** C ------------------------------------------------------------- & ( IDBIA0 , IDBRA0 , & NDIM , NCELET , NCEL , NFAC , NFABOR , NFML , NPRFML , & NNOD , LNDNOD , LNDFAC , LNDFBR , NCELBR , & NVAR , NSCAL , NPHAS , & NBPMAX , NVP , NVP1 , NVEP , NIVEP , & NTERSL , NVLSTA , NVISBR , & KFACE , NBPT , ISUIVI , & NIDEVE , NRDEVE , NITUSE , NRTUSE , & IFACEL , IFABOR , IFMFBR , IFMCEL , IPRFML , & IPNFAC , NODFAC , IPNFBR , NODFBR , & ITYPFB , ITRIFB , IFRLAG , ITEPA , INDEP , & IDEVEL , ITUSER , IA , & XYZCEN , SURFAC , SURFBO , CDGFAC , CDGFBO , XYZNOD , VOLUME , & SURFBN , DT , RTPA , RTP , PROPCE , PROPFA , PROPFB , & COEFA , COEFB , & ETTP , ETTPA , TEPA , PARBOR , VITPAR , VITFLU , AUXL , & RDEVEL , RTUSER , RA ) C ------------------------------------------------------------- C*********************************************************************** C FONCTION : C ---------- c@foncb CFONC CFONC SOUS-PROGRAMME DU MODULE LAGRANGIEN : CFONC ------------------------------------- CFONC CFONC SOUS-PROGRAMME UTILISATEUR (INTERVENTION NON OBLIGATOIRE) CFONC CFONC SOUS-PROGRAMME UTILISATEUR DE GESTION DU COMPORTEMENT CFONC DES PARTICULES A UNE INTERACTION PARTICULE/FRONTIERE CFONC ET D'ENREGISTREMENT DES STATISTIQUES AUX FRONTIERES. CFONC CFONC L'UTILISATEUR N'A PAS A MODIFIER LE PRESENT SOUS-PROGRAMME DANS CFONC LES CONDITIONS D'UTILISATION STANDARD. DANS LE CAS OU IL CFONC SOUHAITE AVOIR DES INTERACTIONS NON PREVU, IL DOIT INTERVENIR CFONC DANS LES RUBRIQUES NUMERO 8 et 10. CFONC CFONC Gestion de l'interaction particule/face de frontiere CFONC en fonction des informations donnees par l'utilisateur CFONC (valeur de IUSCLB par zone) dans la routine USLAG2. CFONC CFONC En fonction du nom stocke dans IUSCLB et affecte CFONC a la face de bord KFACE, on donne le type de comportement CFONC des particules. CFONC CFONC En standard, on a : CFONC CFONC IUSCLB conditions au bord pour les particules CFONC = IENTRL -> zone d'injection de particules CFONC = ISORTL -> sortie du domaine CFONC = IREBOL -> rebond des particules CFONC = IDEPO1 -> deposition definitive CFONC = IDEPO2 -> deposition definitive mais la particule reste en CFONC memoire (utile si IENSI2 = 1 uniquement) CFONC = IDEPO3 -> deposition et remise en suspension possible CFONC suivant les condition de l'ecoulement CFONC = IENCRL -> encrassement (Charbon uniquement IPHYLA = 2) CFONC CFONC CFONC De plus, si on souhaite un autre type d'interaction non prevue CFONC pour zone de faces de bord, dans USLAG2 on affecte dans CFONC IUSCLB(KZONE) un des noms suivants : CFONC CFONC JBORD1, JBORD2, JBORD3, JBORD4, JBORD5 CFONC CFONC et dans cette routine on code le comportement CFONC des particules pour cette zone frontiere. CFONC CFONC ATTENTION : en debut de routine la variable ISUIVI est CFONC initialisee a une valeur aberrante et DOIT etre modifiee CFONC avant la fin de la routine. CFONC Les vitesses de la particule et du fluide vu doivent etre CFONC modifiees en fonction de l'interaction via les tableaux CFONC VITPAR et VITFLU, elles NE DOIVENT PAS etre modifiees CFONC via ETTP et ETTPA dans ce sous-programme. CFONC CFONC Regle de modification de ISUIVI : CFONC ================================= CFONC 1) Mettre ISUIVI = 0 si la particule ne doit pas etre CFONC suivi dans le maillage apres l'interaction de sa trajectoire CFONC et de la face de bord (ex : IENTRL, ISORTL, IDEPO1, IDEPO2). CFONC CFONC 2) mettre ISUIVI = 1 pour continuer a suivre la particule CFONC dans le maillage apres son interaction (ex : IREBOL, IDEPO3). CFONC On peut tres bien avoir ISUIVI = 0 ou ISUIVI = 1 pour un type CFONC d'interaction en fonction de comportement de la particule CFONC (ex : IENCRL). CFONC CFONC REMARQUE : Lorsqu'il y a interaction, les calculs de la vitesse CFONC ^^^^^^^^ de la particule et de la vitesse du fluide vu sont CFONC forcement a l'ordre 1 (meme si on est a l'ordre 2 CFONC par ailleurs). 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 ! LNDNOD ! E ! -> ! DIM. CONNECTIVITE CELLULES->FACES ! CARGU ! LNDFAC ! E ! -> ! LONGUEUR DU TABLEAU NODFAC ! CARGU ! LNDFBR ! E ! -> ! LONGUEUR DU TABLEAU NODFBR ! 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 ! NBPMAX ! E ! -> ! NOMBRE MAX DE PARTICULIES AUTORISE ! CARGU ! NVP ! E ! -> ! NOMBRE DE VARIABLES PARTICULAIRES ! CARGU ! NVP1 ! E ! -> ! NVP SANS POSITION, VFLUIDE, VPART ! CARGU ! NVEP ! E ! -> ! NOMBRE INFO PARTICULAIRES (REELS) ! CARGU ! NIVEP ! E ! -> ! NOMBRE INFO PARTICULAIRES (ENTIERS) ! CARGU ! NTERSL ! E ! -> ! NBR TERMES SOURCES DE COUPLAGE RETOUR! CARGU ! NVLSTA ! E ! -> ! NOMBRE DE VAR STATISTIQUES LAGRANGIEN! CARGU ! NVISBR ! E ! -> ! NOMBRE DE STATISTIQUES AUX FRONTIERES! CARGU ! KFACE ! E ! -> ! NUMERO DE LA FACE D'INTERACTION ! CARGU ! NBPT ! E ! -> ! NUMERO DE LA PARTICULE TRAITEE ! CARGU ! ISUIVI ! E ! <- ! INDICATEUR DE SUIVI DE LA PARTICULE ! 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 ! CARGU ! NODFAC ! TE ! -> ! CONNECTIVITE FACES INTERNES/NOEUDS ! CARGU ! (NFAC+1) ! ! ! ! CARGU ! IPNFBR ! TE ! -> ! POSITION DU PREMIER NOEUD DE CHAQUE ! CARGU ! (LNDFBR) ! ! ! FACE DE BORD DANS NODFBR ! CARGU ! NODFBR ! TE ! -> ! CONNECTIVITE FACES DE BORD/NOEUDS ! CARGU ! (NFABOR+1) ! ! ! ! CARGU ! ITYPFB(NFABOR! TE ! -> ! TYPE DES FACES DE BORD ! CARGU ! NPHAS) ! ! ! ! CARGU ! ITRIFB(NFABOR! TE ! -> ! TAB D'INDIRECTION POUR TRI DES FACES ! CARGU ! NPHAS) ! ! ! ! CARGU ! IFRLAG ! TE ! -> ! NUMERO DE ZONE DE LA FACE DE BORD ! CARGU ! (NFABOR) ! ! ! POUR LE MODULE LAGRANGIEN ! CARGU ! ITEPA ! TE ! <- ! INFO PARTICULAIRES (ENTIERS) ! CARGU ! (NBPMAX,NIVEP! ! ! (CELLULE DE LA PARTICULE,...) ! CARGU ! INDEP ! TE ! <- ! POUR CHAQUE PARTICULE : ! CARGU ! (NBPMAX) ! ! ! NUMERO DE LA CELLULE DE DEPART ! 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 ! CARGU ! (NDIM,NNOD) ! ! ! ! CARGU ! VOLUME(NCELET! TR ! -> ! VOLUME D'UN DES NCELET ELEMENTS ! CARGU ! SURFBN(NFABOR! TR ! -> ! SURFACE DES FACES DE BORD ! CARGU ! DT(NCELET) ! TR ! -> ! PAS DE TEMPS ! CARGU ! RTP, RTPA ! TR ! -> ! VARIABLES DE CALCUL AU CENTRE DES ! CARGU ! (NCELET,*) ! ! ! CELLULES (INSTANT COURANT ET 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 ! ETTP ! TR ! <-> ! TABLEAUX DES VARIABLES LIEES ! CARGU ! (NBPMAX,NVP)! ! ! AUX PARTICULES ETAPE COURANTE ! CARGU ! ETTPA ! TR ! -> ! TABLEAUX DES VARIABLES LIEES ! CARGU ! (NBPMAX,NVP)! ! ! AUX PARTICULES ETAPE PRECEDENTE ! CARGU ! TEPA ! TR ! -> ! INFO PARTICULAIRES (REELS) ! CARGU ! (NBPMAX,NVEP)! ! ! (POIDS STATISTIQUES,...) ! CARGU ! PARBOR(NFABOR! TR ! <-> ! CUMUL DES STATISTIQUES AUX FRONTIERES! CARGU ! NVISBR) ! ! ! ! CARGU ! VITPAR ! TR ! <-> ! VITESSE PARTICULE POUR LE TRAITEMENT ! CARGU ! (NBPMAX,3) ! ! ! INTERACTIONS PARTICULES/FRONTIERES ! CARGU ! VITFLU ! TR ! <-> ! VITESSE FLUIDE VU POUR LE TRAITEMENT ! CARGU ! (NBPMAX,3) ! ! ! INTERACTIONS PARTICULES/FRONTIERES ! CARGU ! AUXL(NBPMAX,3! TR ! - ! TABLEAU DE TRAVAIL ! 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*********************************************************************** C IMPLICIT NONE C C*********************************************************************** C DONNEES EN COMMON C*********************************************************************** C INCLUDE "paramx.h" INCLUDE "cstnum.h" INCLUDE "numvar.h" INCLUDE "optcal.h" INCLUDE "entsor.h" INCLUDE "cstphy.h" INCLUDE "pointe.h" INCLUDE "period.h" INCLUDE "parall.h" INCLUDE "lagpar.h" INCLUDE "lagran.h" INCLUDE "ppppar.h" INCLUDE "ppthch.h" INCLUDE "cpincl.h" C C*********************************************************************** C C ARGUMENTS C INTEGER IDBIA0 , IDBRA0 INTEGER NDIM , NCELET , NCEL , NFAC , NFABOR INTEGER NFML , NPRFML INTEGER NNOD , LNDNOD , LNDFAC , LNDFBR , NCELBR INTEGER NVAR , NSCAL , NPHAS INTEGER NBPMAX , NVP , NVP1 , NVEP , NIVEP INTEGER NTERSL , NVLSTA , NVISBR INTEGER KFACE , NBPT , ISUIVI INTEGER NIDEVE , NRDEVE , NITUSE , NRTUSE 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 ITYPFB(NFABOR,NPHAS) , ITRIFB(NFABOR,NPHAS) INTEGER IFRLAG(NFABOR) , ITEPA(NBPMAX,NIVEP) INTEGER INDEP(NBPMAX) INTEGER IDEVEL(NIDEVE) , ITUSER(NITUSE) INTEGER 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 SURFBN(NFABOR) 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 ETTP(NBPMAX,NVP) , ETTPA(NBPMAX,NVP) DOUBLE PRECISION TEPA(NBPMAX,NVEP) DOUBLE PRECISION PARBOR(NFABOR,NVISBR) , AUXL(NBPMAX,3) DOUBLE PRECISION VITPAR(NBPMAX,3) , VITFLU(NBPMAX,3) DOUBLE PRECISION RDEVEL(NRDEVE) , RTUSER(NRTUSE) DOUBLE PRECISION RA(*) C C VARIABLES LOCALES C INTEGER IDEBIA, IDEBRA INTEGER IP , NFIN , KZONE , N1 , ICHA, IOK C DOUBLE PRECISION AA DOUBLE PRECISION XP , YP , ZP DOUBLE PRECISION XQ , YQ , ZQ DOUBLE PRECISION XK , YK , ZK DOUBLE PRECISION XPQ , YPQ , ZPQ DOUBLE PRECISION XNN , YNN , ZNN DOUBLE PRECISION VNORL(1) , ENC3 , VISCP , MASSE DOUBLE PRECISION DPINIT , DP03 , MP0 , TRAP , VNORM , ANG DOUBLE PRECISION ENERG , ENERGT , VV1 DOUBLE PRECISION UXN , VYN , WZN C C*********************************************************************** C C======================================================================= C 0. GESTION MEMOIRE C======================================================================= C IDEBIA = IDBIA0 IDEBRA = IDBRA0 C C======================================================================= C 1. Traitement en fonction du type de frontiere C======================================================================= C IOK = 0 C C--> numero de la particule traitee C IP = NBPT C C--> Zone de la face de bord a traitee C KZONE = IFRLAG(KFACE) C C-->normale normee sortante de la face KFACE C AA = 1.D0 / SURFBN(KFACE) XNN = SURFBO(1,KFACE) * AA YNN = SURFBO(2,KFACE) * AA ZNN = SURFBO(3,KFACE) * AA C C======================================================================= C 2. Recherche du point d'intersection entre la face de bord et C le rayon. Les coordonnees sont stockees dans XK YK ZK C======================================================================= C C Petit rappel de geometrie 3D : C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C 1) equation d'un plan de normal (a,b,c) : C a*x + b*y + c*z + d = 0 C 2) equation d'une droite qui passe par P et Q : C x = XP + (XQ-XP) * AA C y = YP + (YQ-YP) * AA C z = ZP + (ZQ-ZP) * AA C ou AA est un parametre qui varie dans l'ensemble des reels C C-->on determine le vecteur PQ C XP = ETTPA(IP,JXP) YP = ETTPA(IP,JYP) ZP = ETTPA(IP,JZP) C XQ = ETTP(IP,JXP) YQ = ETTP(IP,JYP) ZQ = ETTP(IP,JZP) C XPQ = XQ - XP YPQ = YQ - YP ZPQ = ZQ - ZP C C-->si la particule n'a pas bougee ( si elle est deposee C sur la face de bord), elle n'est pas traitee. C IF (XPQ.EQ.0.D0 .AND. YPQ.EQ.0.D0 .AND. ZPQ.EQ.0.D0) RETURN C C-->a partir de l'equation du plan de la face et de l'equation C parametre du rayon, on determine le parametre du point d'intersection C AA = XPQ * SURFBO(1,KFACE) & + YPQ * SURFBO(2,KFACE) & + ZPQ * SURFBO(3,KFACE) C IF ( AA.EQ.0.D0 ) THEN WRITE (NFECRA,9010) IP NBPERR = NBPERR + 1 DNBPER = DNBPER + TEPA(IP,JRPOI) ISUIVI = 0 ITEPA(IP,JISOR) = 0 RETURN ENDIF C AA = & ( SURFBO(1,KFACE) * CDGFBO(1,KFACE) & + SURFBO(2,KFACE) * CDGFBO(2,KFACE) & + SURFBO(3,KFACE) * CDGFBO(3,KFACE) & - SURFBO(1,KFACE) * XP & - SURFBO(2,KFACE) * YP & - SURFBO(3,KFACE) * ZP ) & / AA C C-->on reporte ce parametre dans l'equation de la droite du rayon pour C avoir le point d'intersection (XK YK ZK) C XK = XP + XPQ * AA YK = YP + YPQ * AA ZK = ZP + ZPQ * AA C C======================================================================= C 3. Depart de la particule du domaine de calcul, C ou deposition de la particule sur la frontiere C======================================================================= C IF ( IUSCLB(KZONE).EQ.ISORTL .OR. & IUSCLB(KZONE).EQ.IENTRL .OR. & IUSCLB(KZONE).EQ.IDEPO1 ) THEN C ISUIVI = 0 ITEPA(IP,JISOR) = 0 C C mise a jour du debit C DEBLAG(KZONE) = DEBLAG(KZONE)-TEPA(IP,JRPOI)*ETTP(IP,JMP) C C--> La particule sort mais pour la visualisation ensight, on la place C correctement au point d'intersection. C Possibilite que le numero IP soit reutilise pour une autre C particule a entrer. C ETTP(IP,JXP) = XK ETTP(IP,JYP) = YK ETTP(IP,JZP) = ZK C C======================================================================= C 4. Deposition de la particule, celle-ci reste en memoire C======================================================================= C ELSE IF (IUSCLB(KZONE).EQ.IDEPO2) THEN C C--> La particule ne sort pas, elle n'est plus traitee, C mais toujours visualisee. Le numero IP n'est pas reutilisable. C ISUIVI = 0 ITEPA(IP,JISOR) = -ITEPA(IP,JISOR) ETTP(IP,JXP) = XK ETTP(IP,JYP) = YK ETTP(IP,JZP) = ZK C DO N1 = 1,3 VITPAR(IP,N1) = 0.D0 VITFLU(IP,N1) = 0.D0 ENDDO C C======================================================================= C 5. Deposition de la particule, la remise en suspension est possible C======================================================================= C ELSE IF (IUSCLB(KZONE).EQ.IDEPO3) THEN C ISUIVI = 0 ITEPA(IP,JISOR) = IFABOR(KFACE) ETTP(IP,JXP) = XK ETTP(IP,JYP) = YK ETTP(IP,JZP) = ZK C DO N1 = 1,3 VITPAR(IP,N1) = 0.D0 VITFLU(IP,N1) = 0.D0 ENDDO DO N1 = JUP,JWF ETTPA(IP,N1) = 0.D0 ENDDO C C======================================================================= C 6. Deposition de la particule avec force d'attachement, C vitesse conservee et re-entrainement possible C======================================================================= C ELSE IF (IUSCLB(KZONE).EQ.IDEPFA) THEN C C C Calcul du critere C UXN = ETTP(IP,JUP)*XNN VYN = ETTP(IP,JVP)*YNN WZN = ETTP(IP,JWP)*ZNN C ENERG = 0.5D0*ETTP(IP,JMP)*(UXN*UXN+VYN*VYN+WZN*WZN) C ENERGT = 3.34D-12*ETTP(IP,JDP) C IF ( ENERG .GE. ENERGT )THEN C ISUIVI = 0 ITEPA(IP,JISOR) = IFABOR(KFACE) ETTP(IP,JXP) = XK ETTP(IP,JYP) = YK ETTP(IP,JZP) = ZK C C-->changement de la vitesse de la particule au point d'arrive C (comme pour une rebond elastique) C AA = ABS(( VITPAR(IP,1)*XNN & +VITPAR(IP,2)*YNN & +VITPAR(IP,3)*ZNN) )*2.D0 C VITPAR(IP,1) = VITPAR(IP,1) - AA*XNN VITPAR(IP,2) = VITPAR(IP,2) - AA*YNN VITPAR(IP,3) = VITPAR(IP,3) - AA*ZNN C C-->Annule la vitesse du fluide vu au point d'arrive C (comme pour une rebond elastique) C AA = ABS( (VITFLU(IP,1)*XNN & + VITFLU(IP,2)*YNN & + VITFLU(IP,3)*ZNN) ) * 2.D0 C VITFLU(IP,1) = 0.D0 VITFLU(IP,2) = 0.D0 VITFLU(IP,3) = 0.D0 ELSE C ISUIVI = 0 ITEPA(IP,JISOR) = INDEP(IP) ETTP(IP,JXP) = ETTPA(IP,JXP) ETTP(IP,JYP) = ETTPA(IP,JYP) ETTP(IP,JZP) = ETTPA(IP,JZP) C C On suppose que la paroi est en X,Y (conduite en Z) C VV1 = SQRT((ENERGT-ENERG)/(0.5D0*ETTP(IP,JMP))) VITPAR(IP,1) = VV1*XNN VITPAR(IP,2) = VV1*YNN VITPAR(IP,3) = ETTP(IP,JWP) ENDIF C C======================================================================= C 7. Rebond elastique de la particule sur la frontiere C======================================================================= C ELSE IF (IUSCLB(KZONE).EQ.IREBOL) THEN C ISUIVI = 1 ITEPA(IP,JISOR) = IFABOR(KFACE) C C-->changement du point de depart C ETTPA(IP,JXP) = XK ETTPA(IP,JYP) = YK ETTPA(IP,JZP) = ZK C IF (IENSI1.EQ.1) THEN NFIN = 0 CALL ENSLAG C =========== & ( IDEBIA, IDEBRA , & NBPMAX , NVP , NVP1 , NVEP , NIVEP , & NFIN , IP , & ITEPA , & ETTPA , TEPA , RA) ENDIF C C-->changement du point d'arrive C (la valeur absolue est pour eviter les p. scalaires C negatifs impossibles qui interviennent avec des erreurs d'arrondi C machines) C AA = 2.D0 * ABS( (XQ-XK)*XNN + (YQ-YK)*YNN + (ZQ-ZK)*ZNN ) C ETTP(IP,JXP) = XQ - AA*XNN ETTP(IP,JYP) = YQ - AA*YNN ETTP(IP,JZP) = ZQ - AA*ZNN C C-->changement de la vitesse de la particule au point d'arrive C C AA = ABS( (VITPAR(IP,1)*XNN & + VITPAR(IP,2)*YNN & + VITPAR(IP,3)*ZNN) ) * 2.D0 C VITPAR(IP,1) = VITPAR(IP,1) - AA*XNN VITPAR(IP,2) = VITPAR(IP,2) - AA*YNN VITPAR(IP,3) = VITPAR(IP,3) - AA*ZNN C C-->changement de la vitesse du fluide vu au point d'arrive C AA = ABS( (VITFLU(IP,1)*XNN & + VITFLU(IP,2)*YNN & + VITFLU(IP,3)*ZNN) ) * 2.D0 C VITFLU(IP,1) = VITFLU(IP,1) - AA*XNN VITFLU(IP,2) = VITFLU(IP,2) - AA*YNN VITFLU(IP,3) = VITFLU(IP,3) - AA*ZNN C C======================================================================= C 8. Encrassement des grains de charbon C======================================================================= C ELSE IF (IUSCLB(KZONE).EQ.IENCRL) THEN C C--> Encrassement de la particules si ses proprietes le permettent C et en fonction d'une probabilite C ICI c'est si Tp > TPENC C si VISCP > VISCREF C IF ( ETTP(IP,JHP).GT.TPRENC ) THEN C ENC3 = ( (1.D+7 * ENC1)/((ETTP(IP,JHP)-150.D0)**2) ) + ENC2 VISCP = EXP( LOG(10.D0)*ENC3 ) * 0.1D0 C IF ( VISCP.GT.VISREF ) THEN N1 = 1 CALL ZUFALL(N1,VNORL(1)) TRAP = 1.D0-VISREF / VISCP ENDIF C C Si VISCP <= VISREF ===> 100% de chance d'encrasser C Si VISCP > VISREF ===> probabilte TRAP = 1-VISREF/VISCP C d'encrasser C ===> On encrasse si VNORL est compris C entre TRAP et 1. C IF ( VISCP.LE.VISREF .OR. & (VISCP.GT.VISREF .AND. VNORL(1).GE.TRAP) ) THEN C C La calcul de la masse de grains de charbon encrasses est faite plus bas. C NPENCR = NPENCR + 1 ISUIVI = 0 ITEPA(IP,JISOR) = 0 ETTP(IP,JXP) = XK ETTP(IP,JYP) = YK ETTP(IP,JZP) = ZK C ENDIF ENDIF C C--> Si pas encrassement alors rebond elastique C IF ( ITEPA(IP,JISOR).NE.0 ) THEN C ISUIVI = 1 ITEPA(IP,JISOR) = IFABOR(KFACE) C C-->changement du point de depart C ETTPA(IP,JXP) = XK ETTPA(IP,JYP) = YK ETTPA(IP,JZP) = ZK C IF (IENSI1.EQ.1) THEN NFIN = 0 CALL ENSLAG C =========== & ( IDEBIA, IDEBRA , & NBPMAX , NVP , NVP1 , NVEP , NIVEP , & NFIN , IP , & ITEPA , & ETTPA , TEPA , RA) ENDIF C C-->changement du point d'arrive C AA = 2.D0 * ABS((XQ-XK)*XNN + (YQ-YK)*YNN + (ZQ-ZK)*ZNN) C ETTP(IP,JXP) = XQ - AA*XNN ETTP(IP,JYP) = YQ - AA*YNN ETTP(IP,JZP) = ZQ - AA*ZNN C ENDIF C IF (ITEPA(IP,JISOR).GT.0) THEN C C-->changement de la vitesse de la particule au point d'arrive C C AA = ABS( (VITPAR(IP,1)*XNN & + VITPAR(IP,2)*YNN & + VITPAR(IP,3)*ZNN) ) * 2.D0 C VITPAR(IP,1) = VITPAR(IP,1) - AA*XNN VITPAR(IP,2) = VITPAR(IP,2) - AA*YNN VITPAR(IP,3) = VITPAR(IP,3) - AA*ZNN C C-->changement de la vitesse du fluide vu au point d'arrive C AA = ABS( (VITFLU(IP,1)*XNN & + VITFLU(IP,2)*YNN & + VITFLU(IP,3)*ZNN) ) * 2.D0 C VITFLU(IP,1) = VITFLU(IP,1) - AA*XNN VITFLU(IP,2) = VITFLU(IP,2) - AA*YNN VITFLU(IP,3) = VITFLU(IP,3) - AA*ZNN C ENDIF C C======================================================================= C 9. Interaction utilisateur numero 1 : JBORD1 C======================================================================= C C ON PEUT FAIRE DE MEME AVEC JBORD2, JBORD3, JBORD4 et JBORD5 C On ne donnne l'exemple que pour JBORD1 C C On regarde si on est dans la zone cherchee C ELSE IF (IUSCLB(KZONE).EQ.JBORD1) THEN C C si on doit continuer a suivre la particule C ISUIVI = 0 OU 1 C C l'element de maillage concerne C ITEPA(IP,JISOR) = C C on change la localisation de depart de la particule C ETTP(IP,JXP) = C ETTP(IP,JYP) = C ETTP(IP,JZP) = C C la vitesse de la particule au point d'arrivee C VITPAR(IP,1) = C VITPAR(IP,2) = C VITPAR(IP,3) = C C la vitesse du fluide vu au point d'arrivee C VITFLU(IP,1) = C VITFLU(IP,2) = C VITFLU(IP,3) = C C ELSE IF (IUSCLB(KZONE).EQ.JBORD1 C TEST_A_ENLEVER_POUR_UTILISER_LE_SOUS_PROGRAMME_DEBUT C en standard, sans intervention de l'utilisateur, C on ne souhaite pas passer ici C mais on desire C que le test IUSCLB(KZONE).EQ.JBORD1 soit dans le us* exemple C que le source ci-dessous soit compile pour verifier C qu'il n y a pas d'erreur & .AND.(0.EQ.1) C TEST_A_ENLEVER_POUR_UTILISER_LE_SOUS_PROGRAMME_FIN & ) THEN C C ---------------------------------------------------- C EXEMPLE 1 : LA PARTICULE A UNE CHANCE SUR DEUX (TRAP) C DE SE DEPOSER DEFINITIVEMENT A LA PAROI C OU DE REBONDIR VERS L'ECOULEMENT. C ---------------------------------------------------- C C N1 = 1 CALL ZUFALL(N1,VNORL(1)) TRAP = 0.5D0 C IF (VNORL(1).GE.TRAP) THEN C ISUIVI = 0 ITEPA(IP,JISOR) = 0 ETTP(IP,JXP) = XK ETTP(IP,JYP) = YK ETTP(IP,JZP) = ZK C ELSE C ISUIVI = 1 ITEPA(IP,JISOR) = IFABOR(KFACE) C C-->changement du point de depart C ETTPA(IP,JXP) = XK ETTPA(IP,JYP) = YK ETTPA(IP,JZP) = ZK C C-->changement du point d'arrive C AA = 2.D0 * ABS((XQ-XK)*XNN + (YQ-YK)*YNN + (ZQ-ZK)*ZNN) C ETTP(IP,JXP) = XQ - AA*XNN ETTP(IP,JYP) = YQ - AA*YNN ETTP(IP,JZP) = ZQ - AA*ZNN C ENDIF C C-->Inutile de traiter les particule ITEPA(IP,JISOR)=0 car C Elle vont etre eliminees de la listes des particules. C IF (ITEPA(IP,JISOR).GT.0) THEN C C-->changement de la vitesse de la particule au point d'arrive C C AA = ABS( (VITPAR(IP,1)*XNN & + VITPAR(IP,2)*YNN & + VITPAR(IP,3)*ZNN) ) * 2.D0 C VITPAR(IP,1) = VITPAR(IP,1) - AA*XNN VITPAR(IP,2) = VITPAR(IP,2) - AA*YNN VITPAR(IP,3) = VITPAR(IP,3) - AA*ZNN C C-->changement de la vitesse du fluide vu au point d'arrive C AA = ABS( (VITFLU(IP,1)*XNN & + VITFLU(IP,2)*YNN & + VITFLU(IP,3)*ZNN) ) * 2.D0 C VITFLU(IP,1) = VITFLU(IP,1) - AA*XNN VITFLU(IP,2) = VITFLU(IP,2) - AA*YNN VITFLU(IP,3) = VITFLU(IP,3) - AA*ZNN C ENDIF C C C======================================================================= C 10. Verification et sortie si erreur C======================================================================= C ELSE WRITE (NFECRA,9020) KZONE IOK = IOK + 1 ENDIF C IF (IOK.NE.0) THEN CALL CSEXIT (1) C =========== ENDIF C C======================================================================= C 11. Enregistrement de l'interaction particule-frontiere s'il y a lieu C======================================================================= C C L'enregistrement des statistiques parietales debute des que C l'indicateur est a IENSI3 = 1. Cependant tant que le numero C de l'iteration Lagrangienne absolue est inferieur a NSTBOR, C ou que l'ecoulement est instationnaire (ISTTIO = 0), alors C le tableau PARBOR est remis a zero avant d'entrer dans ce C sous-programme. C C NPSTF : nombre d'iterations de calcul de stat aux frontieres C stationnaires C C NPSTFT : nombre d'iterations total des stats iaux frontieres C depuis le du calcul, partie instationnaire comprise C (a utiliser que pour les affichages listing) C C TSTATP : Temps physique d'enregistrement des statistiques C des interactions particules/frontiere stationnaires, C si instatinnaire alors contient DTP le dernier C pas de temps Lagrangien C C Avant impression dans le listing, ou sortie pour le C post-processing, on applique aux statistiques sur les frontieres C une moyenne en fonction des infos fournies dans USLAG1 de la C maniere suivante : C C C C CES LIGNES NE SONT DONNEES QU'A TITRE INDICATIF C C C DO IVAR = 1,NVISBR C C IF (IMOYBR(IVAR).EQ.2) THEN C C DO IFAC = 1,NFABOR C IF (PARBOR(IFAC,INBR).GT.SEUILF) THEN C PARBOR(IFAC,IVAR) = PARBOR(IFAC,IVAR) /PARBOR(IFAC,INBR) C ELSE C PARBOR(IFAC,IVAR) = 0.D0 C ENDIF C ENDDO C C ELSE IF (IMOYBR(IVAR).EQ.1) THEN C C DO IFAC = 1,NFABOR C IF (PARBOR(IFAC,INBR).GT.SEUILF) THEN C PARBOR(IFAC,IVAR) = PARBOR(IFAC,IVAR) / TSTATP C ELSE C PARBOR(IFAC,IVAR) = 0.D0 C ENDIF C ENDDO C ENDIF C ENDDO C C C C IF ( IENSI3.EQ.1 ) THEN C C--> exemple de types d'interaction sur lesquels on souhaite C enregistrer des informations C IF ( IUSCLB(KZONE).EQ.IREBOL .OR. & IUSCLB(KZONE).EQ.IDEPO1 .OR. & IUSCLB(KZONE).EQ.IDEPO2 .OR. & IUSCLB(KZONE).EQ.IDEPO3 ) THEN C IF (INBRBD.EQ.1) THEN PARBOR(KFACE,INBR) = PARBOR(KFACE,INBR) + TEPA(IP,JRPOI) ENDIF C IF (IFLMBD.EQ.1) THEN PARBOR(KFACE,IFLM) = PARBOR(KFACE,IFLM) & + ( TEPA(IP,JRPOI) * ETTP(IP,JMP) /SURFBN(KFACE) ) ENDIF C IF (IANGBD.EQ.1) THEN VNORM = ETTP(IP,JUP) * ETTP(IP,JUP) & + ETTP(IP,JVP) * ETTP(IP,JVP) & + ETTP(IP,JWP) * ETTP(IP,JWP) VNORM = SQRT( VNORM ) ANG = ETTP(IP,JUP) * SURFBO(1,KFACE) & + ETTP(IP,JVP) * SURFBO(2,KFACE) & + ETTP(IP,JWP) * SURFBO(3,KFACE) & / SURFBN(KFACE) & / VNORM ANG = ACOS(ANG) PARBOR(KFACE,IANG) = PARBOR(KFACE,IANG) + ANG*TEPA(IP,JRPOI) ENDIF C IF (IVITBD.EQ.1) THEN VNORM = ETTP(IP,JUP) * ETTP(IP,JUP) & + ETTP(IP,JVP) * ETTP(IP,JVP) & + ETTP(IP,JWP) * ETTP(IP,JWP) VNORM = SQRT( VNORM ) PARBOR(KFACE,IVIT) =PARBOR(KFACE,IVIT) +VNORM*TEPA(IP,JRPOI) ENDIF C IF (NUSBOR.GT.0) THEN DO N1 = 1,NUSBOR PARBOR(KFACE,IUSB(N1)) = 0.D0 ENDDO ENDIF C C--> Cas particulier de la masse de grains de charbon encrasses C ELSE IF ( IUSCLB(KZONE).EQ.IENCRL .AND. ISUIVI.EQ.0 ) THEN C PARBOR(KFACE,INBR) = PARBOR(KFACE,INBR) + TEPA(IP,JRPOI) C IF (IENCBD.EQ.1) THEN C ICHA = ITEPA(IP,JINCH) DPINIT = TEPA(IP,JRD0P) C DP03 = DPINIT * DPINIT * DPINIT MP0 = PI * DP03 * RHO0CH(ICHA) / 6.D0 C MASSE = ETTP(IP,JMCH) + ETTP(IP,JMCK) & + XASHCH(ICHA) * MP0 C PARBOR(KFACE,IENC) = PARBOR(KFACE,IENC) & + TEPA(IP,JRPOI)*MASSE C ENDIF C ENDIF C ENDIF C C C======================================================================= C Archives. Je laisse cette partie au cas ou... C Creation d'un repere local associe a la face frontiere C======================================================================= C C Le repere local (T1,T2,N) est construit de facon a ce que N soit C la normale normee de la face et que T1 et T2 appartiennent a la face C C-->1. Je connais les vecteurs N et PK, je definis T1 tel que C T1 = PK vectoriel N C c XPK = XK - ETTPA(IP,JXP) c YPK = YK - ETTPA(IP,JYP) c ZPK = ZK - ETTPA(IP,JZP) C c XT1 = YPK*ZNN - ZPK*YNN c YT1 = ZPK*XNN - XPK*ZNN c ZT1 = XPK*YNN - YPK*XNN C c AA = SQRT(XT1*XT1 + YT1*YT1 + ZT1*ZT1) c XT1 = XT1 / AA c YT1 = YT1 / AA c ZT1 = ZT1 / AA C C-->2. Maintenant je peux construire T2 = - T1 vectoriel N C c XT2 = YT1*ZNN - ZT1*YNN c YT2 = ZT1*XNN - XT1*ZNN c ZT2 = XT1*YNN - YT1*XNN C c AA = SQRT(XT2*XT2 + YT2*YT2 + ZT2*ZT2) c XT2 = -XT2 / AA c YT2 = -YT2 / AA c ZT2 = -ZT2 / AA C C-------- C FORMATS C-------- C 9010 FORMAT( &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/, &'@ @@ ATTENTION : ARRET A L''EXECUTION DU MODULE LAGRANGIEN ',/, &'@ ********* (USLABO) ',/, &'@ ',/, &'@ La normale d''une face de frontiere est perpendiculaire ',/, &'@ a un rayon PQ : impossible ',/, &'@ ',/, &'@ La particle ',I10,' est ELIMINEE ',/, &'@ ',/, &'@ Contacter l''equipe de developpement. ',/, &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/) C 9020 FORMAT( &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/, &'@ @@ ATTENTION : ARRET A L''EXECUTION DU MODULE LAGRANGIEN ',/, &'@ ********* (USLABO) ',/, &'@ ',/, &'@ Le type de condition aux limites IUSCLB ',/, &'@ n''est pas renseigne pour la frontiere NB = ',I10 ,/, &'@ ',/, &'@ Le calcul ne peut etre execute. ',/, &'@ ',/, &'@ Verifier USLAG2 et USLABO. ',/, &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/) C C---- C FIN C---- C RETURN END c@z