C CLMI SOURCE CB215821 23/01/25 21:15:07 11573 SUBROUTINE CLMI PARAMETER (NTB=1) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 ctype,ferm C*********************************************************************** C C C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMELEME -INC SMCOORD -INC SIZFFB -INC SMMATRIK POINTEUR IMAT1.IMATRI,IMAT2.IMATRI,IMAT3.IMATRI -INC SMCHPOI POINTEUR IZTU1.MPOVAL,MCHPO9.MCHPOI,MSOUP9.MSOUPO,MPOVA9.MPOVAL POINTEUR MCHPO8.MCHPOI,MSOUP8.MSOUPO,MPOVA8.MPOVAL POINTEUR MCHPOA.MCHPOI,MSOUPA.MSOUPO,MPOVAA.MPOVAL POINTEUR MCHPOB.MCHPOI,MSOUPB.MSOUPO,MPOVAB.MPOVAL SEGMENT TRAV REAL*8 UN(NPTI,IDIM),TN(NPTI,1) REAL*8 WT(NP,NPG),WS(NP,NPG) REAL*8 UIL(IDIM,NPG),DUIL(IDIM,NPG),ANUK(1),AIRE(1) REAL*8 HNM(NPTI),BNM(NPTI),Coef1(NPTI),Coef2(NPTI),DNM(NPTI) REAL*8 H32NM(NPTI),CFN(NPTI),CEN(NPTI),DeN(NPTI),D1N(NPTI) REAL*8 HHN(NPTI) REAL*8 UCL(NPTI),Delta(NPTI),SM1(NPTI),SM2(NPTI) REAL*8 FP(NPTI),PHIP(NPTI),VT(NPTI) ENDSEGMENT -INC SMLENTI -INC SMLMOTS POINTEUR LINCO.MLMOTS real*8 nu CHARACTER*8 MTYP,NOME,NOMI,TYPE,NOMZ,TYPC CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB),IXV(9) DATA LTAB/'KIZX '/ C***************************************************************************** C C write(6,*)' Debut CLMI ' CALL LITABS(LTAB,KTAB,NTB,1,IRET) IF (IERR.NE.0) RETURN MTABX=KTAB(1) C C- Récupération de la table EQEX (pointeur MTAB1) C CALL LEKTAB(MTABX,'EQEX',MTAB1) IF(MTAB1.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' EQEX ' MOTERR( 9:16) = ' EQEX ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF CALL ACME(MTAB1,'NAVISTOK',NASTOK) C C- Récupération de la table INCO (pointeur KINC) C CALL LEKTAB(MTAB1,'INCO',KINC) IF(KINC.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' INCO ' MOTERR( 9:16) = ' INCO ' MOTERR(17:24) = ' EQEX ' CALL ERREUR(786) RETURN ENDIF C***************************************************************************** C OPTIONS C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG IAXI=0 IF(IFOMOD.EQ.0)IAXI=2 C C- Récupération de la table des options KOPT (pointeur KOPTI) C CALL LEKTAB(MTABX,'KOPT',KOPTI) IF (KOPTI.EQ.0) THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' KOPT ' MOTERR( 9:16) = ' KOPT ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF C write(6,*)' Avant les options ' CALL ACME(KOPTI,'IDCEN',IDCEN) KDIM=1 IF(IDCEN.EQ.2)KDIM=IDIM CALL ACME(KOPTI,'IKOMP',IKOMP) CALL ACME(KOPTI,'KIMPL',KIMPL) CALL ACME(KOPTI,'KPOIN',KPRE) CALL ACME(KOPTI,'KFORM',KFORM) CALL ACMF(KOPTI,'CMD',CMD) IF(KFORM.NE.0.AND.KFORM.NE.1)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = 'EF/EFM1 ' CALL ERREUR(803) RETURN ENDIF CALL ACMF(KOPTI,'AIMPL',AIMPL) AG=AIMPL AD=AIMPL-1.D0 IF (IERR.NE.0) RETURN C write(6,*)' Apres les options ' C***************************************************************************** C C- Récupération de la table DOMAINE associée au domaine local C c write(6,*)' Recuperation du modele ' CALL ACMM(MTABX,'NOMZONE',NOMZ) TYPE=' ' CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) IF(TYPE.NE.'MMODEL')THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' DOMZ ' MOTERR( 9:16) = ' DOMZ ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF C E/ MMODEL : Pointeur de la table contenant l'information cherchée C /S IPOINT : Pointeur sur la table DOMAINE C /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE C INEFMD=4 LINB CALL LEKMOD(MMDZ,MTABZ,INEFMD) C verifier que la formulation est bonne INEFMD = 1 CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) CALL LEKTAB(MTABZ,'SOMMET',MELEMS) C CALL LEKTAB(MTABZ,'CENTRE',MELEMC) CALL KRIPAD(MELEMS,MLENTI) IF (IERR.NE.0) RETURN C************************************************************************* C VERIFICATIONS SUR LES INCONNUES C C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI C c write(6,*)' Recuperation inconnues ' TYPE='LISTMOTS' CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) IF (IERR.NE.0) RETURN SEGACT LINCO NBINC=LINCO.MOTS(/2) IF(NBINC.NE.1)THEN C Indice %m1:8 : contient plus de %i1 %m9:16 MOTERR( 1:8) = 'LISTINCO' INTERR(1) = 1 MOTERR(9:16) = ' MOTS ' CALL ERREUR(799) RETURN ENDIF NOMI=LINCO.MOTS(1) c write(6,*) 'Nb inco=',NBINC c write(6,*) 'NOMI=',NOMI C C- Récupération de l'inconnue c write(6,*) 'Récupération de l''inconnue' TYPE=' ' CALL ACMO(KINC,NOMI,TYPE,MCHPO1) IF(TYPE.NE.'CHPOINT ')THEN C Indice %m1:8 : ne contient pas un objet de type %m9:16 MOTERR( 1: 8) = 'INC '//NOMI MOTERR( 9:16) = 'CHPOINT ' CALL ERREUR(800) RETURN ELSE CALL LICHT(MCHPO1,IZTU1,TYPC,MELEMI) NINKO = IZTU1.VPOCHA(/2) NPTI = IZTU1.VPOCHA(/1) c write (6,*) 'NINKO =',NINKO c write(6,*)' MCHPOI,MELEMI=',MCHPOI,MELEMI C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes ENDIF C***************************************************************************** C Le domaine de definition est donne par le SPG de la premiere inconnue C Les inconnues suivantes devront posseder ce meme pointeur C On verifie que les points de la zone sont tous inclus dans ce SPG CALL KRIPAD(MELEMI,IPADI) CALL VERPAD(IPADI,MELEME,IRET) IF(IRET.NE.0)THEN C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique MOTERR(1: 8) = 'INC '//NOMI MOTERR(9:16) = 'CHPOINT ' CALL ERREUR(788) RETURN ENDIF C***************************************************************************** NU=1.d-6 C************************************************************************* C Lecture des coefficients C Type du coefficient : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur c write(6,*)' lecture coef ' CALL ACME(MTABX,'IARG',IARG) IF(IARG.NE.6)THEN WRITE(6,*)' Operateur DETO : ' WRITE(6,*)' Nombre d''arguments incorrect : ',IARG WRITE(6,*)' On attend 6 ' C Indice %m1:8 : nombre d'arguments incorrect MOTERR(1:8) = 'IARG ' CALL ERREUR(804) RETURN ENDIF C--Cas incompréssible c-------1er coefficient/type de fermeture IXV(1)=0 IXV(2)=1 IXV(3)=0 CALL LEKCOF('Opérateur DETO :', & MTABX,KINC,1,IXV,IZTG1,MPOVA1,NPT1,NC1,IK1,IRET) IF(IRET.EQ.0)RETURN ferm = MPOVA1.VPOCHA(1,1) c------ 2eme coef/n°de l'équation traitée IXV(1)=0 IXV(2)=1 IXV(3)=0 CALL LEKCOF('Opérateur DETO :', & MTABX,KINC,2,IXV,IZTG2,MPOVA2,NPT2,NC2,IK2,IRET) IF(IRET.EQ.0)RETURN ctype = MPOVA2.VPOCHA(1,1) c------ Ue, vitesse extérieure IXV(1)=MELEMS IXV(2)=0 IXV(3)=0 CALL LEKCOF('Opérateur DETO :', & MTABX,KINC,3,IXV,IZTG3,MPOVA3,NPT3,NC3,IK3,IRET) IF(IRET.EQ.0)RETURN c------ dUe/dX, gradient de vitesse extérieure IXV(1)=MELEMS IXV(2)=0 IXV(3)=0 CALL LEKCOF('Opérateur DETO :', & MTABX,KINC,4,IXV,IZTG4,MPOVA4,NPT4,NC4,IK4,IRET) IF(IRET.EQ.0)RETURN c------ lecture de l'inconnue D2 au temps n-1 IXV(1)=MELEMS IXV(2)=0 IXV(3)=0 CALL LEKCOF('Opérateur DETO :', & MTABX,KINC,5,IXV,IZTG5,MPOVA5,NPT5,NC5,IK5,IRET) IF(IRET.EQ.0)RETURN c------ lecture de l'inconnue D3 au temps n-1(non obligatoire) IXV(1)=MELEMS IXV(2)=0 IXV(3)=0 CALL LEKCOF('Opérateur DETO :', & MTABX,KINC,6,IXV,IZTG6,MPOVA6,NPT6,NC6,IK6,IRET) IF(IRET.EQ.0)RETURN C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C On cree le champoint second membre NAT=2 NSOUPO=1 SEGINI MCHPOI IFOPOI = IFOUR MTYPOI = ' ' MOCHDE = 'cree par CLIM' JATTRI(1)=2 NC=1 SEGINI MSOUPO IPCHP(1)=MSOUPO NOCOMP(1)=NOMI IGEOC=MELEMS N=NPTI SEGINI MPOVAL IPOVAL=MPOVAL c write (6,*) 'longueur1 chpoint = ',VPOCHA(/1) c write (6,*) 'longueur2 chpoint = ',VPOCHA(/2) c write (6,*) 'NPTI = ',NPTI C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C On crée le champoint qui reçoit les valeurs du facteur de forme SEGINI MCHPO9 MCHPO9.JATTRI(1)=2 NC=1 SEGINI MSOUP9 MCHPO9.IPCHP(1)=MSOUP9 MSOUP9.NOCOMP(1)='SCAL' MSOUP9.IGEOC=MELEMS SEGINI MPOVA9 MSOUP9.IPOVAL=MPOVA9 C On crée le champoint qui reçoit les valeurs du coefficient de frottement SEGINI MCHPO8 MCHPO8.JATTRI(1)=2 NC=1 SEGINI MSOUP8 MCHPO8.IPCHP(1)=MSOUP8 MSOUP8.NOCOMP(1)='SCAL' MSOUP8.IGEOC=MELEMS SEGINI MPOVA8 MSOUP8.IPOVAL=MPOVA8 C On crée le champoint qui reçoit les valeurs del'épaisseur de déplacement SEGINI MCHPOA MCHPOA.JATTRI(1)=2 NC=1 SEGINI MSOUPA MCHPOA.IPCHP(1)=MSOUPA MSOUPA.NOCOMP(1)='SCAL' MSOUPA.IGEOC=MELEMS SEGINI MPOVAA MSOUPA.IPOVAL=MPOVAA C On crée le champoint qui reçoit les valeur de d(UE*delta1)/dx SEGINI MCHPOB MCHPOB.JATTRI(1)=2 NC=1 SEGINI MSOUPB MCHPOB.IPCHP(1)=MSOUPB MSOUPB.NOCOMP(1)='SCAL' MSOUPB.IGEOC=MELEMS SEGINI MPOVAB MSOUPB.IPOVAL=MPOVAB C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C MTYP='MAILLAGE' C CALL LIROBJ(MTYP,MELEME,1,IRETOU) C CALL ACTOBJ(MTYP,MELEME,1) C IF (IRETOU.NE.1)THEN C MOTERR( 1: 8) ='MAILLAGE' C CALL ERREUR(471) C RETURN C ENDIF c write(6,*)' On peut commencer a travailler ' c write(6,*)' La dimension de l espace IDIM=',idim SEGACT MELEME NBSOUS=LISOUS(/1) c write(6,*)' Nb de sous maillage NBSOUS=',NBSOUS IF (NBSOUS.NE.0) RETURN IF (NBSOUS.EQ.0)NBSOUS=1 c write(6,*)' ITYPEL=',ITYPEL,NOMS(ITYPEL) NOME=NOMS(ITYPEL)//' ' CALL KALPBG(NOME,'FONFORM ',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=GR(/1) NPG=GR(/3) IAXI=0 NBELEM=NUM(/2) NBNN =NUM(/1) write(6,*)' NBELEM=',nbelem,' NBNN=',nbnn NBEL = NBELEM NP=NBNN MP=NBNN c SEGINI MTEL SEGINI TRAV CALL INITD(UN,(NPTI*IDIM),1.D0) CALL INITD(TN,NPTI,1.D0) C creation de l'objet MATRIK NRIGE=7 NKID =9 NKMT =7 NMATRI=2 SEGINI MATRIK NBME=1 c matrice masse élémentaire SEGINI IMAT1 IRIGEL(1,1)=MELEME IRIGEL(2,1)=MELEME IRIGEL(4,1)=IMAT1 IRIGEL(7,1)=0 IMAT1.LISPRI(1)=NOMI IMAT1.LISDUA(1)=NOMI IMAT1.KSPGP=MELEMS IMAT1.KSPGD=MELEMS SEGINI IPM1 IMAT1.LIZAFM(1,1)=IPM1 c matrice de convection SEGINI IMAT2 IRIGEL(1,2)=MELEME IRIGEL(2,2)=MELEME IRIGEL(4,2)=IMAT2 IRIGEL(7,2)=2 IMAT2.LISPRI(1)=NOMI IMAT2.LISDUA(1)=NOMI IMAT2.KSPGP=MELEMS IMAT2.KSPGD=MELEMS SEGINI IPM2 IMAT2.LIZAFM(1,1)=IPM2 nd1=ipm2.am(/1) nd2=ipm2.am(/2) nd3=ipm2.am(/3) C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C>>>>>>>>>>Relations de fermeture >>>>>>>>>>>>>>>>>>>>>>>>>> C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF (ferm.EQ.1.) THEN write(6,*) 'CAS LAMINAIRE: Méthode simplifiée' IF(ctype.EQ.1.) THEN C----Cas laminaire, méthode simplifiée write(6,*) 'EQUATION DE QDM' DO I=1,VPOCHA(/1) HNM(I) = 2.5911 BNM(I) = 0.22052 MPOVA9.VPOCHA(I,1) = HNM(I) ENDDO ENDIF ENDIF IF(ferm.EQ.2.) THEN C----Cas laminaire, méthode de Karman Pohlhausen write(6,*) 'CAS LAMINAIRE: Méthode de Von Karman-Pohlhausen' IF(ctype.EQ.1.) THEN write(6,*) 'EQUATION DE QDM' CALL Karpohl (MPOVA5.VPOCHA,MPOVA4.VPOCHA,MPOVA3.VPOCHA, & VPOCHA(/1),HNM,BNM) DO I=1,VPOCHA(/1) Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) & *((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=(BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1) & *MPOVA5.VPOCHA(I,1))**(-1)) MPOVA9.VPOCHA(I,1) = HNM(I) ENDDO ENDIF ENDIF IF(ferm.EQ.3.) THEN C----Cas laminaire, méthode à deux équations write(6,*) 'CAS LAMINAIRE: Méthode à 2 équations' IF(ctype.EQ.1.) THEN C--------Equation de Qdm write(6,*) 'EQUATION DE QDM' CALL LAM2(MPOVA5.VPOCHA,MPOVA6.VPOCHA,H32NM,HNM,BNM,DNM,D1N, & VPOCHA(/1)) DO I=1,VPOCHA(/1) Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) & *((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=((BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1) & *MPOVA5.VPOCHA(I,1))**(-1))) MPOVA9.VPOCHA(I,1)=HNM(I) MPOVA8.VPOCHA(I,1)=2*Coef2(I) MPOVAA.VPOCHA(I,1)=D1N(I) ENDDO ENDIF IF(ctype.EQ.2.) THEN c--------Equation d'énegie cinétique write(6,*) 'EQUATION D''ENERGIE CINETIQUE' CALL LAM2(MPOVA6.VPOCHA,MPOVA5.VPOCHA,H32NM,HNM,BNM,DNM,D1N, & VPOCHA(/1)) DO I=1,VPOCHA(/1) Coef1(I)=3*MPOVA4.VPOCHA(I,1)* & ((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=(2*DNM(I)*1e-6*((MPOVA3.VPOCHA(I,1) & *MPOVA5.VPOCHA(I,1))**(-1))*H32NM(I)) MPOVA9.VPOCHA(I,1)=HNM(I) MPOVAA.VPOCHA(I,1)=D1N(I) ENDDO ENDIF ENDIF IF(ferm.EQ.4.) THEN c----Cas Turbulent, méthode à deux équations (Cousteix) write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Cousteix)' IF(ctype.EQ.1.) THEN c--------Equation de qdm write(6,*) 'EQUATION DE QDM' CALL TURB5(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, & HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1)) DO I=1,VPOCHA(/1) Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) & *((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=CFN(I)/2 MPOVA9.VPOCHA(I,1)=HNM(I) MPOVA8.VPOCHA(I,1)=CFN(I) MPOVAA.VPOCHA(I,1)=D1N(I) ENDDO ENDIF IF(ctype.EQ.3.) THEN c--------Equation d'entrainement write(6,*) 'EQUATION D''ENTRAINEMENT' CALL TURB5(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA, & HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1)) DO I=1,VPOCHA(/1) Coef1(I)=(MPOVA4.VPOCHA(I,1)) & *((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=CEN(I) MPOVA9.VPOCHA(I,1)=HNM(I) MPOVA8.VPOCHA(I,1)=CFN(I) MPOVAA.VPOCHA(I,1)=D1N(I) ENDDO ENDIF ENDIF IF(ferm.EQ.5.) THEN c----Cas Turbulent, méthode à deux équations relations de Head write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Head)' IF(ctype.EQ.1.) THEN c--------Equation de qdm write(6,*) 'EQUATION DE QDM' c write(6,*)MPOVA3.VPOCHA(/1),MPOVAL.VPOCHA(/1) CALL TURB6(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, & HHN,HNM,CFN,CEN,D1N,VPOCHA(/1)) DO I=1,VPOCHA(/1) Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) & *((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=CFN(I)/2 MPOVA9.VPOCHA(I,1)=HNM(I) MPOVA8.VPOCHA(I,1)=CFN(I) MPOVAA.VPOCHA(I,1)=D1N(I) ENDDO ENDIF IF(ctype.EQ.3.) THEN c--------Equation d'entrainement write(6,*) 'EQUATION D''ENTRAINEMENT' CALL TURB6(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA, & HHN,HNM,CFN,CEN,D1N,VPOCHA(/1)) DO I=1,VPOCHA(/1) Coef1(I)=(MPOVA4.VPOCHA(I,1)) & *((MPOVA3.VPOCHA(I,1))**(-1)) Coef2(I)=CEN(I) MPOVA9.VPOCHA(I,1)=HNM(I) MPOVA8.VPOCHA(I,1)=CFN(I) MPOVAA.VPOCHA(I,1)=D1N(I) ENDDO ENDIF ENDIF c IF(ferm.EQ.6) THEN C---- Convection naturelle laminaire c write(6,*) 'CONVECTION NATURELLE' c IF(ctype.EQ.1) THEN C--------Equation de quantité de mouvement c write(6,*) 'EQUATION DE QDM' c CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, c & MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1)) c DO I=1,VPOCHA(/1) c Coef1(I) = 0 c Coef2(I) = SM1(I) c MPOVA9.VPOCHA(I,1)=FP(I) c MPOVA8.VPOCHA(I,1)=PHIP(I) c MPOVAA.VPOCHA(I,1)=Delta(I) c ENDDO c ENDIF c IF(ctype.EQ.4.) THEN C--------Equation de la chaleur c write(6,*) 'EQUATION DE LA CHALEUR' c CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, c & MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1)) c DO I=1,VPOCHA(/1) c Coef1(I) = 0 c Coef2(I) = SM2(I) c MPOVA9.VPOCHA(I,1)=FP(I) c MPOVA8.VPOCHA(I,1)=PHIP(I) c MPOVAA.VPOCHA(I,1)=Delta(I) c ENDDO c ENDIF c ENDIF C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c write(6,*)' nd1,nd2,nd3=',nd1,nd2,nd3 c write(6,*)' NBELEM=',NBELEM DO 101 K=1,NBELEM c write(6,*)' Element n ',k c write(6,1001)(num(ii,k),ii=1,nbnn) C? DO 102 I=1,NBNN C? IP1=num(I,K) C? XYZ(1,I)=XCOOR((IP1-1)*(IDIM+1) +1) C? XYZ(2,I)=XCOOR((IP1-1)*(IDIM+1) +2) c write(6,*)'Coordonnees element ',K,' noeud ',ip1 c write(6,1002)(XYZ(M,I),M=1,IDIM) C? 102 CONTINUE C? CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, C? & NES,IDIM,NBNN,NPG,IAXI,AIRE,AJ) KDEB=1 KFIN=1 LRV =1 IKOMP=0 ANUK(1)=1.D0 C CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN DT=0.D0 SEGACT,MCOORD CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, & NES,NP,NPG,IAXI,XCOOR, & NUM(1,K),KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, & TN,LECT,UN,LECT,NPTI,ANUK, & WT,WS,HR,PGSQ,RPG,AJ,AIRE,UIL,DUIL, & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) SEGDES,MCOORD C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C>>>>>>>>>>Calcul des matrices élémentaires>>>>>>>>>>>>>>>>> C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c calcul des matrices masses élémentaires c write(6,*)' NBNN=',nbnn DO 103 I=1,NBNN DO 104 J=1,NBNN U=0.D0 DO 105 L=1,NPG V=0.D0 DO 1051 L1=1,NBNN J1=lect(num(L1,K)) V=V+FN(L1,L)*Coef1(J1) 1051 CONTINUE U=U+(WT(J,L)*FN(I,L)*PGSQ(L))*V 105 CONTINUE IPM1.AM(K,I,J)=U c write(6,*) 'M(',I,',',J,',',K,')=',IPM1.AM(K,I,J) 104 CONTINUE 103 CONTINUE ax=1.d0 ay=1.d0 c calcul des matrices de convection élémentaires DO 106 I=1,NBNN DO 107 J=1,NBNN c J1 = lect(num(J,K)) U=0.D0 DO 108 L=1,NPG V=0.D0 DO 1081 L1=1,NBNN J1 = lect(num(L1,K)) V = V+FN(L1,L)*1. 1081 CONTINUE U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L))*V c write(6,*) 'hr(',J,',',L,')=',HR(1,J,L) c write(6,*) 'FN(',I,',',L,',',K,')=',FN(I,L) c write(6,*) 'WT(',I,',',L,',',K,')=',WT(I,L) c U=U+(WT(I,L)*((ax*HR(1,J,L)) + (ay*HR(1,J,L)) )*PGSQ(L)) c U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L)) c &*MPOVA2.VPOCHA(J1,1) 108 CONTINUE C write(6,*)'K,i,j=', K,i,j IPM2.AM(K,I,J) = U 107 CONTINUE 106 CONTINUE C calcul du terme source DO 109 I=1,NBNN I1= lect(num(I,K)) U=0.D0 DO 110 L=1,NPG V=0.D0 v2=0.D0 DO 1101 J=1,NBNN J1 =lect(num(J,K)) V=V+FN(J,L)*Coef2(J1) 1101 CONTINUE U=U+V*WT(I,L)*PGSQ(L) 110 CONTINUE VPOCHA(I1,1)=VPOCHA(I1,1)+ U 109 CONTINUE C calcul de la vitesse de transpiration DO 111 I=1,NBNN I1 = lect(num(I,K)) U=0.D0 DO 112 J=1,NBNN J1 = lect(num(J,K)) U = U+HR(1,J,1)*(D1N(J1)*MPOVA3.VPOCHA(J1,1)) 112 CONTINUE VT(I1) = U c write(6,*) 'VT(',I1,')=',VT(I1) MPOVAB.VPOCHA(I1,1)= VT(I1) 111 CONTINUE 101 CONTINUE C Sortie des valeurs des coefficients H, Cf, Ce et VT CALL ECMO(KINC,'H','CHPOINT',MCHPO9) CALL ECMO(KINC,'CF','CHPOINT',MCHPO8) CALL ECMO(KINC,'D1','CHPOINT',MCHPOA) CALL ECMO(KINC,'Vtra','CHPOINT',MCHPOB) SEGDES IPM2,IMAT2,IPM1,IMAT1 SEGDES MATRIK SEGDES TRAV SEGDES IZFFM,IZHR SEGDES LINCO C SEGDES MTEL CALL ECROBJ('MATRIK ',MATRIK) CALL ACTOBJ('CHPOINT ',MCHPOI,1) CALL ECROBJ('CHPOINT ',MCHPOI) RETURN 1003 FORMAT(6(1X,1PE11.4)) 1002 FORMAT(10(1X,1PE11.4)) 1001 FORMAT(20(1X,I5)) END