Télécharger la plaquette Cast3M


fr en

Télécharger allee.dgibi

Retour à la liste

Numérotation des lignes : 

  1. *****************************************************
  2. ************************************************************************
  3. ************************************************************************
  4. * fichier :  allee.dgibi                            *
  5. ** modifie le 15/06/2014 passage EQPR -> EQEX       *
  6. *****************************************************
  7.  
  8.  'OPTION' echo 0;
  9. ************************************************************************
  10. * NOM         : allee.dgibi
  11. * DESCRIPTION :
  12. *
  13. * Ce jeu de données est dédié à l'écoulement d'un fluide
  14. * incompressible autour d'un cylindre de section
  15. * circulaire. L'écoulement est supposé laminaire. Le but est de
  16. * retrouver les lâchers tourbillonnaires de von Karman et les
  17. * fréquences de ces derniers.
  18. *
  19. * Références :
  20. *
  21. * U. Frisch. Turbulence. Cambridge Universuty Press, 1995.
  22. * H. Schlichting. Boundary-layer theory. McGraw-Hill book company,
  23. * 1968.
  24. *
  25. * LANGAGE     : GIBIANE
  26. * AUTEUR      : F. PIGEONNEAU (CEA/DRN/DMT/SEMT/LTMF)
  27. *               e-mail : frankp@semt2.smts.cea.fr
  28. ************************************************************************
  29. * APPELES          : 
  30. ************************************************************************
  31. * SYNTAXE GIBIANE    : Uniquement.
  32. ************************************************************************
  33. *             Description des paramètres et variables locaux:
  34. ************************************************************************
  35. *
  36. *
  37. ************************************************************************
  38. * VERSION    : v1, 15/12/98, version initiale.
  39. * HISTORIQUE : v1, 15/12/98, création.
  40. * HISTORIQUE : v2, 15/06/2014 passage EQPR -> EQEX
  41. * HISTORIQUE :
  42. * HISTORIQUE :
  43. ************************************************************************
  44. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  45. * en cas de modification de cette procédure afin de faciliter
  46. * la maintenance !
  47. ************************************************************************
  48.  
  49. DISCR  = 'MACRO'  ;
  50. KPRESS = 'CENTRE' ;
  51. BETAP  = 1.       ;
  52.  
  53. GRAPH=FAUX;
  54. COMPLET = FAUX ;
  55. 'SI' (COMPLET) ;
  56. *  Temps final de calcul :
  57.  
  58.    tfin=1. ;
  59.  
  60. *  Définitions des nombres d'éléments:
  61.  
  62.    nx1 = 10 ;
  63.    nx2 = 20 ;
  64.    nx3 = 40 ;
  65.  
  66. *  Définitions des densites:
  67.  
  68.    dbeg = 0.1 ;
  69.    dend = 0.5 ;
  70.  
  71. 'SINON' ;
  72. *  Temps final de calcul :
  73.  
  74.    tfin = 0.1 ;
  75.  
  76. *  Définitions des nombres d'éléments:
  77.  
  78.    nx1 = 2 ;
  79.    nx2 = 4 ;
  80.    nx3 = 5 ;
  81.  
  82. *  Définitions des densites:
  83.  
  84.    dbeg = 0.4 ;
  85.    dend = 0.5 ;
  86.  
  87. 'FINSI' ;
  88.  
  89. *=====================
  90. * Domaine de calcul 2D:
  91. *=====================
  92.  
  93. 'OPTION' 'DIME' 2 'ELEM' 'QUA8' ;
  94.  
  95. *****************************************
  96. *                                       *
  97. *           Procedure calcul            *
  98. *                                       *
  99. *****************************************
  100.  
  101. 'DEBP' CALCUL ;
  102. *  Calcul de la convergence.
  103. *  Evolution de la vitesse en un point quelconque.
  104.  
  105. 'ARGUMENT' RVX*'TABLE' ;
  106. RV = RVX.'EQEX' ;
  107.  
  108. dd = RV.PASDETPS.'NUPASDT' ;
  109. nn = dd/2 ;
  110.  
  111. 'SI' ((dd-(2*nn))'EGA'0) ;
  112. *  Champs de vitesse aux instants n et n-1:
  113.    un = RV.INCO.'UN' ;
  114.    unm1 = RV.INCO.'UNM1' ;
  115.  
  116. *  Extraction des composantes de la vitesse au pas de temps n et n-1:
  117.    unx   = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' un) ;
  118.    unm1x = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' unm1) ;
  119.    uny   = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' un) ;
  120.    unm1y = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' unm1) ;
  121.  
  122. *  Sauvegarde de la vitesse en un point fixe du domaine à l'instant n:
  123.    z = 'REDU' unx RV.INCO.'PoiObs' ;
  124.    zz = 'PROG' ('EXTR' z 'SCAL' RV.INCO.'PoiObs') ;
  125.    RV.INCO.'Uvit'= RV.INCO.'Uvit' et zz ;
  126.  
  127. *  Calcul de la convergence:
  128.    errx = 'KOPS' unx '-' unm1x ;
  129.    erry = 'KOPS' uny '-' unm1y ;
  130.    elix = 'MAXIMUM' errx 'ABS' ;
  131.    eliy = 'MAXIMUM' erry 'ABS' ;
  132.    elix = 'LOG'(elix + 1.0E-10)/('LOG' 10.) ;
  133.    eliy = 'LOG'(eliy + 1.0E-10)/('LOG' 10.) ;
  134.    erx  = 'PROG' elix ;
  135.    ery  = 'PROG' eliy ;
  136.    RV.INCO.'ERX'  = RV.INCO.'ERX' et erx ;
  137.    RV.INCO.'ERY'  = RV.INCO.'ERY' et ery ;
  138.    it             = 'PROG' RV.PASDETPS.'NUPASDT' ;
  139.    temps          = 'PROG' RV.PASDETPS.'TPS' ;
  140.    RV.INCO.'IT'   = RV.INCO.'IT' et it ;
  141.    RV.INCO.'tps'  = RV.INCO.'tps' et temps ;
  142.    RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET'
  143.                     (RV.INCO.'UN') ;
  144. 'FINS' ;
  145.  
  146. as2 ama1 = 'KOPS' 'MATRIK' ;
  147. 'RESPRO' as2 ama1 ;
  148.  
  149. 'FINP' ;
  150.  
  151.  
  152. *============================
  153. *  Construction du maillage:
  154. *============================
  155.  
  156. *  Construction des points:
  157.  
  158. p0 = 0. 0. ;
  159. p1 = -6. 0. ;
  160. p2 = -0.5 0. ;
  161. p3 = 0. 0.5 ;
  162. p4 = 0.5 0. ;
  163. p5 = 15. 0. ;
  164. p6 = 15. 6. ;
  165. p7 = 3. 6. ;
  166. p8 = -3. 6. ;
  167. p9 = -6. 6. ;
  168.  
  169. *  construction des points de la partie supérieure:
  170.  
  171. p1p2 = p1 d nx1 p2 ;
  172. p2p3 = 'CERCLE' (nx2/2) p2 p0 p3 ;
  173. p3p4 = 'CERCLE' (nx2/2) p3 p0 p4 ;
  174. p4p5 = p4 d nx3 p5 ;
  175. p5p6 = p5 d dini dbeg dfin dend p6 ;
  176. p6p7 = p6 d nx3 p7 ;
  177. p7p8 = p7 d nx2 p8 ;
  178. p8p9 = p8 d nx1 p9 ;
  179. p9p1 = p9 d dini dend dfin dbeg p1 ;
  180.  
  181. *  Création du contour:
  182. peri1 = p1p2 et p2p3 et p3p4 et p4p5 ;
  183. peri2 = p6p7 et p7p8 et p8p9 ;
  184. s1 = 'DALL' peri1 p5p6 peri2 p9p1 'PLAN' ;
  185.  
  186. *  construction des points de la partie inférieure:
  187. p10 = -6. -6. ;
  188. p11 = -3. -6. ;
  189. p12 = 3. -6. ;
  190. p13 = 15. -6. ;
  191. p14 = 0. -0.5 ;
  192. p1p10 = p1 d dini dbeg dfin dend p10 ;
  193. p10p11 = p10 d nx1 p11 ;
  194. p11p12 = p11 d nx2 p12 ;
  195. p12p13 = p12 d nx3 p13 ;
  196. p13p5 = p13 d dini dend dfin dbeg p5 ;
  197. p5p4 = 'INVERSE' p4p5 ;
  198. p4p14 = 'CERCLE' (nx2/2) p4 p0 p14 ;
  199. p14p2 = 'CERCLE' (nx2/2) p14 p0 p2 ;
  200. p2p1 = 'INVERSE' p1p2 ;
  201.  
  202. *  Création du contour:
  203. peri3 = p10p11 et p11p12 et p12p13 ;
  204. peri4 = p5p4 et p4p14 et p14p2 et p2p1 ;
  205. s2 = 'DALL' peri3 p13p5 peri4 p1p10 'PLAN' ;
  206. * Création du domaine total:
  207.  
  208. 'ELIMINER' 1.E-3 s1 s2 ;
  209. domtot = s1 et s2 ;
  210.  
  211. *==============================================
  212. *  Changement des éléments du maillage en QUAF:
  213. *==============================================
  214.  
  215. £domtot = 'CHANGER' domtot 'QUAF' ;
  216. £cercle = 'CHANGER' (p2p3 et p3p4 et p4p14 et p14p2) 'QUAF' ;
  217. £entree = 'CHANGER' (peri2 et p9p1 et p1p10 et peri3) 'QUAF' ;
  218. 'ELIM' 1.E-3 (£domtot et £cercle et £entree) ;
  219.  
  220. *=======================================
  221. *  Formulation du domaine Navier Stokes:
  222. *=======================================
  223.  
  224. $domtot = 'MODE' £domtot 'NAVIER_STOKES' 'MACRO' ;
  225. mdomtot = 'DOMA' $domtot maillage ;
  226.  
  227. *====================================
  228. *  Création des tables de résolution:
  229. *====================================
  230.  
  231. maxit = 10000000 ;
  232. rey = 20. ;
  233. nu = 1./rey ;
  234. arelax = 0.7 ;
  235.  
  236. RV = 'EQEX' $domtot 'ITMA' maxit 'ALFA' arelax
  237.      'ZONE' $domtot
  238.      'OPER' CALCUL ;
  239.  
  240. RV = 'EQEX' RV
  241.      'OPTI' 'EFM1' 'CENTREE' 'EXPL'
  242.      'ZONE' $domtot
  243.      'OPER' 'NS' 'INV_RE' 'INCO' 'UN' ;
  244.  
  245. RV = 'EQEX' RV
  246.      'OPTI' 'EFM1' 'CENTREE'
  247.      'ZONE' $domtot
  248.      'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN' ;
  249.  
  250. RVP = EQEX 'OPTI' 'EF' KPRESS
  251.  'ZONE' $domtot  OPER  KBBT  -1. betap  INCO 'UN' 'PRES'
  252. ;
  253.  
  254.     rvp.'METHINV'.TYPINV=1 ;
  255.     rvp.'METHINV'.IMPINV=0 ;
  256.     rvp.'METHINV'.NITMAX=300;
  257.     rvp.'METHINV'.PRECOND=3 ;
  258.     rvp.'METHINV'.RESID  =1.e-8 ;
  259.     rvp.'METHINV' . 'FCPRECT'=100 ;
  260.     rvp.'METHINV' . 'FCPRECI'=100 ;
  261.  
  262.   RV.'PROJ'     =RVP ;
  263.  
  264. *  Implantation des conditions aux limites:
  265.  
  266. RV = 'EQEX' RV
  267.      'CLIM' 'UN' 'UIMP' £cercle 0. 'UN' 'VIMP' £cercle 0.
  268.             'UN' 'UIMP' £entree 1. 'UN' 'VIMP' £entree 0. ;
  269.  
  270. *  Création de la table des inconnues et initialisation:
  271.  
  272. RV.INCO = 'TABLE' INCO ;
  273. RV.INCO.'UN' = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3) ;
  274. RV.INCO.'PRES' = 'KCHT' $domtot 'SCAL' KPRESS 0. ;
  275. RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3) ;
  276. RV.INCO.'INV_RE'= nu ;
  277. RV.TFINAL = tfin ;
  278.  
  279. *  Détermination d'un point d'observation:
  280.  
  281. obs = mdomtot 'POINT' 'PROC' (0.8 0.5) ;
  282. RV.INCO.'PoiObs'= obs ;
  283.  
  284. * Initialisation des inconnues complémentaires:
  285.  
  286. RV.INCO.'Uvit'= 'PROG' 0. ;
  287. RV.INCO.'ERX' = 'PROG' 0. ;
  288. RV.INCO.'ERY' = 'PROG' 0. ;
  289. RV.INCO.'IT' =  'PROG' 0 ;
  290. RV.INCO.'tps' = 'PROG' 0. ;
  291. RV.INCO.'entree'= £entree ;
  292. RV.INCO.'cercle'= £cercle ;
  293.  
  294. *========================
  295. * Résolution du problème:
  296. *========================
  297.  
  298. EXEC RV ;
  299.  
  300. * Récupération de la vitesse au point PoiObs :
  301.  
  302. unx = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' RV.'INCO'.'UN') ;
  303. z = 'REDU' unx RV.INCO.'PoiObs' ;
  304. Uvit ='EXTR' z 'SCAL' RV.INCO.'PoiObs' ;
  305. ecart = 'ABS' (Uvit - 0.98310);
  306.  
  307. mess 'Uvit=' uvit ;
  308.  
  309. Si GRAPH;
  310. un=rv.inco.'UN';
  311. ung=vect un 0.2 ux uy jaune ;
  312. trace ung mdomtot;
  313. Finsi;
  314.  
  315. 'SI' ('NON' COMPLET) ;
  316.    'SI' (ecart > 1.E-4) ;
  317.       'ERREUR' 5 ;
  318.    'FINSI' ;
  319. 'FINSI' ;
  320.  
  321. * 'FIN' du jeu de données.
  322. 'FIN' ;
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  

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