Télécharger clmi.eso

Retour à la liste

Numérotation des lignes :

  1. C CLMI SOURCE CB215821 19/08/01 21:15:21 10279
  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.  
  397. NBME=1
  398.  
  399. c matrice masse élémentaire
  400. SEGINI IMAT1
  401. IRIGEL(1,1)=MELEME
  402. IRIGEL(2,1)=MELEME
  403. IRIGEL(4,1)=IMAT1
  404. IRIGEL(7,1)=0
  405.  
  406. IMAT1.LISPRI(1)=NOMI
  407. IMAT1.LISDUA(1)=NOMI
  408. IMAT1.KSPGP=MELEMS
  409. IMAT1.KSPGD=MELEMS
  410. SEGINI IPM1
  411. IMAT1.LIZAFM(1,1)=IPM1
  412. c matrice de convection
  413. SEGINI IMAT2
  414. IRIGEL(1,2)=MELEME
  415. IRIGEL(2,2)=MELEME
  416. IRIGEL(4,2)=IMAT2
  417. IRIGEL(7,2)=2
  418.  
  419. IMAT2.LISPRI(1)=NOMI
  420. IMAT2.LISDUA(1)=NOMI
  421. IMAT2.KSPGP=MELEMS
  422. IMAT2.KSPGD=MELEMS
  423. SEGINI IPM2
  424. IMAT2.LIZAFM(1,1)=IPM2
  425. nd1=ipm2.am(/1)
  426. nd2=ipm2.am(/2)
  427. nd3=ipm2.am(/3)
  428.  
  429.  
  430. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  431. C>>>>>>>>>>Relations de fermeture >>>>>>>>>>>>>>>>>>>>>>>>>>
  432. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  433. IF (ferm.EQ.1.) THEN
  434. write(6,*) 'CAS LAMINAIRE: Méthode simplifiée'
  435. IF(ctype.EQ.1.) THEN
  436. C----Cas laminaire, méthode simplifiée
  437. write(6,*) 'EQUATION DE QDM'
  438. DO I=1,VPOCHA(/1)
  439. HNM(I) = 2.5911
  440. BNM(I) = 0.22052
  441. MPOVA9.VPOCHA(I,1) = HNM(I)
  442. ENDDO
  443. ENDIF
  444. ENDIF
  445.  
  446. IF(ferm.EQ.2.) THEN
  447. C----Cas laminaire, méthode de Karman Pohlhausen
  448. write(6,*) 'CAS LAMINAIRE: Méthode de Von Karman-Pohlhausen'
  449. IF(ctype.EQ.1.) THEN
  450. write(6,*) 'EQUATION DE QDM'
  451. CALL Karpohl (MPOVA5.VPOCHA,MPOVA4.VPOCHA,MPOVA3.VPOCHA,
  452. & VPOCHA(/1),HNM,BNM)
  453. DO I=1,VPOCHA(/1)
  454. Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
  455. & *((MPOVA3.VPOCHA(I,1))**(-1))
  456. Coef2(I)=(BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1)
  457. & *MPOVA5.VPOCHA(I,1))**(-1))
  458. MPOVA9.VPOCHA(I,1) = HNM(I)
  459.  
  460. ENDDO
  461. ENDIF
  462. ENDIF
  463.  
  464. IF(ferm.EQ.3.) THEN
  465. C----Cas laminaire, méthode à deux équations
  466. write(6,*) 'CAS LAMINAIRE: Méthode à 2 équations'
  467. IF(ctype.EQ.1.) THEN
  468. C--------Equation de Qdm
  469. write(6,*) 'EQUATION DE QDM'
  470. CALL LAM2(MPOVA5.VPOCHA,MPOVA6.VPOCHA,H32NM,HNM,BNM,DNM,D1N,
  471. & VPOCHA(/1))
  472. DO I=1,VPOCHA(/1)
  473. Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
  474. & *((MPOVA3.VPOCHA(I,1))**(-1))
  475. Coef2(I)=((BNM(I)*1e-6)*((MPOVA3.VPOCHA(I,1)
  476. & *MPOVA5.VPOCHA(I,1))**(-1)))
  477. MPOVA9.VPOCHA(I,1)=HNM(I)
  478. MPOVA8.VPOCHA(I,1)=2*Coef2(I)
  479. MPOVAA.VPOCHA(I,1)=D1N(I)
  480. ENDDO
  481. ENDIF
  482. IF(ctype.EQ.2.) THEN
  483. c--------Equation d'énegie cinétique
  484. write(6,*) 'EQUATION D''ENERGIE CINETIQUE'
  485. CALL LAM2(MPOVA6.VPOCHA,MPOVA5.VPOCHA,H32NM,HNM,BNM,DNM,D1N,
  486. & VPOCHA(/1))
  487. DO I=1,VPOCHA(/1)
  488. Coef1(I)=3*MPOVA4.VPOCHA(I,1)*
  489. & ((MPOVA3.VPOCHA(I,1))**(-1))
  490. Coef2(I)=(2*DNM(I)*1e-6*((MPOVA3.VPOCHA(I,1)
  491. & *MPOVA5.VPOCHA(I,1))**(-1))*H32NM(I))
  492. MPOVA9.VPOCHA(I,1)=HNM(I)
  493. MPOVAA.VPOCHA(I,1)=D1N(I)
  494. ENDDO
  495. ENDIF
  496. ENDIF
  497.  
  498. IF(ferm.EQ.4.) THEN
  499. c----Cas Turbulent, méthode à deux équations (Cousteix)
  500. write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Cousteix)'
  501. IF(ctype.EQ.1.) THEN
  502. c--------Equation de qdm
  503. write(6,*) 'EQUATION DE QDM'
  504. CALL TURB5(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
  505. & HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1))
  506. DO I=1,VPOCHA(/1)
  507. Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
  508. & *((MPOVA3.VPOCHA(I,1))**(-1))
  509. Coef2(I)=CFN(I)/2
  510. MPOVA9.VPOCHA(I,1)=HNM(I)
  511. MPOVA8.VPOCHA(I,1)=CFN(I)
  512. MPOVAA.VPOCHA(I,1)=D1N(I)
  513. ENDDO
  514. ENDIF
  515. IF(ctype.EQ.3.) THEN
  516. c--------Equation d'entrainement
  517. write(6,*) 'EQUATION D''ENTRAINEMENT'
  518. CALL TURB5(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA,
  519. & HHN,HNM,CFN,CEN,D1N,DeN,VPOCHA(/1))
  520. DO I=1,VPOCHA(/1)
  521. Coef1(I)=(MPOVA4.VPOCHA(I,1))
  522. & *((MPOVA3.VPOCHA(I,1))**(-1))
  523. Coef2(I)=CEN(I)
  524. MPOVA9.VPOCHA(I,1)=HNM(I)
  525. MPOVA8.VPOCHA(I,1)=CFN(I)
  526. MPOVAA.VPOCHA(I,1)=D1N(I)
  527. ENDDO
  528. ENDIF
  529. ENDIF
  530. IF(ferm.EQ.5.) THEN
  531. c----Cas Turbulent, méthode à deux équations relations de Head
  532. write(6,*) 'CAS TURBULENT: Méthode à 2 équations(Head)'
  533. IF(ctype.EQ.1.) THEN
  534. c--------Equation de qdm
  535. write(6,*) 'EQUATION DE QDM'
  536. c write(6,*)MPOVA3.VPOCHA(/1),MPOVAL.VPOCHA(/1)
  537. CALL TURB6(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
  538. & HHN,HNM,CFN,CEN,D1N,VPOCHA(/1))
  539. DO I=1,VPOCHA(/1)
  540. Coef1(I)= (HNM(I)+2)*MPOVA4.VPOCHA(I,1)
  541. & *((MPOVA3.VPOCHA(I,1))**(-1))
  542. Coef2(I)=CFN(I)/2
  543. MPOVA9.VPOCHA(I,1)=HNM(I)
  544. MPOVA8.VPOCHA(I,1)=CFN(I)
  545. MPOVAA.VPOCHA(I,1)=D1N(I)
  546. ENDDO
  547. ENDIF
  548. IF(ctype.EQ.3.) THEN
  549. c--------Equation d'entrainement
  550. write(6,*) 'EQUATION D''ENTRAINEMENT'
  551. CALL TURB6(MPOVA6.VPOCHA,MPOVA5.VPOCHA,MPOVA3.VPOCHA,
  552. & HHN,HNM,CFN,CEN,D1N,VPOCHA(/1))
  553. DO I=1,VPOCHA(/1)
  554. Coef1(I)=(MPOVA4.VPOCHA(I,1))
  555. & *((MPOVA3.VPOCHA(I,1))**(-1))
  556. Coef2(I)=CEN(I)
  557. MPOVA9.VPOCHA(I,1)=HNM(I)
  558. MPOVA8.VPOCHA(I,1)=CFN(I)
  559. MPOVAA.VPOCHA(I,1)=D1N(I)
  560. ENDDO
  561. ENDIF
  562. ENDIF
  563. c IF(ferm.EQ.6) THEN
  564. C---- Convection naturelle laminaire
  565. c write(6,*) 'CONVECTION NATURELLE'
  566. c IF(ctype.EQ.1) THEN
  567. C--------Equation de quantité de mouvement
  568. c write(6,*) 'EQUATION DE QDM'
  569. c CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
  570. c & MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1))
  571. c DO I=1,VPOCHA(/1)
  572. c Coef1(I) = 0
  573. c Coef2(I) = SM1(I)
  574. c MPOVA9.VPOCHA(I,1)=FP(I)
  575. c MPOVA8.VPOCHA(I,1)=PHIP(I)
  576. c MPOVAA.VPOCHA(I,1)=Delta(I)
  577. c ENDDO
  578. c ENDIF
  579. c IF(ctype.EQ.4.) THEN
  580. C--------Equation de la chaleur
  581. c write(6,*) 'EQUATION DE LA CHALEUR'
  582. c CALL CNAT(MPOVA5.VPOCHA,MPOVA6.VPOCHA,MPOVA3.VPOCHA,
  583. c & MPOVA4.VPOCHA,UCL,Delta,SM1,SM2,FP,PHIP,VPOCHA(/1))
  584. c DO I=1,VPOCHA(/1)
  585. c Coef1(I) = 0
  586. c Coef2(I) = SM2(I)
  587. c MPOVA9.VPOCHA(I,1)=FP(I)
  588. c MPOVA8.VPOCHA(I,1)=PHIP(I)
  589. c MPOVAA.VPOCHA(I,1)=Delta(I)
  590. c ENDDO
  591. c ENDIF
  592. c ENDIF
  593.  
  594. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  595. c write(6,*)' nd1,nd2,nd3=',nd1,nd2,nd3
  596. c write(6,*)' NBELEM=',NBELEM
  597.  
  598.  
  599. DO 101 K=1,NBELEM
  600. c write(6,*)' Element n ',k
  601. c write(6,1001)(num(ii,k),ii=1,nbnn)
  602.  
  603. C? DO 102 I=1,NBNN
  604. C? IP1=num(I,K)
  605. C? XYZ(1,I)=XCOOR((IP1-1)*(IDIM+1) +1)
  606. C? XYZ(2,I)=XCOOR((IP1-1)*(IDIM+1) +2)
  607. c write(6,*)'Coordonnees element ',K,' noeud ',ip1
  608. c write(6,1002)(XYZ(M,I),M=1,IDIM)
  609. C? 102 CONTINUE
  610.  
  611. C? CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
  612. C? & NES,IDIM,NBNN,NPG,IAXI,AIRE,AJ)
  613.  
  614. KDEB=1
  615. KFIN=1
  616. LRV =1
  617. IKOMP=0
  618. ANUK(1)=1.D0
  619.  
  620. C CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN
  621. DT=0.D0
  622.  
  623. CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  624. & NES,NP,NPG,IAXI,XCOOR,
  625. & NUM(1,K),KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
  626. & TN,LECT,UN,LECT,NPTI,ANUK,
  627. & WT,WS,HR,PGSQ,RPG,AJ,AIRE,UIL,DUIL,
  628. & DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  629.  
  630. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  631. C>>>>>>>>>>Calcul des matrices élémentaires>>>>>>>>>>>>>>>>>
  632. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  633. c calcul des matrices masses élémentaires
  634. c write(6,*)' NBNN=',nbnn
  635. DO 103 I=1,NBNN
  636. DO 104 J=1,NBNN
  637. U=0.D0
  638. DO 105 L=1,NPG
  639. V=0.D0
  640. DO 1051 L1=1,NBNN
  641. J1=lect(num(L1,K))
  642. V=V+FN(L1,L)*Coef1(J1)
  643. 1051 CONTINUE
  644. U=U+(WT(J,L)*FN(I,L)*PGSQ(L))*V
  645. 105 CONTINUE
  646. IPM1.AM(K,I,J)=U
  647. c write(6,*) 'M(',I,',',J,',',K,')=',IPM1.AM(K,I,J)
  648. 104 CONTINUE
  649. 103 CONTINUE
  650.  
  651. ax=1.d0
  652. ay=1.d0
  653. c calcul des matrices de convection élémentaires
  654. DO 106 I=1,NBNN
  655. DO 107 J=1,NBNN
  656. c J1 = lect(num(J,K))
  657. U=0.D0
  658. DO 108 L=1,NPG
  659. V=0.D0
  660. DO 1081 L1=1,NBNN
  661. J1 = lect(num(L1,K))
  662. V = V+FN(L1,L)*1.
  663. 1081 CONTINUE
  664. U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L))*V
  665. c write(6,*) 'hr(',J,',',L,')=',HR(1,J,L)
  666. c write(6,*) 'FN(',I,',',L,',',K,')=',FN(I,L)
  667. c write(6,*) 'WT(',I,',',L,',',K,')=',WT(I,L)
  668. c U=U+(WT(I,L)*((ax*HR(1,J,L)) + (ay*HR(1,J,L)) )*PGSQ(L))
  669. c U=U+(WT(J,L)*((ax*HR(1,I,L)) )*PGSQ(L))
  670. c &*MPOVA2.VPOCHA(J1,1)
  671. 108 CONTINUE
  672. C write(6,*)'K,i,j=', K,i,j
  673. IPM2.AM(K,I,J) = U
  674. 107 CONTINUE
  675. 106 CONTINUE
  676.  
  677. C calcul du terme source
  678. DO 109 I=1,NBNN
  679. I1= lect(num(I,K))
  680. U=0.D0
  681. DO 110 L=1,NPG
  682. V=0.D0
  683. v2=0.D0
  684. DO 1101 J=1,NBNN
  685. J1 =lect(num(J,K))
  686. V=V+FN(J,L)*Coef2(J1)
  687. 1101 CONTINUE
  688. U=U+V*WT(I,L)*PGSQ(L)
  689. 110 CONTINUE
  690. VPOCHA(I1,1)=VPOCHA(I1,1)+ U
  691. 109 CONTINUE
  692.  
  693. C calcul de la vitesse de transpiration
  694. DO 111 I=1,NBNN
  695. I1 = lect(num(I,K))
  696. U=0.D0
  697. DO 112 J=1,NBNN
  698. J1 = lect(num(J,K))
  699. U = U+HR(1,J,1)*(D1N(J1)*MPOVA3.VPOCHA(J1,1))
  700. 112 CONTINUE
  701. VT(I1) = U
  702. c write(6,*) 'VT(',I1,')=',VT(I1)
  703. MPOVAB.VPOCHA(I1,1)= VT(I1)
  704. 111 CONTINUE
  705.  
  706. 101 CONTINUE
  707.  
  708.  
  709. C Sortie des valeurs des coefficients H, Cf, Ce et VT
  710.  
  711. CALL ECMO(KINC,'H','CHPOINT',MCHPO9)
  712. CALL ECMO(KINC,'CF','CHPOINT',MCHPO8)
  713. CALL ECMO(KINC,'D1','CHPOINT',MCHPOA)
  714. CALL ECMO(KINC,'Vtra','CHPOINT',MCHPOB)
  715.  
  716.  
  717. SEGDES IPM2,IMAT2,IPM1,IMAT1
  718. SEGDES MATRIK
  719. SEGDES TRAV
  720. SEGDES IZFFM,IZHR
  721. SEGDES LINCO
  722. C SEGDES MTEL
  723.  
  724. CALL ECROBJ('MATRIK ',MATRIK)
  725. CALL ACTOBJ('CHPOINT ',MCHPOI,1)
  726. CALL ECROBJ('CHPOINT ',MCHPOI)
  727. RETURN
  728. 1003 FORMAT(6(1X,1PE11.4))
  729. 1002 FORMAT(10(1X,1PE11.4))
  730. 1001 FORMAT(20(1X,I5))
  731. END
  732.  
  733.  
  734.  

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