* fichier :  linekmanimp.dgibi
************************************************************************
* Section : Fluides Transport
* Section : Fluides Transitoire
************************************************************************
****************************************************************************
*
*
*         Date : 06/01/99
*
*         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.
*
*                      La solution est permanente.
*
*                      Cet écoulement admet une solution analytique car
*                      les termes de transport sont négligés.
*
*  Type de calcul    : PERMANENT, PAS DE TERME DE TRANSPORT
*
*  Type de méthode   : IMPLICITE
*
*  Type d'éléments   : ELEMENTS QUA9 CENTREP1
*
*  Opérateurs utilisés : KBBT LAPN MDIA
*
*
*
*         Auteurs    : Gilles Bernard-Michel
*                    : Daniel Guilbaud pour la présente version implicite
*
*****************************************************************************
*

*
********************* Les options utilisateur ******************************
*

*
** Déformation du maillage.
*

DEFORM = FAUX ;

*
** Sortie graphique.
*

GRAPH = FAUX ;

*
** Test long ou court.
*

COMPLET = FAUX ;

**************************************************************************

option dime  2 elem qua8 mode axis;

*
************************ 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 ;

*
********************* 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 = 10;
  nptz = 20;
SINON;
  nptx = 30;
  nptz = 60;
FINSI;

*
** Constantes caractérisant les différentes couches limites.
*

ekma = nu / (abs 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 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);

** Amplitude des déformation (nombre arbitraire, 1 donne une déformation
** acceptable pour le maillage).

   ampdef = 1. ;

   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 trace le maillage obtenu.
*

SI GRAPH ;
  titre 'tracé du maillage' ;
  trace mt ;
FINSI ;

*
** On crée le domaine associé au maillage.
*

* transformation des éléments en qua9

£mt    = CHANGER mt    QUAF ;
£haut  = CHANGER haut  QUAF ;
£cnt   = CHANGER cnt   QUAF ;
£cgau  = CHANGER cgau  QUAF ;
£bas   = CHANGER bas   QUAF ;
£cdro  = CHANGER cdro  QUAF ;




ELIM 1.e-6 (£mt et £cnt et £haut et £cgau et £bas et £cdro) ;

* formulation du modèle NAVIER_STOKES

$mt    = MODE £mt    'NAVIER_STOKES' QUAF ;

mmt   = doma $mt    maillage ;

*
************************ On défini les conditions sur les vitesses **************
*

*
** On se donne l'échelle de grandeur pour le tracage.
*

uref= 0.2 ;

*
** Calcul de la solution exacte.
*
corx = coor 1 mmt ;
cory = coor 2 mmt ;
zba = cory / ddek;
dzba = zba * 180. / PI ;
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. / PI ;
mzb = -1 * zb;
pdelt = 2. * ddom * ekma * (EXP mzb) * (SIN dzb) ;
pdom = pdom + pdelt;

*
****************** Mise en place du calcul numérique ********************
*


* equations de qdm sur ur et uz + div V = 0

RV = EQEX $mt
OPTI 'EF'  'IMPL' 'CENTREP1'  'SUPG'
* il faut opti CENTREP1 pour KBBT
ZONE $mt   OPER KBBT (+1.)                INCO 'UN' 'PN'
ZONE $mt    OPER LAPN nu                  INCO 'UN'
ZONE $mt   OPER MDIA 'cour'               INCO 'VT' 'UN'
;

* equation scalaire sur ut

RV = EQEX RV
OPTI 'EF'  'IMPL' 'SUPG'
ZONE $mt    OPER LAPN nu                  INCO 'VT'
ZONE $mt    OPER MDIA 'nusr'              INCO 'VT'
ZONE $mt    OPER MDIA 'cout'              INCO 'UN' 'VT'
;

* conditions aux limites

RV = EQEX RV
   'CLIM' 'UN' 'UIMP' £cnt  urdo
          'UN' 'VIMP' (£cgau et £bas et £cdro)  wzdo
          'VT' 'TIMP' £cnt  vdom
;


*                      ________________
*
*                       INITIALISATION
*                      ________________

RV.INCO = TABLE INCO ;

RV.INCO.'UN'   = KCHT $mt   VECT SOMMET (1.e-8 1.e-8) ;
RV.INCO.'VT'   = KCHT $mt   SCAL SOMMET 1.e-8 ;
RV.INCO.'PN'   = KCHT $mt   SCAL CENTREP1 0. ;
RV.INCO.'cour' = KCHT $mt   VECT SOMMET ( (-2.*omeg)  0. ) ;
RV.INCO.'cout' = KCHT $mt   VECT SOMMET ( (+2.*omeg) 0. ) ;

* terme -nu/r2

coorx = 'KCHT' $mt   scal SOMMET corx ;
nusr2 = KOPS ( +1. * nu ) '/' ( KOPS coorx '*' coorx ) ;

RV.INCO.'nusr' = KCHT $mt   SCAL SOMMET nusr2 ;

*                         __________
*
*                           CALCUL
*                         __________

EXEC RV ;


*
** Calcul de l'erreur par rapport à la solution analytique.
*

urdo = KCHT $mt   SCAL SOMMET urdo ;
wzdo = KCHT $mt   SCAL SOMMET wzdo ;
vdom = KCHT $mt   SCAL SOMMET vdom ;
pdom = KCHT $mt   SCAL SOMMET pdom ;


UNR = EXCO UX RV.INCO.'UN' ;
UNR = KCHT $mt   SCAL SOMMET UNR ;
ERUNR = KOPS ( KOPS UNR '-' urdo) '/' ( (MAXI (ABS UNR) ) + 1.e-13) ;
ERRUNR = MAXI ( ABS ERUNR ) ;
MESS 'ERRUNR = ' ERRUNR ;

UNZ = EXCO UY RV.INCO.'UN' ;
UNZ = KCHT $mt   SCAL SOMMET UNZ ;
ERUNZ = KOPS ( KOPS UNZ '-' wzdo) '/' ( (MAXI (ABS UNZ) ) + 1.e-13) ;
ERRUNZ = MAXI ( ABS ERUNZ ) ;
MESS 'ERRUNZ = ' ERRUNZ ;

VT = RV.INCO.'VT' ;
ERVT = KOPS ( KOPS VT '-' vdom) '/' ( (MAXI (ABS VT) ) + 1.e-13) ;
ERRVT = MAXI ( ABS ERVT ) ;
MESS 'ERRVT = ' ERRVT ;

PNN = ELNO $mt (RV.'INCO'.'PN') CENTREP1 ;
pnn = kcht $mt SCAL SOMMET PNN;
ERPN = KOPS ( KOPS pnn '-' pdom) '/' ( (MAXI (ABS pnn) ) + 1.e-13) ;
ERRPN = MAXI ( ABS ERPN ) ;
MESS 'ERRPN = ' ERRPN ;

SI ( DEFORM ) ;
  erv = 2.e-3 ;
  erp = 1.e-1 ;
SINON ;
  erv = 1.e-5 ;
  erp = 5.e-2 ;
FINSI ;

BOOL1 = ERRUNR <EG erv ;
BOOL2 = ERRUNZ <EG erv ;
BOOL3 = ERRVT  <EG erv ;
BOOL4 = ERRPN  <EG erp ;

BOOLT = BOOL1 et BOOL2 et BOOL3 et BOOL4 ;

SI BOOLT ;
  MESS 'TEST CORRECT' ;
SINON ;
  MESS 'TEST FAUX' ;
  ERREUR 5 ;
FINSI ;

*
**************** Traitement graphique ***************************************
*

SI GRAPH ;

*
** Tracé pour les vitesses dans le plan (rz).
*
   UN = RV.INCO.'UN' ;
   ung1=vect un uref ux uy jaune ;
titre 'vitesse VR VZ du calcul numérique' ;
   trac ung1 £mt titre 'vitesse VR VZ du calcul numérique' ;
titre 'vitesse VT du calcul numérique' ;
   trac vt £mt ;
titre 'pression numérique' ;
   trac pnn £mt ;

   upla1=vect upla uref ux uy roug ;
titre 'vitesse VR VZ du calcul analytique' ;
   trac upla1 £mt titre 'vitesse VR VZ du calcul analytique' ;
titre 'vitesse VT du calcul analytique' ;
   trac vdom £mt ;
titre 'pression analytique' ;
   trac pdom £mt ;

titre 'difference relative de vitesse UR' ;
   trac ERUNR £mt ;

titre 'difference relative de vitesse UZ' ;
   trac ERUNZ £mt ;

titre 'difference relative de vitesse VT' ;
   trac ERVT  £mt ;

titre 'difference relative de pression' ;
   trac ERPN  £mt ;

*
** Tracé pour la vitesse ortho-radiale en 3D.
*
    VT = RV.INCO.'VT' ;
 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;
 titre 'vitesse azimutale 3D' ;
 trac ung2 mt titre 'vitesse azimutale 3D' ;

FINSI ;

fin ;