! ÉCOLE THÉMATIQUE INCENDIE – 30 MAI AU 04 JUIN 2015 ! ! CALCUL DU RAYONNEMENT ÉMIS PAR UNE SOURCE PONCTUELLE DANS UNE CAVITÉ 2D ! PAR LA MÉTHODE DE MONTE CARLO ! CONTACTS: MATTHIEU.DEGENNARO@NOVELTIS.FR ! MAXIME.MENSE@GMAIL.COM ! BERNARD.PORTERIE@UNIV-AMU.FR ! AIX-MARSEILLE UNIVERSITÉ, CNRS, IUSTI UMR 7343, 13453 MARSEILLE, FRANCE ! PROGRAM MMC2D ! RAYONNEMENT DANS UNE CAVITÉ 2D ! ON PLACE UNE SOURCE RADIATIVE À L’INTÉRIEUR D’UNE CAVITÉ 2D CONTENANT UN ! MILIEU SEMI-TRANSPARENT (MST) ABSORBANT-DIFFUSANT, D’ALBÉDO OMEGA ET DE ! COEFFICIENT D’EXTINCTION BETA. ! L’OBJECTIF EST DE CALCULER LE CHAMP DE SOURCES RADIATIVES DANS LE MST (W/M3), ! TERME SOURCE DE L’ÉQUATION DE BILAN D’ÉNERGIE, ET LES FLUX RADIATIFS (W/M2) ! AUX PAROIS. ! LES PAROIS SONT GRISES. ON NOTERA EPS1,EPS2,EPS3,EPS4 LES ! ÉMISSIVITÉS DES PAROIS OUEST (NPAROI=1), EST (NPAROI=2), SUD (NPAROI=3) ! ET NORD (NPAROI=4). ON DÉCOUPE LE DOMAINE DE CALCUL BIDIMENSIONNEL EN ! IMAX*JMAX CELLULES ÉLÉMENTAIRES DE DIMENSIONS ∆X*∆Y. ! IMPLICIT NONE !- MAILLAGE DE LA CAVITÉ : IMAX, JMAX INTEGER, PARAMETER :: IMAX=41, JMAX=41 !- NOMBRE DE QUANTA: NQUANTA INTEGER, PARAMETER :: NMAX=10000000 INTEGER I,J,IS,JS,N,IARR,JARR,NPAROI,NQ(IMAX,JMAX),NQ1,NQ2,NQ3,NQ4 DOUBLE PRECISION XL,YL,DX,DY,X(IMAX),Y(JMAX),EPS1,EPS2,EPS3,EPS4,OMEGA,BETA,PI DOUBLE PRECISION XDEP,YDEP,XARR,YARR,TETA,RS,DS,PHI,ALPHA,DSINTETA,PE,Q0 DOUBLE PRECISION RPHI,RTETA,REPS,ROMEGA !=============== XL=4.D0 ! LONGUEUR DE LA CAVITÉ YL=4.D0 ! LARGEUR DE LA CAVITÉ EPS1=0.6D0;EPS2=0.6D0;EPS3=0.6D0;EPS4=0.6D0 ! EMISSIVITÉ DES PAROIS OMEGA=0.5D0 ! ALBÉDO BETA=0.5D0 ! COEFFICIENT D’EXTINCTION PE=1.D0 ! POUVOIR ÉMISSIF DE LA SOURCE IS=11; JS=21 ! INDICES DE POSITION DE LA SOURCE PI=DACOS(-1.D0) Q0=PE/DFLOAT(NMAX)! PUISSANCE TRANSPORTÉE PAR UN QUANTUM !====================================================================== ! MAILLAGE DE LA CAVITÉ ! ON NOTERA X(I) ET Y(J) LES COORDONNÉES DES CENTRES DES MAILLES !====================================================================== DX=XL/DFLOAT(IMAX) DO I=1,IMAX X(I)=(DFLOAT(I)-0.5D0)*DX ENDDO DY=YL/DFLOAT(JMAX) DO J=1,JMAX Y(J)=(DFLOAT(J)-0.5D0)*DY ENDDO !====================================================================== ! INITIALISATION DU GÉNÉRATEUR DE NOMBRES ALÉATOIRES. !====================================================================== CALL RANDOM_SEED !====================================================================== ! MISE A ZERO DES COMPTEURS DE QUANTA AUX PAROIS !====================================================================== NQ=0; NQ1=0; NQ2=0; NQ3=0; NQ4=0 !====================================================================== ! BOUCLE QUANTA !====================================================================== DO 2 N=1,NMAX IF(MOD(N,NMAX/10).EQ.0) PRINT *,'N=',N XDEP=X(IS) ! COORDONNÉES DE DÉPART (SOURCE) YDEP=Y(JS) NPAROI=0 CALL RANDOM_NUMBER(RTETA) TETA=DACOS(1.D0-2.D0*RTETA) ! ANGLE D'ÉMISSION CALL RANDOM_NUMBER(RPHI) PHI=2.D0*PI*RPHI ! ANGLE D'ÉMISSION CALL RANDOM_NUMBER(RS) DS=-1.D0/BETA*DLOG(RS) ! LONGUEUR PARCOURUE AVANT EXTINCTION XARR=XDEP+DCOS(TETA)*DS ! COORDONNÉES D'ARRIVÉE YARR=YDEP+DSIN(TETA)*DCOS(PHI)*DS 3 CONTINUE ! ! LE QUANTUM ATTEINT LA PAROI OUEST (NPAROI=1), EST (NPAROI=2), ! SUD (NPAROI=3), NORD (NPAROI=4). DANS LE CAS CONTRAIRE, LE ! QUANTUM N'ATTEINT AUCUNE PAROI (NPAROI=0). CALL NUMERO_PAROI(PHI,XDEP,YDEP,XARR,YARR,XL,YL,NPAROI) ! IF(NPAROI.EQ.1) THEN ! LE QUANTUM ATTEINT LA PAROI OUEST CALL RANDOM_NUMBER(REPS) IF(REPS.LE.EPS1) THEN ! ABSORPTION NQ1=NQ1+1 GOTO 2 ELSE ! RÉFLEXION YDEP=YDEP-(XDEP-0.D0)*(YDEP-YARR)/(XDEP-XARR) XDEP=0.D0 CALL RANDOM_NUMBER(RTETA) TETA=DASIN(DSQRT(RTETA)) CALL RANDOM_NUMBER(RPHI) PHI=2.D0*PI*RPHI CALL RANDOM_NUMBER(RS) DS=-1.D0/BETA*DLOG(RS) ! LONGUEUR PARCOURUE AVANT EXTINCTION XARR=XDEP+DCOS(TETA)*DS ! COORDONNÉES D'ARRIVÉE YARR=YDEP+DSIN(TETA)*DCOS(PHI)*DS GOTO 3 END IF ENDIF ! IF(NPAROI.EQ.2) THEN ! LE QUANTUM ATTEINT LA PAROI EST CALL RANDOM_NUMBER(REPS) IF(REPS.LE.EPS2) THEN ! ABSORPTION NQ2=NQ2+1 GOTO 2 ELSE ! RÉFLEXION YDEP=YDEP-(XDEP-XL)*(YDEP-YARR)/(XDEP-XARR) XDEP=XL CALL RANDOM_NUMBER(RTETA) TETA=DASIN(DSQRT(RTETA)) CALL RANDOM_NUMBER(RPHI) PHI=2.D0*PI*RPHI CALL RANDOM_NUMBER(RS) DS=-1.D0/BETA*DLOG(RS) ! LONGUEUR PARCOURUE AVANT EXTINCTION XARR=XDEP-DCOS(TETA)*DS ! COORDONNÉES D'ARRIVÉE YARR=YDEP+DSIN(TETA)*DSIN(PHI)*DS GOTO 3 END IF ENDIF ! IF(NPAROI.EQ.3) THEN ! LE QUANTUM ATTEINT LA PAROI SUD CALL RANDOM_NUMBER(REPS) IF(REPS.LE.EPS3) THEN ! ABSORPTION NQ3=NQ3+1 GOTO 2 ELSE ! RÉFLEXION XDEP=XDEP-(YDEP-0.D0)*(XDEP-XARR)/(YDEP-YARR) YDEP=0.D0 CALL RANDOM_NUMBER(RTETA) TETA=DASIN(DSQRT(RTETA)) CALL RANDOM_NUMBER(RPHI) PHI=2.D0*PI*RPHI CALL RANDOM_NUMBER(RS) DS=-1.D0/BETA*DLOG(RS) ! LONGUEUR PARCOURUE AVANT EXTINCTION XARR=XDEP+DSIN(TETA)*DSIN(PHI)*DS ! COORDONNÉES D'ARRIVÉE YARR=YDEP+DCOS(TETA)*DS GOTO 3 END IF ENDIF ! IF(NPAROI.EQ.4) THEN ! LE QUANTUM ATTEINT LA PAROI NORD CALL RANDOM_NUMBER(REPS) IF(REPS.LE.EPS4) THEN ! ABSORPTION NQ4=NQ4+1 GOTO 2 ELSE ! RÉFLEXION XDEP=XDEP-(YDEP-YL)*(XDEP-XARR)/(YDEP-YARR) YDEP=YL CALL RANDOM_NUMBER(RTETA) TETA=DASIN(DSQRT(RTETA))!-PI/2.D0 CALL RANDOM_NUMBER(RPHI) PHI=2.D0*PI*RPHI CALL RANDOM_NUMBER(RS) DS=-1.D0/BETA*DLOG(RS) ! LONGUEUR PARCOURUE AVANT EXTINCTION XARR=XDEP+DSIN(TETA)*DSIN(PHI)*DS ! COORDONNÉES D'ARRIVÉE YARR=YDEP-DCOS(TETA)*DS GOTO 3 END IF ENDIF ! ! LE QUANTUM N'A RENCONTRE AUCUNE PAROI CALL RANDOM_NUMBER(ROMEGA) ! DIFFUSION PAR LE MILIEU IF (ROMEGA.LT.OMEGA) THEN XDEP=XARR YDEP=YARR CALL RANDOM_NUMBER(RTETA) TETA=DACOS(1.D0-2.D0*RTETA) CALL RANDOM_NUMBER(RPHI) PHI=2.D0*PI*RPHI CALL RANDOM_NUMBER(RS) DS=-1.D0/BETA*DLOG(RS) ! LONGUEUR PARCOURUE AVANT EXTINCTION XARR=XDEP+DCOS(TETA)*DS ! COORDONNÉES D'ARRIVÉE YARR=YDEP+DSIN(TETA)*DCOS(PHI)*DS GOTO 3 END IF ! ABSORPTION PAR LE MILIEU IARR=INT(XARR/DX)+1 JARR=INT(YARR/DY)+1 NQ(IARR,JARR)=NQ(IARR,JARR)+1 ! 2 CONTINUE !====================================================================== ! IMPRESSION DU NOMBRE DE QUANTA ABSORBÉS PAR LES PAROIS !====================================================================== PRINT *,'NQ1 OUEST=',NQ1*Q0/DY PRINT *,'NQ2 EST =',NQ2*Q0/DY PRINT *,'NQ3 SUD =',NQ3*Q0/DX PRINT *,'NQ4 NORD =',NQ4*Q0/DX !====================================================================== ! ENREGISTREMENT SUR FICHIER DU NOMBRE DE QUANTA ABSORBÉS PAR LE MILIEU !====================================================================== OPEN (2,FILE='MC2D.DAT') WRITE(2,*)'ZONE I=',IMAX,' J=',JMAX,' F=POINT' WRITE(2,9) (((I-0.5D0)*DX,(J-0.5D0)*DY,NQ(I,J)*Q0/(DX*DY),I=1,IMAX),J=1,JMAX) ! 9 FORMAT (3(E11.4,1X)) ! !====================================================================== STOP END ! SUBROUTINE NUMERO_PAROI(PHI,XDEP,YDEP,XARR,YARR,XL,YL,NPAROI) IMPLICIT NONE INTEGER NPAROI DOUBLE PRECISION PHI,XL,YL,XDEP,YDEP,XARR,YARR,PI DOUBLE PRECISION PHI1,PHI2,PHI3,PHI4,X1,Y1 !=============== PI=DACOS(-1.D0) ! PHI1=DATAN((YL-YDEP)/(XL-XDEP)) PHI2=PI/2.D0+DATAN(XDEP/(YL-YDEP)) PHI3=PI+DATAN(YDEP/XDEP) PHI4=3.D0*PI/2.D0+DATAN((XL-XDEP)/YDEP) X1=XARR-XDEP Y1=YARR-YDEP PHI=ATAN2(Y1,X1) IF(phi.LT.0.D0) phi=2.d0*pi+phi IF(PHI.GE.0.D0.AND.PHI.LT.PHI1.OR.PHI.GT.PHI4.AND.PHI.LT.2.D0*PI)THEN NPAROI=2 IF(XARR.LT.XL) NPAROI=0 ENDIF IF(PHI.GE.PHI1.AND.PHI.LT.PHI2) THEN NPAROI=4 IF(YARR.LT.YL) NPAROI=0 ENDIF IF(PHI.GE.PHI2.AND.PHI.LT.PHI3) THEN NPAROI=1 IF(XARR.GT.0.D0) NPAROI=0 ENDIF IF(PHI.GE.PHI3.AND.PHI.LT.PHI4) THEN NPAROI=3 IF(YARR.GT.0.D0) NPAROI=0 ENDIF ! RETURN END