Télécharger clmi.eso

Retour à la liste

Numérotation des lignes :

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

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