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.  
  162. -INC PPARAM
  163. -INC CCOPTIO
  164. -INC SMLMOTS
  165. -INC SMCHPOI
  166. -INC SMELEME
  167. POINTEUR MLMVIT.MLMOTS,MLMOTY.MLMOTS,MLMODT.MLMOTS
  168. C
  169. INTEGER IFLU, NBMET, INDMET, IRET, INDIC, NBCOMP
  170. & ,IDOMA, MELEMC, MELEMF, MELEFE, MELLIM, ICHPSU, ICHPDI
  171. & ,ICHPVO, INORM, JGN, JGM
  172. & , NC, NESP, ICACCA, IPGAZ
  173. & ,ILIINC
  174. & ,IROF, IVITF, IPF, IGAMF, IFRMAF, IUINF, IUPRI, IFLIM
  175. & ,ICHFLU, ICHRES, INEFMD, ICOND, MMODEL
  176. PARAMETER (NBMET=16)
  177. REAL*8 DT
  178. CHARACTER*8 LMETO(NBMET), TYPE
  179. CHARACTER*4 MOT
  180. CHARACTER*(40) MESERR
  181. LOGICAL LOGNC, LOGAN, LOGME
  182. C
  183. DATA LMETO/'GODUNOV ','VANLEER','VLH ',
  184. & 'HUSVL ','HUSVLH ','AUSMPLUS','ROE ',
  185. & 'SS ','AUSMPLM','RUSANOV', 'RUSANOLM',
  186. & 'CENTERED','ROELM ','HLLC ','HLLCLM ',
  187. & 'AUSMPUP ' /
  188. C
  189. C**** Initialisation des variables pour la gestion des erreurs.
  190. C
  191. LOGNC = .FALSE.
  192. LOGAN = .FALSE.
  193. MESERR = ' '
  194. C
  195. C**** Metode utilisée
  196. C
  197. CALL LIRMOT(LMETO,NBMET,INDMET,1)
  198. IF(IERR .NE. 0)GOTO 9999
  199. IF(INDMET .EQ. 0)THEN
  200. C
  201. C******** Message d'erreur standard
  202. C 251 2
  203. C Tentative d'utilisation d'une option non implémentée
  204. C
  205. CALL ERREUR(251)
  206. ENDIF
  207. C
  208. C**********************************
  209. C**** Lecture de l'objet MODELE ***
  210. C**********************************
  211. C
  212. ICOND = 1
  213. CALL QUETYP(TYPE,ICOND,IRET)
  214. C
  215. IF(IRET.EQ.0.AND.TYPE.NE.'MMODEL')THEN
  216. WRITE(IOIMP,*)' On attend un objet MMODEL'
  217. GOTO 9999
  218. ENDIF
  219. CALL LIROBJ('MMODEL',MMODEL,ICOND,IRET)
  220. IF(IERR.NE.0)GOTO 9999
  221. CALL LEKMOD(MMODEL,IDOMA,INEFMD)
  222. C INEFMD inutilisé
  223. IF(IERR.NE.0)GOTO 9999
  224. C
  225. C**** Centre, FACE, FACEL
  226. C
  227. CALL LEKTAB(IDOMA,'CENTRE',MELEMC)
  228. IF(IERR .NE. 0) GOTO 9999
  229. C
  230. CALL LEKTAB(IDOMA,'FACE',MELEMF)
  231. IF(IERR .NE. 0) GOTO 9999
  232. C
  233. CALL LEKTAB(IDOMA,'FACEL',MELEFE)
  234. IF(IERR .NE. 0) GOTO 9999
  235. C
  236. C**** Lecture du CHPOINT contenant les surfaces des faces.
  237. C
  238. CALL LEKTAB(IDOMA,'XXSURFAC',ICHPSU)
  239. IF(IERR .NE. 0) GOTO 9999
  240. INDIC = 1
  241. NBCOMP = 1
  242. MOT = 'SCAL'
  243. CALL QUEPOI(ICHPSU, MELEMF, INDIC, NBCOMP, MOT)
  244. IF(IERR .NE. 0) GOTO 9999
  245. C
  246. C**** Lecture du CHPOINT contenant les diametres minimums.
  247. C
  248. CALL LEKTAB(IDOMA,'XXDIEMIN',ICHPDI)
  249. IF(IERR .NE. 0) GOTO 9999
  250. INDIC = 1
  251. NBCOMP = 1
  252. MOT = 'SCAL'
  253. CALL QUEPOI(ICHPDI, MELEMC, INDIC, NBCOMP, MOT)
  254. IF(IERR .NE. 0) GOTO 9999
  255. C
  256. C**** Lecture du CHPOINT contenant les volumes
  257. C
  258. CALL LEKTAB(IDOMA,'XXVOLUM',ICHPVO)
  259. IF(IERR .NE. 0) GOTO 9999
  260. INDIC = 1
  261. NBCOMP = 1
  262. MOT = 'SCAL'
  263. CALL QUEPOI(ICHPVO, MELEMC, INDIC, NBCOMP, MOT)
  264. IF(IERR .NE. 0) GOTO 9999
  265. C
  266. C**** Les normales aux faces
  267. C
  268. IF(IDIM .EQ. 2)THEN
  269. C Que les normales
  270. CALL LEKTAB(IDOMA,'XXNORMAF',INORM)
  271. IF(IERR .NE. 0) GOTO 9999
  272. JGN = 4
  273. JGM = 2
  274. SEGINI MLMVIT
  275. MLMVIT.MOTS(1) = 'UX '
  276. MLMVIT.MOTS(2) = 'UY '
  277. CALL QUEPO1(INORM, MELEMF, MLMVIT)
  278. SEGSUP MLMVIT
  279. IF(IERR .NE. 0) GOTO 9999
  280. ELSE
  281. C Les normales et les tangentes
  282. TYPE = ' '
  283. CALL ACMO(IDOMA,'MATROT',TYPE,INORM)
  284. IF (TYPE .NE. 'CHPOINT ') THEN
  285. CALL MATRAN(IDOMA,INORM)
  286. IF(IERR .NE. 0) GOTO 9999
  287. ENDIF
  288. JGN = 4
  289. JGM = 9
  290. SEGINI MLMVIT
  291. MLMVIT.MOTS(1) = 'UX '
  292. MLMVIT.MOTS(2) = 'UY '
  293. MLMVIT.MOTS(3) = 'UZ '
  294. MLMVIT.MOTS(4) = 'RX '
  295. MLMVIT.MOTS(5) = 'RY '
  296. MLMVIT.MOTS(6) = 'RZ '
  297. MLMVIT.MOTS(7) = 'MX '
  298. MLMVIT.MOTS(8) = 'MY '
  299. MLMVIT.MOTS(9) = 'MZ '
  300. CALL QUEPO1(INORM, MELEMF, MLMVIT)
  301. SEGSUP MLMVIT
  302. ENDIF
  303. C
  304. C********************************
  305. C**** Fin table domaine *********
  306. C********************************
  307. C
  308. C**** On va lire les pointeurs des MCHAMLs
  309. C Lecture du MCHAML 'FACEL' densité
  310. C
  311. TYPE='MCHAML '
  312. CALL LIROBJ(TYPE,IROF,1,IRET)
  313. IF(IERR.NE.0) GOTO 9999
  314. C
  315. C**** Lecture du MCHAML 'FACEL' vitesse
  316. C
  317. TYPE='MCHAML '
  318. CALL LIROBJ(TYPE,IVITF,1,IRET)
  319. IF(IERR .NE. 0) GOTO 9999
  320. C
  321. C**** Lecture du MCHAML 'FACEL' contenant la pression
  322. C
  323. TYPE='MCHAML '
  324. CALL LIROBJ(TYPE,IPF,1,IRET)
  325. IF(IERR .NE. 0) GOTO 9999
  326. C
  327. C**** Lecture du MCHAML 'FACEL' contenant les gamma
  328. C
  329. TYPE='MCHAML '
  330. CALL LIROBJ(TYPE,IGAMF,1,IRET)
  331. IF(IERR .NE. 0) GOTO 9999
  332. C
  333. C**** Si LOGME -> MULTIESPECES
  334. C
  335. IF(LOGME)THEN
  336. C
  337. C******** Lecture du MCHAML 'FACEL' contenant les fractiones massiques
  338. C
  339. TYPE='MCHAML '
  340. CALL LIROBJ(TYPE,IFRMAF,1,IRET)
  341. IF(IERR .NE. 0) GOTO 9999
  342. C
  343. C********** Lecture de la table qui contient le proprieté du gaz
  344. C
  345. TYPE='TABLE '
  346. CALL LIROBJ(TYPE,IPGAZ,1,IRET)
  347. IF(IERR .NE. 0) GOTO 9999
  348. TYPE=' '
  349. CALL ACMO(IPGAZ,'ESPEULE ',TYPE,MLMOTY)
  350. IF(IERR .NE. 0) GOTO 9999
  351. IF(TYPE .NE. 'LISTMOTS')THEN
  352. C
  353. C********** Message d'erreur standard
  354. C -301 0 %m1:40
  355. C
  356. MOTERR(1:40) = 'TABGAS = ??? '
  357. WRITE(IOIMP,*) MOTERR
  358. C
  359. C********** Message d'erreur standard
  360. C 21 2
  361. C Données incompatibles
  362. C
  363. CALL ERREUR(21)
  364. GOTO 9999
  365. ENDIF
  366. C
  367. SEGACT MLMOTY
  368. NESP = MLMOTY.MOTS(/2)
  369. C
  370. ELSE
  371. IFRMAF = 0
  372. NESP=0
  373. ENDIF
  374. C
  375. C**** La list des inconnues
  376. C
  377. TYPE='LISTMOTS'
  378. CALL LIROBJ(TYPE,ILIINC,1,IRET)
  379. IF(IERR .NE. 0) GOTO 9999
  380. MLMOTS = ILIINC
  381. SEGACT MLMOTS
  382. NC = MLMOTS.MOTS(/2)
  383. SEGDES MLMOTS
  384. IF(NC .NE. (IDIM+2+NESP))THEN
  385. MOTERR(1:40) = 'LISTINCO = ???'
  386. WRITE(IOIMP,*) MOTERR
  387. C
  388. C******* Message d'erreur standard
  389. C 21 2
  390. C Données incompatibles
  391. C
  392. CALL ERREUR(21)
  393. GOTO 9999
  394. ENDIF
  395. C
  396. C**** Boundary condition
  397. C
  398. IRET=0
  399. TYPE='MAILLAGE'
  400. CALL LIROBJ(TYPE,IFLIM,0,IRET)
  401. IF(IERR.NE.0)GOTO 9999
  402. IF(IRET .EQ. 0)THEN
  403. MELLIM = 0
  404. ELSE
  405. MELEME=IFLIM
  406. SEGACT MELEME
  407. ICACCA=MELEME.NUM(/2)
  408. IF(ICACCA .EQ. 0)THEN
  409. MELLIM = 0
  410. ELSE
  411. MELLIM = IFLIM
  412. ENDIF
  413. SEGDES MELEME
  414. ENDIF
  415. C
  416. C**** Bas Mach (AUSMPLM, RUSANOLM, ROELM, HLLCLM, AUSMPUP)
  417. C
  418. IF((INDMET .EQ. 9) .OR. (INDMET .EQ. 11) .OR.
  419. & (INDMET .EQ. 13) .OR. (INDMET .EQ. 15) .OR.
  420. & (INDMET .EQ. 16)) THEN
  421. TYPE = 'CHPOINT '
  422. C
  423. C******* Reference speed
  424. C
  425. CALL LIROBJ(TYPE,IUINF,1,IRET)
  426. IF(IERR .NE. 0) GOTO 9999
  427. INDIC = 1
  428. NBCOMP = 1
  429. MOT = 'SCAL'
  430. CALL QUEPOI(IUINF, MELEMC, INDIC, NBCOMP, MOT)
  431. IF(IERR .NE. 0) GOTO 9999
  432. C
  433. C******* Minimal cutoff
  434. C
  435. TYPE = 'CHPOINT '
  436. CALL LIROBJ(TYPE,IUPRI,1,IRET)
  437. IF(IERR .NE. 0) GOTO 9999
  438. INDIC = 1
  439. NBCOMP = 1
  440. MOT = 'SCAL'
  441. CALL QUEPOI(IUPRI, MELEMC, INDIC, NBCOMP, MOT)
  442. IF(IERR .NE. 0) GOTO 9999
  443. C
  444. ELSE
  445. IUINF=0
  446. IUPRI=0
  447. ENDIF
  448. C
  449. C**** Creation des flux aux interfaces
  450. C
  451. TYPE = 'CHPOINT '
  452. CALL KRCHP1(TYPE, MELEMF, ICHFLU, MLMOTS)
  453. C
  454. C**** Calcul des flux et du pas du temps.
  455. C
  456. IF(IDIM .EQ. 2)THEN
  457. CALL KONFL1(LOGME,INDMET,
  458. & IROF,IVITF,IPF,IGAMF,IFRMAF,
  459. & ICHPSU,ICHPDI,
  460. & MELEMC,MELEMF,MELEFE,MELLIM,
  461. & IUINF,IUPRI,
  462. & ICHFLU,
  463. & DT,
  464. & LOGNC,LOGAN,MESERR)
  465. ELSEIF(IDIM .EQ. 3)THEN
  466. CALL KONFL3(LOGME,INDMET,
  467. & IROF,IVITF,IPF,IGAMF,IFRMAF,
  468. & ICHPSU,ICHPDI,
  469. & MELEMC,MELEMF,MELEFE,MELLIM,
  470. & IUINF,IUPRI,
  471. & ICHFLU,
  472. & DT,
  473. & LOGNC,LOGAN,MESERR)
  474. ENDIF
  475. C
  476. IF(LOGAN)THEN
  477. C
  478. C******* Anomalie detectée
  479. C
  480. C
  481. C******* Message d'erreur standard
  482. C -301 0
  483. C %m1:40
  484. C
  485. MOTERR(1:40) = MESERR(1:40)
  486. WRITE(IOIMP,*) MOTERR(1:40)
  487. C
  488. C******* Message d'erreur standard
  489. C 5 3
  490. C Erreur anormale.contactez votre support
  491. C
  492. CALL ERREUR(5)
  493. GOTO 9999
  494. ENDIF
  495. IF(LOGNC)THEN
  496. C
  497. C******* Message d'erreur standard
  498. C -301 0
  499. C %m1:40
  500. C
  501. MOTERR(1:40) = MESERR(1:40)
  502. WRITE(IOIMP,*) MOTERR(1:40)
  503. C
  504. C******* Message d'erreur standard
  505. C 460 2
  506. C Pas de convergence dans les itérations internes
  507. C
  508. CALL ERREUR(460)
  509. GOTO 9999
  510. ENDIF
  511. C
  512. C**** Calcul de residu (IFLU = 2)
  513. C
  514. IF(IFLU .EQ. 2)THEN
  515. TYPE = 'CHPOINT '
  516. CALL KRCHP1(TYPE, MELEMC, ICHRES, MLMOTS)
  517. C
  518. CALL KONRE1(MELEMC,MELEMF,MELEFE,ICHPVO,
  519. & ICHFLU, ICHRES,
  520. & LOGAN,MESERR)
  521. IF(LOGAN)THEN
  522. C
  523. C******* Anomalie detectée
  524. C
  525. C
  526. C******* Message d'erreur standard
  527. C -301 0
  528. C %m1:40
  529. C
  530. MOTERR(1:40) = MESERR(1:40)
  531. WRITE(IOIMP,*) MOTERR(1:40)
  532. C
  533. C******* Message d'erreur standard
  534. C 5 3
  535. C Erreur anormale.contactez votre support
  536. C
  537. CALL ERREUR(5)
  538. GOTO 9999
  539. ENDIF
  540. ELSE
  541. ICHRES = 0
  542. ENDIF
  543. C
  544. C**** Ecriture des resultats
  545. C
  546. CALL ECRREE(DT)
  547. TYPE = 'CHPOINT '
  548. IF(ICHRES .NE. 0) CALL ECROBJ(TYPE,ICHRES)
  549. IF(ICHFLU .NE. 0) CALL ECROBJ(TYPE,ICHFLU)
  550. 9999 CONTINUE
  551. RETURN
  552. END
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  
  565.  
  566.  

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