Télécharger kon11.eso

Retour à la liste

Numérotation des lignes :

  1. C KON11 SOURCE KK2000 14/04/10 21:15:15 8032
  2. SUBROUTINE KON11(IFLU,LOGME)
  3. C************************************************************************
  4. C
  5. C PROJET : CASTEM 2000
  6. C
  7. C NOM : KON11
  8. C
  9. C DESCRIPTION : Subroutine appellée par KONV1
  10. C
  11. C Modelisation 2D/3D des equations d'Euler
  12. C
  13. C Calcul du flux/residu, DELTAT
  14. C
  15. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  16. C
  17. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF
  18. C
  19. C************************************************************************
  20. C
  21. C APPELES (Calcul) : KONFL1 (2D, gaz "calorically perfect")
  22. C KONFL3 (3D, gaz "calorically perfect")
  23. C
  24. C************************************************************************
  25. C
  26. C ENTRÉES : IFLU : si EGAL à 1 on veut calculer le flux
  27. C IFLU : si EGAL à 2 on veut calculer le residu
  28. C
  29. C LOGME : .TRUE./.FALSE. Multi-espece/mono-espece
  30. C
  31. C*** SYNTAXE
  32. C
  33. C Discrétisation en VF "cell-centered" des équations d'Euler pour
  34. C un gaz parfait polytropique
  35. C Inconnues: densités, quantité de mouvement, énergie totale par
  36. C unité de volumes (variables conservatives)
  37. C
  38. C RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMONO' MOT1 MOT2 LMOT1
  39. C MOD1 MCHAM1 MCHAM2 MCHAM3 MCHAM4 (MAIL1) ;
  40. C
  41. C or (Bas Mach mono-espece)
  42. C
  43. C RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMONO' MOT1 MOT2 LMOT1
  44. C MOD1 MCHAM1 MCHAM2 MCHAM3 MCHAM4
  45. C CHPO1 CHPO2 (MAIL1) ;
  46. C
  47. C or (multiespece)
  48. C
  49. C RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMULT' MOT1 MOT2 LMOT1
  50. C MOD1 MCHAM1 MCHAM2 MCHAM3 MCHAM4 MCHAM5 TABP ;
  51. C
  52. C ENTREES
  53. C
  54. C MOT1 : objet de type MOT
  55. C Il vaut 'RESI' si on veut calculer le résidu
  56. C Il vaut 'FLUX' si on veut calculer le flux
  57. C
  58. C MOT2 : objet de type MOT
  59. C Il indique la méthode de décentrement:
  60. C 'GODUNOV' = solveur exacte
  61. C 'VANLEER' = solveur de van Leer
  62. C 'VLH' = solveur de van Leer Hanel
  63. C 'HUSVL' = HUS (van Leer + Osher)
  64. C 'HUSVLH' = HUS (van Leer Hanel + Osher)
  65. C 'AUSMPLUS' = AUSM+
  66. C 'ROE' = solveur de Roe
  67. C 'SS' = solveur choc-choc
  68. C 'RUSANOV' = solveur de Rusanov
  69. C 'AUSMPLM' = AUSM+ low Mach
  70. C 'RUSANOLM' = solveur de Rusanov low Mach
  71. C 'CENTERED' = solveur centré
  72. C 'ROELM' = solveur de Roe-Turkel
  73. C 'HLLC' = solveur 'HLLC'
  74. C 'HLLCLM' = solveur 'HLLC'-Turkel low Mach
  75. C 'AUSMPUP' = AUSM+up low Mach
  76. C
  77. C
  78. C LMOT1 : objet de type LISTMOTS
  79. C Noms de composantes du résultat (RCHPO1)
  80. C Il contient dans l'ordre suivant: le noms de la densité,
  81. C du momentum, de l'énergie totale par unité de volume
  82. C (le densité des constituants qui sont dans TABP.'ESPEULE'
  83. C dans le cas multi-espece)
  84. C
  85. C MOD1 : objet modele de type Navier_Stokes
  86. C
  87. C MCHAM1 : MCHAML contenant la masse volumique, qui a comme
  88. C SPG (support géométrique) l'indice 'FACEL' de
  89. C MOD1 (une composante, 'SCAL')
  90. C
  91. C MCHAM2 : MCHAML contenant la vitesse et les cosinus
  92. C directeurs du repère locale (n,t) dans le repère
  93. C global (x,y) (dans le cas 2D 6 composantes:
  94. C * 'UN' = vitesse normale (SPG =('DOMA' MOD1 'FACEL'))
  95. C * 'UT' = vitesse tangentielle (SPG =('DOMA' MOD1 'FACEL'))
  96. C * 'NX' = n.x (SPG = 'FACE')
  97. C * 'NY' = n.y (SPG = 'FACE')
  98. C * 'TX' = t.x (SPG = 'FACE')
  99. C * 'TY' = t.y (SPG = 'FACE')).
  100. C
  101. C MCHAM3 : MCHAML (SPG =('DOMA' MOD1 'FACEL')) contenant la pression du
  102. C gaz (une seule composante, 'SCAL').
  103. C
  104. C MCHAM4 : MCHAML (SPG =('DOMA' MOD1 'FACEL')) contenant le "gamma" du
  105. C gaz (une seule composante, 'SCAL').
  106. C
  107. C MCHAM5 : MCHAML (SPG =('DOMA' MOD1 'FACEL')) contenant les fractions
  108. C massiques du gaz dans le cas multi-espece
  109. C (noms de composantes dans TABP . 'ESPEULE')
  110. C
  111. C TABP : la table des proprietés du gaz (voir operateur 'PRIM',
  112. C modele 'PERFMULT')
  113. C
  114. C CHPO1 : vitesse de cut-off (SPG =('DOMA' MOD1 'CENTRE'), une seule
  115. C composante, 'SCAL')
  116. C A donner que pour les solveurs Bas Mach
  117. C
  118. C CHPO2 : deuxieme vitesse de cut-off (SPG = ('DOMA' MOD1 'CENTRE'),
  119. C une seule composante, 'SCAL')
  120. C A donner que pour les solveurs Bas Mach
  121. C
  122. C (MAIL1) : maillage de 'POI1' qui contient de point 'FACE' ou on ne
  123. C veut pas calculer les flux
  124. C
  125. C SORTIES
  126. C
  127. C
  128. C RCHPO1 : objet de type CHPOINT (composantes = LMOT1)
  129. C Residu si MOT2 = 'RESI' (SPG = ('DOMA' MOD1 'CENTRE'))
  130. C Flux si MOT2 = 'FLUX' (SPG = ('DOMA' MOD1 'FACE'))
  131. C
  132. C RFLOT1 : objet de type FLOTTANT
  133. C Il est le temps caracteristique associé à l'onde la plus
  134. C rapide.
  135. C
  136. C Remarque
  137. C --------
  138. C
  139. C RCHPO1 est égal à:
  140. C * la derivé temporelle des inconnues si l'option 'RESI' est utilisée
  141. C * la projection du flux convectif sur (DOMA MOD1 'XXNORMAF') si
  142. C l'option 'FLUX' est utilisée
  143. C
  144. C
  145. C************************************************************************
  146. C
  147. C HISTORIQUE (Anomalies et modifications éventuelles)
  148. C
  149. C HISTORIQUE :
  150. C 01/05/2005 par T. KLOCZKO (DEN/DM2S/SFME/LTMF)
  151. C Ajout d'une méthode de décentrement pour les écoulements
  152. C bas-Mach.
  153. C
  154. C MOT2 : 'ROELM' = schéma de ROE avec préconditionnement Turkel
  155. C
  156. C************************************************************************
  157.  
  158. IMPLICIT INTEGER(I-N)
  159. IMPLICIT REAL*8(A-H,O-Z)
  160.  
  161. -INC CCOPTIO
  162. -INC SMLMOTS
  163. -INC SMCHPOI
  164. -INC SMELEME
  165. POINTEUR MLMVIT.MLMOTS,MLMOTY.MLMOTS,MLMODT.MLMOTS
  166. C
  167. INTEGER IFLU, NBMET, INDMET, IRET, INDIC, NBCOMP
  168. & ,IDOMA, MELEMC, MELEMF, MELEFE, MELLIM, ICHPSU, ICHPDI
  169. & ,ICHPVO, INORM, JGN, JGM
  170. & , NC, NESP, ICACCA, IPGAZ
  171. & ,ILIINC
  172. & ,IROF, IVITF, IPF, IGAMF, IFRMAF, IUINF, IUPRI, IFLIM
  173. & ,ICHFLU, ICHRES, INEFMD, ICOND, MMODEL
  174. PARAMETER (NBMET=16)
  175. REAL*8 DT
  176. CHARACTER*8 LMETO(NBMET), TYPE
  177. CHARACTER*4 MOT
  178. CHARACTER*(40) MESERR
  179. LOGICAL LOGNC, LOGAN, LOGME
  180. C
  181. DATA LMETO/'GODUNOV ','VANLEER','VLH ',
  182. & 'HUSVL ','HUSVLH ','AUSMPLUS','ROE ',
  183. & 'SS ','AUSMPLM','RUSANOV', 'RUSANOLM',
  184. & 'CENTERED','ROELM ','HLLC ','HLLCLM ',
  185. & 'AUSMPUP ' /
  186. C
  187. C**** Initialisation des variables pour la gestion des erreurs.
  188. C
  189. LOGNC = .FALSE.
  190. LOGAN = .FALSE.
  191. MESERR = ' '
  192. C
  193. C**** Metode utilisée
  194. C
  195. CALL LIRMOT(LMETO,NBMET,INDMET,1)
  196. IF(IERR .NE. 0)GOTO 9999
  197. IF(INDMET .EQ. 0)THEN
  198. C
  199. C******** Message d'erreur standard
  200. C 251 2
  201. C Tentative d'utilisation d'une option non implémentée
  202. C
  203. CALL ERREUR(251)
  204. ENDIF
  205. C
  206. C**********************************
  207. C**** Lecture de l'objet MODELE ***
  208. C**********************************
  209. C
  210. ICOND = 1
  211. CALL QUETYP(TYPE,ICOND,IRET)
  212. C
  213. IF(IRET.EQ.0.AND.TYPE.NE.'MMODEL')THEN
  214. WRITE(IOIMP,*)' On attend un objet MMODEL'
  215. GOTO 9999
  216. ENDIF
  217. CALL LIROBJ('MMODEL',MMODEL,ICOND,IRET)
  218. IF(IERR.NE.0)GOTO 9999
  219. CALL LEKMOD(MMODEL,IDOMA,INEFMD)
  220. C INEFMD inutilisé
  221. IF(IERR.NE.0)GOTO 9999
  222. C
  223. C**** Centre, FACE, FACEL
  224. C
  225. CALL LEKTAB(IDOMA,'CENTRE',MELEMC)
  226. IF(IERR .NE. 0) GOTO 9999
  227. C
  228. CALL LEKTAB(IDOMA,'FACE',MELEMF)
  229. IF(IERR .NE. 0) GOTO 9999
  230. C
  231. CALL LEKTAB(IDOMA,'FACEL',MELEFE)
  232. IF(IERR .NE. 0) GOTO 9999
  233. C
  234. C**** Lecture du CHPOINT contenant les surfaces des faces.
  235. C
  236. CALL LEKTAB(IDOMA,'XXSURFAC',ICHPSU)
  237. IF(IERR .NE. 0) GOTO 9999
  238. INDIC = 1
  239. NBCOMP = 1
  240. MOT = 'SCAL'
  241. CALL QUEPOI(ICHPSU, MELEMF, INDIC, NBCOMP, MOT)
  242. IF(IERR .NE. 0) GOTO 9999
  243. C
  244. C**** Lecture du CHPOINT contenant les diametres minimums.
  245. C
  246. CALL LEKTAB(IDOMA,'XXDIEMIN',ICHPDI)
  247. IF(IERR .NE. 0) GOTO 9999
  248. INDIC = 1
  249. NBCOMP = 1
  250. MOT = 'SCAL'
  251. CALL QUEPOI(ICHPDI, MELEMC, INDIC, NBCOMP, MOT)
  252. IF(IERR .NE. 0) GOTO 9999
  253. C
  254. C**** Lecture du CHPOINT contenant les volumes
  255. C
  256. CALL LEKTAB(IDOMA,'XXVOLUM',ICHPVO)
  257. IF(IERR .NE. 0) GOTO 9999
  258. INDIC = 1
  259. NBCOMP = 1
  260. MOT = 'SCAL'
  261. CALL QUEPOI(ICHPVO, MELEMC, INDIC, NBCOMP, MOT)
  262. IF(IERR .NE. 0) GOTO 9999
  263. C
  264. C**** Les normales aux faces
  265. C
  266. IF(IDIM .EQ. 2)THEN
  267. C Que les normales
  268. CALL LEKTAB(IDOMA,'XXNORMAF',INORM)
  269. IF(IERR .NE. 0) GOTO 9999
  270. JGN = 4
  271. JGM = 2
  272. SEGINI MLMVIT
  273. MLMVIT.MOTS(1) = 'UX '
  274. MLMVIT.MOTS(2) = 'UY '
  275. CALL QUEPO1(INORM, MELEMF, MLMVIT)
  276. SEGSUP MLMVIT
  277. IF(IERR .NE. 0) GOTO 9999
  278. ELSE
  279. C Les normales et les tangentes
  280. TYPE = ' '
  281. CALL ACMO(IDOMA,'MATROT',TYPE,INORM)
  282. IF (TYPE .NE. 'CHPOINT ') THEN
  283. CALL MATRAN(IDOMA,INORM)
  284. IF(IERR .NE. 0) GOTO 9999
  285. ENDIF
  286. JGN = 4
  287. JGM = 9
  288. SEGINI MLMVIT
  289. MLMVIT.MOTS(1) = 'UX '
  290. MLMVIT.MOTS(2) = 'UY '
  291. MLMVIT.MOTS(3) = 'UZ '
  292. MLMVIT.MOTS(4) = 'RX '
  293. MLMVIT.MOTS(5) = 'RY '
  294. MLMVIT.MOTS(6) = 'RZ '
  295. MLMVIT.MOTS(7) = 'MX '
  296. MLMVIT.MOTS(8) = 'MY '
  297. MLMVIT.MOTS(9) = 'MZ '
  298. CALL QUEPO1(INORM, MELEMF, MLMVIT)
  299. SEGSUP MLMVIT
  300. ENDIF
  301. C
  302. C********************************
  303. C**** Fin table domaine *********
  304. C********************************
  305. C
  306. C**** On va lire les pointeurs des MCHAMLs
  307. C Lecture du MCHAML 'FACEL' densité
  308. C
  309. TYPE='MCHAML '
  310. CALL LIROBJ(TYPE,IROF,1,IRET)
  311. IF(IERR.NE.0) GOTO 9999
  312. C
  313. C**** Lecture du MCHAML 'FACEL' vitesse
  314. C
  315. TYPE='MCHAML '
  316. CALL LIROBJ(TYPE,IVITF,1,IRET)
  317. IF(IERR .NE. 0) GOTO 9999
  318. C
  319. C**** Lecture du MCHAML 'FACEL' contenant la pression
  320. C
  321. TYPE='MCHAML '
  322. CALL LIROBJ(TYPE,IPF,1,IRET)
  323. IF(IERR .NE. 0) GOTO 9999
  324. C
  325. C**** Lecture du MCHAML 'FACEL' contenant les gamma
  326. C
  327. TYPE='MCHAML '
  328. CALL LIROBJ(TYPE,IGAMF,1,IRET)
  329. IF(IERR .NE. 0) GOTO 9999
  330. C
  331. C**** Si LOGME -> MULTIESPECES
  332. C
  333. IF(LOGME)THEN
  334. C
  335. C******** Lecture du MCHAML 'FACEL' contenant les fractiones massiques
  336. C
  337. TYPE='MCHAML '
  338. CALL LIROBJ(TYPE,IFRMAF,1,IRET)
  339. IF(IERR .NE. 0) GOTO 9999
  340. C
  341. C********** Lecture de la table qui contient le proprieté du gaz
  342. C
  343. TYPE='TABLE '
  344. CALL LIROBJ(TYPE,IPGAZ,1,IRET)
  345. IF(IERR .NE. 0) GOTO 9999
  346. TYPE=' '
  347. CALL ACMO(IPGAZ,'ESPEULE ',TYPE,MLMOTY)
  348. IF(IERR .NE. 0) GOTO 9999
  349. IF(TYPE .NE. 'LISTMOTS')THEN
  350. C
  351. C********** Message d'erreur standard
  352. C -301 0 %m1:40
  353. C
  354. MOTERR(1:40) = 'TABGAS = ??? '
  355. WRITE(IOIMP,*) MOTERR
  356. C
  357. C********** Message d'erreur standard
  358. C 21 2
  359. C Données incompatibles
  360. C
  361. CALL ERREUR(21)
  362. GOTO 9999
  363. ENDIF
  364. C
  365. SEGACT MLMOTY
  366. NESP = MLMOTY.MOTS(/2)
  367. C
  368. ELSE
  369. IFRMAF = 0
  370. NESP=0
  371. ENDIF
  372. C
  373. C**** La list des inconnues
  374. C
  375. TYPE='LISTMOTS'
  376. CALL LIROBJ(TYPE,ILIINC,1,IRET)
  377. IF(IERR .NE. 0) GOTO 9999
  378. MLMOTS = ILIINC
  379. SEGACT MLMOTS
  380. NC = MLMOTS.MOTS(/2)
  381. SEGDES MLMOTS
  382. IF(NC .NE. (IDIM+2+NESP))THEN
  383. MOTERR(1:40) = 'LISTINCO = ???'
  384. WRITE(IOIMP,*) MOTERR
  385. C
  386. C******* Message d'erreur standard
  387. C 21 2
  388. C Données incompatibles
  389. C
  390. CALL ERREUR(21)
  391. GOTO 9999
  392. ENDIF
  393. C
  394. C**** Boundary condition
  395. C
  396. IRET=0
  397. TYPE='MAILLAGE'
  398. CALL LIROBJ(TYPE,IFLIM,0,IRET)
  399. IF(IERR.NE.0)GOTO 9999
  400. IF(IRET .EQ. 0)THEN
  401. MELLIM = 0
  402. ELSE
  403. MELEME=IFLIM
  404. SEGACT MELEME
  405. ICACCA=MELEME.NUM(/2)
  406. IF(ICACCA .EQ. 0)THEN
  407. MELLIM = 0
  408. ELSE
  409. MELLIM = IFLIM
  410. ENDIF
  411. SEGDES MELEME
  412. ENDIF
  413. C
  414. C**** Bas Mach (AUSMPLM, RUSANOLM, ROELM, HLLCLM, AUSMPUP)
  415. C
  416. IF((INDMET .EQ. 9) .OR. (INDMET .EQ. 11) .OR.
  417. & (INDMET .EQ. 13) .OR. (INDMET .EQ. 15) .OR.
  418. & (INDMET .EQ. 16)) THEN
  419. TYPE = 'CHPOINT '
  420. C
  421. C******* Reference speed
  422. C
  423. CALL LIROBJ(TYPE,IUINF,1,IRET)
  424. IF(IERR .NE. 0) GOTO 9999
  425. INDIC = 1
  426. NBCOMP = 1
  427. MOT = 'SCAL'
  428. CALL QUEPOI(IUINF, MELEMC, INDIC, NBCOMP, MOT)
  429. IF(IERR .NE. 0) GOTO 9999
  430. C
  431. C******* Minimal cutoff
  432. C
  433. TYPE = 'CHPOINT '
  434. CALL LIROBJ(TYPE,IUPRI,1,IRET)
  435. IF(IERR .NE. 0) GOTO 9999
  436. INDIC = 1
  437. NBCOMP = 1
  438. MOT = 'SCAL'
  439. CALL QUEPOI(IUPRI, MELEMC, INDIC, NBCOMP, MOT)
  440. IF(IERR .NE. 0) GOTO 9999
  441. C
  442. ELSE
  443. IUINF=0
  444. IUPRI=0
  445. ENDIF
  446. C
  447. C**** Creation des flux aux interfaces
  448. C
  449. TYPE = 'CHPOINT '
  450. CALL KRCHP1(TYPE, MELEMF, ICHFLU, MLMOTS)
  451. C
  452. C**** Calcul des flux et du pas du temps.
  453. C
  454. IF(IDIM .EQ. 2)THEN
  455. CALL KONFL1(LOGME,INDMET,
  456. & IROF,IVITF,IPF,IGAMF,IFRMAF,
  457. & ICHPSU,ICHPDI,
  458. & MELEMC,MELEMF,MELEFE,MELLIM,
  459. & IUINF,IUPRI,
  460. & ICHFLU,
  461. & DT,
  462. & LOGNC,LOGAN,MESERR)
  463. ELSEIF(IDIM .EQ. 3)THEN
  464. CALL KONFL3(LOGME,INDMET,
  465. & IROF,IVITF,IPF,IGAMF,IFRMAF,
  466. & ICHPSU,ICHPDI,
  467. & MELEMC,MELEMF,MELEFE,MELLIM,
  468. & IUINF,IUPRI,
  469. & ICHFLU,
  470. & DT,
  471. & LOGNC,LOGAN,MESERR)
  472. ENDIF
  473. C
  474. IF(LOGAN)THEN
  475. C
  476. C******* Anomalie detectée
  477. C
  478. C
  479. C******* Message d'erreur standard
  480. C -301 0
  481. C %m1:40
  482. C
  483. MOTERR(1:40) = MESERR(1:40)
  484. WRITE(IOIMP,*) MOTERR(1:40)
  485. C
  486. C******* Message d'erreur standard
  487. C 5 3
  488. C Erreur anormale.contactez votre support
  489. C
  490. CALL ERREUR(5)
  491. GOTO 9999
  492. ENDIF
  493. IF(LOGNC)THEN
  494. C
  495. C******* Message d'erreur standard
  496. C -301 0
  497. C %m1:40
  498. C
  499. MOTERR(1:40) = MESERR(1:40)
  500. WRITE(IOIMP,*) MOTERR(1:40)
  501. C
  502. C******* Message d'erreur standard
  503. C 460 2
  504. C Pas de convergence dans les itérations internes
  505. C
  506. CALL ERREUR(460)
  507. GOTO 9999
  508. ENDIF
  509. C
  510. C**** Calcul de residu (IFLU = 2)
  511. C
  512. IF(IFLU .EQ. 2)THEN
  513. TYPE = 'CHPOINT '
  514. CALL KRCHP1(TYPE, MELEMC, ICHRES, MLMOTS)
  515. C
  516. CALL KONRE1(MELEMC,MELEMF,MELEFE,ICHPVO,
  517. & ICHFLU, ICHRES,
  518. & LOGAN,MESERR)
  519. IF(LOGAN)THEN
  520. C
  521. C******* Anomalie detectée
  522. C
  523. C
  524. C******* Message d'erreur standard
  525. C -301 0
  526. C %m1:40
  527. C
  528. MOTERR(1:40) = MESERR(1:40)
  529. WRITE(IOIMP,*) MOTERR(1:40)
  530. C
  531. C******* Message d'erreur standard
  532. C 5 3
  533. C Erreur anormale.contactez votre support
  534. C
  535. CALL ERREUR(5)
  536. GOTO 9999
  537. ENDIF
  538. ELSE
  539. ICHRES = 0
  540. ENDIF
  541. C
  542. C**** Ecriture des resultats
  543. C
  544. CALL ECRREE(DT)
  545. TYPE = 'CHPOINT '
  546. IF(ICHRES .NE. 0) CALL ECROBJ(TYPE,ICHRES)
  547. IF(ICHFLU .NE. 0) CALL ECROBJ(TYPE,ICHFLU)
  548. 9999 CONTINUE
  549. RETURN
  550. END
  551.  
  552.  
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  

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