Télécharger cyltest6.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : cyltest6.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. 'OPTI' 'ECHO' 0 ;
  5. ************************************************************************
  6. * *
  7. * Ecoulement autour d'un cylindre circulaire *
  8. * Résolution des équations de Navier Stokes laminaires instationnaires *
  9. * *
  10. * Pompé sur cyltest.dgibi *
  11. * Ajouts Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) *
  12. * mél : gounand@semt2.smts.cea.fr *
  13. * en date du 20/12/2007 *
  14. * teste une méthode couplée et une méthode de projection *
  15. * algébrique incrémentale : vérification du bilan des forces *
  16. * locales et de la force s'exercant sur le cylindre *
  17. * *
  18. ************************************************************************
  19. 'SAUTER' 2 'LIGNE' ;
  20. 'MESSAGE' ' Execution de cyltest6.dgibi' ;
  21. 'SAUTER' 2 'LIGNE' ;
  22. *
  23. interact = FAUX ;
  24. graph = FAUX ;
  25. *
  26. *****************************************
  27. * *
  28. * Procedure tracvit *
  29. * *
  30. *****************************************
  31. 'DEBPROC' TRACVIT;
  32. 'ARGUMENT' RVX*'TABLE';
  33. rv = rvx . 'EQEX' ;
  34. umax = 2. ;
  35. CHUN = 'VECTEUR' (rv.inco.un) (1. '/' (1. * umax)) UX UY JAUN;
  36. 'TRACER' CHUN mt cmt 'TITRE' 'Champ de vitesse'
  37. 'NCLK' ;
  38. as2 ama1 = 'KOPS' 'MATRIK';
  39. 'RESPRO' as2 ama1;
  40. 'FINP';
  41. *****************************************
  42. * *
  43. * Proc Bloc de pression nul *
  44. * *
  45. *****************************************
  46. 'DEBPROC' BLOCPNUL ;
  47. 'ARGUMENT' rvx*'TABLE' ;
  48. RV = RVX . 'EQEX';
  49. *
  50. 'SI' ('NON' ('EXISTE' rv 'MATPNULL')) ;
  51. *
  52. * Creation d'un bloc de pression nulle
  53. *
  54. 'MESSAGE' ('CHAINE' 'Calcul du bloc de pression') ;
  55. $mt = rvx . 'DOMZ' ;
  56. rt = 'EQEX' 'OPTI' 'EF' 'IMPL' 'CENTREE' 'CENTREP1'
  57. 'ZONE' $mt
  58. 'OPER' 'KMAB' 0. 'INCO' 'UN' 'PN' ;
  59. Rt.'INCO' = RV . 'INCO' ;
  60. smbb matb = kmab rt . '1KMAB' ;
  61. chpv = 'KCHT' $mt 'VECT' 'SOMMET' 'COMP' '1UN' '2UN' (0. 0.) ;
  62. matpres = 'KOPS' 'CMCT' matb chpv matb ;
  63. rv . 'MATPNULL' = matpres ;
  64. 'SINON' ;
  65. matpres =rv . 'MATPNULL' ;
  66. 'FINSI' ;
  67. smbvide matvide = 'KOPS' 'MATRIK' ;
  68. 'FINPROC' smbvide matpres ;
  69. *
  70. * Options globales
  71. *
  72. 'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'ECHO' 0 ;
  73. 'SI' interact ;
  74. 'OPTION' 'TRAC' 'X' ;
  75. 'SINON' ;
  76. 'OPTION' 'TRAC' 'PSC';
  77. 'FINSI' ;
  78. 'OPTION' isov suli;
  79. disv = 'QUAF' ;
  80. DISP = 'CENTREP1';
  81. * iraff : raffinement du maillage
  82. * mproj : méthode de projection ou couplée
  83. * miter : méthode itérative ou directe
  84. iraff = 1 ;
  85. impkres = 0 ;
  86. *
  87. lok = VRAI ;
  88. *
  89. 'REPETER' imet 6 ;
  90. 'SAUTER' 1 'LIGN' ;
  91. iimet = &imet ;
  92. 'SI' ('EGA' iimet 1) ;
  93. mproj = FAUX ;
  94. miter = FAUX ;
  95. 'FINSI' ;
  96. 'SI' ('EGA' iimet 2) ;
  97. mproj = FAUX ;
  98. miter = VRAI ;
  99. 'FINSI' ;
  100. 'SI' ('EGA' iimet 3) ;
  101. mproj = VRAI ; lprec = FAUX ;
  102. miter = FAUX ;
  103. 'FINSI' ;
  104. 'SI' ('EGA' iimet 4) ;
  105. mproj = VRAI ; lprec = FAUX ;
  106. miter = VRAI ;
  107. 'FINSI' ;
  108. 'SI' ('EGA' iimet 5) ;
  109. mproj = VRAI ; lprec = VRAI ;
  110. miter = FAUX ;
  111. 'FINSI' ;
  112. 'SI' ('EGA' iimet 6) ;
  113. mproj = VRAI ; lprec = VRAI ;
  114. miter = VRAI ;
  115. 'FINSI' ;
  116. 'SI' mproj ;
  117. 'SI' lprec ;
  118. 'MESSAGE' 'Methode projection preconditionnee' ;
  119. 'SINON' ;
  120. 'MESSAGE' 'Methode projection' ;
  121. 'FINSI' ;
  122. 'SINON' ;
  123. 'MESSAGE' 'Methode directe' ;
  124. 'FINSI' ;
  125. 'SI' miter ;
  126. 'MESSAGE' 'Solveur itératif' ;
  127. 'SINON' ;
  128. 'MESSAGE' 'Solveur direct' ;
  129. 'FINSI' ;
  130. 'SAUTER' 1 'LIGN' ;
  131. *============================
  132. * Construction du maillage:
  133. *============================
  134. * Construction des points:
  135. nc = 8 ;
  136. nx = 22 ; ny = 12 ;
  137. nc = '*' nc iraff ;
  138. nx = '*' nx iraff ;
  139. ny = '*' ny iraff ;
  140. *
  141. p0 = 0. 0. ;
  142. p1 = 0.5 0. ;
  143. p2 = 0. 0.5 ;
  144. p3 = -0.5 0. ;
  145. p4 = 0. -0.5 ;
  146. *
  147. c1 = 'CERCLE' nc p1 p0 p2 ;
  148. c2 = 'CERCLE' nc p2 p0 p3 ;
  149. c3 = 'CERCLE' nc p3 p0 p4 ;
  150. c4 = 'CERCLE' nc p4 p0 p1 ;
  151. *
  152. ligc = c1 'ET' c2 'ET' c3 'ET' c4 ;
  153. ligc = 'INVERSE' ligc ;
  154. *
  155. *
  156. pA = -6. -6. ;
  157. pB = 15. -6. ;
  158. pC = 15. 6. ;
  159. pD = -6. 6. ;
  160. bas = pA 'DROIT' nx pB ;
  161. dro = pB 'DROIT' ny pC ;
  162. hau = pC 'DROIT' nx pD ;
  163. gau = pD 'DROIT' ny pA ;
  164. *
  165. lige = bas 'ET' dro 'ET' hau 'ET' gau ;
  166. *
  167. cntt = ligc 'ET' lige ;
  168. *
  169. mt = tria cntt ;
  170. xmt ymt = 'COORDONNEE' mt ;
  171. rmt = '**' ('+' ('**' xmt 2) ('**' ymt 2)) 0.5D0 ;
  172. *rmt = '*' rmt 3. ;
  173. rmt = '/' rmt 3. ;
  174. rmt = '/' rmt iraff ;
  175. mtr = 'RAFT' rmt mt ;
  176. 'SI' graph ;
  177. 'TRACER' mtr 'TITR' ('CHAINE' 'nbl=' ('NBEL' mtr)) ;
  178. 'FINSI' ;
  179. *'OPTION' 'DONN' 5 ;
  180. cercle = ligc ;
  181. * Création du domaine total:
  182. mt = mtr;
  183. cmt = 'CONTOUR' mt;
  184. *==============================================
  185. * Changement des éléments du maillage en QUAF:
  186. *==============================================
  187.  
  188. _mt = 'CHANGER' mt 'QUAF';
  189. _cmt = 'CHANGER' cmt 'QUAF';
  190. _cercle = 'CHANGER' cercle 'QUAF';
  191. _entree = 'CHANGER' gau 'QUAF';
  192. _sortie = 'CHANGER' dro 'QUAF';
  193. _hau = 'CHANGER' hau 'QUAF';
  194. _bas = 'CHANGER' bas 'QUAF';
  195. 'ELIM' 1.E-5 (_mt et _cercle et _entree 'ET' _sortie
  196. 'ET' _cmt 'ET' _hau 'ET' _bas);
  197.  
  198. NN = 'NBEL' mt;
  199. 'MESSAGE' 'Nombre d éléments du maillage' NN;
  200.  
  201. *=======================================
  202. * Formulation du domaine Navier Stokes:
  203. *=======================================
  204.  
  205. $mt = 'MODE' _mt 'NAVIER_STOKES' disv ;
  206. $cmt = 'MODE' _cmt 'NAVIER_STOKES' disv ;
  207. $cercle = 'MODE' _cercle 'NAVIER_STOKES' disv ;
  208. $entree = 'MODE' _entree 'NAVIER_STOKES' disv ;
  209. $sortie = 'MODE' _sortie 'NAVIER_STOKES' disv ;
  210. $hau = 'MODE' _hau 'NAVIER_STOKES' disv ;
  211. $bas = 'MODE' _bas 'NAVIER_STOKES' disv ;
  212. *
  213. mt = 'DOMA' $mt 'MAILLAGE' ;
  214. cmt = 'DOMA' $cmt 'MAILLAGE' ;
  215. cercle = 'DOMA' $cercle 'MAILLAGE' ;
  216. entree = 'DOMA' $entree 'MAILLAGE' ;
  217. sortie = 'DOMA' $sortie 'MAILLAGE' ;
  218. hau = 'DOMA' $hau 'MAILLAGE' ;
  219. bas = 'DOMA' $bas 'MAILLAGE' ;
  220.  
  221. *====================================
  222. * Création des tables de résolution:
  223. *====================================
  224.  
  225. dt = .5D0 ;
  226. *rey = 30. ;
  227. rey = 100.;
  228. nu = 1./rey;
  229. arelax = 1.;
  230. *
  231. niter = 6 ;
  232. nitime = 3 ;
  233. *
  234. RV = 'EQEX' 'OMEGA' 1.0 'NITER' niter 'ITMA' 1 'FIDT' 1 ;
  235. *
  236. RV = 'EQEX' RV
  237. 'OPTI' 'EF' 'CENTREE' 'IMPL'
  238. * 'ZONE' $mt 'OPER' 'NS' 'INV_RE' 'INCO' 'UN'
  239. 'ZONE' $mt 'OPER' 'KONV' 1. 'UN' 'INV_RE' 'INCO' 'UN'
  240. 'OPTI' 'EF' 'CENTREE' 'IMPL'
  241. 'ZONE' $mt 'OPER' 'LAPN' 'INV_RE' 'INCO' 'UN'
  242. 'OPTI' 'EF' 'CENTREE' 'IMPL'
  243. 'ZONE' $mt 'OPER' 'DFDT' 1. 'UNM' 'DT' 'INCO' 'UN'
  244. 'OPTI' 'EF' 'IMPL' 'CENTREE' DISP
  245. 'ZONE' $mt
  246. 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PN' ;
  247.  
  248. 'SI' graph ;
  249. rv = 'EQEX' rv
  250. 'ZONE' $mt 'OPER' TRACVIT ;
  251. 'FINSI' ;
  252.  
  253. * Implantation des conditions aux limites:
  254.  
  255. RV = 'EQEX' RV
  256. 'CLIM' 'UN' 'UIMP' cercle 0. 'UN' 'VIMP' cercle 0.
  257. 'UN' 'UIMP' entree 1. 'UN' 'VIMP' entree 0.
  258. 'UN' 'VIMP' hau 0.
  259. 'UN' 'VIMP' bas 0. ;
  260. *
  261. * Choix de la méthode de résolution (couplée, projection)
  262. * et des solveurs (direct, itératif)
  263. *
  264. 'SI' ('NON' mproj) ;
  265. * Solveur pour le système vitesse-pression
  266. rvm = rv . 'METHINV' ;
  267. rvm . 'SCALING' = 1 ;
  268. rvm . 'IMPINV' = impkres ;
  269. 'SI' ('NON' miter) ;
  270. rvm . 'TYPINV' = 1 ;
  271. 'SINON' ;
  272. * On rajoute le bloc de pression nul
  273. rv = 'EQEX' rv 'ZONE' $mt 'OPER' BLOCPNUL ;
  274. *
  275. rvm . 'TYPINV' = 4 ;
  276. rvm . 'LBCG' = 2 ;
  277. rvm . 'PRECOND' = 5 ;
  278. rvm . 'ILUTLFIL' = 3.D0 ;
  279. 'FINSI' ;
  280. 'SINON' ;
  281. rv . 'GPROJ' = 'TABLE' ;
  282. rv . 'GPROJ' . 'NOMVIT' = 'UN' ;
  283. rv . 'GPROJ' . 'NOMPRES' = 'PN' ;
  284. 'SI' ('NON' lprec) ;
  285. rv . 'GPROJ' . 'NOPREC' = VRAI ;
  286. 'FINSI' ;
  287. * Solveur pour les composantes de la vitesse
  288. rvm = rv . 'METHINV' ;
  289. rvm . 'SCALING' = 1 ;
  290. rvm . 'IMPINV' = impkres ;
  291. 'SI' ('NON' miter) ;
  292. rvm . 'TYPINV' = 1 ;
  293. 'SINON' ;
  294. rvm . 'TYPINV' = 4 ;
  295. rvm . 'LBCG' = 2 ;
  296. rvm . 'PRECOND' = 3 ;
  297. 'FINSI' ;
  298. * Solveur pour l'équation de pression
  299. rv . 'GPROJ'. 'METHINV' = 'TABLE' 'METHINV' ;
  300. rvgm = rv . 'GPROJ' . 'METHINV' ;
  301. rvgm . 'SCALING' = 1 ;
  302. rvgm . 'IMPINV' = impkres ;
  303. 'SI' ('NON' miter) ;
  304. rvgm . 'TYPINV' = 1 ;
  305. 'SINON' ;
  306. rvgm . 'TYPINV' = 4 ;
  307. rvgm . 'LBCG' = 2 ;
  308. rvgm . 'PRECOND' = 3 ;
  309. * rvgm . 'PRECOND' = 5 ;
  310. * rvgm . 'ILUTLFIL' = 1.D0 ;
  311. 'FINSI' ;
  312. 'FINSI' ;
  313. * Création de la table des inconnues et initialisation:
  314. RV.INCO = 'TABLE' INCO;
  315. *
  316. RV.INCO.'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (1.E-3 1.E-3);
  317. RV.INCO.'UNM' = 'KCHT' $mt 'VECT' 'SOMMET' (1.E-3 1.E-3);
  318. rv . 'INCO' . 'G' = (0. -9.8) ;
  319. RV.INCO.'DT' = dt ;
  320. RV.INCO.'INV_RE'= nu;
  321. RV.INCO.'PN' = 'KCHT' $mt 'SCAL' DISP 0.;
  322. *
  323. * Table avec les résultats en vitesse
  324. *
  325. res = 'TABLE' ;
  326. itres = 0 ;
  327. res . itres = 'COPIER' (RV.INCO.'UN') ;
  328. ***************************************************
  329. *
  330. * Boucle en temps faite main
  331. *
  332. 'TEMPS' 'ZERO' ;
  333. 'REPETER' iitime nitime ;
  334. rv . 'INCO' . 'UNM' = 'COPIER' (res . itres) ;
  335. EXEC RV;
  336. itres = '+' itres 1 ;
  337. res . itres = 'COPIER' (RV.INCO.'UN') ;
  338. 'FIN' iitime ;
  339. TABTPS = TEMP 'NOEC';
  340. tcpu = TABTPS.'TEMPS_CPU'.'INITIAL' ;
  341. *
  342. ***************************************************
  343. lastime = '-' ('DIME' res) 1 ;
  344. vit = res . lastime ;
  345. vitm1 = res . ('-' lastime 1) ;
  346. umax = 2. ;
  347. *
  348. *
  349. * Créer les matrices
  350. *
  351. bidon matkon = 'KONV' (rv . '1KONV') ;
  352. bidon matlap = 'LAPN' (rv . '2LAPN') ;
  353. chdfdt matdfdt = 'DFDT' (rv . '3DFDT') ;
  354. *
  355. rt = 'EQEX' 'ZONE' $mt
  356. 'OPTI' 'EF' 'IMPL' 'CENTREE' DISP
  357. 'OPER' 'KMBT' -1. 'INCO' 'PN' 'UN'
  358. ;
  359. rt.'INCO' = rv . 'INCO' ;
  360. bidon matgrp = 'KMBT' (rt . '1KMBT') ;
  361. *
  362. pre = rv . 'INCO' . 'PN' ;
  363. nomavv = 'MOTS' 'UX' 'UY' ;
  364. nomapv = 'MOTS' '1UN' '2UN' ;
  365. nomavp = 'MOTS' 'SCAL' ;
  366. nomapp = 'MOTS' 'PN' ;
  367. vit = 'NOMC' nomavv nomapv vit ;
  368. vitm1 = 'NOMC' nomavv nomapv vitm1 ;
  369. pre = 'NOMC' nomavp nomapp pre ;
  370. forlap = 'KOPS' matlap '*' vit ;
  371. forkon = 'KOPS' matkon '*' vit ;
  372. forgrp = 'KOPS' matgrp '*' pre ;
  373. fordfdt = 'KOPS' matdfdt '*' ('-' vit vitm1) ;
  374. *
  375. forlap = ('NOMC' nomapv nomavv forlap) ;
  376. forkon = ('NOMC' nomapv nomavv forkon) ;
  377. forgrp = ('NOMC' nomapv nomavv forgrp) ;
  378. fordfdt = ('NOMC' nomapv nomavv fordfdt) ;
  379. *
  380. fort = forlap '+' forkon '+' forgrp '+' fordfdt ;
  381. *
  382. fmax = 'MAXIMUM' fort 'ABS' ;
  383. echf = '/' 0.5 fmax ;
  384. *
  385. vlap = 'VECTEUR' forlap echf UX UY VERT ;
  386. vkon = 'VECTEUR' forkon echf UX UY ROUG ;
  387. vgrp = 'VECTEUR' forgrp echf UX UY JAUN ;
  388. vdfdt = 'VECTEUR' fordfdt echf UX UY TURQ ;
  389. vf = 'VECTEUR' fort echf UX UY BLAN ;
  390. vtot = vlap 'ET' vkon 'ET' vgrp 'ET' vdfdt 'ET' vf ;
  391. *
  392. 'SI' graph ;
  393. tit = 'CHAINE' 'vert=lap;roug=kon;jaun=grp;turq=dfdt;blan=reac;fmax='
  394. fmax ;
  395. 'TRACER' vtot _mt ('CONTOUR' _mt) 'TITRE' tit ;
  396. 'FINSI' ;
  397. *
  398. pmt = 'CHANGER' mt 'POI1' ;
  399. pcmt = 'CHANGER' cmt 'POI1' ;
  400. pint = 'DIFF' pmt pcmt ;
  401. fortint = 'REDU' fort pint ;
  402. fmaxint = 'MAXIMUM' fortint 'ABS' ;
  403. 'MESSAGE' ('CHAINE' 'fmax=' fmax) ;
  404. 'MESSAGE' ('CHAINE' 'fmaxint=' fmaxint) ;
  405. errbf = '*' ('/' fmaxint fmax) 100. ;
  406. 'MESSAGE' ('CHAINE' 'erreur bilan forces =' errbf ' %') ;
  407. forcer = 'REDU' fort cercle ;
  408. rfor = 'RESULT' forcer ;
  409. 'MESSAGE' 'Resultante des forces sur le cylindre' ;
  410. rfx = 'MAXIMUM' ('EXCO' 'UX' rfor) ;
  411. rfy = 'MAXIMUM' ('EXCO' 'UY' rfor) ;
  412. 'MESSAGE' ('CHAINE' ' en x =' rfx) ;
  413. 'MESSAGE' ('CHAINE' ' en y =' rfy) ;
  414. 'MESSAGE' ('CHAINE' 'tcpu=' tcpu) ;
  415. *
  416. * Tests
  417. * Erreur sur le bilan des forces au 16/04/2010
  418. * 2.e-3 % pour la méthode directe
  419. * 7 % pour la méthode projection algébrique
  420. * 1 % pour la méthode double projection algébrique
  421. * Valeur de la force de trainée (fx)
  422. * -6.94e-1 pour la méthode directe
  423. * -6.76e-1 pour la méthode projection algébrique
  424. * -6.93e-1 pour la méthode double projection algébrique
  425. *
  426. 'SI' mproj ;
  427. * vrbf = 8. ;
  428. vrbf = 2. ;
  429. 'SINON' ;
  430. vrbf = 5.D-3 ;
  431. 'FINSI' ;
  432. *vrft = -6.9e-1 ;
  433. *tol = 0.3e-1 ;
  434. vrft = -6.94e-1 ;
  435. tol = 0.02e-1 ;
  436. *
  437. tst = '<' errbf vrbf ;
  438. 'SI' ('NON' tst) ;
  439. cherr = 'CHAINE' '!!! Erreur bilan force > ' vrbf ' %' ;
  440. 'MESSAGE' cherr ;
  441. 'FINSI' ;
  442. lok = lok 'ET' tst ;
  443. *
  444. tst = '<' ('ABS' ('-' rfx vrft)) tol ;
  445. 'SI' ('NON' tst) ;
  446. cherr = 'CHAINE' '!!! Erreur force trainee <> ' vrft ' a ' tol
  447. ' pres' ;
  448. 'MESSAGE' cherr ;
  449. 'FINSI' ;
  450. lok = lok 'ET' tst ;
  451. 'FIN' imet ;
  452. 'SAUTER' 2 'LIGNE' ;
  453. 'SI' lok ;
  454. 'MESSAGE' 'Tout sest bien passe' ;
  455. 'SINON' ;
  456. 'MESSAGE' 'Il y a eu des erreurs' ;
  457. 'FINSI' ;
  458. 'SAUTER' 2 'LIGNE' ;
  459. *
  460. 'SI' interact ;
  461. 'OPTION' 'DONN' 5 'ECHO' 1 ;
  462. 'FINSI' ;
  463. 'SI' ('NON' lok) ;
  464. 'ERREUR' 5 ;
  465. 'FINSI' ;
  466. * 'FIN' du jeu de données.
  467. fin;
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  

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