! É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 DES ORDONNEES DISCRETES A L'EQUILIBRE RADIATIF ! 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 MOD_EQR ! *** CE PROGRAMME CALCULE LE CHAMP DE TEMPERATURE T ET LES FLUX ! *** PARIETAUX QP DANS UNE CAVITE RECTANGULAIRE A PAROIS NOIRES ! *** CONTENANT UN MILIEU GRIS, NON DIFFUSANT, A L'EQUILIBRE RA- ! *** DIATIF. LE RAYONNEMENT INCIDENT G EST CALCULE PAR LA METHODE ! *** DES ORDONNEES DISCRETES IMPLICIT NONE INTEGER I,J,IMAX,ITER,ITMAX,JMAX,N,MSH(3),MMX(3) DOUBLE PRECISION ALPHA,DX,DY,KAPPA,ECART,TM,TOL,XMAX,YMAX DOUBLE PRECISION TP(1:4),QP(1:100,1:4),G(1:100,1:100),T(0:101,0:101),QDI(40,3) ! ORDRE DE LA QUADRATURE N=6 ! PARAMETRES DU PROCESSUS DE CONVERGENCE: NOMBRE D'ITÉRATIONS ET TOLÉRANCE ITMAX=100 TOL=0.0001D0 ! COEFFICIENT D'INTERPOLATION (STEP: ALPHA=1; DIAMANT: ALPHA=0.5D0) ALPHA=0.5D0 ! LONGUEUR L (SUIVANT X) XMAX=1.D0 ! HAUTEUR H (SUIVANT Y) YMAX=1.D0 ! TEMPERATURE DES PAROIS X=0, X=L, Y=0 ET Y=H TP(1)=0.D0;TP(2)=0.D0;TP(3)=1.D0;TP(4)=0.D0 ! COEFFICIENT D ABSORPTION DU MILIEU KAPPA=1.D0 ! NOMBRE DE MAILLES SUIVANT X IMAX=40 ! NOMBRE DE MAILLES SUIVANT Y JMAX=40 DX=XMAX/IMAX DY=YMAX/JMAX !=========================================================================== ! DONNEES DE LA QUADRATURE (LIMITEES A S6) ! COSINUS DIRECTEURS ET POIDS DE LA DIRECTION M: MU=QDI(MM,1); ETA=QDI(MM,2); WM=QDI(MM,3) !=========================================================================== DATA MMX/4,12,24/,MSH/0,4,16/ DATA QDI/ & -.5773503D0,-.5773503D0,0.5773503D0,0.5773503D0,-.9082483D0,-.9082483D0, & -.2958759D0,-.2958759D0,-.2958759D0,-.2958759D0,0.2958759D0,0.2958759D0, & 0.2958759D0,0.2958759D0,0.9082483D0,0.9082483D0,-.9656013D0,-.9656013D0, & -.6950514D0,-.6950514D0,-.6950514D0,-.6950514D0,-.1838670D0,-.1838670D0, & -.1838670D0,-.1838670D0,-.1838670D0,-.1838670D0,0.1838670D0,0.1838670D0, & 0.1838670D0,0.1838670D0,0.1838670D0,0.1838670D0,0.6950514D0,0.6950514D0, & 0.6950514D0,0.6950514D0,0.9656013D0,0.9656013D0,-.5773503D0,0.5773503D0, & -.5773503D0,0.5773503D0,-.2958759D0,0.2958759D0,-.9082483D0,-.2958759D0, & 0.2958759D0,0.9082483D0,-.9082483D0,-.2958759D0,0.2958759D0,0.9082483D0, & -.2958759D0,0.2958759D0,-.1838670D0,0.1838670D0,-.6950514D0,-.1838670D0, & 0.1838670D0,0.6950514D0,-.9656013D0,-.6950514D0,-.1838670D0,0.1838670D0, & 0.6950514D0,0.9656013D0,-.9656013D0,-.6950514D0,-.1838670D0,0.1838670D0, & 0.6950514D0,0.9656013D0,-.6950514D0,-.1838670D0,0.1838670D0,0.6950514D0, & -.1838670D0,0.1838670D0,3.1415927D0,3.1415927D0,3.1415927D0,3.1415927D0, & 1.0471976D0,1.0471976D0,1.0471976D0,1.0471976D0,1.0471976D0,1.0471976D0, & 1.0471976D0,1.0471976D0,1.0471976D0,1.0471976D0,1.0471976D0,1.0471976D0, & 0.3219034D0,0.3219034D0,0.7252938D0,0.7252938D0,0.7252938D0,0.7252938D0, & 0.3219034D0,0.7252938D0,0.3219034D0,0.3219034D0,0.7252938D0,0.3219034D0, & 0.3219034D0,0.7252938D0,0.3219034D0,0.3219034D0,0.7252938D0,0.3219034D0, & 0.7252938D0,0.7252938D0,0.7252938D0,0.7252938D0,0.3219034D0,0.3219034D0/ !=========================================================================== ! INITIALISATION DE LA TEMPERATURE !=========================================================================== T=(TP(1)+TP(2)+TP(3)+TP(4))/4.D0 T(1:IMAX,0)=TP(3) T(1:IMAX,JMAX+1)=TP(4) T(0,1:JMAX)=TP(1) T(IMAX+1,1:JMAX)=TP(2) !=========================================================================== ! PROCESSUS ITERATIF POUR OBTENIR L'EQUILIBRE RADIATIF !=========================================================================== DO 4 ITER=1,ITMAX CALL MOD2D(N,IMAX,JMAX,DX,DY,KAPPA,ALPHA,T,G,QP,QDI,MSH,MMX,ITER) IF(ITER.EQ.1) PRINT *,' PROCESSUS ITERATIF' ECART=0.D0 DO 41 I=1,IMAX DO 41 J=1,JMAX TM=(G(I,J)/4.d0)**0.25D0 ECART =DMAX1(ECART,DABS(TM-T(I,J))) T(I,J)=TM 41 CONTINUE IF (ECART.LE.TOL) GOTO 5 PRINT *,' ITER=',ITER,' ECART=',ECART !PAUSE 4 CONTINUE WRITE (*,*) ' *** NON CONVERGENCE : MAX.ITER =',ITMAX 5 CONTINUE !=========================================================================== ! RESULTATS !=========================================================================== ! CHAMP DE TEMPERATURE (NORMALISE PAR T0) OPEN (UNIT=11,FILE='CHAMP_T.DAT') WRITE(11,*) 'TITLE="CHAMPS_T"' WRITE(11,*) 'VARIABLES="X","Y","T"' WRITE(11,*)'ZONE I=',IMAX,' J=',JMAX,'F=POINT' WRITE(11,2603) (( (I-0.5D0)*DX, (J-0.5D0)*DY, T(I,J),I=1,IMAX),J=1,JMAX) ! FLUX INCIDENT SUR LES PAROIS VERTICALES'//' (NORMALISE PAR SIGMA.T0**4) OPEN (UNIT=12,FILE='QINC_VERTICAL.DAT') WRITE(12,*) 'TITLE="QINC_VERTICAL"' WRITE(12,*) 'VARIABLES="Y","QPW","QPE"' WRITE(12,2603) ((J-0.5D0)*DY,QP(J,1),QP(J,2),J=1,JMAX) ! FLUX INCIDENT SUR LES PAROIS VERTICALES'//' (NORMALISE PAR SIGMA.T0**4) OPEN (UNIT=13,FILE='QINC_HORIZONTAL.DAT') WRITE(12,*) 'TITLE="QINC_HORIZONTAL"' WRITE(12,*) 'VARIABLES="X","QPS","QPN"' WRITE(13,2603) ((I-0.5D0)*DX,QP(I,3),QP(I,4),I=1,IMAX) !=========================================================================== ! FORMATS !=========================================================================== 2603 FORMAT (3(E11.4,1X)) !=========================================================================== STOP END ! SUBROUTINE MOD2D(N,IMAX,JMAX,DX,DY,KAPPA,ALPHA,T,G,QP,QDI,MSH,MMX,ITER) ! *** RESOLUTION DE L'EQUATION DE TRANSFERT RADIATIF PAR LA MOD *** ! *** N : ORDRE N LA QUADRATURE SN (2,4 OU 6) (ENTREE) ! *** IMAX : NOMBRE DE VALEURS DISCRETES DE X (ENTREE) ! *** JMAX : NOMBRE DE VALEURS DISCRETES DE Y (ENTREE) ! *** DX : LARGEUR DE MAILLE SUIVANT X (ENTREE) ! *** DY : LARGEUR DE MAILLE SUIVANT Y (ENTREE) ! *** KAPPA : COEFFICIENT D'ABSORPTION DU MILIEU (ENTREE) ! *** ALPHA : COEFFICIENT D'INTERPOLATION (ENTREE) ! *** T : CHAMP DE TEMPERATURE : T(I,J) (ENTREE) ! *** G : CHAMP DE RAYONNEMENT INCIDENT : G(I,J) (SORTIE) ! *** QP : DENSITE DE FLUX RADIATIF INCIDENT AUX PAROIS ! *** QP(I-OU-J,1-2-3-OU-4) (SORTIE) IMPLICIT NONE INTEGER KI,KJ,I,I1,I2,IMAX,J,J1,J2,JMAX,M,MM,M0,MMAX,PX,PY,N,ITER INTEGER MSH(3),MMX(3) DOUBLE PRECISION AETA,AMU,AX,AY,DV,ETA,KAPPA,L0,LP,LS,LW,LX,LY,MU,PI,SP DOUBLE PRECISION DX,DY,ALPHA,WM,XI DOUBLE PRECISION LJM(1:100),G(1:100,1:100),QP(1:100,1:4),T(0:101,0:101) DOUBLE PRECISION QDI(40,3) ! PI=DACOS(-1.D0) ! INITIALISATION G=0.D0 QP=0.D0 ! ORDRE DE LA QUADRATURE IF (N.GT.6.OR.N.LT.2) N=6 M0=MSH(N/2) MMAX=MMX(N/2) !=========================================================================== ! CONTROLE !=========================================================================== IF(ITER.EQ.1) THEN IF (N.EQ.2) PRINT *,' QUADRATURE S2' IF (N.EQ.4) PRINT *,' QUADRATURE S4' IF (N.EQ.6) PRINT *,' QUADRATURE S6' WRITE(*,11)'M',' MU ',' ETA ',' XI ',' POIDS' DO 1 M=1,MMAX MM=M0+M MU=QDI(MM,1) ETA=QDI(MM,2) XI=SQRT(1.-MU*MU-ETA*ETA) WM=QDI(MM,3) WRITE(*,12) M,MU,ETA,XI,WM 1 CONTINUE 11 FORMAT (70('-')/A10,4(5X,A10)/70('-')) 12 FORMAT (I10,4(5X,F10.7)) ENDIF !=========================================================================== ! BOUCLE SUR L'ENSEMBLE DES DIRECTIONS DE LA QUADRATURE !=========================================================================== DO 2 M=1,MMAX ! COSINUS DIRECTEURS ET POIDS DE LA DIRECTION M MM=M0+M MU=QDI(MM,1) ETA=QDI(MM,2) WM=QDI(MM,3) AMU=dABS(MU) AETA=dABS(ETA) ! SENS DE L'INTEGRATION SPATIALE EN FONCTION DES SIGNES DE MU ET ETA IF (MU.GT.0.d0) THEN I1=1 I2=IMAX KI=1 PX=2 ELSE I1=IMAX I2=1 KI=-1 PX=1 END IF IF (ETA.GT.0.d0) THEN J1=1 J2=JMAX KJ=1 PY=4 ELSE J1=JMAX J2=1 KJ=-1 PY=3 END IF ! VALEUR DE L SUR LA PAROI HORIZONTALE DE DEPART DO 21 I=I1,I2,KI 21 LJM(I)=T(I,J1-KJ)**4/PI ! INTEGRATIONS SUCCESSIVES LE LONG DES LIGNES HORIZONTALES DO 22 J=J1,J2,KJ ! VALEUR DE L SUR LA PAROI VERTICALE DE DEPART LW = T(I1-KI,J)**4/PI ! INTEGRATION LE LONG DE LA LIGNE J DO 221 I=I1,I2,KI ! DISCRETISATION DE L'ETR ! - LP EST CALCULE D'APRES LES VALEURS CONNUES DE LS ET LW ! - LES VALEURS DE LN ET LS (EN FAIT, LES NOUVELLES VALEURS ! DE LS ET LW) SONT OBTENUES PAR EXTRAPOLATION AX=DY AY=DX DV=AX*AY LX=AX*AMU/ALPHA LY=AY*AETA/ALPHA L0=DV*KAPPA SP=KAPPA*T(I,J)**4/PI*DV LS=LJM(I) LP=(LX*LW+LY*LS+SP)/(LX+LY+L0) LW=LW+(LP-LW)/ALPHA LS=LS+(LP-LS)/ALPHA LJM(I)=LS ! LA CONTRIBUTION DE LA DIRECTION M EST AJOUTEE : ! - A LA VALEUR LOCALE DU RAYONNEMENT INCIDENT G ! - AU FLUX INCIDENT SUR LA PAROI VERTICALE D'ARRIVEE G(I,J)=G(I,J)+WM*LP IF(J.EQ.J2) QP(I,PY)=QP(I,PY)+AETA*WM*LS 221 CONTINUE ! FLUX INCIDENT SUR LA PAROI VERTICALE D'ARRIVEE QP(J,PX)=QP(J,PX)+AMU*WM*LW 22 CONTINUE 2 CONTINUE !=========================================================================== RETURN END