Télécharger la plaquette Cast3M


fr en

Télécharger slotevol.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  slotevol.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *  Jeu de données GIBIANE:
  5. * Description :
  6. * -------------
  7. * Ce jeu de données permet la discrétisation éléments finis d'un
  8. * écoulement dans une hypothèse de bas Mach. Il s'agit ici de comparer
  9. * les résultats obtenus avec CASTEM et ceux issus de l'article
  10. * Chenoweth et Paolucci (1986).
  11. *
  12. * Auteur : Franck PIGEONNEAU.
  13. * --------
  14. *
  15. * Date : 22/06/99
  16. * ------
  17.  
  18. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  19. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  20. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  21. *<><><><><><><><><><><>  DEBUT DU JEUX DE DONNEES  <><><><><><><><><><><
  22. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  23. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  24. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  25.  
  26. *=======================================================================
  27. *                   DEFINITION DE L'ESPACE PHYSIQUE:
  28. *=======================================================================
  29.  
  30. 'OPTION' 'DIME' 2 'ELEM' 'QUA8';
  31. 'OPTION' ECHO 0;
  32.  
  33. COMPLET = FAUX;
  34. 'SI' (COMPLET);
  35.  
  36. *  Constantes du maillage :
  37. *  ------------------------
  38.  
  39.    nx = 40;
  40.    ny = 40;
  41.  
  42. *  Nombre d'itérations :
  43.  
  44.    nbit = 100;
  45.  
  46. *  Nombre d'itérations dans la boucle de température :
  47.  
  48.    itmax = 5;
  49. 'SINON';
  50. *  Constantes du maillage :
  51. *  ------------------------
  52.  
  53.    nx = 7;
  54.    ny = 7;
  55.  
  56. *  Nombre d'itérations :
  57.  
  58.    nbit = 1;
  59.  
  60. *  Nombre d'itérations dans la boucle de température :
  61.  
  62.    itmax = 1;
  63. 'FINSI';
  64.  
  65. *=======================================================================
  66. *               DEFINITION DE 'TYPE' DE DISCRETISATION :
  67. *=======================================================================
  68.  
  69. DISVIT = 'QUAF';
  70. DISPRE = 'CENTREP1';
  71.  
  72. *=======================================================================
  73. *                   Définition des constantes du problème
  74. *=======================================================================
  75.  
  76. * Eccart de température adimensionné :
  77. eps   = 0.01;
  78. * Precision spatiale pour les 'ELIMINER' :
  79. epsi  = 1.D-7;
  80. * Rapport d'aspect :
  81. A     = 1.;
  82. * Nombre de Prandtl :
  83. pr    = 0.71;
  84. * Nombre de Rayleigh :
  85. ra    = 1.E+6;
  86. * Inverse du nombre de Froude :
  87. invfr = (ra*pr)/(2.*eps);
  88. * Coefficients des lois de Sutherland pour la conductibilité thermique
  89. * et la viscosité dynamique :
  90. sk    = 0.648;
  91. skp1  = 1.+sk;
  92. smu   = 0.368;
  93. smup1 = 1.+smu;
  94. * Rapport des chaleurs spécifiques :
  95. gamma = 1.4;
  96. * Résilience du fluide :
  97. reli  = (gamma-1.)/gamma;
  98. * Pas de temps :
  99. dt = 1.D-3;
  100.  
  101. * Valeur de la précision sur la pression et la température :
  102. ecarpre = 1.D-4;
  103. * Valeur de la précision sur la vitesse, la pression et la température :
  104. ecarmax = 1.D-6;
  105.  
  106. *=======================================================================
  107. *                                MAILLAGE
  108. *=======================================================================
  109.  
  110.  
  111. * Construction des points :
  112. * -------------------------
  113.  
  114. p1 = 0. 0.;
  115. p2 = 1. 0.;
  116. p3 = 1. (0.5*A);
  117. p4 = 1. A;
  118. p5 = 0. A;
  119. p6 = 0. (0.5*A);
  120. p7 = (0.5*A) 0.;
  121. p8 = (0.5*A) 1.;
  122. p9 = (0.5*A) (0.5*A);
  123.  
  124. * Densités d'éléments :
  125. * ---------------------
  126.  
  127. d1 = 0.8;
  128. d2 = 0.83;
  129.  
  130. * Construction des lignes :
  131. * -------------------------
  132.  
  133. p1p7 = p1 d (-1*nx) dini d1 dfin d2 p7;
  134. p7p2 = p7 d (-1*nx) dini d2 dfin d1 p2;
  135. p1p2 = p1p7 et p7p2;
  136. p2p3 = p2 d (-1*ny) dini d1 dfin d2 p3;
  137. p3p4 = p3 d (-1*ny) dini d2 dfin d1 p4;
  138. p4p8 = p4 d (-1*nx) dini d1 dfin d2 p8;
  139. p8p5 = p8 d (-1*nx) dini d2 dfin d1 p5;
  140. p4p5 = p4p8 et p8p5;
  141. p5p6 = p5 d (-1*ny) dini d1 dfin d2 p6;
  142. p6p1 = p6 d (-1*ny) dini d2 dfin d1 p1;
  143. p3p9 = p3 d (-1*nx) dini d1 dfin d2 p9;
  144. p9p6 = p9 d (-1*nx) dini d2 dfin d1 p6;
  145. p3p6 = p3p9 et p9p6;
  146. p6p3 = 'INVERSE' p3p6;
  147.  
  148. * Construction des contours servant à l'établissement des conditions
  149. * ------------------------------------------------------------------
  150. * aux limites :
  151. * -------------
  152.  
  153. lgau = 'INVERSE' (p5p6 et p6p1);
  154. ldro = 'INVERSE' (p2p3 et p3p4);
  155. lbas = p1p2;
  156. lhau = p4p5;
  157.  
  158. * Construction des surfaces :
  159. * ---------------------------
  160.  
  161. s1 = 'DALLER' p1p2 p2p3 p3p6 p6p1;
  162. s2 = 'DALLER' p6p3 p3p4 p4p5 p5p6;
  163. 'ELIMINER' epsi s1 s2;
  164. domtot = s1 et s2;
  165.  
  166. * Détermination de la surface extérieure :
  167. * ----------------------------------------
  168.  
  169. surtot = 'CONTOUR' domtot;
  170. surtot = 'ENVELOPPE' surtot;
  171.  
  172. *=======================================================================
  173. *                          Création du modèle
  174. *=======================================================================
  175.  
  176. * Transformations des éléments du maillage en QUAF :
  177. * --------------------------------------------------
  178.  
  179. £domtot = 'CHANGER' domtot 'QUAF';
  180. £gau    = 'CHANGER' lgau   'QUAF';
  181. £dro    = 'CHANGER' ldro   'QUAF';
  182. £bas    = 'CHANGER' lbas   'QUAF';
  183. £hau    = 'CHANGER' lhau   'QUAF';
  184. £surtot = 'CHANGER' surtot 'QUAF';
  185.  
  186. * Récupération du nombre de noeuds :
  187. * ----------------------------------
  188.  
  189. ntot = 'NBNO' £domtot;
  190.  
  191. * Elimination des eventuels points en doubles :
  192. * ---------------------------------------------
  193.  
  194. 'ELIMINER' epsi (£domtot et £gau et £dro et £bas et £hau);
  195.  
  196. * Création des objets MMODEL :
  197. * ----------------------------
  198.  
  199. * Modèles Navier Stokes :
  200. * -----------------------
  201.  
  202. $domtot = 'MODELE' £domtot 'NAVIER_STOKES' DISVIT;
  203. $hau    = 'MODELE' £hau    'NAVIER_STOKES' DISVIT;
  204. $bas    = 'MODELE' £bas    'NAVIER_STOKES' DISVIT;
  205. $gau    = 'MODELE' £gau    'NAVIER_STOKES' DISVIT;
  206. $dro    = 'MODELE' £dro    'NAVIER_STOKES' DISVIT;
  207. $surtot = 'MODELE' £surtot 'NAVIER_STOKES' DISVIT;
  208.  
  209. mdomtot = 'DOMA' £domtot;
  210.  
  211. * Modèles Thermique :
  212. * ------------------
  213.  
  214. modelt  = 'MODELE' domtot 'THERMIQUE';
  215.  
  216. *                             ====================
  217. *                             Calcul des volumes :
  218. *                             ====================
  219.  
  220. vol  = 'DOMA' $domtot 'VOLUME';
  221. volt = 'SOMT' vol;
  222. vol  = 'KCHT' $domtot 'SCAL' 'CENTRE' vol;
  223.  
  224. *=======================================================================
  225. *                 CREATION DE LA TABLE DE RESOLUTION
  226. *
  227. * A l'instant tn-1, on connait TNM1 (la température), PTNM1 (la pression
  228. * thermodynamique), MUNM1 (la viscosité dynamique), KNM1 (la conductivité
  229. * thermique), RHONM1 (la masse volumique). L'algorithme de résolution
  230. * est le suivant:
  231. *
  232. * 1- Résolution de l'équation d'énergie permettant le calcul de TN;
  233. * 2- Détermination de MUN et KN à partir des lois de Sutherland;
  234. * 3- Résolution de l'équation portant sur PTN;
  235. * 4- Détermination de RHON à l'aide de l'équation d'état;
  236. * 5- Calcul de RHOM = 0.5(RHONM1+RHON) et de ZN = Ln(RHOM);
  237. * 6- Résolution des équation de Navier-Stokes avec masse volumique
  238. *    variable.
  239. *
  240. *=======================================================================
  241.  
  242. * Discrétisation des équations Navier-Stokes à masse volumique variable:
  243. * ----------------------------------------------------------------------
  244.  
  245. * Discrétisation de la dérivée temporelle :
  246.  
  247. RVU = 'EQEX' 'NITER' 1 'OMEGA' 1.0  'ITMA' 1
  248.       'OPTI' 'EF' 'CENTREE'
  249.       'ZONE' $domtot
  250.       'OPER' 'DFDT' 'RHOM' 'UN' 'DT'
  251.       'INCO' 'UN';
  252.  
  253. * Discrétisation du terme de convection :
  254.  
  255. RVU = 'EQEX' RVU
  256.       'OPTI' 'EF' 'CENTREE' 'IMPL'
  257.       'ZONE' $domtot
  258.       'OPER' 'KONV' 'RHOM' 'UN' 'MUN'
  259.       'INCO' 'UN';
  260.  
  261. * Discrétisation du terme de diffusion visqueuse :
  262.  
  263. RVU = 'EQEX' RVU
  264.       'OPTI' 'EF' 'CENTREE' 'IMPL'
  265.       'ZONE' $domtot
  266.       'OPER' 'LAPN' 'MUN'
  267.       'INCO' 'UN';
  268.  
  269. * Discrétisation du terme source de l'équation de QDM :
  270.  
  271. RVU = 'EQEX' RVU
  272.       'OPTI' 'EF' 'CENTREE'
  273.       'ZONE' $domtot
  274.       'OPER' 'TOIM' 'RHOG'
  275.       'INCO' 'UN';
  276.  
  277. * Discrétisation de la pression :
  278.  
  279. RVU = 'EQEX' RVU
  280.       'OPTI' 'EF' 'CENTREE' DISPRE
  281.       'ZONE' $domtot
  282.       'OPER' 'KBBT' (+1.)
  283.       'INCO' 'UN' 'PN';
  284.  
  285. * Ajout du terme source dans l'équation de continuité :
  286.  
  287. RVU = 'EQEX' RVU
  288.       'OPTI' 'EF' 'CENTREE' 'INCOD' DISPRE
  289.       'ZONE' $domtot
  290.       'OPER' 'FIMP' 'SN'
  291.       'INCO' 'PN';
  292.  
  293. * Conditions aux limites :
  294. * ------------------------
  295.  
  296. RVU = 'EQEX' RVU
  297.       'CLIM' 'UN' 'UIMP' (£bas et £hau et £dro et £gau) 0.
  298.              'UN' 'VIMP' (£bas et £hau et £dro et £gau) 0.
  299.              ;
  300. bcp='ELEM' ('DOMA' $domtot DISPRE)  POI1 (lect 1) ;
  301. * Suppression des modes de pression constante
  302. RVU = EQEX  RVU
  303.    'CLIM'
  304.           'PN' 'TIMP' bcp 0.D0  ;
  305.  
  306. * Etablissement de la table des inconues et conditions initiales:
  307. * ---------------------------------------------------------------
  308.  
  309. RVU.INCO        = 'TABLE' INCO;
  310.  
  311. * Inconnues liées à la table RVU :
  312.  
  313. RVU.INCO.'UN'   = 'KCHT' $domtot 'VECT' 'SOMMET' (0. 0.);
  314. RVU.INCO.'PN'   = 'KCHT' $domtot 'SCAL' DISPRE 0.;
  315.  
  316. * Variables associées à la discrétisation des équations de
  317. * Navier-Stokes :
  318.  
  319. RVU.INCO.'RHON' = 'KCHT' $domtot 'SCAL' 'CENTRE'   1.;
  320. RVU.INCO.'RHOM' = 'KCHT' $domtot 'SCAL' 'CENTRE'   1.;
  321. RVU.INCO.'MUN'  = 'KCHT' $domtot 'SCAL' 'CENTRE'   pr;
  322. RVU.INCO.'RHOG' = 'KCHT' $domtot 'VECT' 'CENTRE'  (0. 0.);
  323.  
  324. * Variables complémentaires :
  325.  
  326. RVU.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' (0. 0.);
  327. RVU.INCO.'UNM'  = 'KCHT' $domtot 'VECT' 'SOMMET' (0. 0.);
  328. RVU.INCO.'PTN'  = 1.;
  329. RVU.INCO.'TN'   = 'KCHT' $domtot 'SCAL' 'SOMMET' 0.;
  330. RVU.INCO.'KN'   = 'KCHT' $domtot 'SCAL' 'CENTRE' 1.;
  331. RVU.INCO.'PRES' = 'PROG' 1.;
  332. RVU.INCO.'TPS'  = 'PROG' 0.;
  333. RVU.INCO.'STN'  = 0.;
  334.  
  335. * Paramètres de la méthode d'inversion :
  336.  
  337. RVU.'METHINV'.TYPINV = 1;
  338. RVU.'METHINV'.NITMAX = 600;
  339. RVU.'METHINV'.RESID  = 1.E-7;
  340.  
  341. * Pas de temps :
  342. * --------------
  343.  
  344. RVU.INCO.'DT'= dt;
  345.  
  346. * Numero du pas de temps :
  347.  
  348. nupadt = 0;
  349.  
  350. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  351. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  352. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  353. *<><><><><><><><><><><><><  PROCEDURE BASMACH  <><><><><><><><><><><><><
  354. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  355. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  356. *<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><
  357.  
  358. 'DEBPROC' BASMACH;
  359. *   ====================================================================
  360. *                       DEBUT DE LA PROCEDURE BASMACH
  361. *   ====================================================================
  362.  
  363. 'ARGUMENT' nbit*'ENTIER';
  364.  
  365. tps = 0.;
  366.  
  367. 'REPETER' BCLTPS nbit;
  368.  
  369.     nupadt= nupadt + 1;
  370.     mess 'Nombre de pas de temps' nupadt;
  371.  
  372. *   Initialisation des champs nécessaire à la résolution de l'équation
  373. *   sur la température :
  374.  
  375.     PTNM1 = RVU.INCO.'PTN';
  376.     PT    = PTNM1;
  377.     TNM1  = RVU.INCO.'TN';
  378.     TT    = TNM1;
  379.     RHON  = RVU.INCO.'RHON';
  380.     UM    = RVU.INCO.'UNM';
  381.     KN    = RVU.INCO.'KN';
  382.     STN   = RVU.INCO.'STN';
  383.  
  384.     'REPETER' BCLTEMP itmax;
  385. *      =================================================================
  386. *                       ETAPE 1 : DETERMINATION DE TN
  387. *      =================================================================
  388.  
  389.        TN = CAL_TEMP TNM1 RHON UM KN STN;
  390.  
  391. *      Détermination de la convergence sur la température :
  392. *      ----------------------------------------------------
  393.  
  394.        errt = 'KOPS' TN '-' TT;
  395.        errt = errt*errt;
  396.        errt = 'CHANGER' 'CHAM' errt modelt 'MASSE';
  397.        errtl2 = 'INTG' modelt errt;
  398.        errtl2 = errtl2**0.5;
  399.        TT = TN;
  400.  
  401. *      =================================================================
  402. *                      ETAPE 2 : DETERMINATION DE PTN
  403. *      =================================================================
  404.  
  405. *      Détermination de l'inverse de la température :
  406. *      ----------------------------------------------
  407.  
  408.        TNC = 'NOEL' $domtot TN;
  409.        TNC = 'KOPS' TNC '+' 1.;
  410.        invtn = 'KOPS' 1. '/' TNC;
  411.        moyinvt = 'SOMT' (vol*invtn);
  412.        moyinvt = moyinvt/volt;
  413.        PTN = 1./moyinvt;
  414.  
  415. *      Détermination de la convergence sur la pression :
  416. *      -------------------------------------------------
  417.  
  418.        errpre = 'ABS' (PTN-PT);
  419.        PT = PTN;
  420.  
  421. *      ==============================================================
  422. *                         ETAPE 3 : CALCUL DE KN
  423. *      ==============================================================
  424.  
  425. *       KN  = 'KOPS' ('KOPS' skp1 '*' ('KOPS' TNC '**' 1.5)) '/' ('KOPS'
  426. *              TNC '+' sk);
  427. *       KN  = 'KCHT' $domtot 'SCAL' 'CENTRE' KN;
  428.        KN  = 'KCHT' $domtot 'SCAL' 'CENTRE' 1.;
  429.  
  430. *      ==============================================================
  431. *                  ETAPE 4 : DETERMINATIONS DE dpdt ET STN
  432. *      ==============================================================
  433.  
  434. *      Calcul de dpdt :
  435. *      ----------------
  436.  
  437. *      Calcul du débit volumique :
  438. *      ---------------------------
  439.  
  440.        usurf = 'REDU' UM £surtot;
  441.        phi = 'DBIT' usurf $surtot;
  442.        phi = phi/volt;
  443.  
  444. *      Calcul du flux de chaleur :
  445. *      ---------------------------
  446.  
  447.        TM1 = 'NOMC' 'T' TN;
  448.        gradt = 'GRAD' modelt TM1;
  449.  
  450. *      NOTA BENE : Attention ici, gradt est un champ par élément.
  451. *      Transformation du champ par élément en champ point :
  452.  
  453.        gradt = 'CHANGER' 'CHPO' modelt gradt;
  454.  
  455. *      Multiplication par KN :
  456. *      -----------------------
  457.  
  458. *      Extraction des composantes :
  459.  
  460.        gradtx = 'EXCO' 'T,X' gradt 'GX' 'NATURE' 'DISCRET';
  461.        gradty = 'EXCO' 'T,Y' gradt 'GY' 'NATURE' 'DISCRET';
  462.        gradtx = 'KCHT' $domtot 'SCAL' 'SOMMET' 'COMP' 'GX' gradtx;
  463.        gradty = 'KCHT' $domtot 'SCAL' 'SOMMET' 'COMP' 'GY' gradty;
  464.  
  465. *      Multiplication de chaque composante par KN :
  466. *      --------------------------------------------
  467.  
  468.        KN1 = 'ELNO' $domtot KN;
  469.        KN1 = 'KCHT' $domtot 'SCAL' 'SOMMET' KN1;
  470.        kngratx = 'KOPS' KN1 '*' gradtx;
  471.        kngraty = 'KOPS' KN1 '*' gradty;
  472.        kngratx = 'NOMC' 'GX' kngratx 'NATURE' 'DISCRET';
  473.        kngraty = 'NOMC' 'GY' kngraty 'NATURE' 'DISCRET';
  474.        kngratx = 'KCHT' $domtot 'SCAL' 'SOMMET' 'COMP' 'GX' kngratx;
  475.        kngraty = 'KCHT' $domtot 'SCAL' 'SOMMET' 'COMP' 'GY' kngraty;
  476.  
  477. *      Calcul de la div(KN GRAD T) :
  478. *      -----------------------------
  479.  
  480. *      1. Composante suivant x :
  481. *      -------------------------
  482.  
  483.        UX = 'NOMC' 'T' kngratx 'NATURE' 'DISCRET';
  484.        gradUX = 'GRAD' modelt UX;
  485.  
  486. *      Extraction de la composante suivant x du gradient :
  487.        UX,X = 'EXCO' 'T,X' gradUX 'div';
  488.  
  489. *      2. Composante suivant y :
  490. *      -------------------------
  491.  
  492.        UY = 'NOMC' 'T' kngraty 'NATURE' 'DISCRET';
  493.        gradUY = 'GRAD' modelt UY;
  494.  
  495. *      Extraction de la composante suivant y du gradient :
  496.        UY,Y = 'EXCO' 'T,Y' gradUY 'div';
  497.  
  498. *      3. Détermination de la divergence par sommation :
  499. *      -------------------------------------------------
  500.  
  501.        div = UX,X+UY,Y;
  502.  
  503. *      Calcul de dpdt/gamma par intégration de div :
  504. *      ---------------------------------------------
  505.  
  506.        dpdtgm1 = 'INTG' modelt div;
  507.        dpdtgm1 = dpdtgm1/volt;
  508.        dpdt = (dpdtgm1*gamma) + ((gamma*PTN)*phi);
  509.  
  510. *      Détermination du terme source de l'équation de la chaleur :
  511. *      -----------------------------------------------------------
  512.  
  513.        STN = (reli*dpdt);
  514.  
  515. *      ==============================================================
  516. *                    ETAPE 5 : DETERMINATION DE RHON
  517. *      ==============================================================
  518.  
  519.        RHON = 'KOPS' PTN '/' TNC;
  520.        RHON = 'KCHT' $domtot 'SCAL' 'CENTRE' RHON;
  521.  
  522. *      Test de la convegence dans la boucle de température :
  523.  
  524.        erregl = 'PROG' errtl2 errpre;
  525.        errmax = 'MAXIMUM' erregl;
  526.        mess 'L erreur est de ' errmax;
  527.        'SI' (errmax < ecarpre);
  528.          'QUITTER' BCLTEMP;
  529.        'FINSI';
  530.    'FIN' BCLTEMP;
  531.  
  532. *   Détermination de la convergence sur la température :
  533. *   ----------------------------------------------------
  534.  
  535.     errt   = 'KOPS' TN '-' TNM1;
  536.     errt   = errt*errt;
  537.     errt   = 'CHANGER' 'CHAM' errt modelt 'MASSE';
  538.     errtl2 = 'INTG' modelt errt;
  539.     errtl2 = errtl2**0.5;
  540.  
  541. *   Détermination de la convergence sur la pression :
  542. *   -------------------------------------------------
  543.  
  544.     errpre = 'ABS' (PTN-PTNM1);
  545.  
  546. *   Sauvegarde pour l'itération temporelle suivante :
  547.  
  548.     RVU.INCO.'TN' = TN;
  549.     RVU.INCO.'KN' = KN;
  550.     RVU.INCO.'PTN'  = PTN;
  551.     RVU.INCO.'STN'  = STN;
  552.     mess 'La pression est égale à ' PTN;
  553.     PRES = 'PROG' PTN;
  554.     RVU.INCO.'PRES'= RVU.INCO.'PRES' et PRES;
  555.  
  556. *   ====================================================================
  557. *                  ETAPE 6 : DETERMINATIONS DE RHOM et RHOG
  558. *   ====================================================================
  559.  
  560.     RHONM1 = RVU.INCO.'RHON';
  561.     RVU.INCO.'RHON' = RHON;
  562.  
  563. *   Calcul de la masse volumique moyenne entre les instants tn et tn-1 :
  564.  
  565.     RHOM  = 0.5*(RHON+RHONM1);
  566.     RVU.INCO.'RHOM' = 'KCHT' $domtot 'SCAL' 'CENTRE' RHOM;
  567.  
  568. *   Détermination du terme source de l'équation de QDM :
  569.  
  570.     rhog = 'KOPS' invfr '*' ('KOPS' RHOM '-' 1.);
  571.     rgy  = 'NOMC' 'GY' rhog 'NATURE' 'DISCRET';
  572.     RHGX = 'KCHT' $domtot 'SCAL' 'CENTRE' 'COMP' 'GX' 0.;
  573.     RHGY = 'KCHT' $domtot 'SCAL' 'CENTRE' 'COMP' 'GY' rgy;
  574.     RVU.INCO.'RHOG' = 'KCHT' $domtot 'VECT' 'CENTRE' 'COMP' 'GX' 'GY'
  575.                       (RHGX ET RHGY);
  576.  
  577. *   ====================================================================
  578. *            ETAPE 7 : DETERMINATION DU TERME SOURCE SN DE
  579. *                           L'EQUATION DE CONTINUITE
  580. *   ====================================================================
  581.  
  582. *   Calcul de la dérivée logarithmique de PTN divivée par gamma :
  583. *   -------------------------------------------------------------
  584.  
  585.     dlnpdt = dpdtgm1/PTN;
  586.  
  587. *   Calcul de la -div(KN GRAD T)/PTN :
  588. *   ----------------------------------
  589.  
  590.     div = div/(-1.*PTN);
  591.  
  592. *   Création du CHAML dlnpdt :
  593.  
  594.     dlnpdt = 'KCHT' $domtot 'SCAL' 'SOMMET' 'COMP' 'div' dlnpdt;
  595.     dlnpdt = 'CHANGER' 'CHAM' dlnpdt modelt 'MASSE';
  596.  
  597. *   Calcul de la dlnpdt - div(KN GRAD T)/PTN :
  598. *   ------------------------------------------
  599.  
  600.     div = div + dlnpdt;
  601.  
  602. *   Changement de div de CHAML en 'CHPO' :
  603.  
  604.     div = 'CHANGER' 'CHPO' modelt div;
  605.  
  606. *   Application de div sur le support du modèle Navier Stokes :
  607.  
  608.     RVU.INCO.'SN' = 'KCHT' $domtot 'SCAL' 'SOMMET' 'COMP' 'div' div;
  609.  
  610. *   ====================================================================
  611. *                      ETAPE 8 : CALCUL DE MUN
  612. *   ====================================================================
  613.  
  614. *    MUN = 'KOPS' ('KOPS' (pr*smup1) '*' ('KOPS' TNC '**' 1.5)) '/'
  615. *          ('KOPS' TNC '+' smu);
  616. *    RVU.INCO.'MUN' = 'KCHT' $domtot 'SCAL' 'CENTRE' MUN;
  617.     RVU.INCO.'MUN' = 'KCHT' $domtot 'SCAL' 'CENTRE' pr;
  618.  
  619. *   ====================================================================
  620. *              ETAPE 9 : RESOLUTION DES EQUATIONS DE NAVIER-STOKES
  621. *   ====================================================================
  622.  
  623.     EXEC RVU;
  624.  
  625. *   ====================================================================
  626. *                ETAPE 10 : CALCUL DE UNM ET SAUVEGARDE DE UNM1
  627. *   ====================================================================
  628.  
  629.     UN   = RVU.INCO.'UN';
  630.     UNM1 = RVU.INCO.'UNM1';
  631.     UNM  = 0.5*((3.*UN)-UNM1);
  632.  
  633. *   Calcul de la convergence sur la vitesse :
  634.  
  635.     DUN    = 'KOPS' UN '-' UNM1;
  636.  
  637.     DUNS   =  DUN (mots 'UX' 'UY') 'PSCA' DUN (mots 'UX' 'UY');
  638.     DUNS   = 'CHANGER' 'CHAM' DUNS modelt 'MASSE';
  639.     errvl2 = 'INTG' modelt DUNS;
  640.     errvl2 = errvl2**0.5;
  641.  
  642. *   Sauvegarde des champs par points dans le liste d'inconnues de RVU :
  643.  
  644.     RVU.INCO.'UNM'  = 'KCHT' $domtot 'VECT' 'SOMMET' UNM;
  645.     RVU.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' UN;
  646.  
  647. *   Incrémentation du pas de temps :
  648.  
  649.     tps = tps + RVU.INCO.'DT' ;
  650.     tt  = 'PROG' tps;
  651.     RVU.INCO.'TPS' = RVU.INCO.'TPS' et tt;
  652.  
  653. *   Création de la liste de réels regroupant l'ensemble des erreurs
  654. *   ---------------------------------------------------------------
  655. *   sur la température et sur la pression :
  656. *   ---------------------------------------
  657.  
  658.     erregl = 'PROG' errtl2 errpre errvl2;
  659.     errmax = 'MAXIMUM' erregl;
  660.  
  661.     'SI' (errmax < ecarmax);
  662.        mess 'L état stationnaire est atteint.';
  663.        'QUITTER' BCLTPS;
  664.     'FINSI';
  665.  
  666. 'FIN' BCLTPS;
  667.  
  668. *   ====================================================================
  669. *                          FIN DE LA PROCEDURE BASMACH
  670. *   ====================================================================
  671. 'FINPROC';
  672.  
  673. 'DEBPROC' CAL_TEMP;
  674. *   ====================================================================
  675. *                       DEBUT DE LA PROCEDURE TEMPERATURE
  676. *   ====================================================================
  677.  
  678. 'ARGUMENT' TNM1*'CHPOINT' RHON*'CHPOINT' UM*'CHPOINT' KN*'CHPOINT'
  679.            STN*'FLOTTANT';
  680.  
  681. *                =============================================
  682. *                Discrétisation de l'équation sur la chaleur :
  683. *                =============================================
  684.  
  685. * Discrétisation de la dérivée temporelle :
  686. * -----------------------------------------
  687.  
  688. RVT = 'EQEX' 'NITER' 1 'OMEGA' 1.0  'ITMA' 1
  689.       'OPTI' 'EF' 'CENTREE'
  690.       'ZONE' $domtot
  691.       'OPER' 'DFDT' 'RHON' 'TNM1' 'DT'
  692.       'INCO' 'TN';
  693.  
  694. * Discrétisation du terme de convection :
  695. * ---------------------------------------
  696.  
  697. RVT = 'EQEX' RVT
  698.       'OPTI' 'EF' 'CENTREE' 'IMPL'
  699.       'ZONE' $domtot
  700.       'OPER' 'KONV' 'RHON' 'UM' 'KN'
  701.       'INCO' 'TN';
  702.  
  703. * Discrétisation du terme de diffusion :
  704. * --------------------------------------
  705.  
  706. RVT = 'EQEX' RVT
  707.       'OPTI' 'EF' 'CENTREE' 'IMPL'
  708.       'ZONE' $domtot
  709.       'OPER' 'LAPN' 'KN'
  710.       'INCO' 'TN';
  711.  
  712. * Discrétisation du terme source :
  713. * --------------------------------
  714.  
  715. RVT = 'EQEX' RVT
  716.       'OPTI' 'EF'
  717.       'ZONE' $domtot
  718.       'OPER' 'FIMP' 'STN'
  719.       'INCO' 'TN';
  720.  
  721. * Discrétisation de la condition au limite de type Neumann sur les
  722. * ----------------------------------------------------------------
  723. * parois supérieure et inférieure :
  724. * ---------------------------------
  725.  
  726. RVT = 'EQEX' RVT
  727.       'OPTI' 'EF'
  728.       'ZONE' ($bas et $hau)
  729.       'OPER' 'FIMP' 0.
  730.       'INCO' 'TN';
  731.  
  732. * Conditions aux limites sur l'équation de la chaleur :
  733. * -----------------------------------------------------
  734.  
  735. RVT = 'EQEX' RVT
  736.       'CLIM' 'TN' 'TIMP' £gau eps
  737.              'TN' 'TIMP' £dro (-1.*eps);
  738.  
  739. * Etablissement de la table inco relative à RVT :
  740. * -----------------------------------------------
  741.  
  742. RVT.INCO        = 'TABLE' INCO;
  743. RVT.INCO.'TN'   = 'KCHT' $domtot 'SCAL' 'SOMMET' 0.;
  744. RVT.INCO.'TNM1' = TNM1;
  745. RVT.INCO.'KN'   = KN;
  746. RVT.INCO.'STN'  = STN;
  747. RVT.INCO.'RHON' = RHON;
  748. RVT.INCO.'UM'   = UM;
  749. RVT.INCO.'DT'   = dt;
  750.  
  751. *=======================================================================
  752. *                    RESOLUTION DE L'EQUATION SUR LA TEMPERATURE
  753. *=======================================================================
  754.  
  755. EXEC RVT;
  756.  
  757. TN = RVT.INCO.'TN';
  758.  
  759. *=======================================================================
  760. *                       FIN DE LA PROCEDURE TEMPERATURE
  761. *=======================================================================
  762. 'FINPROC' TN;
  763.  
  764. *=======================================================================
  765. *                      EXECUTION DE LA RESOLUTION
  766. *=======================================================================
  767.  
  768. BASMACH nbit;
  769.  
  770. * 'FIN' de jeux de données:
  771. fin;
  772.  
  773.  
  774.  
  775.  
  776.  
  777.  
  778.  
  779.  
  780.  

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