1 : C CLMI SOURCE MAGN 10/05/18 21:15:13 6675 2 : SUBROUTINE CLMI 3 : PARAMETER (NTB=1) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : REAL*8 ctype,ferm 7 : C*********************************************************************** 8 : C 9 : C 10 : C*********************************************************************** 11 : 12 : -INC CCOPTIO 13 : -INC CCGEOME 14 : -INC SMELEME 15 : -INC SMCOORD 16 : -INC SIZFFB 17 : -INC SMMATRIK 18 : POINTEUR IMAT1.IMATRI,IMAT2.IMATRI,IMAT3.IMATRI 19 : -INC SMCHPOI 20 : POINTEUR IZTU1.MPOVAL,MCHPO9.MCHPOI,MSOUP9.MSOUPO,MPOVA9.MPOVAL 21 : POINTEUR MCHPO8.MCHPOI,MSOUP8.MSOUPO,MPOVA8.MPOVAL 22 : POINTEUR MCHPOA.MCHPOI,MSOUPA.MSOUPO,MPOVAA.MPOVAL 23 : POINTEUR MCHPOB.MCHPOI,MSOUPB.MSOUPO,MPOVAB.MPOVAL 24 : SEGMENT TRAV 25 : REAL*8 UN(NPTI,IDIM),TN(NPTI,1) 26 : REAL*8 WT(NP,NPG),WS(NP,NPG) 27 : REAL*8 UIL(IDIM,NPG),DUIL(IDIM,NPG),ANUK(1),AIRE(1) 28 : REAL*8 HNM(NPTI),BNM(NPTI),Coef1(NPTI),Coef2(NPTI),DNM(NPTI) 29 : REAL*8 H32NM(NPTI),CFN(NPTI),CEN(NPTI),DeN(NPTI),D1N(NPTI) 30 : REAL*8 HHN(NPTI) 31 : REAL*8 UCL(NPTI),Delta(NPTI),SM1(NPTI),SM2(NPTI) 32 : REAL*8 FP(NPTI),PHIP(NPTI),VT(NPTI) 33 : ENDSEGMENT 34 : 35 : -INC SMLENTI 36 : -INC SMLMOTS 37 : POINTEUR LINCO.MLMOTS 38 : real*8 nu 39 : CHARACTER*8 MTYP,NOME,NOMI,TYPE,NOMZ,TYPC 40 : CHARACTER*8 LTAB(NTB) 41 : DIMENSION KTAB(NTB),IXV(9) 42 : DATA LTAB/'KIZX '/ 43 : 44 : C***************************************************************************** 45 : C 46 : C write(6,*)' Debut CLMI ' 47 : 47 : CALL LITABS(LTAB,KTAB,NTB,1,IRET) 48 : IF (IERR.NE.0) RETURN 49 : MTABX=KTAB(1) 50 : C 51 : C- Récupération de la table EQEX (pointeur MTAB1) 52 : C 52 : CALL LEKTAB(MTABX,'EQEX',MTAB1) 53 : IF(MTAB1.EQ.0)THEN 54 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 55 : MOTERR( 1: 8) = ' EQEX ' 56 : MOTERR( 9:16) = ' EQEX ' 57 : MOTERR(17:24) = ' KIZX ' 57 : CALL ERREUR(786) 58 : RETURN 59 : ENDIF 59 : CALL ACME(MTAB1,'NAVISTOK',NASTOK) 60 : C 61 : C- Récupération de la table INCO (pointeur KINC) 62 : C 62 : CALL LEKTAB(MTAB1,'INCO',KINC) 63 : IF(KINC.EQ.0)THEN 64 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 65 : MOTERR( 1: 8) = ' INCO ' 66 : MOTERR( 9:16) = ' INCO ' 67 : MOTERR(17:24) = ' EQEX ' 67 : CALL ERREUR(786) 68 : RETURN 69 : ENDIF 70 : 71 : C***************************************************************************** 72 : C OPTIONS 73 : C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI 74 : C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC 75 : C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG 76 : 77 : IAXI=0 78 : IF(IFOMOD.EQ.0)IAXI=2 79 : C 80 : C- Récupération de la table des options KOPT (pointeur KOPTI) 81 : C 81 : CALL LEKTAB(MTABX,'KOPT',KOPTI) 82 : IF (KOPTI.EQ.0) THEN 83 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 84 : MOTERR( 1: 8) = ' KOPT ' 85 : MOTERR( 9:16) = ' KOPT ' 86 : MOTERR(17:24) = ' KIZX ' 86 : CALL ERREUR(786) 87 : RETURN 88 : ENDIF 89 : 90 : C write(6,*)' Avant les options ' 90 : CALL ACME(KOPTI,'IDCEN',IDCEN) 91 : KDIM=1 92 : IF(IDCEN.EQ.2)KDIM=IDIM 92 : CALL ACME(KOPTI,'IKOMP',IKOMP) 92 : CALL ACME(KOPTI,'KIMPL',KIMPL) 92 : CALL ACME(KOPTI,'KPOIN',KPRE) 92 : CALL ACME(KOPTI,'KFORM',KFORM) 92 : CALL ACMF(KOPTI,'CMD',CMD) 93 : 94 : IF(KFORM.NE.0.AND.KFORM.NE.1)THEN 95 : C Option %m1:8 incompatible avec les données 96 : MOTERR( 1: 8) = 'EF/EFM1 ' 96 : CALL ERREUR(803) 97 : RETURN 98 : ENDIF 98 : CALL ACMF(KOPTI,'AIMPL',AIMPL) 99 : AG=AIMPL 100 : AD=AIMPL-1.D0 101 : IF (IERR.NE.0) RETURN 102 : 103 : C write(6,*)' Apres les options ' 104 : C***************************************************************************** 105 : C 106 : C- Récupération de la table DOMAINE associée au domaine local 107 : C 108 : c write(6,*)' Recuperation du modele ' 108 : CALL ACMM(MTABX,'NOMZONE',NOMZ) 109 : TYPE=' ' 109 : CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ) 110 : IF(TYPE.NE.'MMODEL')THEN 111 : C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 112 : MOTERR( 1: 8) = ' DOMZ ' 113 : MOTERR( 9:16) = ' DOMZ ' 114 : MOTERR(17:24) = ' KIZX ' 114 : CALL ERREUR(786) 115 : RETURN 116 : ENDIF 117 : 118 : C E/ MMODEL : Pointeur de la table contenant l'information cherchée 119 : C /S IPOINT : Pointeur sur la table DOMAINE 120 : C /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE 121 : C INEFMD=4 LINB 122 : 122 : CALL LEKMOD(MMDZ,MTABZ,INEFMD) 123 : C verifier que la formulation est bonne INEFMD = 1 124 : 124 : CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) 124 : CALL LEKTAB(MTABZ,'SOMMET',MELEMS) 124 : C CALL LEKTAB(MTABZ,'CENTRE',MELEMC) 124 : CALL KRIPAD(MELEMS,MLENTI) 125 : IF (IERR.NE.0) RETURN 126 : 127 : 128 : C************************************************************************* 129 : C VERIFICATIONS SUR LES INCONNUES 130 : C 131 : C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI 132 : C 133 : c write(6,*)' Recuperation inconnues ' 134 : TYPE='LISTMOTS' 134 : CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) 135 : IF (IERR.NE.0) RETURN 136 : SEGACT LINCO 137 : NBINC=LINCO.MOTS(/2) 138 : IF(NBINC.NE.1)THEN 139 : C Indice %m1:8 : contient plus de %i1 %m9:16 140 : MOTERR( 1:8) = 'LISTINCO' 141 : INTERR(1) = 1 142 : MOTERR(9:16) = ' MOTS ' 142 : CALL ERREUR(799) 143 : RETURN 144 : ENDIF 145 : 146 : NOMI=LINCO.MOTS(1) 147 : c write(6,*) 'Nb inco=',NBINC 148 : c write(6,*) 'NOMI=',NOMI 149 : C 150 : C- Récupération de l'inconnue 151 : c write(6,*) 'Récupération de l''inconnue' 152 : TYPE=' ' 152 : CALL ACMO(KINC,NOMI,TYPE,MCHPO1) 153 : IF(TYPE.NE.'CHPOINT ')THEN 154 : C Indice %m1:8 : ne contient pas un objet de type %m9:16 155 : MOTERR( 1: 8) = 'INC '//NOMI 156 : MOTERR( 9:16) = 'CHPOINT ' 156 : CALL ERREUR(800) 157 : RETURN 158 : ELSE 158 : CALL LICHT(MCHPO1,IZTU1,TYPC,MELEMI) 159 : NINKO = IZTU1.VPOCHA(/2) 160 : NPTI = IZTU1.VPOCHA(/1) 161 : c write (6,*) 'NINKO =',NINKO 162 : c write(6,*)' MCHPOI,MELEMI=',MCHPOI,MELEMI 163 : C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour 164 : C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes 165 : ENDIF 166 : 167 : C***************************************************************************** 168 : C Le domaine de definition est donne par le SPG de la premiere inconnue 169 : C Les inconnues suivantes devront posseder ce meme pointeur 170 : C On verifie que les points de la zone sont tous inclus dans ce SPG 171 : 171 : CALL KRIPAD(MELEMI,IPADI) 172 : 172 : CALL VERPAD(IPADI,MELEME,IRET) 173 : IF(IRET.NE.0)THEN 174 : C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique 175 : MOTERR(1: 8) = 'INC '//NOMI 176 : MOTERR(9:16) = 'CHPOINT ' 176 : CALL ERREUR(788) 177 : RETURN 178 : ENDIF 179 : 180 : C***************************************************************************** 181 : 182 : NU=1.d-6 183 : 184 : C************************************************************************* 185 : C Lecture des coefficients 186 : C Type du coefficient : 187 : C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur 188 : 189 : c write(6,*)' lecture coef ' 189 : CALL ACME(MTABX,'IARG',IARG) 190 : IF(IARG.NE.6)THEN 191 : WRITE(6,*)' Operateur DETO : ' 192 : WRITE(6,*)' Nombre d''arguments incorrect : ',IARG 193 : WRITE(6,*)' On attend 6 ' 194 : C Indice %m1:8 : nombre d'arguments incorrect 195 : MOTERR(1:8) = 'IARG ' 195 : CALL ERREUR(804) 196 : RETURN 197 : ENDIF 198 : 199 : 200 : C--Cas incompréssible 201 : c-------1er coefficient/type de fermeture 202 : IXV(1)=0 203 : IXV(2)=1 204 : IXV(3)=0 204 : CALL LEKCOF('Opérateur DETO :', 205 : & MTABX,KINC,1,IXV,IZTG1,MPOVA1,NPT1,NC1,IK1,IRET) 206 : IF(IRET.EQ.0)RETURN 207 : ferm = MPOVA1.VPOCHA(1,1) 208 : c------ 2eme coef/n°de l'équation traitée 209 : IXV(1)=0 210 : IXV(2)=1 211 : IXV(3)=0 211 : CALL LEKCOF('Opérateur DETO :', 212 : & MTABX,KINC,2,IXV,IZTG2,MPOVA2,NPT2,NC2,IK2,IRET) 213 : IF(IRET.EQ.0)RETURN 214 : ctype = MPOVA2.VPOCHA(1,1) 215 : c------ Ue, vitesse extérieure 216 : IXV(1)=MELEMS 217 : IXV(2)=0 218 : IXV(3)=0 218 : CALL LEKCOF('Opérateur DETO :', 219 : & MTABX,KINC,3,IXV,IZTG3,MPOVA3,NPT3,NC3,IK3,IRET) 220 : IF(IRET.EQ.0)RETURN 221 : c------ dUe/dX, gradient de vitesse extérieure 222 : IXV(1)=MELEMS 223 : IXV(2)=0 224 : IXV(3)=0 224 : CALL LEKCOF('Opérateur DETO :', 225 : & MTABX,KINC,4,IXV,IZTG4,MPOVA4,NPT4,NC4,IK4,IRET) 226 : IF(IRET.EQ.0)RETURN 227 : c------ lecture de l'inconnue D2 au temps n-1 228 : IXV(1)=MELEMS 229 : IXV(2)=0 230 : IXV(3)=0 230 : CALL LEKCOF('Opérateur DETO :', 231 : & MTABX,KINC,5,IXV,IZTG5,MPOVA5,NPT5,NC5,IK5,IRET) 232 : IF(IRET.EQ.0)RETURN 233 : c------ lecture de l'inconnue D3 au temps n-1(non obligatoire) 234 : IXV(1)=MELEMS 235 : IXV(2)=0 236 : IXV(3)=0 236 : CALL LEKCOF('Opérateur DETO :', 237 : & MTABX,KINC,6,IXV,IZTG6,MPOVA6,NPT6,NC6,IK6,IRET) 238 : IF(IRET.EQ.0)RETURN 239 : 240 : 241 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 242 : C On cree le champoint second membre 243 : NAT=2 244 : NSOUPO=1 245 : SEGINI MCHPOI 246 : JATTRI(1)=2 247 : NC=1 248 : SEGINI MSOUPO 249 : IPCHP(1)=MSOUPO 250 : NOCOMP(1)=NOMI(1:4) 251 : NOCONS(1)='QDM ' 252 : IGEOC=MELEMS 253 : N=NPTI 254 : SEGINI MPOVAL 255 : IPOVAL=MPOVAL 256 : 257 : c write (6,*) 'longueur1 chpoint = ',VPOCHA(/1) 258 : c write (6,*) 'longueur2 chpoint = ',VPOCHA(/2) 259 : c write (6,*) 'NPTI = ',NPTI 260 : 261 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 262 : C On crée le champoint qui reçoit les valeurs du facteur de forme 263 : SEGINI MCHPO9 264 : MCHPO9.JATTRI(1)=2 265 : NC=1 266 : SEGINI MSOUP9 267 : MCHPO9.IPCHP(1)=MSOUP9 268 : MSOUP9.NOCOMP(1)='SCAL' 269 : MSOUP9.NOCONS(1)='H ' 270 : MSOUP9.IGEOC=MELEMS 271 : SEGINI MPOVA9 272 : MSOUP9.IPOVAL=MPOVA9 273 : C On crée le champoint qui reçoit les valeurs du coefficient de frottement 274 : SEGINI MCHPO8 275 : MCHPO8.JATTRI(1)=2 276 : NC=1 277 : SEGINI MSOUP8 278 : MCHPO8.IPCHP(1)=MSOUP8 279 : MSOUP8.NOCOMP(1)='SCAL' 280 : MSOUP8.NOCONS(1)='Cf ' 281 : MSOUP8.IGEOC=MELEMS 282 : SEGINI MPOVA8 283 : MSOUP8.IPOVAL=MPOVA8 284 : C On crée le champoint qui reçoit les valeurs del'épaisseur de déplacement 285 : SEGINI MCHPOA 286 : MCHPOA.JATTRI(1)=2 287 : NC=1 288 : SEGINI MSOUPA 289 : MCHPOA.IPCHP(1)=MSOUPA 290 : MSOUPA.NOCOMP(1)='SCAL' 291 : MSOUPA.NOCONS(1)='D1 ' 292 : MSOUPA.IGEOC=MELEMS 293 : SEGINI MPOVAA 294 : MSOUPA.IPOVAL=MPOVAA 295 : C On crée le champoint qui reçoit les valeur de d(UE*delta1)/dx 296 : SEGINI MCHPOB 297 : MCHPOB.JATTRI(1)=2 298 : NC=1 299 : SEGINI MSOUPB 300 : MCHPOB.IPCHP(1)=MSOUPB 301 : MSOUPB.NOCOMP(1)='SCAL' 302 : MSOUPB.NOCONS(1)='Vtra' 303 : MSOUPB.IGEOC=MELEMS 304 : SEGINI MPOVAB 305 : MSOUPB.IPOVAL=MPOVAB 306 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 307 : 308 : C MTYP='MAILLAGE' 308 : C CALL LIROBJ(MTYP,MELEME,1,IRETOU) 309 : C IF (IRETOU.NE.1)THEN 310 : C MOTERR( 1: 8) ='MAILLAGE' 310 : C CALL ERREUR(471) 311 : C RETURN 312 : C ENDIF 313 : 314 : c write(6,*)' On peut commencer a travailler ' 315 : c write(6,*)' La dimension de l espace IDIM=',idim 316 : 317 : SEGACT MELEME 318 : NBSOUS=LISOUS(/1) 319 : c write(6,*)' Nb de sous maillage NBSOUS=',NBSOUS 320 : IF (NBSOUS.NE.0)RETURN 321 : IF (NBSOUS.EQ.0)NBSOUS=1 322 : 323 : c write(6,*)' ITYPEL=',ITYPEL,NOMS(ITYPEL) 324 : NOME=NOMS(ITYPEL)//' ' 324 : CALL KALPBG(NOME,'FONFORM ',IZFFM) 325 : 326 : SEGACT IZFFM*MOD 327 : IZHR=KZHR(1) 328 : SEGACT IZHR*MOD 329 : 330 : NES=GR(/1) 331 : NPG=GR(/3) 332 : IAXI=0 333 : 334 : NBELEM=NUM(/2) 335 : NBNN =NUM(/1) 336 : write(6,*)' NBELEM=',nbelem,' NBNN=',nbnn 337 : NBEL = NBELEM 338 : NP=NBNN 339 : MP=NBNN 340 : 341 : c SEGINI MTEL 342 : SEGINI TRAV 342 : CALL INITD(UN,(NPTI*IDIM),1.D0) 342 : CALL INITD(TN,NPTI,1.D0) 343 : 344 : C creation de l'objet MATRIK 345 : 346 : NRIGE=7 347 : NKID =9 348 : NKMT =7 349 : NMATRI=2 350 : 351 : SEGINI MATRIK 351 : CALL ECROBJ('MATRIK',MATRIK) 351 : CALL ECROBJ('CHPOINT',MCHPOI) 352 : 353 : NBME=1 354 : 355 : c matrice masse élémentaire 356 : SEGINI IMAT1 357 : IRIGEL(1,1)=MELEME 358 : IRIGEL(2,1)=MELEME 359 : IRIGEL(4,1)=IMAT1 360 : IRIGEL(7,1)=0 361 : 362 : IMAT1.LISPRI(1)=NOMI 363 : IMAT1.LISDUA(1)=NOMI 364 : IMAT1.KSPGP=MELEMS 365 : IMAT1.KSPGD=MELEMS 366 : SEGINI IPM1 367 : IMAT1.LIZAFM(1,1)=IPM1 368 : c matrice de convection 369 : SEGINI IMAT2 370 : IRIGEL(1,2)=MELEME 371 : IRIGEL(2,2)=MELEME 372 : IRIGEL(4,2)=IMAT2 373 : IRIGEL(7,2)=2 374 : 375 : IMAT2.LISPRI(1)=NOMI 376 : IMAT2.LISDUA(1)=NOMI 377 : IMAT2.KSPGP=MELEMS 378 : IMAT2.KSPGD=MELEMS 379 : SEGINI IPM2 380 : IMAT2.LIZAFM(1,1)=IPM2 381 : nd1=ipm2.am(/1) 382 : nd2=ipm2.am(/2) 383 : nd3=ipm2.am(/3) 384 : 385 : 386 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 387 : C>>>>>>>>>>Relations de fermeture >>>>>>>>>>>>>>>>>>>>>>>>>> 388 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 389 : IF (ferm.EQ.1.) THEN 390 : write(6,*) 'CAS LAMINAIRE: Méthode simplifiée' 391 : IF(ctype.EQ.1.) THEN 392 : C----Cas laminaire, méthode simplifiée 393 : write(6,*) 'EQUATION DE QDM' 394 : DO I=1,VPOCHA(/1) 395 : HNM(I) = 2.5911 396 : BNM(I) = 0.22052 397 : MPOVA9.VPOCHA(I,1) = HNM(I) 398 : ENDDO 399 : ENDIF 400 : ENDIF 401 : 402 : IF(ferm.EQ.2.) THEN 403 : C----Cas laminaire, méthode de Karman Pohlhausen 404 : write(6,*) 'CAS LAMINAIRE: Méthode de Von Karman-Pohlhausen' 405 : IF(ctype.EQ.1.) THEN 406 : write(6,*) 'EQUATION DE QDM' 406 : CALL Karpohl (MPOVA5.VPOCHA,MPOVA4.VPOCHA,MPOVA3.VPOCHA, 407 : & VPOCHA(/1),HNM,BNM) 408 : DO I=1,VPOCHA(/1) 409 : Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) 410 : & *((MPOVA3.VPOCHA(I,1))**(-1)) 411 : Coef2(I)=(BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1) 412 : & *MPOVA5.VPOCHA(I,1))**(-1)) 413 : MPOVA9.VPOCHA(I,1) = HNM(I) 414 : 415 : ENDDO 416 : ENDIF 417 : ENDIF 418 : 419 : IF(ferm.EQ.3.) THEN 420 : C----Cas laminaire, méthode à deux équations 421 : write(6,*) 'CAS LAMINAIRE: Méthode à 2 équations' 422 : IF(ctype.EQ.1.) THEN 423 : C--------Equation de Qdm 424 : write(6,*) 'EQUATION DE QDM' 424 : CALL LAM2(MPOVA5.VPOCHA,MPOVA6.VPOCHA,H32NM,HNM,BNM,DNM,D1N, 425 : & VPOCHA(/1)) 426 : DO I=1,VPOCHA(/1) 427 : Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) 428 : & *((MPOVA3.VPOCHA(I,1))**(-1)) 429 : Coef2(I)=((BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1) 430 : & *MPOVA5.VPOCHA(I,1))**(-1))) 431 : MPOVA9.VPOCHA(I,1)=HNM(I) 432 : MPOVA8.VPOCHA(I,1)=2*Coef2(I) 433 : MPOVAA.VPOCHA(I,1)=D1N(I) 434 : ENDDO 435 : ENDIF 436 : IF(ctype.EQ.2.) THEN 437 : c--------Equation d'énegie cinétique 438 : write(6,*) 'EQUATION D''ENERGIE CINETIQUE' 438 : CALL LAM2(MPOVA6.VPOCHA,MPOVA5.VPOCHA,H32NM,HNM,BNM,DNM,D1N, 439 : & VPOCHA(/1)) 440 : DO I=1,VPOCHA(/1) 441 : Coef1(I)=3*MPOVA4.VPOCHA(I,1)* 442 : & ((MPOVA3.VPOCHA(I,1))**(-1)) 443 : Coef2(I)=(2*DNM(I)*1e-6*((MPOVA3.VPOCHA(I,1) 444 : & *MPOVA5.VPOCHA(I,1))**(-1))*H32NM(I)) 445 : MPOVA9.VPOCHA(I,1)=HNM(I) 446 : MPOVAA.VPOCHA(I,1)=D1N(I) 447 : ENDDO 448 : ENDIF 449 : ENDIF 450 : 451 : IF(ferm.EQ.4.) THEN 452 : c----Cas Turbulent, méthode à deux équations (Cousteix) 453 : write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Cousteix)' 454 : IF(ctype.EQ.1.) THEN 455 : c--------Equation de qdm 456 : write(6,*) 'EQUATION DE QDM' 456 : CALL TURB5(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, 457 : & HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1)) 458 : DO I=1,VPOCHA(/1) 459 : Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) 460 : & *((MPOVA3.VPOCHA(I,1))**(-1)) 461 : Coef2(I)=CFN(I)/2 462 : MPOVA9.VPOCHA(I,1)=HNM(I) 463 : MPOVA8.VPOCHA(I,1)=CFN(I) 464 : MPOVAA.VPOCHA(I,1)=D1N(I) 465 : ENDDO 466 : ENDIF 467 : IF(ctype.EQ.3.) THEN 468 : c--------Equation d'entrainement 469 : write(6,*) 'EQUATION D''ENTRAINEMENT' 469 : CALL TURB5(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA, 470 : & HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1)) 471 : DO I=1,VPOCHA(/1) 472 : Coef1(I)=(MPOVA4.VPOCHA(I,1)) 473 : & *((MPOVA3.VPOCHA(I,1))**(-1)) 474 : Coef2(I)=CEN(I) 475 : MPOVA9.VPOCHA(I,1)=HNM(I) 476 : MPOVA8.VPOCHA(I,1)=CFN(I) 477 : MPOVAA.VPOCHA(I,1)=D1N(I) 478 : ENDDO 479 : ENDIF 480 : ENDIF 481 : IF(ferm.EQ.5.) THEN 482 : c----Cas Turbulent, méthode à deux équations relations de Head 483 : write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Head)' 484 : IF(ctype.EQ.1.) THEN 485 : c--------Equation de qdm 486 : write(6,*) 'EQUATION DE QDM' 487 : c write(6,*)MPOVA3.VPOCHA(/1),MPOVAL.VPOCHA(/1) 487 : CALL TURB6(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, 488 : & HHN,HNM,CFN,CEN,D1N,VPOCHA(/1)) 489 : DO I=1,VPOCHA(/1) 490 : Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1) 491 : & *((MPOVA3.VPOCHA(I,1))**(-1)) 492 : Coef2(I)=CFN(I)/2 493 : MPOVA9.VPOCHA(I,1)=HNM(I) 494 : MPOVA8.VPOCHA(I,1)=CFN(I) 495 : MPOVAA.VPOCHA(I,1)=D1N(I) 496 : ENDDO 497 : ENDIF 498 : IF(ctype.EQ.3.) THEN 499 : c--------Equation d'entrainement 500 : write(6,*) 'EQUATION D''ENTRAINEMENT' 500 : CALL TURB6(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA, 501 : & HHN,HNM,CFN,CEN,D1N,VPOCHA(/1)) 502 : DO I=1,VPOCHA(/1) 503 : Coef1(I)=(MPOVA4.VPOCHA(I,1)) 504 : & *((MPOVA3.VPOCHA(I,1))**(-1)) 505 : Coef2(I)=CEN(I) 506 : MPOVA9.VPOCHA(I,1)=HNM(I) 507 : MPOVA8.VPOCHA(I,1)=CFN(I) 508 : MPOVAA.VPOCHA(I,1)=D1N(I) 509 : ENDDO 510 : ENDIF 511 : ENDIF 512 : c IF(ferm.EQ.6) THEN 513 : C---- Convection naturelle laminaire 514 : c write(6,*) 'CONVECTION NATURELLE' 515 : c IF(ctype.EQ.1) THEN 516 : C--------Equation de quantité de mouvement 517 : c write(6,*) 'EQUATION DE QDM' 517 : c CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, 518 : c & MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1)) 519 : c DO I=1,VPOCHA(/1) 520 : c Coef1(I) = 0 521 : c Coef2(I) = SM1(I) 522 : c MPOVA9.VPOCHA(I,1)=FP(I) 523 : c MPOVA8.VPOCHA(I,1)=PHIP(I) 524 : c MPOVAA.VPOCHA(I,1)=Delta(I) 525 : c ENDDO 526 : c ENDIF 527 : c IF(ctype.EQ.4.) THEN 528 : C--------Equation de la chaleur 529 : c write(6,*) 'EQUATION DE LA CHALEUR' 529 : c CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA, 530 : c & MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1)) 531 : c DO I=1,VPOCHA(/1) 532 : c Coef1(I) = 0 533 : c Coef2(I) = SM2(I) 534 : c MPOVA9.VPOCHA(I,1)=FP(I) 535 : c MPOVA8.VPOCHA(I,1)=PHIP(I) 536 : c MPOVAA.VPOCHA(I,1)=Delta(I) 537 : c ENDDO 538 : c ENDIF 539 : c ENDIF 540 : 541 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 542 : c write(6,*)' nd1,nd2,nd3=',nd1,nd2,nd3 543 : c write(6,*)' NBELEM=',NBELEM 544 : 545 : 546 : DO 101 K=1,NBELEM 547 : c write(6,*)' Element n ',k 548 : c write(6,1001)(num(ii,k),ii=1,nbnn) 549 : 550 : C? DO 102 I=1,NBNN 551 : C? IP1=num(I,K) 552 : C? XYZ(1,I)=XCOOR((IP1-1)*(IDIM+1) +1) 553 : C? XYZ(2,I)=XCOOR((IP1-1)*(IDIM+1) +2) 554 : c write(6,*)'Coordonnees element ',K,' noeud ',ip1 555 : c write(6,1002)(XYZ(M,I),M=1,IDIM) 556 : C? 102 CONTINUE 557 : 557 : C? CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG, 558 : C? & NES,IDIM,NBNN,NPG,IAXI,AIRE,AJ) 559 : 560 : KDEB=1 561 : KFIN=1 562 : LRV =1 563 : IKOMP=0 564 : ANUK(1)=1.D0 565 : 565 : CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, 566 : & NES,NP,NPG,IAXI,XCOOR, 567 : & NUM(1,K),KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, 568 : & TN,LECT,UN,LECT,NPTI,ANUK, 569 : & WT,WS,HR,PGSQ,RPG,AJ,AIRE,UIL,DUIL, 570 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 571 : 572 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 573 : C>>>>>>>>>>Calcul des matrices élémentaires>>>>>>>>>>>>>>>>> 574 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 575 : c calcul des matrices masses élémentaires 576 : c write(6,*)' NBNN=',nbnn 577 : DO 103 I=1,NBNN 578 : DO 104 J=1,NBNN 579 : U=0.D0 580 : DO 105 L=1,NPG 581 : V=0.D0 582 : DO 1051 L1=1,NBNN 583 : J1=lect(num(L1,K)) 584 : V=V+FN(L1,L)*Coef1(J1) 585 : 1051 CONTINUE 586 : U=U+(WT(J,L)*FN(I,L)*PGSQ(L))*V 587 : 105 CONTINUE 588 : IPM1.AM(K,I,J)=U 589 : c write(6,*) 'M(',I,',',J,',',K,')=',IPM1.AM(K,I,J) 590 : 104 CONTINUE 591 : 103 CONTINUE 592 : 593 : ax=1.d0 594 : ay=1.d0 595 : c calcul des matrices de convection élémentaires 596 : DO 106 I=1,NBNN 597 : DO 107 J=1,NBNN 598 : c J1 = lect(num(J,K)) 599 : U=0.D0 600 : DO 108 L=1,NPG 601 : V=0.D0 602 : DO 1081 L1=1,NBNN 603 : J1 = lect(num(L1,K)) 604 : V = V+FN(L1,L)*1. 605 : 1081 CONTINUE 606 : U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L))*V 607 : c write(6,*) 'hr(',J,',',L,')=',HR(1,J,L) 608 : c write(6,*) 'FN(',I,',',L,',',K,')=',FN(I,L) 609 : c write(6,*) 'WT(',I,',',L,',',K,')=',WT(I,L) 610 : c U=U+(WT(I,L)*((ax*HR(1,J,L)) + (ay*HR(1,J,L)) )*PGSQ(L)) 611 : c U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L)) 612 : c &*MPOVA2.VPOCHA(J1,1) 613 : 108 CONTINUE 614 : C write(6,*)'K,i,j=', K,i,j 615 : IPM2.AM(K,I,J) = U 616 : 107 CONTINUE 617 : 106 CONTINUE 618 : 619 : C calcul du terme source 620 : DO 109 I=1,NBNN 621 : I1= lect(num(I,K)) 622 : U=0.D0 623 : DO 110 L=1,NPG 624 : V=0.D0 625 : v2=0.D0 626 : DO 1101 J=1,NBNN 627 : J1 =lect(num(J,K)) 628 : V=V+FN(J,L)*Coef2(J1) 629 : 1101 CONTINUE 630 : U=U+V*WT(I,L)*PGSQ(L) 631 : 110 CONTINUE 632 : VPOCHA(I1,1)=VPOCHA(I1,1)+ U 633 : 109 CONTINUE 634 : 635 : C calcul de la vitesse de transpiration 636 : DO 111 I=1,NBNN 637 : I1 = lect(num(I,K)) 638 : U=0.D0 639 : DO 112 J=1,NBNN 640 : J1 = lect(num(J,K)) 641 : U = U+HR(1,J,1)*(D1N(J1)*MPOVA3.VPOCHA(J1,1)) 642 : 112 CONTINUE 643 : VT(I1) = U 644 : c write(6,*) 'VT(',I1,')=',VT(I1) 645 : MPOVAB.VPOCHA(I1,1)= VT(I1) 646 : 111 CONTINUE 647 : 648 : 101 CONTINUE 649 : 650 : 651 : C Sortie des valeurs des coefficients H, Cf, Ce et VT 652 : 652 : CALL ECMO(KINC,'H','CHPOINT',MCHPO9) 652 : CALL ECMO(KINC,'CF','CHPOINT',MCHPO8) 652 : CALL ECMO(KINC,'D1','CHPOINT',MCHPOA) 652 : CALL ECMO(KINC,'Vtra','CHPOINT',MCHPOB) 653 : 654 : 655 : SEGDES IPM2,IMAT2,IPM1,IMAT1 656 : SEGDES MATRIK 657 : SEGDES TRAV 658 : SEGDES IZFFM,IZHR 659 : SEGDES MELEME 660 : SEGDES MPOVAL,MPOVA8,MPOVA9,MPOVAA,MPOVAB 661 : SEGDES MSOUPO,MSOUP8,MSOUP9,MSOUPA,MSOUPB 662 : SEGDES MCHPOI,MCHPO9,MCHPO8,MCHPOA,MCHPOB 663 : SEGDES LINCO 664 : 665 : C SEGDES MTEL 666 : 667 : RETURN 668 : 1003 FORMAT(6(1X,1PE11.4)) 669 : 1002 FORMAT(10(1X,1PE11.4)) 670 : 1001 FORMAT(20(1X,I5)) 671 : END 672 : 673 : 674 : 675 : 676 : 677 : 678 : 679 : 680 : 681 :
© Cast3M 2003 - All rights reserved.
Disclaimer