Télécharger clmi.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales