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 MTTSSC C ***************** C ------------------------------------------------------------- & ( IDBIA0 , IDBRA0 , & NDIM , NCELET , NCEL , NFAC , NFABOR , NFML , NPRFML , & NNOD , LNDFAC , LNDFBR , NCELBR , & NVAR , NSCAL , NPHAS , NCEPDP , NCKPDP , NCESMP , & NIDEVE , NRDEVE , NITUSE , NRTUSE , ISCAL , & IFACEL , IFABOR , IFMFBR , IFMCEL , IPRFML , & IPNFAC , NODFAC , IPNFBR , NODFBR , ICEPDC , ICETSM , ITYPSM , & ICONRA , & IDEVEL , ITUSER , IA , & XYZCEN , SURFAC , SURFBO , CDGFAC , CDGFBO , XYZNOD , VOLUME , & DT , RTPA , RTP , PROPCE , PROPFA , PROPFB , & COEFA , COEFB , CKUPDC , SMACEL , & CRVEXP , CRVIMP , & VISCF , VISCB , XAM , & W1 , W2 , W3 , W4 , W5 , & W6 , W7 , W8 , W9 , W10 , W11 , & RDEVEL , RTUSER , RA ) C ------------------------------------------------------------- C*********************************************************************** C FONCTION : C ---------- c@foncb CFONC CFONC CFONC CALCUL DES TERMES SOURCES POUR LE SCALAIRE ISCAL CFONC CFONC SOUS-PROGRAMME SPECIFIQUE A MATISSE (COPIE DE USTSSC) CFONC CFONC LES SEULS SCALAIRES TRAITES SONT CFONC ISCA(ITAAMT) TEMPERATURE DE L'AIR CFONC ISCA(ITPCMT) TEMPERATURE DE PEAU DES COLIS CFONC ISCA(ITPPMT) TEMPERATURE DE PEAU DES MURS CFONC CFONC CFONC ON RESOUT RHO*VOLUME*D(VAR)/DT = CRVIMP*VAR + CRVEXP CFONC CFONC ON FOURNIT ICI CRVIMP ET CRVEXP (ILS CONTIENNENT RHO*VOLUME) CFONC CRVEXP en kg variable/s : CFONC ex : pour les temperatures kg degres/s CFONC pour les enthalpies Joules/s CFONC CRVIMP en kg /s : CFONC CFONC VEILLER A UTILISER UN CRVIMP NEGATIF CFONC (ON IMPLICITERA CRVIMP CFONC IE SUR LA DIAGONALE DE LA MATRICE, LE CODE AJOUTERA : CFONC MAX(-CRVIMP,0) EN SCHEMA STANDARD EN TEMPS CFONC -CRVIMP SI LES TERMES SOURCES SONT A L'ORDRE 2 CFONC CFONC CES TABLEAUX SONT INITIALISES A ZERO AVANT APPEL A CE SOUS CFONC PROGRAMME ET AJOUTES ENSUITE AUX TABLEAUX PRIS EN COMPTE CFONC POUR LA RESOLUTION CFONC CFONC EN CAS D'ORDRE 2 DEMANDE SUR LES TERMES SOURCES, ON DOIT CFONC FOURNIR CRVEXP A L'INSTANT N (IL SERA EXTRAPOLE) ET CFONC CRVIMP A L'INSTANT N+1/2 (IL EST DANS LA MATRICE, CFONC ON LE SUPPOSE NEGATIF) CFONC c@fonce C----------------------------------------------------------------------- c@argub CARGU ARGUMENTS 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 ! NCEPDP ! E ! -> ! NOMBRE DE CELLULES AVEC PDC ! CARGU ! NCKPDP ! E ! -> ! NBR DE COEF DU TENSEUR DE PDC (3 OU 6! CARGU ! NCESMP ! E ! -> ! NOMBRE DE CELLULES A SOURCE DE MASSE ! CARGU ! NIDEVE NRDEVE! E ! -> ! LONGUEUR DE IDEVEL RDEVEL ! CARGU ! NITUSE NRTUSE! E ! -> ! LONGUEUR DE ITUSER RTUSER ! CARGU ! ISCAL ! E ! -> ! NUMERO DU SCALAIRE ! 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 ! ICEPDC(NCELET! TE ! -> ! NUMERO DES NCEPDP CELLULES AVEC PDC ! CARGU ! ICETSM(NCESMP! TE ! -> ! NUMERO DES CELLULES A SOURCE DE MASSE! CARGU ! ITYPSM ! TE ! -> ! TYPE DE SOURCE DE MASSE POUR LES ! CARGU ! (NCESMP,NVAR)! ! ! VARIABLES (cf. USTSMA) ! CARGU ! ICONRA ! TE ! <-> ! TAB DE CONNECTIVITE POUR ! CARGU ! (NCELET+1) ! ! ! LE RAYONNEMENT ET LES PANACHES ! 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 ! CKUPDC(NCEPDP! TR ! -> ! TABLEAU DE TRAVAIL POUR PDC ! CARGU ! , NCKPDP)! ! ! ! CARGU ! SMACEL ! TR ! -> ! VALEUR DES VARIABLES ASSOCIEE A LA ! CARGU ! (NCESMP,* )! ! ! SOURCE DE MASSE ! CARGU ! ! ! ! POUR IVAR=IPR, SMACEL=FLUX DE MASSE ! CARGU ! CRVEXP(NCELET! TR ! <- ! TABLEAU DE TRAVAIL POUR PART EXPLICIT! CARGU ! CRVIMP(NCELET! TR ! <- ! TABLEAU DE TRAVAIL POUR TERME INSTAT ! CARGU ! VISCF(NFAC) ! TR ! - ! TABLEAU DE TRAVAIL FACES INTERNES ! CARGU ! VISCB(NFABOR ! TR ! - ! TABLEAU DE TRAVAIL FACES DE BORD ! CARGU ! XAM(NFAC,2) ! TR ! - ! TABLEAU DE TRAVAIL FACES DE BORD ! CARGU ! W1..11(NCELET! TR ! - ! TABLEAU DE TRAVAIL CELLULES ! 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 "cstnum.h" 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 "matiss.h" C 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 NCEPDP , NCKPDP , NCESMP INTEGER NIDEVE , NRDEVE , NITUSE , NRTUSE INTEGER ISCAL 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 ICEPDC(NCEPDP) INTEGER ICETSM(NCESMP), ITYPSM(NCESMP,NVAR) INTEGER ICONRA(NCELET) 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 CKUPDC(NCEPDP,NCKPDP), SMACEL(NCESMP,NVAR) DOUBLE PRECISION CRVEXP(NCELET), CRVIMP(NCELET) DOUBLE PRECISION VISCF(NFAC), VISCB(NFABOR) DOUBLE PRECISION XAM(NFAC,2) DOUBLE PRECISION W1(NCELET), W2(NCELET), W3(NCELET) DOUBLE PRECISION W4(NCELET), W5(NCELET), W6(NCELET) DOUBLE PRECISION W7(NCELET), W8(NCELET), W9(NCELET) DOUBLE PRECISION W10(NCELET), W11(NCELET) DOUBLE PRECISION RDEVEL(NRDEVE), RTUSER(NRTUSE), RA(*) C C VARIABLES LOCALES C CHARACTER*80 CHAINE INTEGER IDEBIA, IDEBRA INTEGER IVAR , IISCVR, IPCROM, IPHAS INTEGER ICOUL , IFML , IEL , JEL INTEGER MODNTL INTEGER NZONES, IZONE DOUBLE PRECISION XNU0 , XLAMB0, PRAND0, PHI0 DOUBLE PRECISION GMAT , UNTIER DOUBLE PRECISION FFORM , EMIS DOUBLE PRECISION HPLUS , HMOIN , HMAX , DIFX , DIFY , DHMUR DOUBLE PRECISION FRAGH , XXAV , TBULK , UN0 , PTOT DOUBLE PRECISION RALEIZ, REYNOL DOUBLE PRECISION XNSRAZ, XNSREY, XNSALV DOUBLE PRECISION HRAZ , HREY , H0 , HRAY DOUBLE PRECISION HMREY , HMRAZ , HM0 , HMRAY , HNB DOUBLE PRECISION CPCONT, CPPARO DOUBLE PRECISION XLIMIN, XLIMAX, YRGMIN, YRGMAX, ZALMIN, ZALMAX DOUBLE PRECISION SQZ DOUBLE PRECISION DHE , DHS , PTCN DOUBLE PRECISION FPTSTO, PVCCSC, PCCSC C C*********************************************************************** C======================================================================= C 1. INITIALISATION C======================================================================= C C --- Gestion memoire IDEBIA = IDBIA0 IDEBRA = IDBRA0 C C C --- Numero du scalaire a traiter : ISCAL (argument) C C --- Numero de la variable associee au scalaire a traiter ISCAL IVAR = ISCA(ISCAL) C C --- Nom de la variable associee au scalaire a traiter ISCAL CHAINE = NOMVAR(IPPRTP(IVAR)) C C --- Indicateur de variance C Si ISCAVR = 0 : C le scalaire ISCAL n'est pas une variance C Si ISCAVR > 0 et ISCAVR < NSCAL + 1 : C le scalaire ISCAL est une variance associee C au scalaire ISCAVR IISCVR = ISCAVR(ISCAL) C C --- Numero de phase associee au scalaire ISCAL IPHAS = IPHSCA(ISCAL) C C --- Masse volumique IPCROM = IPPROC(IROM(IPHAS)) C C C --- Reperage du pas de temps pour les impressions listing IF(NTLIST.GT.0) THEN MODNTL = MOD(NTCABS,NTLIST) ELSEIF(NTLIST.EQ.-1.AND.NTCABS.EQ.NTMABS) THEN MODNTL = 0 ELSE MODNTL = 1 ENDIF C C --- Impression IF(IWARNI(IVAR).GE.1) THEN WRITE(NFECRA,1000) CHAINE(1:8) ENDIF IF((IRANGP.LE.0).AND.(MODNTL.EQ.0)) THEN WRITE(NFECRA,2001) CHAINE(1:8) ENDIF C C --- Numerique UNTIER = 1.D0/3.D0 C C======================================================================= C 1. DIMENSIONS GEOMETRIQUES C======================================================================= C C C======================================================================= C 2. PROPRIETES PHYSIQUES C======================================================================= C C --- Proprietes de l'air C C Nombre de Prandtl moleculaire PRAND0 = 0.7D0 C Viscosite cinematique moleculaire (nu en m2.s-1) XNU0 = XMUMAT/RO0(IPHAS) C Condutivite thermique moleculaire (lambda en W.m-1.K-1) XLAMB0 = (XNU0/PRAND0)*CP0(IPHAS)*RO0(IPHAS) C C C --- Capacite thermique fictive pour le solide C C On cherche un etat stationnaire comme limite d'un transitoire C Les equations portant sur la temperature T des solides C comportent un terme en CP dT/dt qui tend vers zero a C convergence. La valeur de CP est utilisee pour accelerer C la convergence : elle correspond a l'utilisation d'un C pas de temps adapte pour la temperature des solides (comme C on pourrait le faire avec CDTVAR) C CPCONT = 1.D0 CPPARO = 1.D0 C C C======================================================================= C 2. EROSION DES PANACHES DE CONVECTION NATURELLE C======================================================================= C C Le but est de reporter dans la zone de ciel une partie de la C puissance transmise a l'air dans la zone de stockage. C Ceci permet de simuler l'effet des panaches, le debit C enthalpique de ceux-ci etant prescrit par l'utilisateur. C C Calcul de DHS, DHE C C --- Debit enthalpique de convection naturelle calcule (initialisation) DHS = 0.D0 C C --- Debit enthalpique de convection naturelle impose pour les panaches C C Nul par defaut (pas de modelisation des panaches) DHE = 0.D0 C Calcul, si modelisation des panaches IF(IMDCNT.EQ.1) THEN DHE = DHPCNT/FRDTRA ENDIF C C======================================================================= C 2. RAYONNEMENT C======================================================================= C C --- Flux d'un conteneur (puissance ramenee a la surface en W.m-2) C C En cas de stockage en alveole, il y a un plenum inferieur libre IF(IALVEO.EQ.1) THEN PHI0 = PUICON/((HRESO-HPLEN)*DMCONT*PI) ELSE PHI0 = PUICON/(HRESO*DMCONT*PI) ENDIF C C --- Facteur de forme (alveole ou mur) / conteneur C C En alveole, le conteneur voit l'alveole entiere et c'est tout IF(IALVEO.EQ.1) THEN C FFORM = 1.D0 C C Sans alveole, le facteur de forme depend de la structure du reseau ELSE C C Pour un reseau en ligne, considerer une rangee suffit C (valable jusqu'a PTRRES/DMCONT=4) FFORM = SQRT( SQRT(DMCONT/PTRRES) + (DMCONT/PTRRES) ) C C Pour un reseau en quinconce (reseau a pas triangulaire), on C complete par la contribution d'une seconde rangee IF(ICONLG.EQ.0) THEN FFORM = MAX( (FFORM + 0.22D0), 1.D0 ) ENDIF C ENDIF C C C --- Emissivite equivalente (conteneur / mur) EMIS = 1.D0/(1.D0/EMICON + 1.D0/EMIMUR) C C --- Densite de surface d'echange (par unite de hauteur) XXAV = PI*DMCONT/(PTRRES*PLGRES) C C C --- Connectivite pour le rayonnement colis / plafond C et pour la modelisation des panaches C Une connectivite est necessaire pour relier la couche C superieure de la zone occupee par des colis et les C mailles situees directement sous le plafond C C Cote du haut des mailles du haut des conteneurs HPLUS = HRESO C Cote du bas des mailles du haut des conteneurs HMOIN = HRESO - EPCHEL C Cote du bas des mailles touchant le plafond HMAX = EPCHEL*NCHEST - EPCHEL C C Calcul de la connectivite C Attention, c'est quadratique ... C mais on ne le fait qu'au premier passage C C Calcul uniquement si pas deja fait et si utile C (utile si les colis ne touchent pas le plafond) IF( (ICNROK.EQ.0).AND.(HRESO.LT.EPCHEL*NCHEST) ) THEN C On repere les mailles de plafond dans la zone de stockage DO IEL = 1, NCEL IFML = IFMCEL(IEL ) ICOUL = IPRFML(IFML,1) IF(ICOUL.EQ.ICMTST.AND.XYZCEN(3,IEL).GT.HMAX)THEN C On repere la maille du haut des colis situee dessous DO JEL = 1, NCEL DIFX = ABS(XYZCEN(1,JEL)-XYZCEN(1,IEL)) DIFY = ABS(XYZCEN(2,JEL)-XYZCEN(2,IEL)) IF(DIFX.LT.PTRRES*1.D-2.AND.DIFY.LT.PLGRES*1.D-2) THEN IF ( XYZCEN(3,JEL).LT.HPLUS.AND. & XYZCEN(3,JEL).GT.HMOIN ) THEN ICONRA(IEL) = JEL ENDIF ENDIF ENDDO ENDIF ENDDO C On indique que la connectivite a ete calculee ICNROK = 1 ENDIF C C======================================================================= C C ON VA COMPLETER CI-DESSOUS LES TERMES SOURCES C - PUREMENT EXPLICITES (PUISSANCE) C - PARTIELLEMENT EXPLICITES ET IMPLICITES (ECHANGES) C C C ON UTILISE LE MODELE DE TERME SOURCE SUIVANT, POUR UN SCALAIRE F : C C S = A * F + B C C APPARAISSANT DANS LES EQUATIONS DE TRANSPORT DE F SOUS LA FORME : C C RHO VOLUME D(F)/Dt = VOLUME*S C C C CE TERME A UNE PARTIE QU'ON VEUT IMPLICITER : A C ET UNE PARTIE QU'ON VA TRAITER EN EXPLICITE : B C C C PAR EXEMPLE SI ON A : C A = - RHO / TAUF C B = RHO * PRODF C AVEC C TAUF = 10.D0 [secondes ] (TEMPS DE DISSIPATION DE F) C PRODF = 100.D0 [variable/s] (PRODUCTION DE F PAR UNITE DE TEMPS) C C ON A ALORS C CRVIMP(IEL) = VOLUME(IEL)* A = - VOLUME(IEL) (RHO / TAUF ) C CRVEXP(IEL) = VOLUME(IEL)* B = VOLUME(IEL) (RHO * PRODF) C C======================================================================= C C === Initialisation C =========================================================== C C --- Compteurs C pour moyenne et affichage, mais pas seulement C C Coefficient d'echange moyen de convection forcee HMREY = 0.D0 C Coefficient d'echange moyen de convection naturelle HMRAZ = 0.D0 C Coefficient d'echange convectif moyen efficace HM0 = 0.D0 C Coefficient d'echange radiatif moyen pour les parois HMRAY = 0.D0 C Compteur pour moyenne HNB = 0.D0 C Puissance totale sans panache PTOT = 0.D0 C Puissance totale avec panache PTCN = 0.D0 C C Gravite GMAT = SQRT(GX**2+GY**2+GZ**2) C C C === Boucle sur les cellules et selection selon couleur C =========================================================== C DO IEL = 1, NCEL C C Vitesse horizontale UN0 = SQRT(RTP(IEL,IU(IPHAS))**2+RTP(IEL,IV(IPHAS))**2) C C Temperature de reference Kelvin TBULK = (TINIT+RTP(IEL,ISCA(ITAAMT)))*0.5D0 + TKELVI C C Raleigh/DeltaT = FRAGH = g*H**3*beta/(nu a) C . g = GMAT C . H = XYZCEN(3,IEL) : on suppose que le sol est a z=0 C . beta = 1/TBULK C . nu = XNU0 C . a = nu/Pr = XNU0/PRAND0 C FRAGH = GMAT*XYZCEN(3,IEL)**3*PRAND0/(TBULK*XNU0**2) C C Couleur pour selection de la zone de stockage IFML = IFMCEL(IEL ) ICOUL = IPRFML(IFML,1) C IF(ICOUL.EQ.ICMTST) THEN C C C === Calcul du coefficient d'echange par rayonnement C =========================================================== C C Coeff d'echange = flux / (S (Tcolis-Tparoi)) C = Sigma Epsilon FactForme (Tcolis+Tparoi)(Tcolis**2+Tparoi**2) C C La paroi consideree est le plafond ou les murs ou les alveoles C C --- Zone de ciel C On utilise la connectivite calculee precedemment IF(ICNROK.EQ.1.AND.XYZCEN(3,IEL).GT.HMAX) THEN C C En presence d'alveoles, C le plafond echange avec les alveoles face a face C (seule l'alveole directement en face est vue, C l'alveole voit un disque de plafond) C La temperature du plafond et des alveoles est representee C par le meme scalaire C Ce coefficient est utilise pour le calcul des termes sources C permettant de determiner la temperature de peau des parois C (plafond, mur, alveoles) C On le verra aussi intervenir pour un calcul de temperature C de colis, mais dans la zone de ciel ou la temperature de colis C n'a pas de signification (et comme il n'y a pas de diffusion C en alveoles, c'est ok) IF(IALVEO.EQ.1) THEN HRAY = STEPHN*EMIMUR & * (0.25D0*PI*DMCONT**2)/(PTRRES*PLGRES) & * ( RTP(ICONRA(IEL),ISCA(ITPPMT))+TKELVI & +RTP( IEL ,ISCA(ITPPMT))+TKELVI) & * ((RTP(ICONRA(IEL),ISCA(ITPPMT))+TKELVI)**2 & +(RTP( IEL ,ISCA(ITPPMT))+TKELVI)**2) C C Sans alveoles, C le plafond echange avec les colis face a face C Ce coefficient sera utilise pour le calcul des termes sources C permettant de determiner la temperature de peau du plafond ELSE HRAY = STEPHN*EMIS & * (0.25D0*PI*DMCONT**2)/(PTRRES*PLGRES) & * ( RTP(ICONRA(IEL),ISCA(ITPCMT))+TKELVI & +RTP( IEL ,ISCA(ITPPMT))+TKELVI) & * ((RTP(ICONRA(IEL),ISCA(ITPCMT))+TKELVI)**2 & +(RTP( IEL ,ISCA(ITPPMT))+TKELVI)**2) ENDIF C C --- Zones de murs ou d'alveoles C les murs et les alveoles echangent avec les colis C le facteur de forme calcule precedemment s'applique C il prend en compte le cas avec ou sans alveoles C Ce coefficient sera utilise pour le calcul des termes sources C permettant de determiner C - la temperature de peau des parois C - la temperature de peau des colis, en alveole (FFORM=1) C ELSE HRAY = STEPHN*EMIS & * FFORM & * ( RTP(IEL,ISCA(ITPCMT))+TKELVI & +RTP(IEL,ISCA(ITPPMT))+TKELVI) & * ((RTP(IEL,ISCA(ITPCMT))+TKELVI)**2 & +(RTP(IEL,ISCA(ITPPMT))+TKELVI)**2) ENDIF C C C === Puissance transmise a l'air et a la peau des colis par conduction C =========================================================== C C Pour l'air : C en stationnaire, toute la puissance des colis se retrouve tot ou C tard dans l'air (elle peut transiter par la peau des colis, C par les alveoles, par les murs, mais elle est extraite par C l'air). C C Pour la peau des colis : C toute la puissance des colis est transmise par conduction a C la peau des colis. C C IF (ISCAL.EQ.ITAAMT.OR.ISCAL.EQ.ITPCMT) THEN C C --- Puissance des colis par unite de surface au sol C On prend en compte la carte de remplissage ensuite C C On doit fournit un CRVEXP en Watt/CP, soit de la dimension de C PUICON/CP. C C Comme la puissance integree sur la hauteur du colis doit etre C la puissance totale du colis, PUICON, on doit imposer dans C chaque maille une fraction de PUICON. C Pour une repartition uniforme en hauteur, cette fraction est C la hauteur de la maille VOLUME/(PTRRES*PLGRES) divisee par C la hauteur totale du colis. C Les repartitions non uniformes en hauteur sont traitees par C les cartes de repartition en altitude, plus bas. C On impose donc CRVEXP = PUICON/CP * VOLUME/(PTRRES*PLGRES) C puis on divise par la hauteur de colis (plus bas, lors du C traitement des cartes de repartition en altitude). C C Noter que comme on multiplie par VOLUME/(PTRRES*PLGRES) ici C et qu'on divise plus bas par EPCHEL, on pourrait simplifier, C puisque, a priori, EPCHEL = VOLUME/(PTRRES*PLGRES) C C C Air ambiant IF (ISCAL.EQ.ITAAMT) THEN CRVEXP(IEL) = PUICON/CP0(IPHAS) & * VOLUME(IEL)/(PTRRES*PLGRES) C Peau des colis ELSEIF(ISCAL.EQ.ITPCMT) THEN CRVEXP(IEL) = PUICON/CPCONT & * VOLUME(IEL)/(PTRRES*PLGRES) ENDIF C C --- Prise en compte des cartes de remplissage (adimensionnelles) C Les entrepots ne sont pas uniformement remplis. C C Pour représenter la carte d'encombrement des lignes de colis, on C repere les lignes occupees (coordonnee X, indicateur ILIGNE) C par leurs bornes XLIMIN et XLIMAX qui sont stockees dans VIZCAR C (avec l'indicateur ICPUIS faisant reference a la puissance) NZONES = NZOCAR(ILIGNE,ICPUIS) DO IZONE = 1, NZONES XLIMIN = PTRRES*VIZCAR(1,IZONE,ILIGNE,ICPUIS) XLIMAX = PTRRES*VIZCAR(2,IZONE,ILIGNE,ICPUIS) IF ((XYZCEN(1,IEL).GE.XLIMIN).AND. & (XYZCEN(1,IEL).LT.XLIMAX)) THEN CRVEXP(IEL) = CRVEXP(IEL) * VCARTH(IZONE,ILIGNE) ENDIF ENDDO C C Pour représenter la carte d'encombrement des rangees de colis, on C repere les rangees occupees (coordonnee Y, indicateur IRANGE) C par leurs bornes YRGMIN et YRGMAX qui sont stockees dans VIZCAR C (avec l'indicateur ICPUIS faisant reference a la puissance) NZONES = NZOCAR(IRANGE,ICPUIS) DO IZONE = 1, NZONES YRGMIN = PLGRES*VIZCAR(1,IZONE,IRANGE,ICPUIS) YRGMAX = PLGRES*VIZCAR(2,IZONE,IRANGE,ICPUIS) IF ((XYZCEN(2,IEL).GE.YRGMIN).AND. & (XYZCEN(2,IEL).LT.YRGMAX)) THEN CRVEXP(IEL) = CRVEXP(IEL) * VCARTH(IZONE,IRANGE) ENDIF ENDDO C C Pour représenter la carte d'encombrement en hauteur de colis, on C repere les couches occupees (coordonnee Z, indicateur IALTIT) C par leurs bornes ZALMIN et ZALMAX qui sont stockees dans VIZCAR C (avec l'indicateur ICPUIS faisant reference a la puissance) C SQZ permet d'integrer VCARTH sur les zones afin que C de normer VCARTH de telle sorte que la multiplication par C VCARTH/SQZ permette de retrouver, en integrant sur la hauteur, C toute la puissance du colis. SQZ = 0.D0 NZONES = NZOCAR(IALTIT,ICPUIS) DO IZONE = 1, NZONES ZALMIN = EPCHEL*VIZCAR(1,IZONE,IALTIT,ICPUIS) ZALMAX = EPCHEL*VIZCAR(2,IZONE,IALTIT,ICPUIS) SQZ = SQZ + VCARTH(IZONE,IALTIT) * & ( VIZCAR(2,IZONE,IALTIT,ICPUIS) - & VIZCAR(1,IZONE,IALTIT,ICPUIS) ) ENDDO DO IZONE = 1, NZONES ZALMIN = EPCHEL*VIZCAR(1,IZONE,IALTIT,ICPUIS) ZALMAX = EPCHEL*VIZCAR(2,IZONE,IALTIT,ICPUIS) IF ((XYZCEN(3,IEL).GE.ZALMIN).AND. & (XYZCEN(3,IEL).LT.ZALMAX)) THEN CRVEXP(IEL) = CRVEXP(IEL)* & VCARTH(IZONE,IALTIT)/(SQZ*EPCHEL) ENDIF ENDDO C C On calcule la puissance totale injectee (Watt) C C Air ambiant IF (ISCAL.EQ.ITAAMT) THEN PTOT = PTOT + CRVEXP(IEL)*CP0(IPHAS) ENDIF C ENDIF C C C === Echange convectif conteneur / air et radiatif conteneur / mur C =========================================================== C C Pour la determination de la temperature de peau des colis C C On calcule le coefficient d'echange convectif efficace selon C la configuration (convection naturelle, forcee, mixte). C On ajoute le terme de rayonnement a partir du coefficient C equivalent HRAY calcule precedemment. C IF (ISCAL.EQ.ITPCMT) THEN C C --- Convection naturelle verticale le long des colis C Van Vliet&Ross et Vliet&Liu C C Rayleigh a partir du flux d'un conteneur RALEIZ = FRAGH*PHI0*XYZCEN(3,IEL)/XLAMB0 C Nusselt associe XNSRAZ = 0.17D0 * RALEIZ**0.25D0 XNSRAZ = MAX(XNSRAZ,(0.6D0 * RALEIZ**0.2D0)) C Coefficient d'echange de convection naturelle conteneur / air C Hypothese que le zero est a la base des colis HRAZ = XNSRAZ*XLAMB0/XYZCEN(3,IEL) C Coefficient d'echange moyen de convection naturelle HMRAZ = HMRAZ + HRAZ C C C --- Convection forcee transversale dans le reseau de colis C Zukauskas C C Reynolds (vitesse horizontale gap, diametre REYNOL = (UN0*PTRRES/(PTRRES-DMCONT)) * DMCONT / XNU0 C Nusselt associe XNSREY = 0.35D0*(PTRRES/PLGRES)**0.2D0 & * 0.89D0*REYNOL**0.6D0 C Coefficient d'echange de convection forcee conteneur / air HREY = XNSREY*XLAMB0/DMCONT C Coefficient d'echange moyen de convection forcee HMREY = HMREY + HREY C C --- Compteur pour les moyennes C HNB = HNB + 1.D0 C C C --- Prise en compte du type de reseau (en ligne ou en quinconce) C et de la presence d'alveoles C C - Alveoles C convection naturelle uniquement C IF(IALVEO.EQ.1) THEN C C Nusselt associe XNSALV = 0.17D0 * RALEIZ**0.25D0 C Coefficient d'echange de convection naturelle conteneur / air C Hypothese que le zero est a la base des colis H0 = XNSALV*XLAMB0/XYZCEN(3,IEL) C C C - Reseau en ligne (pas carre) sans alveoles C 2/3 conv.force (face+cotes) + 1/3 conv.nat (sillage arriere) C ELSEIF(ICONLG.EQ.1)THEN H0 = (HRAZ+2.D0*HREY)/3.D0 C C - Reseau en quinconce (pas triangulaire) sans alveoles C conv.force (pas de protection du colis amont) C ELSE H0 = HREY ENDIF C C - Sans alveoles C le coefficient d'echange est toujours au moins egal a celui C de la convection naturelle IF(IALVEO.NE.1) THEN H0 = MAX(HRAZ,H0) ENDIF C C - Coefficient d'echange convectif moyen efficace HM0 = HM0 + H0 C C C --- Ajout du terme conductif et du terme radiatif C C CRVEXP contient deja la puissance des colis transmise C par conduction C On veut ajouter dans CRVEXP des Watt/CP soit H0*Surface*T/CP. C On a ici XXAV*VOLUME = PI*DMCONT*EPCHEL, qui est la surface C d'échange dans la cellule. C C C - Alveoles C echange par convection naturelle et C les colis rayonnent avec le plafond et les alveoles C IF(IALVEO.EQ.1) THEN C CRVEXP(IEL) = CRVEXP(IEL) & +H0 *(XXAV*VOLUME(IEL)) & *RTP(IEL,ISCA(ITAAMT))/CPCONT & +HRAY*(XXAV*VOLUME(IEL)) & *RTP(IEL,ISCA(ITPPMT))/CPCONT CRVIMP(IEL) = -(H0+HRAY)*(XXAV*VOLUME(IEL))/CPCONT C C C - Sans alveole C echange par convection mixte seul C les colis ne rayonnent qu'entre eux, et ceci est pris en compte C par diffusion C ELSE C CRVEXP(IEL) = CRVEXP(IEL) & +H0 *(XXAV*VOLUME(IEL)) & *RTP(IEL,ISCA(ITAAMT))/CPCONT CRVIMP(IEL) = - H0 *(XXAV*VOLUME(IEL))/CPCONT C ENDIF C ENDIF C C C C === Echange convectif mur / air et radiatif mur / conteneur C =========================================================== C C Pour la determination de la temperature de peau des parois C (murs et alveoles) C IF (ISCAL.EQ.ITPPMT) THEN C C --- Echange par convection paroi / air C C - Alveoles C convection naturelle uniquement C IF(IALVEO.EQ.1) THEN C C Rayleigh a partir du flux d'un conteneur RALEIZ = FRAGH*PHI0*XYZCEN(3,IEL)/XLAMB0 C Nusselt associe XNSALV = 0.17D0 * RALEIZ**0.25D0 C Coefficient d'echange de convection naturelle alveole / air C Hypothese que le zero est a la base des colis HRAZ = XNSALV*XLAMB0/XYZCEN(3,IEL) C Coefficient d'echange efficace alveole / air C Convection naturelle seule H0 = HRAZ C C C - Sans alveole C convection naturelle (Mac Adams) ou convection forcee C ELSE C C .Convection naturelle C C Rayleigh a partir de l'ecart de temperature paroi - air RALEIZ = & FRAGH*(RTP(IEL,ISCA(ITPPMT))-RTP(IEL,ISCA(ITAAMT))) C Nusselt associe (Mac Adams) XNSRAZ = 0.12D0*(ABS(RALEIZ))**UNTIER C Coefficient d'echange de convection naturelle paroi / air HRAZ = XNSRAZ*XLAMB0/XYZCEN(3,IEL) C Coefficient d'echange moyen de convection naturelle HMRAZ = HMRAZ + HRAZ C C C .Convection forcee C C Espace mur reseau DHMUR = PTRRES-DMCONT C Reynolds associé (vitesse horizontale hors reseau, pres des murs) REYNOL = DHMUR*UN0/XNU0 C Nusselt associe (Colburn) XNSREY = 0.023D0*REYNOL**0.8D0*PRAND0**UNTIER C Coefficient d'echange de convection forcee paroi / air HREY = XNSREY*XLAMB0/DHMUR C Coeficient d'echange moyen de convection forcee HMREY = HMREY + HREY C C C .Coefficient d'echange efficace (le plus efficace est dominant) H0 = MAX(HRAZ,HREY) C ENDIF C C C - Coefficient d'echange convectif moyen efficace HM0 = HM0 + H0 C - Coefficient d'echange moyen de convection naturelle HMRAZ = HMRAZ + HRAZ C - Compteur pour les moyennes HNB = HNB + 1.D0 C C C --- Ajout du terme conductif et du terme radiatif C C On veut ajouter dans CRVEXP des Watt/CP soit H0*Surface*T/CP. C On calcule H0*VOLUME*T/CP et on divise par une distance C plus bas. On pourrait faire les choses plus directement. C C Les echanges radiatifs sont les echanges avec les colis. C C C - Coefficient d'echange radiatif moyen pour les parois HMRAY = HMRAY + HRAY C C C - Echanges convectifs CRVEXP(IEL) = H0*VOLUME(IEL)*RTP(IEL,ISCA(ITAAMT))/CPPARO CRVIMP(IEL) =-H0*VOLUME(IEL)/CPPARO C C C - Echange radiatif dans le ciel IF(ICNROK.EQ.1.AND.XYZCEN(3,IEL).GT.HMAX) THEN C C .Alveoles : le plafond echange avec les alveoles IF(IALVEO.EQ.1) THEN CRVEXP(IEL) = CRVEXP(IEL) & +HRAY*VOLUME(IEL) & *RTP(ICONRA(IEL),ISCA(ITPPMT))/CPPARO C C .Sans alveole : le plafond echange avec le dessus des colis ELSE CRVEXP(IEL) = CRVEXP(IEL) & +HRAY*VOLUME(IEL) & *RTP(ICONRA(IEL),ISCA(ITPCMT))/CPPARO ENDIF C C C - Echange radiatif avec les colis pour les autres parois C (murs et alveoles) ELSE CRVEXP(IEL) = CRVEXP(IEL) & +HRAY*VOLUME(IEL)*RTP(IEL,ISCA(ITPCMT))/CPPARO C ENDIF C C C - La partie implicite est la meme pour tous CRVIMP(IEL) = CRVIMP(IEL) - HRAY*VOLUME(IEL)/CPPARO C C C - On divise par une distance C Plus exactement, on multiplie par la surface d'echange et C on divise par le volume (on aurait donc pu eviter de multiplier C par le volume). C C C .Au dessus ou au dessous des colis, le plafond ou le sol echangent C sur une surface de PTRRES*PLGRES C dans des cellules de volume EPCHEL*PTRRES*PLGRES. C Le rapport Surface/Volume est donc 1/EPCHEL C IF(XYZCEN(3,IEL).LT.EPCHEL.OR.XYZCEN(3,IEL).GT.HMAX) THEN CRVEXP(IEL) = CRVEXP(IEL)/EPCHEL CRVIMP(IEL) = CRVIMP(IEL)/EPCHEL C C C .En stockage en alveoles, les alveoles echangent C sur une surface de PI*DHALVE*EPCHEL C dans des cellules de volume EPCHEL*PTRRES*PLGRES. C Le rapport Surface/Volume est donc PI*DHALVE/(PTRRES*PLGRES) C A priori DHALVE est positif, mais on laisse le ABS (pour le C cas, ou il serait defini comme la difference entre le C diametre des colis et des alveoles) C ELSEIF(IALVEO.EQ.1) THEN CRVEXP(IEL) = & CRVEXP(IEL)*PI*ABS(DHALVE)/(PTRRES*PLGRES) CRVIMP(IEL) = & CRVIMP(IEL)*PI*ABS(DHALVE)/(PTRRES*PLGRES) C C C .Sans alveoles, les murs echangent C sur une surface de PLGRES*EPCHEL C dans des cellules de volume EPCHEL*PTRRES*PLGRES. C Le rapport Surface/Volume est donc 1/PTRRES C La ou il n'y a pas de mur, la temperature de mur ne signifie C rien, mais on peut quand meme l'y calculer (il n'y a pas de C phenomene de diffusion, tout est local) C ELSE CRVEXP(IEL) = CRVEXP(IEL)/PTRRES CRVIMP(IEL) = CRVIMP(IEL)/PTRRES C ENDIF C ENDIF C C C Fin du test sur la zone de stockage ENDIF ENDDO C C C === Pour l'air, sans modelisation des panaches de convection naturelle C =========================================================== C IF(ISCAL.EQ.ITAAMT) THEN C IF(IMDCNT.EQ.0)THEN C C --- Calcul du debit enthalpique de convection naturelle C sortant en haut de la zone des colis : rho*W*S*Cp*DeltaT C (il aurait ete plus correct d'utiliser FLUMAS) C DO IEL = 1, NCEL C IFML = IFMCEL(IEL ) ICOUL = IPRFML(IFML,1) IF(ICOUL.EQ.ICMTST) THEN C IF(ICNROK.EQ.1.AND.XYZCEN(3,IEL).GT.HMAX) THEN JEL = ICONRA(IEL) DHS = DHS + & PROPCE(JEL,IPCROM)*RTP(JEL,IW(IPHAS)) & *(VOLUME(JEL)/EPCHEL) & *CP0(IPHAS)*(RTP(JEL,ISCA(ITAAMT))-TINIT) ENDIF C ENDIF ENDDO C C C --- Calcul de la puissance en kW, C avec correction si representation 3D en 2D C PUITOT = PTOT*FRDTRA*1.D-3 C C C --- Calcul du debit enthalpique de convection naturelle en kW, C avec correction si representation 3D en 2D C DEBCON = DHS*FRDTRA*1.D-3 C C --- Impression de PUITOT et DEBCON IF (IRANGP.LE.0.AND.MODNTL.EQ.0) THEN WRITE(NFECRA,2002) PUITOT, DEBCON ENDIF C C C === Pour l'air, modelisation des panaches de convection naturelle C =========================================================== C ELSEIF(IMDCNT.EQ.1) THEN C C C --- Le debit enthalpique des panaches ne peut pas etre plus grand C que la puissance totale disponible ; sinon, c'est qu'il C a ete mal predit ou que la puissance a ete mal imposee : C les cartes ou la puissance d'un conteneur peuvent etre en C cause C On fait ce test ici, dans la mesure ou PTOT est calcule juste C au dessus. C IF(DHPCNT.GE.PTOT*FRDTRA)THEN WRITE(NFECRA,9001) IMDCNT, DHPCNT, PTOT*FRDTRA, & DHPCNT, PTOT*FRDTRA CALL CSEXIT (1) C =========== ENDIF C C --- Modification du terme source pour la modelisation des panaches C C On redistribue la puissance transmise a l'air uniquement : C une fraction de la puissance totale est retiree de la zone C de stockage et injectee directement dans le ciel au dessus C des conteneurs. La repartition se fait a priori de maniere C homogene (le panache de tous les conteneurs presents recoit C la meme fraction de la puissance totale), puis on applique C la carte de repartition thermique XY. C C - Calcul des grandeurs independantes de la maille C C Fraction de la puissance totale a conserver C dans la zone de stockage FPTSTO = 1.D0-DHE/PTOT C C Puissance Volumique a ajouter dans la zone de Ciel au dessus C d'un Conteneur de puissance PUICON, divise par CP PVCCSC = ( (DHE/PTOT)*PUICON / (HERCNT*PTRRES*PLGRES) ) & / CP0(IPHAS) C DO IEL = 1, NCEL C IFML = IFMCEL(IEL ) ICOUL = IPRFML(IFML,1) C IF(ICOUL.EQ.ICMTST) THEN C C - Reduction dans la zone de stockage C C La puissance est reduite pour tous les conteneurs de maniere C homogene, en fonction du rapport entre le debit enthalpique C des panaches impose par l'utilisateur et la puissance totale C presente dans l'entrepot. Les cartes de repartition thermique C ayant deja ete appliquees a CRVEXP, c'est bien le rapport de C reduction (et non pas la valeur absolue de la reduction) qui C est homogene. C CRVEXP(IEL) = CRVEXP(IEL)*FPTSTO C C IF ( XYZCEN(3,IEL).LT.(HERCNT+HRESO).AND. & XYZCEN(3,IEL).GT.HRESO) THEN C C - Report dans la zone des panaches (1/2) C C Calcul de la puissance a ajouter dans la zone de ciel, au dessus C d'un conteneur de puissance PUICON (deja divise par CP) C PCCSC = PVCCSC * VOLUME(IEL) C C Les entrepots ne sont pas uniformement remplis. C Application de la carte thermique horizontale XY pour C selectionner les positions auxquelles un conteneur (ou C une fraction de conteneur est presente, puisque l'on C peut avoir des demi-conteneurs, par exemple sur un plan C de symetrie) C C Pour représenter la carte d'encombrement des lignes de colis, on C repere les lignes occupees (coordonnee X, indicateur ILIGNE) C par leurs bornes XLIMIN et XLIMAX qui sont stockees dans VIZCAR C (avec l'indicateur ICPUIS faisant reference a la puissance) NZONES = NZOCAR(ILIGNE,ICPUIS) DO IZONE = 1, NZONES XLIMIN = PTRRES*VIZCAR(1,IZONE,ILIGNE,ICPUIS) XLIMAX = PTRRES*VIZCAR(2,IZONE,ILIGNE,ICPUIS) IF ((XYZCEN(1,IEL).GE.XLIMIN).AND. & (XYZCEN(1,IEL).LT.XLIMAX)) THEN PCCSC = PCCSC * & VCARTH(IZONE,ILIGNE) ENDIF ENDDO C C Pour représenter la carte d'encombrement des rangees de colis, on C repere les rangees occupees (coordonnee Y, indicateur IRANGE) C par leurs bornes YRGMIN et YRGMAX qui sont stockees dans VIZCAR C (avec l'indicateur ICPUIS faisant reference a la puissance) NZONES = NZOCAR(IRANGE,ICPUIS) DO IZONE = 1, NZONES YRGMIN = PLGRES*VIZCAR(1,IZONE,IRANGE,ICPUIS) YRGMAX = PLGRES*VIZCAR(2,IZONE,IRANGE,ICPUIS) IF ((XYZCEN(2,IEL).GE.YRGMIN).AND. & (XYZCEN(2,IEL).LT.YRGMAX)) THEN PCCSC = PCCSC * & VCARTH(IZONE,IRANGE) ENDIF ENDDO C C Ajout de la puissance des panaches a la puissance deja existante C CRVEXP(IEL) = CRVEXP(IEL) + PCCSC C C Fin si zone de panache ENDIF C C Fin si zone de stockage ENDIF C C C --- Calcul de la puissance totale incluant la modelisation C des panaches de convection natuelle C (normalement la puissance totale n'a pas ete modifiee) C PTCN = PTCN + CRVEXP(IEL)*CP0(IPHAS) C C Fin boucle NCEL ENDDO C C --- Calcul et impression de la puissance totale incluant la modelisation C des panaches de convection naturelle C en kW, avec correction si representation 3D en 2D C PUITOT = PTCN*FRDTRA*1.D-3 IF ((IRANGP.LE.0).AND.(MODNTL.EQ.0)) THEN WRITE(NFECRA,2003) PUITOT ENDIF C C Fin si on modelise les panaches ENDIF C Fin si scalaire = air ambiant ENDIF C C C C === Impressions pour la temperature de peau des colis et des parois C =========================================================== C C --- Valeur fictive : pas de convection forcee en alveole C (valeur bizarre, mais uniquement pour affichage) IF (IALVEO.EQ.1) THEN HMREY = -1.00001D10 ENDIF C C --- Temperature de peau des colis C IF(ISCAL.EQ.ITPCMT) THEN C C Si rien a moyenner, on prend une valeur negative C (mais les numerateurs seront nuls, normalement) IF (HNB.LT.EPZERO) THEN HNB = -EPZERO ENDIF C C Stockage en common pour impression ulterieure CFECCA = HM0/HNB C C Impressions des coefficients d'echange moyens IF ((IRANGP.LE.0).AND.(MODNTL.EQ.0)) THEN WRITE(NFECRA,2004) HMREY/HNB, HMRAZ/HNB, HM0/HNB ENDIF C ENDIF C C --- Temperature de peau des parois C IF(ISCAL.EQ.ITPPMT) THEN C C Si rien a moyenner, on prend une valeur negative C (mais les numerateurs seront nuls, normalement) IF (HNB.LT.EPZERO) THEN HNB = -EPZERO ENDIF C C Stockage en common pour impression ulterieure CFECMA = HM0/HNB C C Impressions des coefficients d'echange moyens IF ((IRANGP.LE.0).AND.(MODNTL.EQ.0)) THEN WRITE(NFECRA,2005) & HMREY/HNB, HMRAZ/HNB, HMRAY/HNB, HM0/HNB ENDIF C ENDIF C C C-------- C FORMATS C-------- C 1000 FORMAT(' TERMES SOURCES MATISSE POUR LA VARIABLE ',A8,/) C 2001 FORMAT (/,3X,'** INFORMATIONS SUR MATISSE, VARIABLE ',A8,/, & 3X,' -------------------------------------------') 2002 FORMAT ( &'mati --------------------------------------------------------',/, &'mati Puissance (sans modelisation panaches) ',/, &'mati -------------------------------------- ',/, &'mati Puissance totale ', E12.5 ,' kW',/, &'mati Debit enthalpique de conv. naturelle ', E12.5 ,' kW',/, &'mati --------------------------------------------------------',/) 2003 FORMAT ( &'mati --------------------------------------------------------',/, &'mati Puissance (avec modelisation panaches) ',/, &'mati -------------------------------------- ',/, &'mati Puissance totale avec erosion panaches', E12.5 ,' kW',/, &'mati --------------------------------------------------------',/) 2004 FORMAT ( &'mati --------------------------------------------------------',/, &'mati Coeff. d echange moyens (calcul T° de peau des colis) ',/, &'mati ----------------------------------------------------- ',/, &'mati Convection forcee ', E12.5 ,' W/m2/°C',/, &'mati Convection naturelle ', E12.5 ,' W/m2/°C',/, &'mati Efficace ', E12.5 ,' W/m2/°C',/, &'mati --------------------------------------------------------',/) 2005 FORMAT ( &'mati --------------------------------------------------------',/, &'mati Coeff. d echange moyens (calcul T° de peau des parois) ',/, &'mati ------------------------------------------------------ ',/, &'mati Convection forcee ', E12.5 ,' W/m2/°C',/, &'mati Convection naturelle ', E12.5 ,' W/m2/°C',/, &'mati Rayonnement ', E12.5 ,' W/m2/°C',/, &'mati Efficace ', E12.5 ,' W/m2/°C',/, &'mati --------------------------------------------------------',/) C 9001 FORMAT( &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/, &'@ @@ ATTENTION : ARRET MATISSE (MTTSSC) ',/, &'@ ********* ',/, &'@ DEBIT ENTHALPIQUE DES PANACHES INADAPTE. ',/, &'@ ',/, &'@ La modelisation des panaches est activee ',/, &'@ avec IMDCNT = ',I10 ,' ',/, &'@ Le debit enthalpique des panaches de convection ',/, &'@ naturelle doit etre strictement inferieur a la ',/, &'@ puissance totale des colis presents dans l''entrepot ',/, &'@ puisqu''il en represente une partie. ',/, &'@ ',/, &'@ Le debit enthalpique prescrit est DHPCNT =',E12.5 ,/, &'@ La puissance totale reelle est PTOT*FRDTRA =',E12.5 ,/, &'@ ',/, &'@ L''inegalite suivante n''est pas verifiee : ',/, &'@ DHPCNT < PTOT*FRDTRA ',/, &'@ ',E12.5 ,' ',E12.5 ,' ',/, &'@ ',/, &'@ Le calcul ne sera pas execute. ',/, &'@ ',/, &'@ Desactiver la modelisation des panaches ou ',/, &'@ modifier la valeur du debit enthalpique prescrit ou ',/, &'@ modifier la puissance totale (carte de remplissage ou ',/, &'@ puissance des colis). ',/, &'@ ',/, &'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/, &'@ ',/) C C---- C FIN C---- C RETURN C END c@z