Télécharger Marangoni2.dgibi

Retour à la liste

* fichier :  Marangoni2.dgibi
*                      _____________________
*
*                       CAS TEST EFFET MARANGONI
*             UTILISATION DE LA LOI DE CHAN, CHEN, MAZUMDER
*           Journal of Heat Transfert, vol 110, 1988
*                      _____________________
*                      OPTIONS TESTEES PLAN QUAF
 
***************************************************************
** 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.35;
'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 'DIRE' (1. 0.) 'DROI' p1 p2) ) coul vert ;
scon11 = (inve (smm1 proj 'DIRE' (1. 0.) 'DROI' p1 p2) ) coul vert ;
scon1 = scon11 et scon10 ;
 
dz0 = (smm0 proj 'DIRE' (1. 0.) 'DROI' p01 p32) coul bleu ;
dz1 = (smm1 proj 'DIRE' (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 'DIRE' (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 'DIRE' (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 ; 
 
 
**************procedure gradient******
DEBPROC GRADS FK*'CHPOINT ' LIGNK*'MAILLAGE' ;
* procédure de calcul de df/ds
* gradient d'une fonction FK suivant l'abscisse curviligne de
* la ligne LIGNK
NELIGNK = NBEL LIGNK ;
REPETER BLIGNK NELIGNK ;
ELEMI = ELEM LSURLI &BLIGNK ;
PO1I = POIN ELEMI INITIAL ;
PO2I = POIN ELEMI FINAL   ;
XK1I = COOR 1 PO1I ;
XK2I = COOR 1 PO2I ;
ZK1I = COOR 2 PO1I ;
ZK2I = COOR 2 PO2I ;
FK1I = EXTR FK 'SCAL' PO1I ;
FK2I = EXTR FK 'SCAL' PO2I ;
DXKI = XK2I - XK1I ;
DZKI = ZK2I - ZK1I ;
DSKI = ( ( DXKI * DXKI ) + ( DZKI * DZKI ) ) ** 0.5 ;
DFKI = FK2I - FK1I ;
SI ( DXKI &lt;EG 0. ) ;
  SIGNI = -1. ;
SINON ;
  SIGNI = 1. ;
FINSI ;
DFDI = 0.5 * SIGNI * DFKI / DSKI ;
CDFDI = 'MANU' 'CHPO' (PO1I ET PO2I ) 1 SCAL ( PROG 2 * DFDI )
'NATURE' DISC ;
SI ( &BLIGNK EGA 1 ) ;
  DFDS = CDFDI ;
SINON ;
  DFDS = DFDS ET CDFDI ;
FINSI ;
FIN BLIGNK ;
FINPROC DFDS ;
*************** fin de la procedure *******
 
 
** calcul de l'effet Marangoni
LSURLI  = CHANGER SURLI SEG2 ;
GTK = GRADS TK LSURLI ;
MARA = KOPS lev '*' GTK ;
MARA = 'NOMC' UX MARA ;
RV.INCO.'tau' = KCHT $SURLI VECT CENTRE MARA ;
MARA = 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 &lt;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' ;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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