* fichier : linekmanimp.dgibi
****************************************************************************
*
*
* 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 ;
© Cast3M 2003 - All rights reserved.
Disclaimer