Accueil | Cast3M

* fichier :  centrif.dgibi
************************************************************************
************************************************************************
****************************************************************************
*
*
*         Date : 19/3/97
*
*         Version 1
*
*         Objet      : Le but de ce test est de valider la mise
*                      en place sous CASTEM 2000 de calculs en
*                      "vrai axi-symétrique", c'est à dire pour
*                      lesquels la composante ortho-radiale de
*                      la vitesse n'est pas nulle mais par contre
*                      toutes les dérivées par rapport à la composante
*                      ortho-radiale du repère cylindrique sont nulles.
*
*                      Les résultats issus du calcul numérique seront
*                      comparés aux résultats prédits par la théorie
*                      (voir M. Ungarish, Hydrodynalics of Suspensions,
*                      p. 47, Springer-Verlag, 1993) et le rapport
*                      CEA DMT 97/202 : étude phéno d'un bol en rotation.
*
*
*         Description: Le cas test est réalisé avec un tore muni
*                      d'une entrée et d'une sortie. Le fluide est injecté
*                      axialement par en dessous et sort librement par
*                      au dessus en un autre endroit. On s'attend à trouver
*                      deux colonnes de fluide au niveau des surfaces d'entrée
*                      et de sortie du fluide ainsi que des couches limites
*                      d'Ekman le long des parois horizontales.
*
*                      On choisit winj = 1 m/s
*                                 omeg = 100 rad/s
*                                 nu = 1/10 SI
*
*                      Ce qui permet d'avoir E = 10^-3  (nombre d'Ekman)
*                                            RO = 1/100 (nombre de Rossby)
*
*          Représentation graphique
*          ------------------------
*
*          Les traits figurants à l'intérieur du rectangle ne sont là
*          que pour mettre en évidence la symétrie du dessin.
*
*          On a laissé des blancs pour représenter les entrées et sortie
*          du fluide.
*
*     z
*     ^  axe de révolution
*     |
*     |--------------------------------------------->
*     |           inj2
*     |
*     |------------------------------------->
*     |           inj1
*     |
*     |                       out1 out2      sym1 sym2
*     |               P4_______|    |___________________________P3
*     |                |                                        |
*     |  couche Ekman  |       |    |        |      |           |
*     |             P44+----------------------------------------+P44
*     |                |       |    |        |      |           |
*     |                |                                        |
*     |                |       |    |        |      |           |
*     |                |                                        |
*     |                |       |    |        |      |           |
*     |                |                                        |
*     |              P8+       |    |        |      |           +P6
*     |                |                                        |
*     |                |       |    |        |      |           |
*     |                |                                        |
*     |                |       |    |        |      |           |
*     |                |                                        |
*     |                |       |    |        |      |           |
*     |             P11+----------------------------------------+P22
*     |   couche Ekman |       |    |        |      |           |
*    0|_ _ _ _ _ _ _ _ |______________________      ____________| _ _ _\ x
*     |0              P1      syo1 syo2      ent1 ent2          P2     /
*     |
*     |-------------->
*     |    rint
*     |
*     |---------------------->
*     |        exi1
*     |
*     |---------------------------->
*     |            exi2
*     |
*     |--------------------------------------------------------->
*     |                rext 
*     |
*
*         Les deux colonnes qui se forment à l'entrée et la sortie du fluide
*         sont entourées de couches limites verticales dites
*         de Stewartson (non représentées).
*
*         
*
*         Auteur     : Gilles Bernard-Michel
*
*         Révision 1 : 16/06/2014 Passage EQPR -> EQEX
*
*****************************************************************************
*
 
 
 
*
** Choix des éléments
*
 
option 'DIME' 2 'ELEM' qua8 'MODELISER' axis ;
 
DISCR  = 'MACRO'  ;
KPRESS = 'CENTRE' ;
BETAP  = 1.       ;
 
*
************************ Quelques procédures ********************************
*
 
*
** Max de l'erreur entre deux Champoints
*
 
DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  vmax = 'MAXIMUM' vit 'ABS' + 1.e-13; 
  dimax = 'MAXIMUM' (vitp1 - vit) 'ABS' ;
  err = (dimax /vmax) ;
FINPROC err ;
 
*
************************ On définit les propriétés physiques ***************
*
 
*
** viscosité dynamique mu/rho.
*
 
nu =1.E-1 ;
 
*
** La vitesse de rotation initiale.
*
 
omeg = 100 ;
 
*
** Ordre de grandeur de la vitesse d'injection.
*
 
winj = 0.2 ;
 
*
********************* Les options utilisateur ******************************
*
 
*
** Critère de convergence implicite.
*
 
 crit = 0.10 ;
 
*
** Test long ou court.
*
 
 COMPLET = FAUX;
 
*
** Nombre max de boucles d'implicitation.
*
 
MAXIMP = 3;
 
*
** Sortie graphique.
*
 
GRAPH = 'N';
 
*
** Calcul des évolutions de vitesses.
*
 
 CALCEV = VRAI;
 
*
** periode de calcul des evolutions.
*
 
period = 10;
 
*
********************* On traite les aspects géométriques ********************
*
 
*
** Les paramètres de la géométrie.
*
 
rint     = 1. ;
rext     = 2.1;
hh       = 1. ;
inj1     = 1.8;
inj2     = 1.9;
exi1     = 1.2;
exi2     = 1.3;
 
*
** Coefficients affectés aux couches d'Ekman et Stewartson.
*
 
 coefekma = 12. ;
 coefste1 = 0.8 ;
 coefste2 = 0.4 ;
 
*
** Densité du maillage.
*
 
npts = 14;
d1 = 1.0 / npts;
d2 = d1 / 1. ;
 
*
** Calcul des épaisseurs des couches d'Ekman et de Stewartson.
** Ces épaisseurs sont des ordres de grandeurs sans dimension.
*
 
ekma = nu / (omeg * (rext**2));
ddek = ekma**0.5;
stew = ekma**0.25;
 
*
** Modifications de ces épaisseurs en tenant compte des coefficents
** introduits.
*
 
dbek = coefekma * ddek * hh ;
fiek = hh * (1. - (coefekma * ddek)) ;
dste = rint * coefste1 * stew ;
fste = rext * coefste2 * stew ;
 
*
** Les sommets du bol
*
 
p1 = rint 0. ;
p2 = rext 0. ;
p3 = rext hh ;
p4 = rint hh ;
 
*
** Les milieux des cotés.
*
 
mil1 = ((rint + rext)/2.0);
mil2 = hh/2.0 ;
p6 = rext mil2 ;
p8 = rint mil2 ;
 
*
** Les points délimitant les couches d'Ekman.
*
 
p11 = rint dbek ;
p22 = rext dbek ;
p33 = rext fiek ;
p44 = rint fiek ;
 
*
** Les points d'entrée du fluide et leurs symétriques.
** Sont rajoutées les couches de Stewartson qui les 
** entourent.
*
 
ent1 = inj1 0. ;
ent2 = inj2 0. ;
sym1 = inj1 hh ;
sym2 = inj2 hh ;
en11 = (inj1 - fste) 0. ;
en22 = (inj2 + fste) 0. ;
sy11 = (inj1 - fste) hh ;
sy22 = (inj2 + fste) hh ;
 
*
** Les points de sortie du fluide et leurs symétriques.
** Sont rajoutées les couches de Stewartson qui les 
** entourent.
*
 
out1 = exi1 hh ;
out2 = exi2 hh ;
syo1 = exi1 0. ;
syo2 = exi2 0. ;
ou11 = (exi1 - dste) hh ;
ou22 = (exi2 + dste) hh ;
so11 = (exi1 - dste) 0. ;
so22 = (exi2 + dste) 0. ;
 
*
** Assemblage des bords du domaine.
*
 
bas1 = p1 'DROIT' dini d1 dfin d1 so11
          'DROIT' dini d1 dfin d2 syo1
          'DROIT' dini d2 dfin d2 syo2
          'DROIT' dini d2 dfin d1 so22
          'DROIT' dini d1 dfin d1 en11
          'DROIT' dini d1 dfin d2 ent1;
entr = ent1 'DROIT' dini d2 dfin d2 ent2 ;
bas2 = ent2 'DROIT' dini d2 dfin d1 en22
            'DROIT' dini d1 dfin d1 p2; 
 cdro = p2 'DROIT' dini d2 dfin d2 p22
           'DROIT' dini d2 dfin d1 p6 
           'DROIT' dini d1 dfin d2 p33
           'DROIT' dini d2 dfin d2 p3; 
hau2 = p3 'DROIT' dini d1 dfin d1 sy22
          'DROIT' dini d1 dfin d2 sym2 
          'DROIT' dini d2 dfin d2 sym1
          'DROIT' dini d2 dfin d1 sy11 
          'DROIT' dini d1 dfin d1 ou22
          'DROIT' dini d1 dfin d2 out2;
sor  = out2 'DROIT' dini d2 dfin d2 out1;
hau1 = out1 'DROIT' dini d2 dfin d1 ou11
            'DROIT' dini d1 dfin d1 p4; 
 cgau = p4 'DROIT' dini d2 dfin d2 p44
           'DROIT' dini d2 dfin d1 p8
           'DROIT' dini d1 dfin d2 p11
           'DROIT' dini d2 dfin d2 p1; 
 cot1 = bas1 'ET' entr 'ET' bas2;
 cot2 = cdro;
 cot3 = hau2 'ET' sor 'ET' hau1;
 cot4 = cgau;
 cnt  = cot1 'ET' cot2 'ET' cot3 'ET' cot4 ;
 
*
** On maille avec des rectangles.
*
 
mt= cot1 cot2 cot3 cot4 'DALLER' ;
'TASSER' mt ;
 
*
** On réoriente les éléments.
*
 
mt = 'ORIENTER' mt ;
 
*
** On crée le domaine associé au maillage.
*
Mmt = 'CHAN' mt 'QUAF' ;
$mt = 'MODE' Mmt 'NAVIER_STOKES' DISCR ;
 'DOMA' $mt 'IMPR';
mt = 'DOMA' $mt 'MAILLAGE' ;
 
*
************************ On défini les conditions sur les vitesses **************
*
 
*
** On se donne l'échelle de grandeur pour le tracage.
*
 
uref= 0.6 ;
 
*
** le terme gb défini la présence de terme source
** pour l'équation radiale et son absence dans
** l'équation axiale.
*
 
gb = (-1.0 0.0) ;
 
*
** Les conditions d'injection pour t>= 0.
** On adopte un profil qui permet de se raccorder
** aux éléments entourant l'entrée.
*
 
 corx = coor 1 mt ;
winj = winj * -4. * (corx - inj1) * (corx - inj2) / ((inj2 - inj1)**2);
 
*
** Les conditions initiales dans le domaine t < 0.
*
 
udom = 0. ;
vdom = 0. ;
wdom = 0. ;
 
*
****************** Mise en place du calcul numérique ********************
*
 
*
** On crée une table RX pour stocker les vitesses.
*
 
RX = TABLE ;
 
*
** On initialise les CHPOINTS utilisés.
*
 
RX.'UN' = 'KCHT' $mt vect SOMMET (udom wdom) ;
RX.'VT' = 'KCHT' $mt scal  SOMMET vdom ;
 
*
** Création de la table RV associé à la résolution des équations
** radiales et axiales avec 'NS' 'ET' azimutale avec 'TSCA'.
** On utilise la version Boussinesq pour avoir
** les termes sources aux sommets et non au centre.
*
 
RV = 'EQEX' $mt 'OPTI' ALE 'ITMA' 1 'ALFA' 0.8
   'ZONE' $mt 'OPER' 'NS' nu gb 'SR' 0. 'WN' 'INCO'  'UN' 
   'ZONE' $mt 'OPER' 'TSCA' nu 'WN' 'ST' 'INCO'  'VT'
      ;
 
  RV=EQEX RV OPTI EFM1 'CENTREE'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
  'ZONE' $MT      'OPER' 'DFDT' 1. 'VT' 'DELTAT' 'INCO' 'VT'
  ;
 
 
RV = 'EQEX' RV
   'CLIM' 'UN' 'UIMP' cnt                                    0.
          'UN' 'VIMP' (cgau et bas1 et bas2 et hau1 et hau2) 0.
          'UN' 'VIMP' entr                                   winj
   'CLIM' 'VT' 'TIMP' (cot1 et cot3)   0.
                    ;
 
*
** On charge la pression (méthode de Choleski). 
*
 
RVP = EQEX 'OPTI' 'EF' KPRESS
 'ZONE' $MT  OPER  KBBT  -1. betap  INCO 'UN' 'PRES'
;
 
    rvp.'METHINV'.TYPINV=1 ;
    rvp.'METHINV'.IMPINV=0 ;
    rvp.'METHINV'.NITMAX=300;
    rvp.'METHINV'.PRECOND=3 ;
    rvp.'METHINV'.RESID  =1.e-8 ;
    rvp.'METHINV' . 'FCPRECT'=100 ;
    rvp.'METHINV' . 'FCPRECI'=100 ;
 
  RV.'PROJ'     =RVP ;
 
RV.'INCO' = TABLE 'INCO';
 
*
** On initialise aussi les champs de vitesse.
*
 
RV.INCO.'PRES'  = 'KCHT' $mt scal KPRESS 0. ;
RV.INCO.'WN'    = 'KCHT' $mt vect SOMMET (0. 0.) ;
RV.INCO.'UN' = 'COPIER' RX.'UN';
RV.INCO.'VT' = 'COPIER' RX.'VT';
 
*
****************** Boucle itérative sur les pas de temps ***************
*
 
indice = 0 ;
nper =1;
 
*
** On définit le nombre d'iterations.
** Pour l'instant seul le cas 'NON' COMPLET est possible.
*
 
SI (NON COMPLET);
  NBITER = 10;
FINSI;
 
REPETER ITER NBITER;
 
indice = indice + 1 ;
 
*
** stocke les vitesses dans des variables intermédiaires.
*
 
un = RX.'UN' ;
vt = RX.'VT' ;
 
 
*
** Bouclage du calcul implicite des termes sources et du
** terme transportant.
*
 
nimpl = 0 ;
 
REPETER IMPLIC MAXIMP ;
 
nimpl = nimpl + 1;
 
*
** On charge les termes sources implicites.
*
 
uimpl = 'KCHT' $mt vect SOMMET RV.INCO.'UN';
vimpl = 'KCHT' $mt SCAL SOMMET RV.INCO.'VT';
sr   = (vimpl * vimpl / corx) + (2. * omeg * vimpl) ;
uscal  = 'EXCO' RV.INCO.'UN' ux scal;
st   = (-1.0 * uscal * vimpl / corx) - (2. * omeg * uscal)
      - (nu*vt/(corx*corx)) ;
stt = 'KCHT' $mt SCAL SOMMET st;
st = 'NOEL' $mt stt;
RV.INCO.'SR' = 'KCHT' $mt SCAL SOMMET sr;
RV.INCO.'ST' = 'KCHT' $mt SCAL CENTRE st;
 
*
** On réinitialise les champs de vitesse transportants.
*
 
RV.INCO.'WN' = 'KCHT' $mt vect SOMMET RV.INCO.'UN';
 
*
** On reinitialise les tables de vitesses non implicites
*
 
RV.INCO.'UN'    = 'KCHT' $mt vect SOMMET un ;
RV.INCO.'VT' = 'KCHT' $mt SCAL SOMMET vt;
 
 
*
** On execute le calcul.
*
 
EXEC RV ;
 
*
** On test la convergence de notre implicitation.
** Ce qui a un sens pour MAXIMP >= 2.
*
 
'SI' (MAXIMP > 1);
 
   err1 = CALCERR uimpl RV.INCO.'UN';
   err2 = CALCERR vimpl RV.INCO.'VT';
 
*
** On charge la premiere evolution.
*
 
   'SI' (nimpl EGA 1);
      errdebu = err1 ;
      errdebv = err2 ;
   'FINSI';
 
*
** On effectue le test.
*
 
   bool1 = nimpl > 1;
   bool2 = err1 &lt;EG (crit*errdebu);
   bool3 = err2 &lt;EG (crit*errdebv);
   bool = bool1 et bool2 et bool3;
   'SI' (bool);
      QUITTER IMPLIC;
   'FINSI';
 
   'SI' ((nimpl EGA MAXIMP) ET (NON bool) ET (nimpl > 1)) ;
***** Le calcul d'erreur n'a pas de sens pour le premier pas de temps. 
      'SI' (indice > 2);
      mess 'implicitation insuffisante ou critère trop severe';
      'FINSI';
   'FINSI';
'FINSI';
 
'FIN' IMPLIC ;
 
*
** Calcul de l'évolution des vitesses.
*
 
'SI' (CALCEV 'ET' (indice 'EGA' (nper*period)));
un1 = 'EXCO' un ux;
un2 = vt;
un3 = 'EXCO' un uy;
un11 = 'EXCO' RV.INCO.'UN' ux;
un22 = RV.INCO.'VT';
un33 = 'EXCO' RV.INCO.'UN' uy;
erru = CALCERR un1 un11 ;
errv = CALCERR un2 un22 ;
errw = CALCERR un3 un33 ;
 
*
** Stockage des erreurs en listes.
*
 
'SI' (nper EGA 1);
   convu = 'PROG' erru;
   convv = 'PROG' errv;
   convw = 'PROG' errw;
   liter = 'PROG' indice;
SINON ; 'SI' (nper > 1);
   convu = convu 'ET' ('PROG' erru);
   convv = convv 'ET' ('PROG' errv);
   convw = convw 'ET' ('PROG' errw); 
   liter = liter 'ET' ('PROG' indice);
   'FINSI';
'FINSI';
nper = nper + 1;
'FINSI';
 
*
** On sauve les expressions obtenues pour les vitesses.
*
 
RX.'UN' = 'COPIER' RV.INCO.'UN' ;
vt = RV.INCO.'VT';
vtn = 'KCHT' $mt SCAL SOMMET vt;
RX.'VT' = 'COPIER' vtn ;
 
'FIN' ITER ;
 
*
******************** Test de l'erreur ************************
*
 
*
** On n'autorise pas de test long donc meme test
** dans les deux cas.
*
si vrai ;
erdif1 = 'ABS' (erru - 2.696e-2);
bool1 = erdif1 &lt;EG 1.e-4;
erdif2 = 'ABS' (errv - 6.285e-2);
bool2 = erdif2 &lt;EG 1.e-4;
erdif3 = 'ABS' (errw - 1.0465e-2);
bool3 = erdif3 &lt;EG 1.e-4;
bool4 = bool1 ET bool2 ET bool3;
'SI' ('NON' bool4);
   'MESSAGE' erdif1 erdif2 erdif3;
   ERREUR 5;
'FINSI';
finsi ;
*
**************** Traitement des résultats ******************************
*
 
SI ('EGA' GRAPH 'O');
 
*
** Tracé pour les vitesses dans le plan (rz).
*
 
   ung1='VECTEUR' un uref ux uy jaune ;
   'TRACER' ung1 mt ;
*
** Tracé des lignes de courant.
*
 
Mcot1 = 'CHAN' cot1 'QUAF';
Mcot2 = 'CHAN' cot2 'QUAF';
Mcot3 = 'CHAN' cot3 'QUAF';
Mcot4 = 'CHAN' cot4 'QUAF';
'ELIM' (Mmt et Mcot1 et Mcot2 et Mcot3 et Mcot4) 1.e-5;
 
$cot1 = 'MODE' Mcot1 'NAVIER_STOKES' DISCR ;
$cot2 = 'MODE' Mcot2 'NAVIER_STOKES' DISCR ;
$cot3 = 'MODE' Mcot3 'NAVIER_STOKES' DISCR ;
$cot4 = 'MODE' Mcot4 'NAVIER_STOKES' DISCR ;
cot1 = 'DOMA' $cot1 'MAILLAGE' ;
cot2 = 'DOMA' $cot2 'MAILLAGE' ;
cot3 = 'DOMA' $cot3 'MAILLAGE' ;
cot4 = 'DOMA' $cot4 'MAILLAGE' ;
 
un1 = 'EXCO' un ux;
un2 = 'EXCO' un uy;
un1 = -1. * corx * un1;
un3 = corx * un2;
un1 = 'KCHT' $mt SCAL SOMMET un1;
un2 = 'KCHT' $mt SCAL SOMMET un2;
un3 = 'KCHT' $mt SCAL SOMMET un3;
uncz = 'NOEL' $mt un2;
uncz = 'KCHT' $mt SCAL CENTRE uncz;
rt2d= 'KOPS' un 'ROT' $mt ;
 corx = 'KCHT' $mt SCAL SOMMET corx;
 corm = 'NOEL' $mt corx;
unn1 = 'KCHT' $cot1 SCAL SOMMET un1;
unn2 = 'KCHT' $cot2 SCAL SOMMET (-1 * un3);
unn3 = 'KCHT' $cot3 SCAL SOMMET (-1 * un1);
unn4 = 'KCHT' $cot4 SCAL SOMMET un3;
unn1 = 'NOEL' $cot1 unn1;
unn1 = 'KCHT' $cot1 SCAL CENTRE unn1;
unn2 = 'NOEL' $cot2 unn2;
unn2 = 'KCHT' $cot2 SCAL CENTRE unn2;
unn3 = 'NOEL' $cot3 unn3;
unn3 = 'KCHT' $cot3 SCAL CENTRE unn3;
unn4 = 'NOEL' $cot4 unn4;
unn4 = 'KCHT' $cot4 SCAL CENTRE unn4;
sw = 2*uncz - (corm * rt2d) ;
rk = 'EQEX' $MT 'OPTI' 'EF' 'IMPL'
ZONE  $mt OPER 'LAPN' 1.  INCO 'PSI'
ZONE  $mt OPER 'FIMP' sw  INCO 'PSI'
ZONE  $cot1 OPER 'FIMP' unn1 INCO 'PSI'
ZONE  $cot2 OPER 'FIMP' unn2 INCO 'PSI'
ZONE  $cot3 OPER 'FIMP' unn3 INCO 'PSI'
ZONE  $cot4 OPER 'FIMP' unn4 INCO 'PSI'
  CLIM 'PSI' 'TIMP' cot4 0. ;
rk.'INCO'.'PSI' = 'KCHT' $mt scal sommet 0. ;
EXEC rk ;
psi = rk.'INCO'.'PSI' ;
'OPTION' ISOV LIGNE;
'TRACER' psi mt ('CONTOUR' (mt)) ;
 
*
** Tracé pour la vitesse ortho-radiale en 3D.
*
 
vtx = 'KCHT'  $mt SCAL SOMMET 0;
vtx = 'NOMC'  'UX' vtx 'NATU' DISCRET;
vty = 'KCHT'  $mt SCAL SOMMET 0.;
vty = 'NOMC'  'UY' vty 'NATU' DISCRET;
 
* 
*  La vitesse est inversée car le maillage
*  est en réalité dans le plan rz et non xy.
*  Pour recreer la véritable orientation, on change
*  le signe de la vitesse ortho-radiale.
*
 
mvt = -1. *  vt;
vtz = 'NOMC' 'UZ' mvt 'NATU' DISCRET;
'OPTION' 'DIME' 3;
uuu = 'KCHT' $mt vect SOMMET (vtx ET vty ET vtz);
ung2 = 'VECTEUR' uuu uref ux uy uz jaune;
'TRACER' ung2 mt;
 
*
** tracé de l'évolution des vitesses en fonction du
** pas de temps. Il n'a rien si le nombre d'iterations
** est insuffisant.
*
 
'SI' (CALCEV);
   EVU    = 'EVOL' 'MANU' 'ITERATION' liter
                     'EVOLUTION EN UR' convu;
   EVV    = 'EVOL' 'MANU' 'ITERATION' liter
                     'EVOLUTION EN UTETA' convv;
   EVW    = 'EVOL' 'MANU' 'ITERATION' liter
                     'EVOLUTION EN UZ' convw;
   'DESSIN' evu logy;
   'DESSIN' evv logy;
   'DESSIN' evw logy;
'FINSI';
'FINSI';
 
 
*
** Sauvegarde des données.
*
 
*REPIX RV;
*fsauv = '  ' ;
*OPTI SAUV fsauv; SAUV;
 
'FIN' ;