Télécharger la plaquette Cast3M


fr en

Télécharger cyltest.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  cyltest.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. *                                                                      *
  6. *             Ecoulement autour d'un cylindre circulaire               *
  7. * Résolution des équations de Navier Stokes laminaires instationnaires *
  8. *                                                                      *
  9. ************************************************************************
  10.  
  11.  
  12. *'OPTION' 'TRAC' 'PSC';
  13.  'OPTION' 'TRAC' 'X'  ;
  14. 'OPTION' isov suli;
  15.  
  16.  GRAPH = VRAI ;
  17.  GRAPH = FAUX ;
  18.  
  19. COMPLET = FAUX;
  20. *COMPLET = VRAI;
  21. 'SI' (COMPLET);
  22. *  Temps final de calcul :
  23.  
  24.    tfin = 2.;
  25.    maxit  = 2000;
  26.    DT0    = 0.05 ;
  27.  
  28. *  Définitions des nombres d'éléments:
  29.  
  30.    nx1 = 10;
  31.    nx2 = 12;
  32.    nx3 = 25;
  33.  
  34. *  Définitions des densites:
  35.  
  36.    dbeg = 0.1;
  37.    dend = 0.5;
  38.  
  39. 'SINON';
  40.  
  41. *  Temps final de calcul :
  42.  
  43.    tfin = 2.;
  44.    maxit  = 200;
  45.    DT0    = 0.5  ;
  46.  
  47. *  Définitions des nombres d'éléments:
  48.  
  49.    nx1 = 6;
  50.    nx1 = 3;
  51.    nx2 = 4;
  52.    nx3 = 14;
  53.    nx3 = 7 ;
  54.  
  55. *  Définitions des densites:
  56.  
  57.    dbeg = 0.4;
  58.    dend = 0.5;
  59.    dend = 1. ;
  60.  
  61. 'FINSI';
  62.  
  63. *=====================
  64. * Domaine de calcul 2D:
  65. *=====================
  66.  
  67. 'OPTION' 'DIME' 2 'ELEM' 'QUA8';
  68.  
  69. DISP  = 'CENTRE';
  70.  
  71. *****************************************
  72. *                                       *
  73. *           Procedure calcul            *
  74. *                                       *
  75. *****************************************
  76.  
  77. 'DEBP' CALCUL;
  78. *  Calcul de la convergence.
  79. *  Evolution de la vitesse en un point quelconque.
  80.  
  81. 'ARGUMENT' RVX*'TABLE';
  82.  
  83. dd = RV.PASDETPS.'NUPASDT';
  84. nn = dd/2;
  85. nf=2 ;
  86. nfg=50 ;
  87. nn = dd/nf;
  88.  
  89. 'SI' ((dd-(nf*nn)) 'EGA' 0);
  90.  
  91. res = 1.E-6;
  92.  
  93. *  Champs de vitesse aux instants n et n-1:
  94.    un   = RV.INCO.'UN';
  95.    unm1 = RV.INCO.'UNM1';
  96.  
  97. *  Extraction des composantes de la vitesse aux pas de temps n et n-1:
  98.    unx   = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' un);
  99.    unm1x = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' unm1);
  100.    uny   = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' un);
  101.    unm1y = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' unm1);
  102.  
  103. *  Sauvegarde de la vitesse en un point fixe du domaine à l'instant n:
  104.    z  = 'REDU' unx RV.INCO.'PoiObs';
  105.    zz = 'PROG' ('EXTR' z 'SCAL' RV.INCO.'PoiObs');
  106.    RV.INCO.'Uvit'= RV.INCO.'Uvit' 'ET' zz;
  107.  
  108. *  Calcul de la convergence:
  109.    errx = 'KOPS' unx '-' unm1x;
  110.    erry = 'KOPS' uny '-' unm1y;
  111.    elix = 'MAXIMUM' errx 'ABS';
  112.    eliy = 'MAXIMUM' erry 'ABS';
  113.    elix = 'LOG'(elix + 1.0E-10)/('LOG' 10.);
  114.    eliy = 'LOG'(eliy + 1.0E-10)/('LOG' 10.);
  115.    erx  = 'PROG' elix;
  116.    ery  = 'PROG' eliy;
  117.    RV.INCO.'ERX' = RV.INCO.'ERX' et erx;
  118.    RV.INCO.'ERY' = RV.INCO.'ERY' et ery;
  119.    it   = 'PROG' RV.PASDETPS.'NUPASDT';
  120.    temps = 'PROG' RV.PASDETPS.'TPS';
  121.    RV.INCO.'IT'   = RV.INCO.'IT' et it;
  122.    RV.INCO.'tps'  = RV.INCO.'tps' et temps;
  123.    RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET'
  124.                     (RV.INCO.'UN');
  125.  
  126.    PRES = RV.INCO.'PRESSION';
  127.    PRESNM1 = RV.INCO.'PRESNM1';
  128.  
  129. 'SI' ((dd-(nfg*nn)) 'EGA' 0);
  130. *  Tracé du champ de pression
  131.    PN = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE'
  132.                 ('NOMC' 'SCAL' (RV.INCO.'PRESSION')));
  133.    'TRACER' PN domtot 'TITRE' 'Champ de pression';
  134.  
  135. *  Tracé du champ de pression lambda
  136.    PT = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE'
  137.                     ('NOMC' 'SCAL' (RV.INCO.'PTILDE')));
  138. *   'TRACER' PT domtot 'TITRE' 'Champ de pression lambda';
  139.  
  140. *  Tracé de la pression sur le cylindre
  141.    presc = 'REDU' PN £cercle;
  142.    precer = 'EVOL' (rouge) 'CHPO' presc cercle ;
  143. *   'DESSIN' precer 'TITRE' 'Valeurs de pression sur le cercle';
  144.  
  145. *  Tracé du champ de vitesse
  146.    VX = 'EXCO' 'UX' (RV.INCO.'UN');
  147.    VY = 'EXCO' 'UY' (RV.INCO.'UN');
  148.    modu = ((VX**2) '+' (VY**2))**(0.5);
  149.    umax = 'MAXIMUM' modu;
  150.  
  151.    CHUN = 'VECTEUR' (RV.INCO.'UN') (1. '/' (1. * umax)) UX UY BLEU;
  152.    'TRACER' CHUN £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de vitesse';
  153.  
  154. *  Tracé de la vitesse au point d'observation
  155.    evv = 'EVOL' 'MANU' 'Temps' (RV.INCO.'tps') 'Valeur de U'
  156.                         (RV.INCO.'Uvit');
  157. *   'DESSIN' evv 'TITRE' 'Evolution de la vitesse au point (0.8 0.5)';
  158.  
  159. *  Tracé de la convergence
  160.    evconv = 'EVOL' 'MANU' 'Nombre d itérations'
  161.               (RV.INCO.'IT') 'log(err)' (RV.INCO.'ERY');
  162. *   'DESSIN' evconv 'TITRE' 'Convergence sur la vitesse';
  163.  
  164. 'FINSI';
  165. *  Convergence en pression : norme l2
  166. *  mess 'Convergence en pression : norme l2';
  167.    'SI' ('EGA' (RV.PASDETPS.'NUPASDT') 1);
  168.         errpl2 = 1.;
  169.    'SINON';
  170.  
  171.         DPDT   =  PRES - PRESNM1;
  172.         DP2    = (nomc dpdt 'SCAL') * (nomc dpdt 'SCAL');
  173. *        DP2    = 'KOPS' DPDT 'PSCA' DPDT;
  174.  
  175.         DP2E   = 'SOMT' DP2;
  176.         errpl2 = DP2E '/' volt;
  177.         errpl2 = errpl2**0.5;
  178.         'MESSAGE' 'convergence sur la pression' errpl2;
  179.    'FINSI';
  180.    errp = 'PROG' errpl2;
  181.    RV.INCO.'errp'=  RV.INCO.'errp' 'ET' errp;
  182.    RV.INCO.'PRESNM1' = 'COPIER' RV.INCO.'PRESSION';
  183.  
  184.    'SI' (errpl2 < res);
  185.      'MESSAGE' 'La pression est résolue à 10E-6';
  186.      'QUITTER' CALCUL;
  187.    'FINSI';
  188.  
  189. 'FINSI';
  190.  
  191.  
  192. as2 ama1 = 'KOPS' 'MATRIK';
  193. 'RESPRO' as2 ama1;
  194.  
  195.  
  196. 'FINP';
  197.  
  198.  
  199.  
  200. *============================
  201. *  Construction du maillage:
  202. *============================
  203.  
  204. *  Construction des points:
  205.  
  206. p0 = 0. 0.;
  207. p1 = -6. 0.;
  208. p2 = -0.5 0.;
  209. p3 = 0. 0.5;
  210. p4 = 0.5 0.;
  211. p5 = 15. 0.;
  212. p6 = 15. 6.;
  213. p7 = 3. 6.;
  214. p8 = -3. 6.;
  215. p9 = -6. 6.;
  216.  
  217. *  construction des points de la partie supérieure:
  218.  
  219. p1p2 = p1 'DROIT' ((-1)*nx1) dini 1. dfin 0.1 p2;
  220. p2p3 = 'CERCLE' (nx2/2) p2 p0 p3;
  221. p3p4 = 'CERCLE' (nx2/2) p3 p0 p4;
  222. p4p5 = p4 'DROIT' ((-1)*nx3) dini 0.1 dfin 1.2 p5;
  223. p5p6 = p5 'DROIT' dini dbeg dfin dend p6;
  224. p6p7 = p6 'DROIT' ((-1)*nx3) dini 1.2 dfin 0.1 p7;
  225. p7p8 = p7 'DROIT' nx2 p8;
  226. p8p9 = p8 'DROIT' ((-1)*nx1) dini 0.1 dfin 1. p9;
  227. p9p1 = p9 'DROIT' dini dend dfin dbeg p1;
  228.  
  229. *  Création du contour:
  230. peri1 = p1p2 et p2p3 et p3p4 et p4p5;
  231. peri2 = p6p7 et p7p8 et p8p9;
  232. s1 = 'DALL' peri1 p5p6 peri2 p9p1 'PLAN';
  233.  
  234. *  construction des points de la partie inférieure:
  235. p10 = -6. -6.;
  236. p11 = -3. -6.;
  237. p12 = 3. -6.;
  238. p13 = 15. -6.;
  239. p14 = 0. -0.5;
  240. p1p10 = p1 'DROIT' dini dbeg dfin dend p10;
  241. p10p11 = p10 'DROIT' ((-1)*nx1) dini 1. dfin 0.1 p11;
  242. p11p12 = p11 'DROIT' nx2 p12;
  243. p12p13 = p12 'DROIT' ((-1)*nx3) dini 0.1 dfin 1.2 p13;
  244. p13p5 = p13 'DROIT' dini dend dfin dbeg p5;
  245. p5p4 = 'INVERSE' p4p5;
  246. p4p14 = 'CERCLE' (nx2/2) p4 p0 p14;
  247. p14p2 = 'CERCLE' (nx2/2) p14 p0 p2;
  248. p2p1 = 'INVERSE' p1p2;
  249.  
  250. *  Création du contour:
  251. peri3 = p10p11 et p11p12 et p12p13;
  252. peri4 = p5p4 et p4p14 et p14p2 et p2p1;
  253. s2 = 'DALL' peri3 p13p5 peri4 p1p10 'PLAN';
  254.  
  255. cercle = (p2p3 et p3p4 et p4p14 et p14p2);
  256. 'ELIMINER' 1.E-3 s1 s2;
  257.  
  258. * Création du domaine total:
  259. domtot = s1 et s2;
  260. DOM    = 'CONTOUR' domtot;
  261.  
  262. *==============================================
  263. *  Changement des éléments du maillage en QUAF:
  264. *==============================================
  265.  
  266. £domtot = 'CHANGER' domtot 'QUAF';
  267. £cercle = 'CHANGER' cercle 'QUAF';
  268. £entree = 'CHANGER' (p9p1 et p1p10) 'QUAF';
  269. £sortie = 'CHANGER' (p5p6 et p13p5) 'QUAF';
  270. £peri2  = 'CHANGER' peri2           'QUAF';
  271. £peri3  = 'CHANGER' peri3           'QUAF';
  272. £gau    = 'CHANGER' (p9p1 et p1p10) 'QUAF';
  273. £droi   = 'CHANGER' (p5p6 et p13p5) 'QUAF';
  274. £limite = 'CHANGER' (peri2 et p9p1 et p1p10 et peri3) 'QUAF';
  275. 'ELIM' 1.E-3 (£domtot et £cercle et £entree);
  276.  
  277. NN   = 'NBEL' domtot;
  278. 'MESSAGE' 'Nombre d éléments du maillage' NN;
  279.  
  280. *=======================================
  281. *  Formulation du domaine Navier Stokes:
  282. *=======================================
  283.  
  284. $domtot = 'MODE' £domtot 'NAVIER_STOKES' 'MACRO';
  285. $entree = 'MODE' £entree 'NAVIER_STOKES' 'MACRO';
  286. $sortie = 'MODE' £sortie 'NAVIER_STOKES' 'MACRO';
  287. mdomtot = 'DOMA' $domtot maillage;
  288. modelt  = 'MODELE' domtot 'THERMIQUE';
  289.  
  290. vol  = 'DOMA' $domtot 'VOLUME';
  291. volt = 'SOMT' vol;
  292.  
  293. *====================================
  294. *  Création des tables de résolution:
  295. *====================================
  296.  
  297. rey    = 300.;
  298. nu     = 1./rey;
  299. arelax = 1.;
  300.  
  301. RV = 'EQEX' 'OMEGA' 1. 'NITER' 1  'ITMA' maxit  'FIDT' 1
  302.      'ZONE' $domtot 'OPER' CALCUL;
  303.  
  304. RV = 'EQEX' RV
  305.    'OPTI' 'EFM1' 'CENTREE' 'IMPL'
  306. *      'OPTI' 'EF' 'CENTREE'
  307.      'ZONE' $domtot
  308.      'OPER' 'DFDT' 1. 'UN' 'DT' 'INCO' 'UN';
  309.  
  310. RV = 'EQEX' RV
  311.      'OPTI' 'EF' 'CENTREE' 'IMPL'
  312.         'ZONE' $domtot 'OPER' 'NS'   1. 'UN' 'INV_RE' 'INCO' 'UN'
  313. *       'ZONE' $domtot 'OPER' 'KONV' 1. 'UN' 'INV_RE' 'INCO' 'UN'
  314. *       'ZONE' $domtot 'OPER' 'LAPN' 'INV_RE' 'INCO' 'UN'
  315.         ;
  316.  
  317.  
  318. *  Implantation des conditions aux limites:
  319.  
  320. RV = 'EQEX' RV
  321.       'CLIM' 'UN' 'UIMP' £cercle 0. 'UN' 'VIMP' £cercle 0.
  322.              'UN' 'UIMP' £entree 1. 'UN' 'VIMP' £entree 0.
  323.              'UN'                        'VIMP' £peri2  0.
  324.              'UN'                        'VIMP' £peri3  0.
  325.       ;
  326.  
  327. *RV = 'EQEX' RV
  328. *     'CLIM' 'UN' 'UIMP' £cercle 0. 'UN' 'VIMP' £cercle 0.
  329. *            'UN' 'UIMP' £limite 1. 'UN' 'VIMP' £limite 0.;
  330.  
  331.  RV.'METHINV'.TYPINV=3 ;
  332. *RV.'METHINV'.TYPINV=1 ;
  333.  RV.'METHINV'.IMPINV=0 ;
  334.  RV.'METHINV'.NITMAX=300 ;
  335.  RV.'METHINV'.PRECOND=3 ;
  336.  RV.'METHINV'.RESID  =1.e-8  ;
  337.  RV. 'METHINV' . 'FCPRECT'=100 ;
  338.  RV. 'METHINV' . 'FCPRECI'=100 ;
  339.  
  340. RVP = 'EQEX'
  341.       'OPTI' 'EF' 'IMPL' INCOD DISP
  342.       'ZONE' $domtot
  343.       'OPER' 'KBBT' -1.   1.    'INCO' 'UN' 'PRES'
  344.       ;
  345.  
  346. RV.'PROJ'= RVP ;
  347.  
  348. RVP.'METHINV'.TYPINV   =1 ;
  349. RVP.'METHINV'.IMPINV   =0 ;
  350. RVP.'METHINV'.NITMAX   =150;
  351. RVP.'METHINV'.PRECOND  =3 ;
  352. RVP.'METHINV'.RESID    =1.e-12;
  353. RVP.'METHINV'.'FCPRECT'=100 ;
  354. RVP.'METHINV'.'FCPRECI'=100 ;
  355.  
  356.  
  357. *  Création de la table des inconnues et initialisation:
  358.  
  359. RV.INCO         = 'TABLE' INCO;
  360.  
  361.  
  362. RV.INCO.'UN'    = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3);
  363. RV.INCO.'UNM1'  = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3);
  364. RV.INCO.'DT'    = DT0  ;
  365. RV.INCO.'INV_RE'= nu;
  366. RV.INCO.'EVOLP' = 'PROG' 0.;
  367. RV.INCO.'PRES'  = 'KCHT' $domtot 'SCAL' DISP 0.;
  368. RV.INCO.'PRESNM1' = 'KCHT' $domtot 'SCAL' DISP 'COMP' 'PRES' 0.;
  369. RV.INCO.'PTILDE'  = 'KCHT' $domtot 'SCAL' DISP 0.;
  370. RVP.INCO        = RV.INCO;
  371.  
  372. *  Détermination d'un point d'observation:
  373. obs = mdomtot 'POINT' 'PROC' (0.8 0.5);
  374. RV.INCO.'PoiObs'= obs;
  375.  
  376. * Initialisation des inconnues complémentaires:
  377.  
  378. RV.INCO.'Uvit'  = 'PROG' 0.;
  379. RV.INCO.'ERX'   = 'PROG' 0.;
  380. RV.INCO.'ERY'   = 'PROG' 0.;
  381. RV.INCO.'errp'  = 'PROG' 0.;
  382. RV.INCO.'IT'    = 'PROG' 0;
  383. RV.INCO.'tps'   = 'PROG' 0.;
  384. RV.INCO.'entree'= £entree;
  385. RV.INCO.'cercle'= £cercle;
  386. RV.'TYPROJ'='VPI2';
  387.  
  388. EXEC RV;
  389.  
  390.  
  391.  
  392. * ------------ Controle calcul ----------------------------
  393.  
  394. Si (NON COMPLET) ;
  395.    PN = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE'
  396.                     ('NOMC' 'SCAL' (RV.INCO.'PRESSION')));
  397.    presc = 'REDU' PN £cercle;
  398.    precer = 'EVOL' (rouge) 'CHPO' presc cercle ;
  399.    tefp= extr precer 'ORDO' ;
  400.    list tefp ;
  401.    refp = PROG
  402.    .38960       .15815      -.34055      -.47886      -.42818
  403.   -.20337     -6.11799E-03  -.12634      -.26610      -.47590
  404.   -.58950      -1.0184      -1.3341      -1.1573      -.73313
  405.   -4.62908E-02   .38960;
  406.    errpc = abs ( tefp - refp )  maxi ;
  407.    mess ' errpc=' errpc ;
  408.    Si (errpc  > 0.005 ); erreur 5 ; Finsi ;
  409. Finsi ;
  410.  
  411. * ------------ Post traitement ----------------------------
  412.  
  413. Si GRAPH ;
  414. VX = 'EXCO' 'UX' (RV.INCO.'UN');
  415. VY = 'EXCO' 'UY' (RV.INCO.'UN');
  416. modu = ((VX**2) '+' (VY**2))**(0.5);
  417. umax = 'MAXIMUM' modu;
  418.  
  419. CHUN = 'VECTEUR' (RV.INCO.'UN') (1. '/' (1. * umax)) UX UY BLEU;
  420. 'TRACER' CHUN £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de vitesse';
  421.  
  422. *  Tracé du champ de pression
  423.    'TRACER' PN £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de pression';
  424.  
  425. *  Tracé du champ de pression lambda
  426.    PT = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE'
  427.                     ('NOMC' 'SCAL' (RV.INCO.'PTILDE')));
  428.    'TRACER' PT £domtot ('CONTOUR' £domtot)
  429.                     'TITRE' 'Champ de pression lambda';
  430.  
  431. *  Tracé de la pression sur le cylindre
  432.  
  433.    'DESSIN' precer 'TITRE' 'Valeurs de pression sur le cercle';
  434.  
  435. *  Tracé de la vitesse au point d'observation
  436.  
  437.    evv = 'EVOL' 'MANU' 't' (RV.INCO.'tps') 'Valeur de U'
  438.                         (RV.INCO.'Uvit');
  439.    'DESSIN' evv 'TITRE' 'Evolution de la vitesse au point (0.8 0.5)';
  440.  
  441. *  Tracé de la convergence
  442.  
  443.    evconv = 'EVOL' 'MANU' 'Iteration' (RV.INCO.'IT') 'log(err)'
  444.                                    (RV.INCO.'ERY');
  445.    'DESSIN' evconv 'TITRE' 'Convergence Vitesse ';
  446.  
  447.    evconvp = 'EVOL' 'MANU' 'Iteration' (RV.INCO.'IT')
  448.                            'errp' (RV.INCO.'errp');
  449.    'DESSIN' evconvp 'TITRE' 'Convergence pression';
  450.  
  451. * Calcul des lignes de courant:
  452. y = 'COOR' 2 (RV.INCO.'entree');
  453. psiim = 'KOPS' 1. '*' y;
  454. psii = 'NOMC' 'psi' psiim 'NATURE' 'DISCRET';
  455.  
  456. vitesse = RV.INCO.'UN';
  457. tourb = 'KOPS' vitesse 'ROT' $domtot;
  458. RK = 'EQEX' $domtot 'OPTI' 'EF' 'SEMI' 0.5
  459. *rk = 'EQEX' $domtot 'OPTI' 'EF' 'IMPL'
  460.      'ZONE' $domtot 'OPER' 'LAPN' 1. 'INCO' 'PSI'
  461.      'ZONE' $domtot 'OPER' 'FIMP' tourb 'INCO' 'PSI'
  462.      'CLIM' 'PSI' 'TIMP' (RV.INCO.'cercle') 0.
  463.      'CLIM' 'PSI' 'TIMP' (RV.INCO.'entree') psii;
  464.  
  465. RK.'INCO' = table 'INCO';
  466. RK.'INCO'.'PSI' = 'KCHT' $domtot 'SCAL' 'SOMMET' 0.;
  467.  
  468. EXEC RK;
  469.  
  470. * Traçage des lignes de courant:
  471. psi =RK.'INCO'.'PSI';
  472.  
  473. *'OPTION' trac ps;
  474. 'OPTION' isov ligne;
  475.  
  476. 'LISTE' ('MAXIMUM' PSI);
  477. 'LISTE' ('MINIMUM' PSI);
  478. 'TRACER' psi £domtot ('CONTOUR' £domtot)
  479. ('TITRE' 'Lignes de courant, Re='rey 'à t='RV.PASDETPS.'TPS');
  480. Finsi ;
  481.  
  482. * 'FIN' du jeu de données.
  483. fin;
  484.  
  485.  
  486.  
  487.  
  488.  

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