* fichier :  Marangoni1.dgibi
*                      _____________________
*
*                       CAS TEST EFFET MARANGONI
*             UTILISATION DE LA LOI DE CHAN, CHEN, MAZUMDER
*           Journal of Heat Transfert, vol 110, 1988
*                      _____________________
*                      OPTIONS TESTEES PLAN MACRO

***************************************************************
** DESCRIPTION DU TEST :
**
** Ce test a pour but d'étudier l'opérateur toimp au travers de l'effet
** Marangoni.
**
** Une poche de metal liquide est soumise au bombardement d'un flux de
** chaleur surfacique ayant une loi parabolique. La loi de CHAN décrit
** ce type d'écoulement. Cette loi  donne le profil de température du
** liquide ainsi que le profil des vitesses a la surface de la poche
** liquide.
**
** Lors de ce test, nous imposons directement sur la surface du bain de
** liquide le profil de température obtenue grâce a la loi de CHAN. Ce
** profil de température permet de calculer les forces de Marangoni, seul
** terme moteur de l'écoulement. Le but est de retomber sur le profil
** théorique des vitesses donné par la loi de CHAN.
**
** Ref : Journal of Heat Transfert, vol 110, 1988
**       "Asymptotic solution for thermocapillary flow at high and low
**        Prandtl numbers due to concentrated surface heating."
**       CHAN, CHEN, MAZUMDER.
**
***************************************************************

***********************
** OPTIONS A TESTER :
***********************
option DIME 2 ELEM qua8 MODE plan ;   MOAXI = FAUX ;
*option DIME 2 ELEM qua8 MODE axis;   MOAXI = VRAI ;
*MACRO = QUAF ;

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

*DIMPR = VRAI ;
DIMPR = FAUX ;
GRAPH = FAUX ;
COMPLET= FAUX ;
option trac psc ;

*                 _________________________________
*
*                  différents parametres du calcul
*                 _________________________________

* Nombre de mailles
NMAIL = 10 ;

* nombre maximal d'itérations
NITER = 8  ;
TERR=0.6   ;
'SI' COMPLET ;NITER = 50 ; TERR=0.05 ; 'FINSI' ;

* coefficient de relaxation pour le champ de vitesse
ALFAU = 0.5 ;

* précision requise
PREC = 1.e-3 ;

*   _______________________________________________________
*
*    géométrie du creuset : tronc de cone (dimension en m)
*   _______________________________________________________
*

* rayon du bas du domaine
rbas  = 0.0199  ;

* rayon du haut du domaine
rhaut = 0.020  ;

* hauteur de liquide
hliq  = 0.005 ;

* épaisseur de la couche limite
eclm = 1.e-4 ;

* rayon moyen du faisceau électonique
rbe   = 0.010  ;

*  __________________________________________________________________
*
*   Caractéristiques du faisceau d'électrons à répartition d'énergie
*   parabolique : puissance PCANON et "rayon" RFOCAL
*  __________________________________________________________________

PCANON = 1.e3 ;
RFOCAL = rbe  ;


*                 ______________________
*
*                  constantes physiques
*                 ______________________

* constante des gaz parfaits
RG = 8.32 ;

* constante de STEPHAN-BOLTZMAN
SIGMA = 5.67E-8 ;



*
*         __________________________________________________
*
*          propriétés physiques du fluide incompressible
*         __________________________________________________

* viscosité dynamique (kg/m/s)
mu = 1.3e-3       ;

* masse volumique (kg/m3)
rho0 = 6600. ;

* conductivité thermique  (W/m/K)
lambda = 5. ;

* chaleur spécifique (J/kg/K)
Cp = 270 ;

* tension de surface (N/m)
sig0 = 0.34 ;

* 1/sig0 dsig0/dT (1/K)
gam0 = 5.88e-4 ;

* viscosité cinématique (m2/s)
nu = mu / rho0 ;

* effet MARANGONI
lev =  1.  * gam0 * sig0 / rho0 ;

*                     ______________________
*
*                      CREATION DU MAILLAGE
*                     ______________________

p0  = 0.    0.          ;
p1  = rhaut  0.          ;
p2  = rbas   (-1. * hliq) ;
p3  = 0.    (-1. * hliq) ;

p01 = rbe    0.          ;
p32 = rbe    (-1. * hliq) ;

p03 = 0.    (-10. * eclm) ;
p12 = rhaut  (-10. * eclm) ;

d0 = rhaut / NMAIL ;
d1 = rhaut / NMAIL ;
d2 = rhaut / NMAIL ;
surli0 = (p01 droi dini d1 dfin d1 p0 ) coul roug ;
surli1 = (p1  droi dini d2 dfin d1 p01) coul roug ;
surli = surli1 et surli0 ;

d0 = hliq / 9 / NMAIL ;
d1 = hliq / 2 / NMAIL ;
d2 = hliq * 3 / NMAIL ;
smm0 = (p0  droi dini d0 dfin d1 p03) coul jaun ;
smm1 = (p03 droi dini d1 dfin d2 p3 ) coul jaun ;
smm = smm1 et smm0 ;

scon10 = (inve (smm0 proj 'CYLI' (1. 0.) 'DROI' p1 p2) ) coul vert ;
scon11 = (inve (smm1 proj 'CYLI' (1. 0.) 'DROI' p1 p2) ) coul vert ;
scon1 = scon11 et scon10 ;

dz0 = (smm0 proj 'CYLI' (1. 0.) 'DROI' p01 p32) coul bleu ;
dz1 = (smm1 proj 'CYLI' (1. 0.) 'DROI' p01 p32) coul bleu ;
dz = dz0 et dz1 ;

zc = (-1. * (rhaut - rbe ) * hliq ) / ( rhaut - rbas) ;

scon21= (inve (surli1 proj 'CONI' (rbe zc) 'DROI' p2 p3) ) coul vert ;
scon20= (inve (surli0 proj 'CYLI' (0. -1.) 'DROI' p2 p3) ) coul vert ;
scon2 = scon20 et scon21 ;

dr1 = (surli1 proj 'CONI' (rbe zc) 'DROI' p12 p03) coul rose ;
dr0 = (surli0 proj 'CYLI' (0. -1.) 'DROI' p12 p03) coul rose ;
dr = dr1 et dr0 ;

scon = scon1 et scon2 ;

elim 1.e-6 ( surli1 et surli0 et scon10 et scon11 et scon20
           et scon21 et dr0 et dr1 et dz0 et dz1 et smm0 et smm1 ) ;

meta1 = daller surli1 dz0 dr1 scon10 ;
meta2 = daller (inve dr1) dz1 scon21 scon11 ;
meta3 = daller (inve dz0) surli0 smm0 dr0 ;
meta4 = daller (inve dz1) (inve dr0) smm1 scon20 ;

meta = meta1 et meta2 et meta3 et meta4 ;


oublier meta1 ;
oublier meta2 ;
oublier meta3 ;
oublier meta4 ;
*trac meta ;
*
** On réoriente les éléments.
*

meta = 'ORIENTER' meta ;
tass meta ;
nbel1 = nbel meta ; list nbel1 ;

* transformation des éléments en qua9

£META  = CHANGER META  MACRO ;
£SURLI = CHANGER SURLI MACRO ;
£SCON  = CHANGER SCON  MACRO ;
£SMM   = CHANGER SMM   MACRO ;

ELIM 1.e-6 (£META ET £SURLI ET £SCON ET £SMM) ;

* formulation du modèle NAVIER_STOKES

$META  = MODE £META  'NAVIER_STOKES' MACRO ;
$SURLI = MODE £SURLI 'NAVIER_STOKES' MACRO ;
$SCON  = MODE £SCON  'NAVIER_STOKES' MACRO ;
$SMM   = MODE £SMM   'NAVIER_STOKES' MACRO ;

mmeta = doma $META  maillage ;
ssurli = doma $SURLI maillage ;



*       _______________________________________________
*
*        SYSTEME D'EQUATIONS ET CONDITIONS AUX LIMITES
*       _______________________________________________



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

RV = EQEX $META
OPTI 'EF'  'IMPL' CENTREP1  'SUPG'
ZONE $META OPER KBBT 1 (-1.)              INCO 'UN' 'PN'
ZONE $META OPER NS 1. 'UN' nu             INCO 'UN'
OPTI 'CENTREE'
ZONE $SURLI OPER TOIMP 'tau'              INCO 'UN'
;


* conditions aux limites
RV = EQEX RV
CLIM UN UIMP  £SMM     0.
     UN VIMP £SURLI    0.
;

*                      ________________
*
*                       INITIALISATION
*                      ________________

RV.INCO = TABLE INCO ;

RV.INCO.'UN'   = KCHT $META VECT SOMMET (1.e-8 1.e-8) ;
RV.INCO.'PN'   = KCHT $META SCAL CENTREP1 0. ;
RV.INCO.'tau'  = KCHT $SURLI VECT CENTRE (0. 0. ) ;


********************
* terme tau
* grad T < 0 ,lev > 0 => Marangoni < 0  donc opposé à la qdm
***********************

* Courbe de temperature selon CHAN, CHEN, MAZUMDER
'SI' ( MOAXI ) ;
bb = 0.70185 ;
cc = 0.35093 ;
f1 = 0.89872 ;
SINON ;
bb = 1. ;
cc = 0.33333 ;
f1 = 1. ;
'FINSI' ;

I0 = PCANON /(PI * (RFOCAL**2.)) ;
q0 = I0 ;
q1 = q0 '/' rbe '/' rbe ;
DTmax = bb * (cc ** -.25) * q0 * (q1 ** -.25) 
        * ((2 * (rho0 ** 2) * (Cp ** 3) * lev * mu) ** -.25) ;
'LISTE' DTmax ;
XK = ssurli coor 1 ;
Ttheo = 'KOPS' DTmax * ('KOPS' 1 '-' ('KOPS' (cc / bb / rbe / rbe)
         * ('KOPS' XK ** 2 ) ) ) ;
Ttheo = 'KOPS' Ttheo '+' 1200. ;
TK = KCHT $SURLI SCAL SOMMET Ttheo ; 

** calcul de l'effet Marangoni
GTK = KOPS GRAD TK $SURLI ;
MARA = KOPS lev '*' GTK ;
RV.INCO.'tau' = KCHT $SURLI VECT CENTRE MARA ;
MARAN = ELNO $SURLI MARA ;

*                     __________________
*
*                      BOUCLE ITERATIVE
*                     __________________


REPETER BLOCKI NITER ;

MESS 'ITERATION INTERNE N°' &BLOCKI ;

UN0 = COPI RV.INCO.'UN' ;



*
* CALCUL D'UNE ITERATION
*
RV.METHINV.FCPRECT = 1 ;
RV.METHINV.FCPRECI = 1 ;
EXEC RV ;

*
*                      TESTS DE CONVERGENCE
*
UNR = EXCO UX RV.INCO.'UN' ;
UNZ = EXCO UY RV.INCO.'UN' ;
URMAX = MAXI UNR ;
UZMAX = MAXI UNZ ;

MESS 'URMAX = ' URMAX ;
MESS 'UZMAX = ' UZMAX ;


UR0 = UN0 EXCO 'UX' scal ;
UR0 = KCHT $META SCAL SOMMET UR0 ;
UNR = KCHT $META SCAL SOMMET UNR ;
ERUNR = KOPS ( KOPS UNR '-' UR0 ) '/' ( (MAXI (ABS UNR) ) + 1.e-13) ;
ERRUNR = MAXI ( ABS ERUNR ) ;
MESS 'ERRUNR = ' ERRUNR ;

UZ0 = EXCO UY UN0 ;
UZ0 = KCHT $META SCAL SOMMET UZ0 ;
UNZ = KCHT $META SCAL SOMMET UNZ ;
ERUNZ = KOPS ( KOPS UNZ '-' UZ0 ) '/' ( (MAXI (ABS UNZ) ) + 1.e-13) ;
ERRUNZ = MAXI ( ABS ERUNZ ) ;
MESS 'ERRUNZ = ' ERRUNZ ;


ERR = ERRUNR ; VAR = MOTS 'UNR' ;

SI ( ERRUNZ > ERR ) ;
  ERR = ERRUNZ ; VAR = MOTS 'UNZ' ;
FINSI ;


VAR1 = EXTR 1 VAR ;
MESS 'erreur relative maximale sur ' VAR1 ERR ;

* RELAXATION

RV.INCO.'UN' = KOPS ( KOPS ALFAU '*' RV.INCO.'UN')
                    '+' ( KOPS (1. - ALFAU) '*' UN0 ) ;


SI (ERR <EG PREC) ;
* convergence des itérations internes
  MESS 'CONVERGENCE EN ' &BLOCKI 'ITERATIONS' ;
  QUITTER BLOCKI ;
FINSI ;

SI (&BLOCKI EGA NITER ) ;
  MESS 'PAS DE CONVERGENCE EN ' NITER 'ITERATIONS' ;
FINSI ;

******************************
** FIN DE LA BOUCLE ITERATIVE
*******************************
FIN BLOCKI ;

*************************************
** COURBE THEORIQUE DE Ur SELON CHAN 
*************************************
XK = ssurli coor 1 ;
URtheo = ((cc * q1 * 2 * lev / mu / Cp) ** .5) * f1 ;
URtheo = 'KOPS' URtheo * XK ;
URtheo = KCHT $SURLI SCAL SOMMET URtheo ; 
evURtheo = evol roug chpo URtheo (inve ssurli ) ;

******************
** IMPRESSIONS
******************
'SI' (DIMPR) ;

titre 'temperature imposee le long de surli.' ;
evtk = evol roug chpo TK (inve ssurli ) ;
'SI' GRAPH ; dess evtk ; 'FINSI' ;

titre 'vitesses calculees (ampl= 0.001) (m/s)' ;
UN = RV.INCO.'UN' ;
VIT1 = VECT UN 0.001  UX UY ROUG ;
'SI' GRAPH ; trac VIT1 mmeta ; 'FINSI' ;

titre 'vitesse radiale le long de surli (m/s) (en rouge, la theorie)' ;
evux = evol chpo UN UX (inve ssurli ) ;
'SI' GRAPH ; dess (evux 'ET' evURtheo) 'LEGE' ; 'FINSI' ;

titre 'vitesse radiale theorique le long de surli (m/s)' ;
'SI' GRAPH ; dess evURtheo ; 'FINSI' ;

titre 'contrainte de Marangoni le long de surli (Pa)' ;
evmaran = evol chpo MARAN UX (inve ssurli ) ;
evmaran = -1. * evmaran ;
evmaran = rho0 * evmaran ;
'SI' GRAPH ;dess evmaran ; 'FINSI' ;

* affichage de la composante imposee tau
nomper = 'EXTRAIRE' 3 (rv . 'LISTOPER') ;
notable= 'MOT' ('TEXTE' ('CHAINE' 3 nomper)) ;
mess 'Operateur ' nomper ;
msi mai=('TEXTE' nomper) (rv . notable) ;
msir = redu msi surli ;
ev1 = evol vert chpo msir 1UN (INVE surli) ;
'SI' ( MOAXI ) ;
XK = ssurli coor 1 ;
XK = KCHT $SURLI SCAL SOMMET XK ;
MARAN = 'EXCO' MARAN UX ;
MARAN = KCHT $SURLI SCAL SOMMET MARAN ;
marx = KOPS (KOPS MARAN * (-1. * rhaut / 2 / NMAIL )) * (2 * pi * XK) ;
SINON ;
MARAN = 'EXCO' MARAN UX ;
marx = KOPS MARAN * (-1. * rhaut / 2 / NMAIL ) ;
'FINSI' ;
titre 'effet marangoni par maille (la theorie en rouge)' ;
evmarx = evol roug chpo marx (inve ssurli ) ;
'SI' GRAPH ;dess (evmarx et ev1) 'LEGE' ; 'FINSI' ;

'FINSI' ;

************
titre 'vitesse radiale le long de surli (m/s) (en rouge, la theorie)' ;
UN = RV.INCO.'UN' ;
evux = evol chpo UN UX (inve ssurli ) ;
'SI' GRAPH ;dess (evux 'ET' evURtheo) 'LEGE' ; 'FINSI' ;

**** TEST DE FONCTIONNEMENT ***********
UNR = EXCO UX RV.INCO.'UN' ;
UNR = 'REDU' UNR SURLI ;
UNR = KCHT $SURLI SCAL SOMMET UNR ;
ERtheo = KOPS ( KOPS UNR - URtheo ) '/' ( (MAXI (ABS UNR) ) + 1.e-13) ;
ERRtheo = MAXI ( ABS ERtheo ) ;
MESS 'erreur de la vitesse radiale par rapport à la théorie =' ERRtheo ;

'SI' (ERRtheo < TERR) ;
   'ERREUR' 0 ;
'SINON' ;
   'ERREUR' 5 ;
'FINSI' ;     

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

'FIN' ;