C FISS SOURCE CB215821 23/11/02 21:15:03 11779 SUBROUTINE FISS IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C__________________________________________________________________ C operateur FUITE: C Calcul du debit air-vapeur dans une fissure plane C C C__________________________________________________________________ -INC PPARAM -INC CCOPTIO -INC SMELEME -INC SMCHAML -INC SMCOORD -INC SMCHPOI -INC SMLENTI -INC SMMODEL -INC SMTABLE C SEGMENT MAVAL REAL*8 VAL(NVAL) ENDSEGMENT C SEGMENT PTRAV1 INTEGER KASN(NNPT),KG(NNPT) REAL*8 XNC(NNPT),TP(NNPT),EP(NNPT),BP(NNPT) ENDSEGMENT SEGMENT PTRAV2 REAL*8 XX(NMX) REAL*8 XP(NMX),XPV(NMX),XT(NMX),XY(NMX),XU(NMX) REAL*8 XQ(NMX),XQA(NMX),XQW(NMX),XHF(NMX),XRE(NMX),XDH(NMX) REAL*8 ZP(NNPT),ZPV(NNPT),ZT(NNPT),ZY(NNPT),ZU(NNPT) REAL*8 ZQ(NNPT),ZQA(NNPT),ZQW(NNPT),ZHF(NNPT),ZRE(NNPT),ZDH(NMX) REAL*8 ZTRA(NNPT) ENDSEGMENT C INTEGER GPARF CHARACTER*8 CMATE CHARACTER*12 CMAT1,CMAT2,CMAT3 CHARACTER*(NCONCH) CONM CHARACTER*8 MCO3,MCO4,NOM CHARACTER*(LOCOMP) MCO1,MCO2,MOT REAL*8 VALCL(5) POINTEUR MCHMAT.MCHELM POINTEUR PCHMAT.MCHAML POINTEUR DCHMAT.MELVAL POINTEUR PCOMP.MAVAL C CHARACTER*10 MOALIR CHARACTER*8 TYPOBJ CHARACTER*28 MOTTAB CHARACTER*72 CHARRE LOGICAL LOGIN,LOGRE REAL*8 XVALIN,XVALRE C TYPOBJ=' ' KIMP=IIMPI QINI = -1. C*********************************************************************** C lecture obligatoire de l objet modele C IPMODL=0 CALL LIROBJ('MMODEL ',IPMODL,1,IRET) CALL ACTOBJ('MMODEL ',IPMODL,1) IF (IERR.NE.0) RETURN C C*********************************************************************** C lecture obligatoire du champ de proprietes materielles C IPCHA1 = 0 CALL LIROBJ('MCHAML ',IPIN,1,IRET1) CALL ACTOBJ('MCHAML ',IPIN,1) IF (IERR.NE.0) RETURN CALL REDUAF(IPIN,IPMODL,IPCHA1,0,IR,KER) IF(IR .NE. 1) CALL ERREUR(KER) IF(IERR .NE. 0) RETURN C verification de l egalite des supports du modele et du ch prop mat call rayn1(IPMODL,IPCHA1) C C*********************************************************************** C lecture obligatoire de la table des conditions aux limites C IPCHE3 = 0 IPCHA2 = 0 IPCHA3 = 0 IPCHA4 = 0 CALL LIROBJ('TABLE ',IPTAB,1,IRETOU) IF(IERR.NE.0) RETURN MTABLE = IPTAB SEGACT MTABLE NDIMTAB = MLOTAB C MOTTAB='PRESSION_TOTALE_AMONT' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN,IOBIN, . TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHE3) PE = XVALRE MOTTAB='PRESSION_VAPEUR_AMONT' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN,IOBIN, . TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHE3) PVE = XVALRE MOTTAB='TEMPERATURE_AMONT' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN,IOBIN, . TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHE3) TE = XVALRE MOTTAB='PRESSION_TOTALE_AVAL' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN,IOBIN, . TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHE3) PS = XVALRE MOTTAB='TEMPERATURE_PAROI' TYPOBJ=' ' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN,IOBIN, . TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHA2) MOTTAB='OUVERTURE' TYPOBJ=' ' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN,IOBIN, . TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHA3) * IF(NDIMTAB.EQ.7) THEN MOTTAB='DEBIT_INITIAL' TYPOBJ=' ' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN, . IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHE3) IF(TYPOBJ.EQ.'FLOTTANT') THEN QINI = XVALRE ELSE MOTTAB='ETENDUE' TYPOBJ=' ' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN, . IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHA4) ENDIF ENDIF * IF(NDIMTAB.EQ.8) THEN MOTTAB='DEBIT_INITIAL' TYPOBJ=' ' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN, . IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHE3) QINI = XVALRE MOTTAB='ETENDUE' TYPOBJ=' ' CALL ACCTAB(IPTAB,'MOT ',IVALIN,XVALIN,MOTTAB,LOGIN, . IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IPCHA4) ENDIF SEGDES MTABLE C*********************************************************************** C ACTIVATION DU MODELE C -------------------- MMODEL = IPMODL SEGACT MMODEL C NSOUS = KMODEL(/1) IMODEL=KMODEL(1) SEGACT IMODEL C C initialisations IPMAIL= IMAMOD NFOR = FORMOD(/2) NMAT = MATMOD(/2) C C verification de la formulation CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CMATE,MATE,INAT) IF (CMATE.EQ.' ')THEN CALL ERREUR(251) IRET=0 RETURN ENDIF C lecture de la formulation CMAT1 = CMATE(1:2) CMAT2 = CMATE(3:4) CMAT3 = CMATE(5:8) IF (CMAT1.EQ.'MA')THEN XW = 0.D0 ELSE XW = 1.D0 ENDIF IF (CMAT2.EQ.'PA')THEN GPARF = 1 ELSE GPARF = 0 ENDIF C C*********************************************************************** C ACTIVATION DU MCHELM ch prop mat C -------------------- MCHMAT = IPCHA1 SEGACT MCHMAT C NBS = MCHMAT.IMACHE(/1) C pas de boucle sur les modeles elementaires car NSOUS=1 C pas de boucle sur les maillages elementaires car NBS=1 C initialilsations RUG = 0. C lecture des prop materielles PCHMAT = MCHMAT.ICHAML(1) SEGACT PCHMAT NBS2 = PCHMAT.IELVAL(/1) NVAL = NBS2 SEGINI PCOMP C boucle sur le nb de composantes du ch prop mat DO 10 I =1,NBS2 MOT = PCHMAT.NOMCHE(I) DCHMAT = PCHMAT.IELVAL(I) SEGACT DCHMAT C N1PTEL = DCHMAT.VELCHE(/1) C N1EL = DCHMAT.VELCHE(/2) PCOMP.VAL(I) = DCHMAT.VELCHE(1,1) IF(MOT.EQ.'RUGO') RUG = PCOMP.VAL(I) IF(MOT.EQ.'REC ') RECU = PCOMP.VAL(I) IF(MOT.EQ.'FK ') XKUL = PCOMP.VAL(I) IF(MOT.EQ.'FA ') XKUT1 = PCOMP.VAL(I) IF(MOT.EQ.'FB ') XKUT2 = PCOMP.VAL(I) IF(MOT.EQ.'FC ') XKUT3 = PCOMP.VAL(I) IF(MOT.EQ.'FD ') XKUT4 = PCOMP.VAL(I) 10 CONTINUE IF (CMAT3.EQ.'BLAS') THEN RECU = -1D5 XKUL = 0. XKUT1 = 0. XKUT2 = 0. XKUT3 = 0. XKUT4 = 0. ELSE IF (CMAT3.EQ.'COLE')THEN RECU = -1D5 XKUL = 0. XKUT1 = 0. XKUT2 = 0. XKUT3 = 0. XKUT4 = 0. * ELSE IF (CMAT3.EQ.'ENT1')THEN * RECU = PCOMP.VAL(2) ELSE IF (CMAT3.EQ.'ENT2')THEN RECU =-1D0*RECU ELSE IF (CMAT3.EQ.'ENT3')THEN RECU = -1D6 XKUT1 = 0. XKUT2 = 0. XKUT3 = 0. XKUT4 = 0. ELSE IF (CMAT3.EQ.'ENT4')THEN RECU = -1D7 XKUT1 = 0. XKUT2 = 0. XKUT3 = 0. XKUT4 = 0. ENDIF SEGSUP PCOMP C C********************************************************************* C TRAITEMENT des 2 champs ep et tp C -------------------------------- IPT2=IPMAIL SEGACT IPT2 NBSO = IPT2.LISOUS(/1) IF (NBSO.NE.0) THEN CALL ERREUR(25) RETURN ENDIF NEL = IPT2.NUM(/2) NPT2 = IPT2.NUM(/1) IF (NPT2.NE.2) THEN CALL ERREUR(16) RETURN ENDIF C segment de travail 1 NNPT = NEL + 1 SEGINI PTRAV1 XNC(1)=0.D0 KG(1)= IPT2.NUM(1,1) C calcul des abscisses curvilignes SEGACT,MCOORD IF (IDIM.EQ.1) THEN DO I =1,NEL KG(I+1) = IPT2.NUM(2,I) IREF1 = (IDIM+1)*(IPT2.NUM(1,I)-1) IREF2 = (IDIM+1)*(IPT2.NUM(2,I)-1) DX1 = XCOOR(IREF2+1)-XCOOR(IREF1+1) XNC(I+1) = XNC(I) + DX1 ENDDO ENDIF IF (IDIM.EQ.2) THEN DO I =1,NEL KG(I+1) = IPT2.NUM(2,I) IREF1 = (IDIM+1)*(IPT2.NUM(1,I)-1) IREF2 = (IDIM+1)*(IPT2.NUM(2,I)-1) DX1 = XCOOR(IREF2+1)-XCOOR(IREF1+1) DX2 = XCOOR(IREF2+2)-XCOOR(IREF1+2) DXC = SQRT(DX1*DX1+DX2*DX2) XNC(I+1) = XNC(I) + DXC ENDDO ENDIF IF (IDIM.EQ.3) THEN DO I =1,NEL KG(I+1) = IPT2.NUM(2,I) IREF1 = (IDIM+1)*(IPT2.NUM(1,I)-1) IREF2 = (IDIM+1)*(IPT2.NUM(2,I)-1) DX1 = XCOOR(IREF2+1)-XCOOR(IREF1+1) DX2 = XCOOR(IREF2+2)-XCOOR(IREF1+2) DX3 = XCOOR(IREF2+3)-XCOOR(IREF1+3) DXC = SQRT(DX1*DX1+DX2*DX2+DX3*DX3) XNC(I+1) = XNC(I) + DXC ENDDO ENDIF SEGDES,MCOORD XL = XNC(NNPT) DO I =1,NNPT XNC(I)=XNC(I)/XL ENDDO * recherche du plus petit DX DXMINAV = 0.01 DO I =1,NNPT-1 DXNC = XNC(I+1)-XNC(I) DXMIN = MIN(DXNC,DXMINAV) DXMINAV = DXMIN ENDDO DX = DXMIN/5.D0 IF(KIMP.GE.1) THEN WRITE(6,*) 'fiss : taille relative de maille fluide= ',DX ENDIF * segment de travail 2 NMX=3*INT(1.D0/DX) SEGINI PTRAV2 IF (KIMP.EQ.2) THEN WRITE(6,*) 'fiss : XL ',XL WRITE(6,*) (XNC(I),I=1,NNPT) WRITE(6,*) 'fiss : KG(I) = ',(KG(I),I=1,NNPT) ENDIF C activation segment MCHPO1 (chpoint Temperature et ouverture) CALL FUCHPO(IPCHA2,IPCHA3,IRETOU) MCHPO1=IRETOU SEGACT MCHPO1 NSOUPO = MCHPO1.IPCHP(/1) IF(NSOUPO.NE.1) THEN CALL ERREUR(25) RETURN ENDIF MSOUP1 = MCHPO1.IPCHP(1) SEGACT MSOUP1 NC = MSOUP1.NOHARM(/1) IF(NC.NE.2) THEN CALL ERREUR(21) RETURN ENDIF IPT1 = MSOUP1.IGEOC CALL KRIPAD(IPT1,MLENT1) SEGACT IPT1, MLENT1 C NTOT = MLENT1.LECT(/1) MPOVA1 = MSOUP1.IPOVAL SEGACT MPOVA1 NNPT1 = MPOVA1.VPOCHA(/1) C IF (NNPT1.NE.NNPT) THEN CALL ERREUR(21) RETURN ENDIF DO I = 1,NNPT KASN(I) = MLENT1.LECT(KG(I)) ENDDO C remplissage des tableaux classes MCO1=MSOUP1.NOCOMP(1) MCO2=MSOUP1.NOCOMP(2) IF (MCO1.EQ.'T ') THEN IF (MCO2.EQ.'OUV ') THEN DO I=1,NNPT TP(I) = MPOVA1.VPOCHA(KASN(I),1) EP(I) = MPOVA1.VPOCHA(KASN(I),2) ENDDO ELSE MOTERR= MCO2 CALL ERREUR(197) RETURN ENDIF ELSEIF (MCO1.EQ.'OUV ') THEN IF (MCO2.EQ.'T ') THEN DO I=1,NNPT TP(I) = MPOVA1.VPOCHA(KASN(I),2) EP(I) = MPOVA1.VPOCHA(KASN(I),1) ENDDO ELSE MOTERR= MCO2 CALL ERREUR(197) RETURN ENDIF ELSE MOTERR= MCO1 CALL ERREUR(197) RETURN ENDIF IF (IIMPI.EQ.2) THEN write(6,*) 'fiss : TP= ',(TP(I),I=1,NNPT) write(6,*) 'fiss : EP= ',(EP(I),I=1,NNPT) ENDIF C cas etendue unite (champ d etendue non defini par utilisateur) IF (IPCHA4.EQ.0) THEN DO I=1,NNPT BP(I) = 1. ENDDO ENDIF C activation segment MCHPO2 (chpoint Temperature et etendue) IF (IPCHA4.GT.0) THEN CALL FUCHPO(IPCHA2,IPCHA4,IRETOU) MCHPO2=IRETOU SEGACT MCHPO2 NSOUPO = MCHPO2.IPCHP(/1) IF(NSOUPO.NE.1) THEN CALL ERREUR(25) RETURN ENDIF MSOUP2 = MCHPO2.IPCHP(1) SEGACT MSOUP2 NC = MSOUP2.NOHARM(/1) IF(NC.NE.2) THEN CALL ERREUR(21) RETURN ENDIF IPT2 = MSOUP2.IGEOC CALL KRIPAD(IPT2,MLENT1) SEGACT IPT2, MLENT1 C NTOT = MLENT1.LECT(/1) MPOVA2 = MSOUP2.IPOVAL SEGACT MPOVA2 NNPT1 = MPOVA2.VPOCHA(/1) C IF (NNPT1.NE.NNPT) THEN CALL ERREUR(21) RETURN ENDIF DO I = 1,NNPT KASN(I) = MLENT1.LECT(KG(I)) ENDDO C remplissage des tableaux classes MCO1=MSOUP2.NOCOMP(1) MCO2=MSOUP2.NOCOMP(2) IF (MCO1.EQ.'T ') THEN IF (MCO2.EQ.'ETEN') THEN DO I=1,NNPT TP(I) = MPOVA2.VPOCHA(KASN(I),1) BP(I) = MPOVA2.VPOCHA(KASN(I),2) ENDDO ELSE MOTERR= MCO2 CALL ERREUR(197) RETURN ENDIF ELSEIF (MCO1.EQ.'ETEN') THEN IF (MCO2.EQ.'T ') THEN DO I=1,NNPT TP(I) = MPOVA2.VPOCHA(KASN(I),2) BP(I) = MPOVA2.VPOCHA(KASN(I),1) ENDDO ELSE MOTERR= MCO2 CALL ERREUR(197) RETURN ENDIF ELSE MOTERR= MCO1 CALL ERREUR(197) RETURN ENDIF IF (IIMPI.EQ.2) THEN write(6,*) 'fiss : BP= ',(BP(I),I=1,NNPT) ENDIF ENDIF C compatibilite de la loi de frottement avec la valeur de rugosite DO I=1,NNPT RUGD = RUG/(2.*EP(I)) IF (RUGD.GE.1e-4.AND.MATMOD(3).EQ.'POISEU_BLASIUS') THEN MOTERR(1:12)= 'RUGO/2e' REAERR(1) = RUGD REAERR(2) = 0.0001 C %m1:18 = %r1 superieur a %r2 CALL ERREUR(43) RETURN ENDIF IF (RUGD.LT.1e-4.AND.MATMOD(3).EQ.'POISEU_COLEBROOK') THEN MOTERR(1:18)= 'RUGO/2e' REAERR(1) = RUGD REAERR(2) = 0.0001 C %m1:18 = %r1 inferieur a %r2 CALL ERREUR(41) RETURN ENDIF ENDDO C C*********************************************************************** C CALCUL C ------ IF (GPARF.EQ.1) THEN CALL BECALC(PE,PVE,TE,PS,XL,DX,RUG,QINI,XW,NNPT,XNC,TP,EP,BP, $ KIMP,NMX,NX,XX,XP,XPV,XT,XY,XU,XHF,XQ,XQW,XQA,XRE,XDH,RECU, $ XKUL,XKUT1,XKUT2,XKUT3,XKUT4) IF (IERR.NE.0) RETURN ELSE CALL BECALC2(PE,PVE,TE,PS,XL,DX,RUG,QINI,XW,NNPT,XNC,TP,EP,BP, $ KIMP,NMX,NX,XX,XP,XPV,XT,XY,XU,XHF,XQ,XQW,XQA,XRE,XDH,RECU, $ XKUL,XKUT1,XKUT2,XKUT3,XKUT4) IF (IERR.NE.0) RETURN ENDIF C creation des champs nodaux resultats CALL BYRETO(XX,XP,NX,ZP,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XPV,NX,ZPV,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XT,NX,ZT,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XY,NX,ZY,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XU,NX,ZU,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XHF,NX,ZHF,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XQ,NX,ZQ,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XQA,NX,ZQA,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XQW,NX,ZQW,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XRE,NX,ZRE,XNC,NNPT,KASN,ZTRA) CALL BYRETO(XX,XDH,NX,ZDH,XNC,NNPT,KASN,ZTRA) C ecriture du champoint multicomposante NC = 11 N = NNPT SEGINI MPOVAL DO I=1,NNPT VPOCHA(I,1) = ZP(I) VPOCHA(I,2) = ZPV(I) VPOCHA(I,3) = ZT(I) VPOCHA(I,4) = ZY(I) VPOCHA(I,5) = ZU(I) VPOCHA(I,6) = ZHF(I) VPOCHA(I,7) = ZQ(I) VPOCHA(I,8) = ZQA(I) VPOCHA(I,9) = ZQW(I) VPOCHA(I,10) = ZRE(I) VPOCHA(I,11) = ZDH(I) ENDDO SEGSUP PTRAV1 SEGSUP PTRAV2 SEGINI MSOUPO NOCOMP(1) = 'P ' NOCOMP(2) = 'PV ' NOCOMP(3) = 'TF ' NOCOMP(4) = 'X ' NOCOMP(5) = 'U ' NOCOMP(6) = 'H ' NOCOMP(7) = 'Q ' NOCOMP(8) = 'QA ' NOCOMP(9) = 'QE ' NOCOMP(10)= 'RE ' NOCOMP(11)= 'F ' IGEOC = IPT1 IPOVAL = MPOVAL DO I=1,NC NOHARM(I) = MSOUP1.NOHARM(1) ENDDO NSOUPO = 1 NAT = MCHPO1.JATTRI(/1) SEGINI MCHPOI MTYPOI = MCHPO1.MTYPOI MOCHDE = MCHPO1.MOCHDE DO I=1,NAT JATTRI(I) = MCHPO1.JATTRI(1) ENDDO IPCHP(1) = MSOUPO IFOPOI = MCHPO1.IFOPOI CALL ACTOBJ('CHPOINT ',MCHPOI,1) CALL ECROBJ('CHPOINT ',MCHPOI) END