* 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 - All rights reserved.
Disclaimer