Télécharger linekman.dgibi

Retour à la liste

* fichier :  linekman.dgibi
************************************************************************
************************************************************************
****************************************************************************
*
*
*         Date : 19/3/97
*
*         Description: Cas test simulant l'écoulement dans une région
*                      limitée par une plaque horizontale infinie en
*                      rotation autour d'un axe perpendiculaire.
*
*                      l'ensemble du fluide tourne à une vitesse Omeg
*                      tandis que la plaque tourne à la vitesse
*                      (Omeg + ddom).
*
*                      Les calculs étant effectués en axi-symétrique,
*                      seule la moitié de la région est représentée.
*                      L'axe de rotation constitue l'axe des z de notre
*                      étude. La plaque inférieure a pour équation z = 0.
*                      
*                      La région étant par ailleurs infinie, on est bien
*                      entendu obligé de la limiter dans notre étude.
*                      On met les conditions adéquates au bord.
*                      
*
*         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): linear ekman layer.
*
*                      Cet écoulement admet une solution analytique.
*
*
*         Auteur     : Gilles Bernard-Michel
*
*         Révision   : Ce fichier diffère du précédent car il n'y a
*                      pas d'implicitation. 
*
*****************************************************************************
*
 
*
** Choix des éléments
*
 
option dime  2 elem qua8 mode axis;
 
*
************************ Quelques procédures ********************************
*
 
*
** Erreur Linfini entre deux Champoints.
*
 
DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  vmax = MAXI vit 'ABS' + 1.e-13; 
  dimax = MAXI (vitp1 - vit) 'ABS' ;
  err = (dimax /vmax) ;
FINPROC err ;
 
*
** Calcul de la norme L2 d'un champ.
*
 
DEBPROC CALCNOL2 vit*'CHPOINT' mod*'MMODEL' 
  mt*'MAILLAGE';
  opti dime 2 mode PLAN;
  vgau = CHANGE 'CHAM' vit mt 'RIGIDITE';
  usol2 = vgau * vgau;
  us2 = INTG mod usol2;
  opti dime 2 mode axis;
FINPROC us2 ;
 
*
** Calcul de l'erreur L2 pour une composante de la vitesse.
*
 
DEBPROC CALCURR vit*'CHPOINT' vitsol*'CHPOINT' mod*'MMODEL' 
  mt*'MAILLAGE';
  opti dime 2 mode PLAN;
  vgau = CHANGE 'CHAM' vit mt 'RIGIDITE';
  vsolgau = CHANGE 'CHAM' vitsol mt 'RIGIDITE';
  udiff = vgau - vsolgau;
  udif2 = udiff*udiff;
  ud2 = INTG mod udif2;
  usol2 = vsolgau * vsolgau;
  us2 = INTG mod usol2;
  us2 = us2 + 1.e-14;
  err = ud2 / us2;
  err = err**0.5;
  opti dime 2 mode axis;
FINPROC err ;
 
*
** Erreur L2 pour l'ensemble des composantes de la vitesse.
*
 
DEBPROC CALCERTO vit1*'CHPOINT' vit2*'CHPOINT' vit3*'CHPOINT'
 vi1*'CHPOINT' vi2*'CHPOINT' vi3*'CHPOINT' mod*'MMODEL'
 mt*'MAILLAGE';
vit1 = vit1 - vi1;
vit2 = vit2 - vi2;
vit3 = vit3 - vi3;
nom1 = CALCNOL2 vit1 mod mt;
nom2 = CALCNOL2 vit2 mod mt;
nom3 = CALCNOL2 vit3 mod mt;
nomi = nom1 + nom2 + nom3;
deno1 = CALCNOL2 vi1 mod mt;
deno2 = CALCNOL2 vi2 mod mt;
deno3 = CALCNOL2 vi3 mod mt;
denom = deno1 + deno2 + deno3 + 1.e-14;
err = nomi/denom;
err = err**0.5;
FINPROC err;
 
*
************************ On définit les propriétés physiques ***************
*
** viscosité dynamique mu/rho.
*
 
nu =1.E-1 ;
 
*
** La vitesse de rotation initiale.
*
 
omeg = 1 ;
 
*
** Le nombre de Rossby.
*
 
Rosb = 0.0 ;
 
*
** L'impulsion.
*
 
ddom = 1 ;
 
*
********************* Les options utilisateur ******************************
*
 
*
** Critère de convergence implicite.
*
 
crit = 0.07 ;
 
*
** Nombre max de boucles d'implicitation.
*
 
MAXIMP = 1 ;
 
*
** Déformation du maillage.
*
 
DEFORM = FAUX;
 
*
** Amplitude des déformation (nombre arbitraire, 1 donne une déformation
** acceptable pour le maillage).
*
 
ampdef = 1. ;
 
*
** Sortie graphique.
*
 
GRAPH = 'N';
 
*
** Test long ou court.
*
 
COMPLET = FAUX;
 
*
** On calcule les erreurs L2.
*
 
ERRCALC = VRAI;
 
*
** Periode de calcul des erreurs L2.
*
 
period = 5;
 
*
********************* On traite les aspects géométriques ********************
*
 
*
** Les paramètres de la géométrie.
*
 
rint     = 0.5 ;
rext     = 1.;
hh       = 1. ;
 
*
** Densité du maillage.
*
 
SI (NON COMPLET);
  nptx = 3;
  nptz = 4;
  NBITER = 60;
SINON;
  nptx = 30;
  npty = 60;
  NBITER = 15000;
FINSI;
 
*
** Constantes caractérisant les différentes couches limites.
*
 
ekma = nu / omeg;
ddek = ekma**0.5;
 
*
** Les sommets du bol
*
 
p1 = rint 0. ;
p2 = rext 0. ;
p3 = rext hh ;
p4 = rint hh ;
 
*
** Assemblage des bords du domaine.
*
 
bas  = d nptx p1 p2;
cdro = d nptz p2 p3;
haut = d nptx p3 p4;
cgau = d nptz p4 p1;
 
*
** On maille.
*
 
cnt = bas et cdro et haut et cgau;
mt= bas cdro haut cgau daller ;
tass mt ;
 
*
** On trace le maillage obtenu.
*
 
*trace mt ;
*opti donn 5;
 
*
** On réoriente les éléments.
*
 
mt = orienter mt ;
 
*
** On déforme le maillage avec un bruit gaussien.
** Sous réserve que l'option soit selectionnée.
*
 
SI (DEFORM);
   mmt = CHANGE mt POI1;
   ccnt= CHANGE cnt POI1;
   inmt= diff mmt ccnt;
   denn = (0.03 * ampdef) / nptx;
   brbl1 = BRUIT BLAN GAUS 0 denn inmt;
   brbl2 = BRUIT BLAN GAUS 0 denn inmt;
   brbl = (NOMC 'UR' brbl1) ET (NOMC 'UZ' brbl2);
   DEPLAC mt PLUS brbl;
FINSI;
 
*
** On crée le domaine associé au maillage.
*
 
$mt   = DOMA mt 'MACRO' impr;
mt = $mt.maillage;
 
*
** On crée un modèle pour pouvoir intégrer des champs.
*
MODL1 = MODE mt 'MECANIQUE';
 
*
************************ On défini les conditions sur les vitesses **************
*
 
*
** On se donne l'échelle de grandeur pour le tracage.
*
 
uref= 0.2 ;
 
*
** 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) ;
 
*
** Calcul de la solution exacte.
*
 
corx = coor 1 mt ;
cory = coor 2 mt ;
zba = cory / ddek;
dzba = zba * 180 / 3.141592654;
mzba = -1*zba;
urdo = ddom * corx * (EXP mzba) * (SIN dzba) ;
udom = NOMC 'UX' (kcht $mt SCAL SOMMET urdo) 'NATU' 'DISCRET';
wzdo = -1 * ddom * ddek * (1 - 
       ((EXP mzba) * ((SIN dzba) + (COS dzba))));
wdom = NOMC 'UY' (kcht $mt SCAL SOMMET wzdo) 'NATU' 'DISCRET';
upla = kcht $mt VECT SOMMET (udom ET wdom);
vdom = ddom * corx * (EXP mzba) * (COS dzba) ;
pdom = -2. * ddom * ekma * (EXP mzba) * (SIN dzba) ;
zb = 1. / ddek;
dzb = zb * 180 / 3.141592654;
mzb = -1 * zb;
pdelt = 2. * ddom * ekma * (EXP mzb) * (SIN dzb) ;
pdom = pdom + pdelt;
 
*
****************** 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 à t=0.
*
 
RX.'UN' = kcht $mt VECT SOMMET (0. 0.) ;
RX.'VT' = kcht $mt SCAL SOMMET 0. ;
 
*
** Création de la table RV associé à la résolution des équations
** radiales et axiales. 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
   'CLIM' 'UN' 'UIMP' cnt  urdo
          'UN' 'VIMP' (cgau et bas et cdro)  wzdo
          'VT' 'TIMP' cnt  vdom
                    ;
 
*
** On charge la pression (méthode de Choleski). 
*
 
RVP = EQPR $mt KTYPI 1
         ZONE $mt  'OPER' 'PRESSION' 0.
*         PIMP 0.
      ;
 
RV.'INCO' = TABLE 'INCO';
RV.PRESSION = RVP;
 
*
** Les champs transportants sont nuls à tout t>=0.
** On initialise aussi les champs de vitesse.
*
 
RV.INCO.'WN' = kcht $mt VECT SOMMET (0. 0.) ;
RV.INCO.'UN' = copi RX.'UN';
RV.INCO.'VT' = copi RX.'VT';
 
*
****************** Boucle itérative sur les pas de temps ***************
*
 
indice = 0 ;
nper = 1;
 
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   = 2. * omeg * vimpl ;
uscal  = exco RV.INCO.'UN' ux scal;
st   = -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 reinitialise les tables de vitesses non implicites
*
 
RV.INCO.'UN'    = kcht $mt VECT SOMMET un ;
RV.INCO.'VT' = kcht $mt SCAL SOMMET vt;
RV.PRESSION =RVP ;
 
*
** On execute le calcul des équations radiales et axiales.
*
 
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 <EG (crit*errdebu);
   bool3 = err2 <EG (crit*errdebv);
   bool = bool1 et bool2 et bool3;
   SI (bool);
      QUITTER IMPLIC;
   FINSI;
 
   SI ((nimpl EGA MAXIMP) ET (NON bool) ET (nimpl > 1)) ;
      mess 'erreur d' 'implicitation' err1 err2;
   FINSI;
FINSI;
 
 
FIN IMPLIC ;
 
*
** Calcul de l'erreur par rapport à la solution analytique.
*
 
SI (ERRCALC ET (indice EGA (nper*period)));
   un1 = exco RV.INCO.'UN' ux;
   un2 = RV.INCO.'VT';   
   un3 = exco RV.INCO.'UN' uy;
   pn = RV.PRESSION.PRESSION;
   pnn = elno $mt pn;
   pnn = kcht $mt SCAL SOMMET pnn;
   errl2u = CALCURR un1 urdo MODL1 mt;
   errl2v = CALCURR un2 vdom MODL1 mt;
   errl2w = CALCURR un3 wzdo MODL1 mt;
   errl2p = CALCURR pnn pdom MODL1 mt;
   errtot = CALCERTO un1 un2 un3 urdo vdom wzdo MODL1 mt;
 
*mess 'pas de temps' indice ;
*mess 'valeur des derivees' errdebu errdebv;
*mess 'erreur l2' errtot errl2u errl2v errl2w errl2p;
 
*
** Stockage des erreurs en listes.
*
 
   SI (nper EGA 1);
      lerl2u = PROG errl2u;
      lerl2v = PROG errl2v;
      lerl2w = PROG errl2w;
      lerl2p = PROG errl2p;
      lerl2t = PROG errtot;
      liter = PROG indice;
   SINON; SI (nper > 1);
      lerl2u = lerl2u ET (PROG errl2u);
      lerl2v = lerl2v ET (PROG errl2v);
      lerl2w = lerl2w ET (PROG errl2w);
      lerl2p = lerl2p ET (PROG errl2p); 
      lerl2t = lerl2t ET (PROG errtot); 
      liter = liter ET (PROG indice);
      FINSI;
   FINSI;
   nper = nper + 1;
FINSI;
 
*
** On sauve les expressions obtenues pour les vitesses.
*
 
RX.'UN' = copi RV.INCO.'UN' ;
vt = RV.INCO.'VT';
vtn = kcht $mt SCAL SOMMET vt;
RX.'VT' = copi vtn ;
 
FIN ITER ;
 
*
******************** Test de l'erreur ************************
*
 
SI (NON COMPLET);
   erdif1 = ABS (errtot - 2.51e-2);
   mess ' erdif1 = ' erdif1 ' errtot = ' errtot;
   bool1 = erdif1 <EG 1.e-4;
   erdif2 = ABS (errl2p - .69732);
   mess ' erdif2 = ' erdif2 ' errl2p = ' errl2p;
   bool2 = erdif2 <EG 1.e-4;
   bool3 = bool1 ET bool2;
   SI (NON bool3);
      mess 'houra';
      ERREUR 5;
   FINSI;
SINON;
   bool1 = errtot < 0.02e-2;
   bool2 = errl2p < 0.01;
   bool3 = bool1 ET bool2;
   SI (NON bool3);
      ERREUR 5;
   FINSI;
FINSI;
 
*
**************** Traitement graphique ***************************************
*
 
SI ('EGA' GRAPH 'O');
 
*
** Tracé pour les vitesses dans le plan (rz).
*
 
   ung1=vect un uref ux uy jaune ;
   trac ung1 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.
 vtz = NOMC 'UZ' (-1.*vt) 'NATU' DISCRET;
 option dime 3;
 uuu = kcht $mt VECT SOMMET (vtx ET vty ET vtz);
 ung2 = vect uuu uref ux uy uz jaune;
 trac ung2 mt;
 
*
** tracé de l'erreur d'implicitation en fonction du
** pas de temps.
*
 
   UL2    = EVOL 'MANU' 'ITERATION' liter 
                     'ERREUR L2 EN UR' lerl2u;
   VL2    = EVOL 'MANU' 'ITERATION' liter 
                     'ERREUR L2 EN UTETA' lerl2v;
   WL2    = EVOL 'MANU' 'ITERATION' liter 
                     'ERREUR L2 EN UZ' lerl2w;
   PL2    = EVOL 'MANU' 'ITERATION' liter 
                     'ERREUR L2 EN P' lerl2p;
   TL2    = EVOL 'MANU' 'ITERATION' liter 
                     'ERREUR L2 VITESSE' lerl2t;
   dessin ul2 logy;
   dessin vl2 logy;
   dessin wl2 logy;
   dessin tl2 logy;
   dessin pl2 logy;
FINSI;
 
*
** Sauvegarde des données.
*
 
*REPIX RV;
*REPIX RA;
*fsauv = '/test4/?????.sauv' ;
*OPTI SAUV fsauv; SAUV;
 
FIN;
 
 
 
 
 
 
 

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