Accueil | Cast3M

* fichier :  aerosol1.dgibi
************************************************************************
************************************************************************
 
GRAPH = FAUX ;
 
**************************************************************************
*                        AEROSOL1.DGIBI                                  *
**************************************************************************
*                                                                        *
*    NATURE DU PROBLEME  :  TRANSPORT DE PARTICULES AVEC DEPOT           *
*                           CONVECTION FORCEE                            *
*                           CALCUL DE L'ECOULEMENT DANS UN PREMIER TEMPS *
*                           CALCUL DU TRANSPORT DANS UN DEUXIEME TEMPS   *
*                           TURBULENT                                    *
*                                                                        *
*    ALGORITHME          :  SEMI-IMPLICITE                               *
*                                                                        *
*    OPERATEURS          :  NSKE , FPU                                   *
*                           TSCAL , ECHIMP , KFPA                        *
*                                                                        *
*    GEOMETRIE           :  LOCAL  2D                                    *
*                                                                        *
*    TEST SUR LA VITESSE DE DEPOT EN 2 POINTS DE LA COUCHE LIMITE        *
*                                                                        *
*P.CORNET-SEMT/TTMF-DECEMBRE 1997                                       *
**************************************************************************
 
*---------------------- CREATION DU MAILLAGE -----------------------------------
 
OPTI ISOV LIGN;
TITRE 'MAILLAGE AEROSOL POUR TURBULENCE'  ;
OPTION DIME 2 ELEM QUA4    ;
 
* EPAISSEUR DE LA COUCHE LIMITE  :
 
E = 0.05 ;
 
* POINTS  :
 
PA1 = 0. 0. ;
 
PB1 = E  E       ;
PB2 = E  6.5     ;
PB3 = E  7.      ;
PB4 = E  9.      ;
PB5 = E  9.5     ;
PB6 = E  25.     ;
PB7 = E  41.     ;
PB8 = E  (50.-E) ;
 
PC1 = (20.-E)  E       ;
PC2 = (20.-E)  6.5     ;
PC3 = (20.-E)  7.      ;
PC4 = (20.-E)  9.      ;
PC5 = (20.-E)  9.5     ;
PC6 = (20.-E)  25.     ;
PC7 = (20.-E)  41.     ;
PC8 = (20.-E)  (50.-E) ;
 
PE8 = (8.+E)  (50.-E)  ;
PE9 = (8.+E)  53.      ;
 
PH8 = (12.-E)  (50.-E) ;
PH9 = (12.-E)  53.     ;
 
* SEGMENTS ELEMENTAIRES :
 
DB12 = PB1 DROIT 13 PB2 ;
DB23 = PB2 DROIT  1 PB3 ;
DB34 = PB3 DROIT  4 PB4 ;
DB45 = PB4 DROIT  1 PB5 ;
DB56 = PB5 DROIT DINI 0.5 DFIN 2.  PB6 ;
DB67 = PB6 DROIT DINI 2.  DFIN 0.5 PB7 ;
DB78 = PB7 DROIT 18 PB8 ;
 
DC12 = PC1 DROIT 13 PC2 ;
DC23 = PC2 DROIT  1 PC3 ;
DC34 = PC3 DROIT  4 PC4 ;
DC45 = PC4 DROIT  1 PC5 ;
DC56 = PC5 DROIT DINI 0.5 DFIN 2.  PC6 ;
DC67 = PC6 DROIT DINI 2.  DFIN 0.5 PC7 ;
DC78 = PC7 DROIT 18 PC8 ;
 
DB1 = DB12 et DB23 et DB34 et DB45 et DB56 ;
DB2 = DB67 et DB78 ;
DB1 = inve DB1 ;
DB2 = inve DB2 ;
 
DD1 = DC12 et DC23 et DC34 et DC45 et DC56 ;
DD2 = DC67 et DC78 ;
 
DE89 = PE8 DROIT  5 PE9 ;
DE98 = INVE DE89 ;
DH89 = PH8 DROIT  5 PH9 ;
DHE9 = PH9 DROIT  8 PE9 ;
 
DBC1 = PB1 DROIT 40 PC1 ;
DBC6 = PB6 DROIT 40 PC6 ;
DCB6 = INVE DBC6 ;
 
DEH8 = PE8 DROIT 8 PH8 ;
DHE8 = inve DEH8 ;
DAE8 = PB8 DROIT 16 PE8 ;
DEA8 = INVE DAE8 ;
DDH8 = PC8 DROIT 16 PH8 ;
 
DDA8 = DDH8 et DHE8 et DEA8 ;
 
DA58 = DB56 et DB67 et DB78 ;
DA85 = inve DA58 ;
DA21 = inve DB12 ;
CLEXT1 = DA21 et DBC1 et DD1 et DD2 et DDH8 et DH89 ;
CLEXT2 = DE98 et DEA8 et DA85 ;
 
* DOMAINE TOTAL :
 
DOMBAS = DBC1 DD1 DCB6 DB1 DALLER PLAN ;
DOMHAU = DBC6 DD2 DDA8 DB2 DALLER PLAN ;
DOMSORT = DEH8 DH89 DHE9 DE98 DALLER PLAN ;
DOMTOT1 = DOMBAS ET DOMHAU ET DOMSORT ;
 
* ON NOMME LES OBJETS MAILLAGE INTERESSANTS :
 
DENTREE = DB34 ;
dentree = inve dentree ;
DSORTIE = DHE9 ;
DSORTI2 = DHE8 ;
DCOUPE  = DBC6 ;
DENTREV = DB23 ET DB34 ET DB45 ;
 
TROUS   = DENTREE ET DSORTIE   ;
BORTOT  = CONTOUR DOMTOT1      ;
BORDOM  = BORTOT DIFF TROUS    ;
CLPRES  = CLEXT1 ET CLEXT2 ;
 
DOMTOT = ORIENTER DOMTOT1 ;
DEPL DOMTOT HOMO 0.1 PA1  ;
 
*-------------------- CREATION DES MODELES NAVIER_STOKES -----------------------
 
MDOMTOT = CHAN DOMTOT QUAF ;
MCLEXT1 = CHAN CLEXT1 QUAF ;
MCLEXT2 = CHAN CLEXT2 QUAF ;
MCLPRES = CHAN CLPRES QUAF ;
MBORDOM = CHAN BORDOM QUAF ;
ELIM (MDOMTOT et MCLEXT1 et MCLEXT2 et MCLPRES et MBORDOM) 1.E-2 ;
 
$DOMTOT = MODE MDOMTOT 'NAVIER_STOKES' LINE ;
$CLEXT1 = MODE MCLEXT1 'NAVIER_STOKES' LINE ;
$CLEXT2 = MODE MCLEXT2 'NAVIER_STOKES' LINE ;
$CLPRES = MODE MCLPRES 'NAVIER_STOKES' LINE ;
$BORDOM = MODE MBORDOM 'NAVIER_STOKES' LINE ;
DOMA $DOMTOT 'IMPR' ;
 
*-------------- PARAMETRES DU PROBLEME : AIR ET PARTICULES ---------------------
 
NU  = 1.5E-5            ;
U0  = 0.75              ;
K0  = (U0/100.)**2.     ;
E0  = K0**(1.5)/0.2     ;
 
C0  = 1.                                             ;
RAP = 0.5E-6                                         ;
ROP = 1000.                                          ;
ROG = 0. -9810.                                      ;
VS0  = 12277.*RAP*RAP*9.81*ROP                        ;
DIF0 = 1./(5.62E21*NU*RAP)                             ;
VAR = EXP(1.667E7*RAP)                                  ;
VS0   = VS0*(1. + (8.296E-8/RAP) + (2.64E-8/(RAP*VAR)))  ;
DIF0 = DIF0*(1. + (8.296E-8/RAP) + (2.64E-8/(RAP*VAR)))  ;
 
YP1  = 0.15       ;
YP2  = 0.30       ;
 
UET0 = U0/100 ;
UET1 = kcht $CLEXT1 'SCAL' 'CENTRE' UET0 ;
UET2 = kcht $CLEXT2 'SCAL' 'CENTRE' UET0 ;
NUT  = kcht $DOMTOT 'SCAL' 'CENTRE' 1.e-10 ;
 
*-------------- DEFINITION DES TABLES POUR LE CALCUL DU REGIME -----------------
*-------------------- PERMANENT ETABLI POUR L'AIR ------------------------------
 
rv = eqex $DOMTOT 'ITMA' 100   'ALFA' 0.8  'TFINAL' 0.1
 'ZONE' $DOMTOT 'OPER' 'NSKE' NU NUT     'INCO' 'UN' 'KN' 'EN'
 'ZONE' $CLEXT1 'OPER' 'FPU' NU UET1 YP1 'INCO' 'UN' 'KN' 'EN'
 'ZONE' $CLEXT2 'OPER' 'FPU' NU UET2 YP2 'INCO' 'UN' 'KN' 'EN'
;
 
rv = eqex rv OPTI EFM1 'CENTREE'
  'ZONE' $DOMTOT  'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
  'ZONE' $DOMTOT  'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
  'ZONE' $DOMTOT  'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN'
;
 
rv = eqex rv
'CLIM' 'KN' 'TIMP' DENTREE K0 'EN' 'TIMP' DENTREE E0
       'UN' 'UIMP' DENTREE U0 'UN' 'VIMP' DENTREE 0. ;
 
rvp= eqpr $DOMTOT
 'ZONE' $DOMTOT   'OPER' 'PRESSION'  0.
 'ZONE' $CLPRES   'OPER' 'VNIMP' 0.;
 
*----------------------- INITIALISATIONS ---------------------------------
 
rv.inco=table 'INCO' ;
rv.inco.'UN' = kcht $DOMTOT 'VECT' 'SOMMET' (1.E-5 1.E-5 ) ;
rv.inco.'KN' = kcht $DOMTOT 'SCAL' 'SOMMET' 1.E-8 ;
rv.inco.'EN' = kcht $DOMTOT 'SCAL' 'SOMMET' 1.E-10 ;
rv.'FIDT'    = 1 ;
rv.'PASDETPS'.'DELTAT' = 0.01 ;
 
rv.pression=rvp ;
 
*------------------ CALCUL DE L'ECOULEMENT D'AIR -------------------------------
 
 EXEC rv ;
 
*------------------ DONNEES POUR L'EQUATION DE CONCENTRATION -------------------
 
CHPC = DOMA $DOMTOT 'NORMALE' ;
NORM1 = KCHT $CLEXT1 vect centre CHPC  ;
NORM2 = KCHT $CLEXT2 vect centre CHPC  ;
 
VSF = kcht $DOMTOT 'VECT' 'SOMMET' ( 0. VS0 ) ;
VSF = CHAN 'ATTRIBUT' VSF 'NATURE' 'DISCRET' ;
VS0 = 0. - VS0 ;
VSP = kcht $BORDOM 'VECT' 'SOMMET' ( 0. VS0 ) ;
VSP = CHAN 'ATTRIBUT' VSP 'NATURE' 'DISCRET' ;
VS = kcht $DOMTOT  'VECT' 'SOMMET' (VSF et VSP) ;
VN = kops (rv.inco.'UN') '-' VS ;
 
KN = abs  (rv.inco.'KN') ;
EN = abs  (rv.inco.'EN') ;
KN = kops (kops KN '*' KN) '/' EN ;
DIF = kops (kops KN '*' 0.09) '+' DIF0 ;
 
AK1 = KFPA NU YP1 UET1 NORM1 ROG RAP ;
AK2 = KFPA NU YP2 UET2 NORM2 ROG RAP ;
 
*------------ TABLE POUR LE CALCUL DES CONCENTRATIONS --------------------------
 
CONC=TABLE 'ARG';
TETA1=KCHT $CLEXT1 'SCAL' 'CENTRE' 0.;
TETA2=KCHT $CLEXT2 'SCAL' 'CENTRE' 0.;
 
RV1 = EQEX $DOMTOT  'ITMA' 100  'ALFA' 0.7  'TFINAL' 0.1
 'ZONE' $DOMTOT 'OPER' TSCAL 'DIFFU' 'VITESSE'  0.     'INCO' 'CN'
 OPTI EFM1 'CENTREE'
 'ZONE' $CLEXT1 'OPER' ECHIMP AK1 TETA1                'INCO' 'CN'
 'ZONE' $CLEXT2 'OPER' ECHIMP AK2 TETA2                'INCO' 'CN'
 'ZONE' $DOMTOT 'OPER' DFDT   1. 'CN' 'DELTAT'         'INCO' 'CN'
 
*---- conditions aux limites : -----
 'CLIM' 'CN' 'TIMP' DENTREE  C0 ;
 
 
*------------------------ INITIALISATIONS --------------------------------
 
RV1.INCO=TABLE 'INCO' ;
RV1.INCO.'CN' = KCHT $DOMTOT 'SCAL' 'SOMMET' 0. ;
RV1.'INCO'.VITESSE= VN ;
RV1.'INCO'.DIFFU = (NOEL $DOMTOT DIF) ;
RV1.'FIDT'= 1   ;
 
*-------------------- EXECUTION : CONCENTRATION --------------------------------
 
EXEC RV1 ;
 
*-------------------- DESSINS --------------------------------------------
 
SI GRAPH ;
 
  UET1NCH = ELNO $CLEXT1 UET1 ;
  UET2NCH = ELNO $CLEXT2 UET2 ;
  AKE1 = ELNO $CLEXT1 AK1 ;
  AKE2 = ELNO $CLEXT2 AK2 ;
 
  TITRE 'VITESSES DES PARTICULES' ;
  UNCH = VECT (RV1.'INCO'.VITESSE) 1. UX UY VERT ;
  TRAC UNCH BORDOM ;
 
  TITRE 'UET LE LONG DE LA COUCHE LIMITE 1 '  ;
  EVUT1CS = EVOL CHPO UET1NCH CLEXT1 ;
  DESS  EVUT1CS ;
 
  TITRE 'AK LE LONG DE LA COUCHE LIMITE 1 '  ;
  EVAK1CS = EVOL CHPO AKE1 CLEXT1 ;
  DESS  EVAK1CS ;
 
  TITRE 'UET LE LONG DE LA COUCHE LIMITE 2 '  ;
  EVUT2CS = EVOL CHPO UET2NCH CLEXT2 ;
  DESS EVUT2CS ;
 
  TITRE 'AK LE LONG DE LA COUCHE LIMITE 2 '  ;
  EVAK2CS = EVOL CHPO AKE2 CLEXT2 ;
  DESS EVAK2CS ;
 
  TITRE 'CONCENTRATIONS' ;
  TRAC DOMTOT (RV1.INCO.'CN') BORDOM ;
 
FINSI ;
 
*-------------------- TEST SUR LA VITESSE DE DEPOT -----------------------
 
P1 =  1.0    (E/10.) ;
P2 = (E/10.)  2.5    ;
PT1 = (DOMA $CLEXT1 CENTRE) 'POIN' 'PROC' P1 ;
PT2 = (DOMA $CLEXT2 CENTRE) 'POIN' 'PROC' P2 ;
VD1 = 'EXTR' AK1 SCAL  PT1 ;
VD2 = 'EXTR' AK2 SCAL  PT2 ;
DV1 = (VD1 - 3.52115E-05)/3.52115E-05 ;
DV2 = (VD2 - 6.66232E-08)/6.66232E-08 ;
 
SI ( (ABS DV1) > 0.05 ) ;
     ERREUR 5 ;
FINSI ;
 
SI ( (ABS DV2) > 0.05 ) ;
     ERREUR 5 ;
FINSI ;
 
FIN ;