Télécharger injairA.dgibi

Retour à la liste

*=======================================================================
*
* fichier :  injairA.dgibi
*
* Air injection inside an adiabatic axisymetric cavity
* Comparaison in average with an analytic solution
*
*=======================================================================
* Analytical solution
*
* Mean pressure, temperature and density are computed according to the
* injection at a constant mass flow rate and temperature of a non
* condensable (nc) gas inside a cavity. At the 
* initial time t=0, the cavity is filled with the SAME nc gas.
*
* Hypothesis : Laplace coefficient Gama=cp/cv is constant 
*
'DEBP' SOREF Ltps*'LISTREEL'  QM*'FLOTTANT'
               Vr*'FLOTTANT'  Sr*'FLOTTANT'  Pnm*'FLOTTANT'
            tetai*'FLOTTANT' INCOND*'MOT'  COUL1*'MOT' ;
*
* E/  : Ltps   : List of time steps (in s)
* E/  : QM     : Mas flow rate at the injection (in kg/s)
* E/  : Vr     : Volume of the cavity (in m3)
* E/  : Sr     : Surface of the wall (in m2)
* E/  : Pnm    : Initial pressure (in Pa)
* E/  : TETAI  : Initial temperature (in C)
* E/  : INCOND : Name of the nc gas
* E/  : COUL1  : Color of the evolutions
*
*  /S : evp1   : EVOLUTION : Mean pressure evolution in time
*  /S : evt1   : EVOLUTION : Mean temperature evolution in time
*  /S : evrho1 : EVOLUTION : Mean nc qas density evolution in time
*
  'SI' ('EGA' INCOND 'AIR' ) ;
     a1 a2 a3 a4 a5 a6 Cp = calcp tetai ;
     Rinc = 287.1 ;
  'SINO' ;
     'SI' ('EGA' INCOND 'N2') ;
        a1 a2 a3 Cp a5 a6 a7 = calcp tetai ;
        Rinc = 296.9 ;
     'SINO' ;
        'MESS' ' Unknown non condensable gas' ' ' incond ;
        'ERRE' 5 ;
     'FINS' ;
  'FINS' ;
  'MESS' 'Reference temperature (C)   :' tetai ; 
  'MESS' 'Thermal capacity cp (J/kg/K):' cp ;
  'MESS' 'Specific constant r (J/kg/K):' Rinc;
*
  ntps  = 'DIME' Ltps ;
  'SI' (ntps '>' 1) ;
    tpsnm = 'EXTR' ltps 1 ;
  'SINO' ;
    'MESS' ' Bad time list :' ;
    'LIST' ltps ;
    'ERRE' 5 ;
  'FINS' ;
*
* Initial conditions
  Rhonm = Pnm '/' Rinc '/' (tetai + 273.15) ;
  Lrho  = 'PROG' Rhonm ;
  Lp    = 'PROG' Pnm   ;
  LT    = 'PROG' tetai ;
*
  gama  = Cp '/' (Cp '-' Rinc) ;
  'REPE' BLOC (ntps '-' 1) ;
    tpsn = 'EXTR' ltps (&BLOC '+' 1) ;
    DT0  = tpsn '-' tpsnm ;
    dpdt = gama '*' Rinc '*' (tetai + 273.15) '*' Qm '/' Vr ;
    Pn   = DT0 '*' dpdt '+' Pnm ;
    Rhon = DT0 '*' Qm '/' Vr '+' Rhonm ;
    Tn   = Pn '/' Rhon '/' Rinc '-' 273.15 ;
*  
*    'SI' ('EGA' INCOND 'AIR' ) ;
*       a1 a2 a3 a4 a5 a6 Cp = calcp tn ;
*    'SINO' ;
*       a1 a2 a3 Cp a5 a6 a7 = calcp tn ;
*    'FINS' ;
*
    Lrho = Lrho 'ET' ('PROG' Rhon) ;
    Lp   = Lp   'ET' ('PROG' Pn)   ;
    LT   = LT   'ET' ('PROG' Tn)   ;
*
    tpsnm = tpsn ;
    Rhonm = Rhon ;
    Pnm   = Pn ;
  'FIN' BLOC ;
  evp1   = 'EVOL' coul1 'MANU' ltps lp    ;
  evt1   = 'EVOL' coul1 'MANU' ltps lt    ;
  evrho1 = 'EVOL' coul1 'MANU' ltps lrho  ;
'FINP' evp1 evt1 evrho1 ;
*=======================================================================
*
'OPTI' 'DIME' 2 'ELEM' 'QUA4' 'MODE' 'AXIS' 'TRAC' 'PSC' 'EPTR' 10 ;
'DENS' 1. ;
*
COMPLET = faux ;
GRAPH   = faux ;
*
'SI' COMPLET ;
  Tmax = 500. ;
  DT0  = 1.   ;
  n1 = 5 ;
  n2 =10 ;
  n3 =20 ;
'SINO' ;
  Tmax = 500. ;
  DT0 = 2. ;
  n1 = 3 ;
  n2 = 3 ;
  n3 = 3 ;
'FINS' ;
nbit  = 'ENTI' (Tmax / DT0) ;
Qinj  = 0.2     ;
PT0   = 1.e5    ;
TF0   = 300. - 273.15 ;
hext  = 0.      ;
Text  = 0.      ;
*
*
* Mesh
*
*
*  R1 = Radius of the vertical cylinder
*  H1 = Half-height of the vertical cylinder
*  H2 = Half-height of the cavity
*  DB = Radius of the injection
* -HB = Elevation of the injection 
*
*  episo = Width of the internal thermal insulator
*
epsi0 = 1.e-5 ;
*
EPISO = 0.044 ;
R1    = 1.200 - episo ;
H1    = 1.605 ;
H2    = 2.342 - episo ;
DB    = 0.1   ;
HB    = 0.835 ;
*
* Points for rotation purposes
CA1  = 0. 0.275 ;
CA2  = (0.33*R1) H1 ;
*
* Points for symetry purposes 
P0   = 0. 0. ;
P2   = DB 0. ;
PA0  = R1 0.;
*
* Specific points
PAXb = 0.      (-1. '*' HB) ;
P2b  = DB      (-1. '*' HB) ;
PAb  = R1      (-1. '*' HB) ;
PAbm = (R1/2.) (-1. '*' HB) ;
PAX2 = 0. H2 ;
PA20 = PAX2 'TOUR'  -5. CA1 ;
PA2  = PA20 'TOUR' -15. CA1 ;
PA1  = R1 H1 ;
PAXM = (PAX2 'PLUS' PAXb) '*' 0.5 ;
*
PbX2 = PAX2 'SYME' 'DROITE' P0 PA0;
Pb1  = R1 (-1. '*' H1) ;
*
* Specific treatment (lips) at the injection
breche = PAXb 'DROI' n1 P2b ;
entb = breche 'PLUS' (0. 0.);
'DEPL' entb 'TOUR' 15. P2b;
P0b  = 'POIN' 1 entb ;
*
* Upper volume
* From the injection elevation to the top
bas    = P2b  'DROI' 'DINI' (DB/n1)  'DFIN' (R1/10.) PABM
              'DROI' 'DINI' (R1/10.) 'DFIN' (db/n1)  PAB ;
n2     = 'NBEL' bas ;
*
plaf0 = PAX2 'CERC' n1 CA1 PA20     ;
plaf1 = PA20 'CERC' (n2/2) CA1 PA2  ;
plaf2 = PA2  'CERC' (n2/2) CA2 PA1  ;
plaf  = 'INVE' (plaf0 'ET' plaf1 'ET' plaf2) ;
*
axeh = PAX2 'DROI' 'DINI' (h2/n3/4.) 'DFIN' (h1/n3)    PAXM
            'DROI' 'DINI' (h1/n3)    'DFIN' (h2/n3/4.) PAXB ;
n3n  = 'NBEL' axeh ;
*
paroih = PAb 'DROI' n3n PA1 ;
*
mth = 'DALL' (breche 'ET' bas) paroih plaf axeh;
*
* Lower volume
* From the bottom to the injection elevation
*
axeb   = PBX2 'DROI' n3 P0b ;
paroib = PAb  'DROI' n3 Pb1 ;
fond   = plaf 'SYME' 'DROITE' P0 PA0 ;
'ELIM' (entb 'ET' bas 'ET' paroib 'ET' fond 'ET' axeb) epsi0 ;
mtb    = 'DALL' (entb 'ET' bas) paroib fond axeb ;
*
* Mesh of the cavity and specific boundaries (axis, containment wall) 
axe    = axeb 'ET' ('INVE' axeh) ;
parext = ('INVE' (plaf 'ET' paroih)) 'ET' paroib 'ET' fond ;
mt     = mth  'ET' mtb  ;
'ELIM' mt epsi0 ;
*
* Control
'SI' GRAPH ;
  'TRAC' (axeh 'ET' breche 'ET' bas 'ET' paroih 'ET' plaf 'ET' ca1) ;
  'TRAC' mt ;
  pla = paroih 'ET' plaf ;
  pla = pla 'ET' ('SYME' pla 'DROIT' P0 pax2) ;
  'TRAC' (pla 'ET' ca1 'ET' ca2 'ET' PA1) ;
  'TRAC' ('ELIM' (mt 'ET' ('SYME' mt 'DROIT' P0 pax2)) epsi0) ;
'FINS' ;
*
*
* Data for execrxt.procedur
*
*
rxt = 'TABLE' ;             
rxt . 'VERSION' = 'V0'                                           ;
rxt . 'vtf'     = mt ;
rxt . 'axe'     = axe ;
rxt . 'epsi'    = epsi0 ;
rxt . 'pi'      = 0.5 0.5 ;
 
rxt . 'DISCR'   = 'LINE';
rxt . 'KPRE'    = 'MSOMMET' ;
rxt . 'DT0'     = DT0 ;
 
rxt . 'MODTURB' = 'LMEL' ;
rxt . 'LMEL'    = 0.01 ;
 
*
rxt . 'TF0'    = TF0 ;
rxt . 'PT0'    = PT0 ;
*
rxt . 'Breches'                    = 'TABLE' ;
rxt . 'Breches' . 'A'              = 'TABLE' ;
rxt . 'Breches' . 'A' . 'Maillage' =  breche ;
rxt . 'Breches' . 'A' . 'diru'     =  (0. 1.) ;
rxt . 'Breches' . 'A' . 'scenario' = 'TABLE' ;
rxt . 'Breches' . 'A' . 'scenario' . 't'    = 'PROG' 0.0  1000. ;
rxt . 'Breches' . 'A' . 'scenario' . 'qair' = 'PROG' Qinj Qinj  ;
rxt . 'Breches' . 'A' . 'scenario' . 'tinj' = 'PROG' TF0  TF0   ;
*
rxt . 'DETMAT'  = vrai ;
rxt . 'RENU'    = 'RIEN' ;
rxt . 'GRAPH'   = faux ;
execrxt  0 rxt ;
rxt . 'GRAPH' = GRAPH ;
*
'MESS' ' Volume :' ('SOMT' rxt . 'GEO' . 'Diag') ;
*
* Transient computation
* (the numerical diffusion is set to zero for the equation devoted to fluid temperature)
'SI' COMPLET ;
   rxt . 'TBT' . 'RTF' . '1TSCA' . 'KOPT' . 'CMD' = 0. ;
'FINS' ;
execrxt  nbit rxt ;
*
*
* Plots devoted to inj... cases
*
*
Ltps = rxt . 'TIC' . 'LTPS' ;
Vr   = 'MAXI' ('RESU' ('DOMA' rxt . 'GEO' . '$vtf'   'VOLUME')) ;
Sr   = 'MAXI' ('RESU' ('DOMA' rxt . 'GEO' . '$menvf' 'VOLUME')) ;
'MESS' vr ; 'MESS' sr ;
*Vr = 17.4 ; Sr = 35.914 ;
evp1 evt1 evrho1 = SOREF Ltps Qinj Vr Sr PT0 TF0 'AIR' 'ROUG' ;
evm1 = evrho1 * Vr ;
 
evp2   = 'EVOL' 'MANU' Ltps rxt . 'TIC' . 'PT'    ;
evt2   = 'EVOL' 'MANU' Ltps rxt . 'TIC' . 'Tfm'   ;
evrho2 = 'EVOL' 'MANU' Ltps rxt . 'TIC' . 'Rhomn' ;
evm2   = evrho2 * Vr ;
'SI' GRAPH ;
  TAB1 = 'TABLE' ;
  TAB1 . 'TITRE' = 'TABLE' ;
  TAB1 . 1 = 'MOT' ' MARQ CROI REGU ';
  TAB1 . 2 = 'MOT' ' MARQ CARR REGU TIRC';
  TAB1 . 3 = 'MOT' ' MARQ LOSA REGU TIRC';
  TAB1 . 4 = 'MOT' ' MARQ TRIU REGU TIRC';
  TAB1 . 5 = 'MOT' ' MARQ TRID REGU TIRC';
  TAB1 . 'TITRE' . 1 = 'Sol Ref';
  TAB1 . 'TITRE' . 2 = 'Nautilus CFD';
*
  'DESS' (evp1 'ET' evp2)
  'TITR' 'Pressure' 'MIMA'
  'GRIL' 'POIN' 'GRIS' 'TITX' 's' 'TITY' ' Pa' 'LEGE' tab1 ;
  'DESS' (evt1 'ET' evt2)
  'TITR' 'Mean gas temperature' 
  'GRIL' 'POIN' 'GRIS' 'TITX' 's' 'TITY' ' C' 'LEGE' tab1 ;
  'DESS' (evrho1 'ET' evrho2)
  'TITR' 'Mean gas density'                                     
  'GRIL' 'POIN' 'GRIS' 'TITX' 's' 'TITY' 'kg/m3' 'LEGE' tab1 ;
  'DESS' (evm1 'ET' evm2)
  'TITR' 'Air mass'                                     
  'GRIL' 'POIN' 'GRIS' 'TITX' 's' 'TITY' 'kg' 'LEGE' tab1 ;
'FINS' ;
*
*
* Tests
*
*
'SI' ('NON' COMPLET) ;
  ERR1 = 0 ;
  devp = evp1 - evp2 'ABS' ;
  erp  = 'MAXI' ('EXTR' devp 'ORDO') ;
  'MESS' ' ERP =' ERP;
  devt = evt1 - evt2 'ABS' ;
  ert  = 'MAXI' ('EXTR' devt 'ORDO') ;
  'MESS' ' ERT =' ERT ;
  t_calc     = 'IPOL' ('PROG' 50. 250. 500.) evt2 ; 'LIST' t_calc ;
  t_calc_ref = 'PROG' 63.032 103.80 116.46 ;
  ert_calc   = 'MAXI' (t_calc '-' t_calc_ref 'ABS') ;
  'MESS' ' ERT_calc =' ert_calc ; 
  'SI' (ERP '>' 0.2D5) ; err1 = err1 '+' 1 ; 'FINS' ;
  'SI' (ERT '>' 10.)   ; err1 = err1 '+' 1 ; 'FINS' ;
  'SI' (ert_calc '>' 0.01) ; err1 = err1 '+' 1 ; 'FINS' ;
  'SI' ('NEG' ERR1 0) ;
    'ERRE' 5 ;
  'FINS' ;
'FINS' ;
*
'FIN' ;
 
 
 

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