* fichier : Marangoni3.dgibi
* _____________________
*
* CAS TEST D'UN ECOULEMENT SOUS UN FLUX D'ENERGIE CONCENTRE
* UTILISATION DE LA LOI DE CHAN, CHEN, MAZUMDER
* Journal of Heat Transfert, vol 110, 1988
* _____________________
***************************************************************
** DESCRIPTION DU TEST :
**
** Une poche de Cerium 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 donc un flux surfacique sur une poche
** de metal liquide. Nous comparons les profils obtenus de température
** et de vitesse a la surface du liquide avec les valeurs théoriques
** données par la loi de CHAN.
***************************************************************
option DIME 2 ELEM qua8 MODE axis; MOAXI = VRAI ;
*option DIME 2 ELEM qua8 MODE plan; MOAXI = FAUX ;
option trac psc ;
* _________________________________
*
* différents parametres du calcul
* _________________________________
COMPLET=FAUX ;
GRAPH =FAUX ;
* Nombre de mailles
NMAIL = 20 ;
* nombre maximal d'itérations
'SI' COMPLET ;NITER = 200 ;
'SINON' ;NITER = 4 ;
'FINSI' ;
* coefficient de relaxation pour le champ de température
ALFAT = 0.2 ;
* coefficient de relaxation pour le champ de vitesse
ALFAU = 0.5 ;
* précision requise
PREC = 5.e-3 ;
* pseudo pas de temps pour la qdm
dt = .01 ;
* _______________________________________________________
*
* géométrie du creuset : tronc de cone (dimension en m)
* _______________________________________________________
*
* rayon du bas du creuset
rbas = 0.0399 ;
* rayon du haut du creuset
rhaut = 0.04 ;
* hauteur de liquide
hliq = 0.020 ;
* épaisseur de la couche limite
eclm = 1.e-4 ;
* on ne maille pas l'axe, à cause des termes en 1/r et 1/r2
ref = 0. ;
* rayon moyen du faisceau électonique
rbe = 0.010 ;
* __________________________________________________________________
*
* Caractéristiques du faisceau d'électrons à répartition d'énergie
* gaussienne : 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 ;
* accélération de la pesanteur
g = -9.81 ;
*
* __________________________________________________
*
* propriétés physiques du cérium évaluées à 2500 K
* __________________________________________________
* viscosité dynamique (kg/m/s)
mu = 6.5e-3 ;
* masse volumique (kg/m3)
rho0 = 6310. ;
* conductivité thermique (W/m/K)
lambda = 25. ;
* chaleur spécifique (J/kg/K)
Cp = 270 ;
* température de fusion (K)
tfus = 1070. ;
* cooefficient de dilatation volumique (1/K)
beta = 4.2e-5 ;
* tension de surface (N/m)
sig0 = 0.34 ;
* 1/sig0 dsig0/dT (1/K)
gam0 = 5.88e-4 ;
* émissivité
emis = 0.34 ;
* coefficient d'absorption des électrons (% d'énergie absorbée)
CABSOR = 1. ;
* masse molaire (kg/mole)
MM = 0.14012 ;
* chaleur latente de vaporisation (J/kg)
LAT = 2.95e6 ;
* coefficients pour le calcul de la pression de vapeur saturante
PA = 6.023 ;
PB = -21278. ;
PC = -0.1127 ;
* coefficient de recondensation de la vapeur
RET = 0. ;
* viscosité cinématique (m2/s)
nu = mu / rho0 ;
* diffusivité thermique (m2/s)
khi = lambda / ( rho0 * Cp ) ;
* effet MARANGONI
lev = 1. * gam0 * sig0 / rho0 ;
* ______________________
*
* CREATION DU MAILLAGE
* ______________________
p0 = ref 0. ;
p1 = rhaut 0. ;
p2 = rbas (-1. * hliq) ;
p3 = ref (-1. * hliq) ;
p01 = rbe 0. ;
p32 = rbe (-1. * hliq) ;
p03 = ref (-10. * eclm) ;
p12 = rhaut (-10. * eclm) ;
d0 = rbe / NMAIL ;
d1 = rbe * 4 / NMAIL ;
d2 = rhaut / 2 / NMAIL ;
d3 = rhaut * 6 / NMAIL ;
surli0 = (p01 droi dini d1 dfin d0 p0 ) coul roug ;
surli1 = (p1 droi dini d3 dfin d2 p01) coul roug ;
surli = surli1 et surli0 ;
d0 = eclm * 20 / NMAIL ;
d1 = eclm * 60 / NMAIL ;
d2 = hliq / NMAIL ;
d3 = hliq * 8 / NMAIL ;
smm0 = (p0 droi dini d0 dfin d1 p03) coul jaun ;
smm1 = (p03 droi dini d2 dfin d3 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 ;
'SI' GRAPH ;trac meta ; 'FINSI' ;
*
** 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 ;
£SURLI0 = CHANGER SURLI0 MACRO ;
£SCON = CHANGER SCON MACRO ;
£SCON2 = CHANGER SCON2 MACRO ;
£SMM = CHANGER SMM MACRO ;
ELIM 1.e-6 (£META ET £SURLI ET SURLI0 ET £SCON ET SCON2 ET £SMM) ;
* formulation du modèle NAVIER_STOKES
$META = MODE £META 'NAVIER_STOKES' MACRO ;
$SURLI = MODE £SURLI 'NAVIER_STOKES' MACRO ;
$SURLI0 = MODE £SURLI0 'NAVIER_STOKES' MACRO ;
$SCON = MODE £SCON 'NAVIER_STOKES' MACRO ;
$SCON2 = MODE £SCON2 'NAVIER_STOKES' MACRO ;
$SMM = MODE £SMM 'NAVIER_STOKES' MACRO ;
mmeta = doma $META maillage ;
ssurli = doma $SURLI maillage ;
ssurli0 = doma $SURLI0 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'
ZONE $META OPER DFDT 1. 'UN' dt 'UN' nu INCO 'UN'
OPTI 'CENTREE'
ZONE $SURLI OPER TOIMP 'tau' INCO 'UN'
;
* equation scalaire sur la température + flux imposé sur surli
RV = EQEX RV
OPTI 'EF' 'IMPL' 'SUPG'
ZONE $META OPER TSCA 1. 'UN' khi 0. INCO 'TN'
ZONE $SURLI OPER FIMP 'F1E' INCO 'TN'
;
* conditions aux limites
RV = EQEX RV
CLIM UN UIMP £SMM 0.
UN UIMP £SCON2 0.
UN VIMP £SURLI 0.
TN TIMP £SCON2 tfus
;
* ________________
*
* INITIALISATION
* ________________
RV.INCO = TABLE INCO ;
RV.INCO.'UN' = KCHT $META VECT SOMMET (1.e-8 1.e-8) ;
RV.INCO.'TN' = KCHT $META SCAL SOMMET tfus ;
RV.INCO.'PN' = KCHT $META SCAL CENTREP1 0. ;
RV.INCO.'tau' = KCHT $SURLI VECT CENTRE (0. 0. ) ;
* terme f1e flux imposé sur la surface libre
I0 = PCANON /(PI * (RFOCAL**2.)) ;
XK = ssurli coor 1 ;
XK = KCHT $SURLI SCAL SOMMET XK ;
FBA = 'KOPS' I0 * ('KOPS' 1. - ( 'KOPS' ('KOPS' XK '/' rbe) ** 2.) );
mask1 = XK MASQUE 'EGINF' rbe ;
mask1 = 'KCHT' $SURLI SCAL SOMMET mask1 ;
FBA = 'KOPS' FBA * mask1 ;
FBA = KCHT $SURLI SCAL SOMMET FBA ;
F1 = 'KOPS' FBA '/' ( rho0 * Cp ) ;
F1E = NOEL $SURLI F1 ;
RV.INCO.'F1E' = KCHT $SURLI SCAL CENTRE F1E ;
* __________________
*
* BOUCLE ITERATIVE
* __________________
REPETER BLOCKI NITER ;
*
** calcul des termes de couplages
*
MESS 'ITERATION INTERNE N°' &BLOCKI ;
UN0 = COPI RV.INCO.'UN' ;
TN0 = COPI RV.INCO.'TN' ;
* terme tau
* grad T < 0 ,lev > 0 => Marangoni < 0 donc opposé à la qdm
*
TK = REDU TN0 ssurli ;
TK = KCHT $SURLI SCAL SOMMET TK ;
GTK = KOPS GRAD TK $SURLI ;
mask1 = GTK MASQUE 'EGINF' 0. ;
mask1 = 'KCHT' $SURLI VECT CENTRE mask1 ;
GTK = 'KOPS' GTK * mask1 ;
MARK = KOPS lev '*' GTK ;
RV.INCO.'tau' = KCHT $SURLI VECT CENTRE MARK ;
MARKN = ELNO $SURLI MARK ;
*
* 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 ;
TMAX = MAXI RV.INCO.'TN' ;
MESS 'TMAX = ' TMAX ;
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 ;
TN = RV.INCO.'TN' ;
ERTN = KOPS ( KOPS TN '-' TN0 ) '/' ( (MAXI (ABS TN) ) + 1.e-13) ;
ERRTN = MAXI ( ABS ERTN ) ;
MESS 'ERRTN = ' ERRTN ;
ERR = ERRUNR ; VAR = MOTS 'UNR' ;
SI ( ERRUNZ > ERR ) ;
ERR = ERRUNZ ; VAR = MOTS 'UNZ' ;
FINSI ;
SI ( ERRTN > ERR ) ;
ERR = ERRTN ; VAR = MOTS 'TN' ;
FINSI ;
VAR1 = EXTR 1 VAR ;
MESS 'erreur relative maximale sur ' VAR1 ERR ;
* RELAXATION
RV.INCO.'TN' = KOPS ( KOPS ALFAT '*' RV.INCO.'TN')
'+' ( KOPS (1. - ALFAT) '*' TN0 ) ;
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
******************************
'SI' ( MOAXI ) ;
bb = 0.70185 ;
cc = 0.35093 ;
f1 = 0.89872 ;
SINON ;
bb = 1. ;
cc = 0.33333 ;
f1 = 1. ;
'FINSI' ;
q0 = I0 ;
q1 = q0 '/' rbe '/' rbe ;
XK = ssurli coor 1 ;
** Courbe theorique de temperature sur la surface selon CHAN
DTmax = bb * (cc ** -.25) * q0 * (q1 ** -.25)
* ((2 * (rho0 ** 2) * (Cp ** 3) * lev * mu) ** -.25) ;
'LISTE' DTmax ;
Ttheo = 'KOPS' DTmax * ('KOPS' 1 '-' ('KOPS' (cc / bb / rbe / rbe)
* ('KOPS' XK ** 2 ) ) ) ;
Ttheo = 'KOPS' Ttheo '+' Tfus ;
Ttheo = KCHT $SURLI SCAL SOMMET Ttheo ;
** Courbe theorique de la vitesse radiale sur la surface selon CHAN
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
******************
titre 'Temperature le long de surli0 (W/m2)' ;
TN = RV.INCO.'TN' ;
evtemp = evol noir chpo TN (inve ssurli0 ) ;
evTtheo = evol roug chpo Ttheo (inve ssurli0 ) ;
'SI' GRAPH ;dess (evtemp 'ET' evTtheo) 'LEGE' ; 'FINSI' ;
titre 'vitesses (ampl= 0.01) (m/s)' ;
UN = RV.INCO.'UN' ;
VIT1 = VECT UN 0.01 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 'pression (Pa)' ;
PNN = ELNO $META (RV.'INCO'.'PN') CENTREP1 ;
'SI' GRAPH ;trace PNN mmeta ; 'FINSI' ;
***************************************
**** TEST DE FONCTIONNEMENT ***********
***************************************
TN = RV.INCO.'TN' ;
TNmax = 'EXTR' TN 'SCAL' (p0) ;
TNtmax = 'EXTR' Ttheo 'SCAL' (p0) ;
'LISTE' TNmax ;
'LISTE' TNtmax ;
ERtheo = ( TNmax - TNtmax ) / ( TNmax - Tfus) ;
ERRtheo = ABS ERtheo ;
MESS 'erreur de la temperature max par rapport à la théorie =' ERRtheo ;
'SI' COMPLET ;
'SI' (ERRtheo < .05) ;
'ERREUR' 0 ;
'SINON' ;
'ERREUR' 5 ;
'FINSI' ;
'SINON' ;
'SI' (ERRtheo < .4 ) ;
'ERREUR' 0 ;
'SINON' ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
***************************************
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer