Télécharger back_proj_8.dgibi

Retour à la liste

* fichier :  back_proj_8.dgibi
****************  CAS TEST : back_proj_8.dgibi  ********************
*
* Ce test permet de vérifier le bon fonctionnement des opérateurs
* utilisés pour la résolution des équations de NAVIER_STOKES en EF
* par un algorithme de projection algébrique incrémental.
*
* Le cas étudié est celui d'un écoulement laminaire dans un canal
* en présence d'une marche descendante. On teste la position du
* point de réattachement par rapport au pied de la marche.
*
* Référence : K. Morgan, J. Périaux anf F. Thomasset, editors,
* Analysis of laminar flow over a backward facing step, Vol9 of
* Notes on numerical Fluid Mechanics, Vieweg, 1984.
*
*-------------------------------------------------------------------
* Auteur : F.Dabbene (LTMF)                                    03/03
*        + S.Gounand (LTMF) (copie de back_impl_4 + ajout des
*                           trois lignes pour la projection) 12/2007
*-------------------------------------------------------------------
*
'SAUT' 'PAGE' ;
'OPTION' 'DIME' 2 'ELEM' 'QUA8' 'ECHO' 0 ;
*
*- Pilotage du calcul
*
*  DISCR : Element en vitesse (LINE/MACRO/QUAF)
*  KPRES : Element en pression (CENTRE/CENTREP1, pas MSOMMET)
*  KSUPG : Méthode de décentrement (CENTREE/SUPG/SUPGDC)
*  KMASS : Matrice masse lumpée ou non (EFM1/EF)
*
*  DT      : Valeur du pas de temps
*  EPS0    : Test d'arret
*  RAF     : Taille de maille définie par 0.1*raf
*
*  FREQ0   : Fréquence d'évaluation du résidu
*  GRAPH   : Booleen pour les tracés
*  COMPLET : Booleen volume de calcul
*
DISCR = 'QUAF' ;
KPRES = 'CENTREP1' ;
'SI' ('EGA' kpres 'MSOMMET') ;
     'ERRE' 'PRESSION MSOMMET interdite en implicite' 5 ;
'FINSI' ;
KSUPG = 'CENTREE' ;
KMASS = 'EF' ;
GRAPH   = FAUX ;
COMPLET = FAUX ;
'SI' COMPLET ;
   DT    = 2. ;
   FREQ0 = 1  ;
   EPS0  = 1.D-5 ;
   RAF   = 1. ;
'SINO' ;
   DT    = 2. ;
   FREQ0 = 2 ;
   EPS0  = 1.D-3 ;
   RAF   = 3. ;
'FINSI' ;
*
*---------------------------------------------------------------------------
* Recherche du point de réattachement (point où dUx/dy=0)
*---------------------------------------------------------------------------
* 1/ Après avoir calculé le gradient de Ux, on ne conserve que les valeurs
*    sur la frontière qui nous intéresse sous la forme d'une évolution.
* 2/ On borne l'évolution au voisinage du point de réattachement afin
*    d'avoir une variation monotone sur l'intervale de dUx/dy.
* 3/ On recherche par interpolation le zero de la fonction (dUx/dy(s),s)
*    La valeur obtenu est l'abscisse curviligne cherchée
*
* Remarques :
*    (i)  1/ permettrait de calculer le coeff de frottement à la paroi :
*         il suffirait de diviser EV1 par le bon coefficient (Re/2 ici)
*    (ii) Il est impératif que dUx/dy soit monotone sur l'intervale 2/ afin
*         qu'il y ait unicité du zero (principe des valeurs intermédiaires)
*---------------------------------------------------------------------------
*
'DEBPROC' attac ;
* 1/
Ux    = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UX' 0. RV.'INCO'.'UN' ;
DUxDY = 'EXCO' 'UY' ('KOPS' Ux 'GRADS' $DOMTOT) 'SCAL' ;
EV1   = 'EVOL' 'CHPO' DUxDY BOTTOM ;
* 2/
EV2   = 'EXTR' EV1 'APRE' 5.  ;
EV3   = 'EXTR' EV2 'AVAN' 10. ;
* 3/
LX1   = 'EXTR' EV3 'ABSC' ;
LY1   = 'EXTR' EV3 'ORDO' ;
LY1 = 'ORDO' LY1 ;
Ymin  = 'MINI' LY1 ;
Ymax  = 'MAXI' LY1 ;
Delta = Ymin * Ymax ;
'SI' (Delta < 0.) ;
   Y0 = 0. ;
'SINON' ;
   'MESS' 'Fonction non monotone' ;
   Y0 = Ymin '+' Ymax '/' 2. ;
*   'ERRE''Fonction non monotone' 5 ;
'FINSI' ;
V1 = 'IPOL' Y0  LY1 LX1 ;
'FINP' V1 ;
*
*==================================================================
* Calcul du résidu basé sur la composante horizontale de la vitesse
* et arrêt suivant un critère transmis
*==================================================================
* E/  : RVX      : TABLE     : TABLE des données créées par EQEX
*                  ARG1      : Fréquence d'impression
*                  ARG2      : Critère d'arrêt
*                  ARG3      : Booleen de tracé
*  /S : MAT1     : MATRIK    : Objet vide
*  /S : CHP1     : CHPO      : Objet vide
*------------------------------------------------------------------
* MAT1 et CHP1 permettent d'assurer la compatibilité des opérateurs
* de discrétisation avec les procédures personnelles
*------------------------------------------------------------------
'DEBPROC' residu rvx*table ;
   RV    = RVX . 'EQEX' ;
   FREQ  = RVX . 'ARG1' ;
   EPS0  = RVX . 'ARG2' ;
   GRAPH = RVX . 'ARG3' ;
   NITER = RV . 'NITER' ;
   DD = RV . 'PASDETPS' . 'NUPASDT' ;
   NN = DD '/' FREQ ;
   L0 = 'EGA' (DD '-' (FREQ*NN)) 0 ;
   'SI' L0 ;
      RANG0 = RV . 'PASDETPS' . 'NUPASDT' ;
      TIME0 = RV . 'PASDETPS' . 'TPS'     ;
      UN0   = 'EXCO' 'UX' RV . 'INCO' . 'UN'  'SCAL' ;
      UNM0  = 'EXCO' 'UX' RV . 'INCO' . 'UN2' 'SCAL' ;
      ERR0  = ('MAXIMUM' ('ABS' (UN0 '-' UNM0))) '+' 1.D-20 ;
      ERR10 = ('LOG' ERR0 ) '/' ('LOG' 10.) ;
      'MESSAGE' 'Résidu en vitesse suivant X au pas'
           RANG0 '(t=' TIME0 ') :' ERR0 ':' ERR10 ;
      RV . 'INCO' . 'IT'  = RV . 'INCO' . 'IT' 'ET' ('PROG' RANG0) ;
      RV . 'INCO' . 'TI'  = RV . 'INCO' . 'TI' 'ET' ('PROG' TIME0) ;
      RV . 'INCO' . 'ER'  = RV . 'INCO' . 'ER' 'ET' ('PROG' ERR10) ;
      V1 = attac ;
      RV . 'INCO' . 'POSI' = RV . 'INCO' . 'POSI' 'ET' ('PROG' V1) ;
      Y1 = ('LOG' EPS0) '/' ('LOG' 10) ; Y2 = 0. ;
     'SI' GRAPH ;
          EV1 = 'EVOL' 'MANU' (RV . 'INCO' . 'IT')(RV . 'INCO' . 'ER') ;
         'DESSIN' EV1 'YBOR' Y1 Y2 'NCLK' ;
      'FINSI' ;
      'SI' ((ERR10 < Y1) 'ET' (DD > ('MAXI' ('LECT' 10 FREQ)))) ;
         RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ;
      'FINSI' ;
   'FINSI' ;
   RV . 'INCO' . 'UN2' = 'COPIER' RV . 'INCO' . 'UN' ;
   mat1 chp1 = 'KOPS' 'MATRIK' ;
'FINP' mat1 chp1 ;
 
*
*==========
* Maillage
*==========
*
*
*-----------------------------------------------------------------------
*
*
*  <    L1     X                L2 - L1                           >
*
*  7___________6__________________________________________________5
*          ^                                     ^
*  INLET   H1
*          v                                     H2               OUTLET
*  1___________2
*              |                                 v
*              3__________________________________________________4
*                                   BOTTOM
*
*
*-----------------------------------------------------------------------
*
*
* L1   : Longueur de la section d'entrée (avant la marche)
* L2   : Longueur de la totalité du dispositif
* H1   : Hauteur de la section d'entrée
* H2   : Hauteur de la section de sortie
* d1   : Dimension caractéristique d'une maille
L1 = 3.0  ;
L2 = 22.0 ;
H1 = 1.0  ;
H2 = 1.5  ;
d1 = 0.1 * raf ;
*
* H2-H1 : Hauteur de la marche servant à l'adimensionnalisation
HDIM = H2 - H1 ;
L1 = L1 / HDIM ;
L2 = L2 / HDIM ;
H1 = H1 / HDIM ;
H2 = H2 / HDIM ;
d1 = d1 / HDIM ;
*
* Points du maillage
p1 = 0. (H2-H1) ;
p2 = L1 (H2-H1) ;
p3 = L1 0. ;
p4 = L2 0. ;
p5 = L2 H2 ;
p6 = L1 H2 ;
p7 = 0. H2 ;
p8 = L2 (H2-H1) ;
*
* Section d'éntrée
p1p2  = p1 'DROIT' 'DINI' d1 'DFIN' d1 p2 ;
p2p6  = p2 'DROIT' 'DINI' d1 'DFIN' d1 p6 ;
p6p7  = p6 'DROIT' 'DINI' d1 'DFIN' d1 p7 ;
p7p1  = p7 'DROIT' 'DINI' d1 'DFIN' d1 p1 ;
mesh1 = 'DALL' p1p2 p2p6 p6p7 p7p1 'PLAN' ;
*
* Section de sortie
p6p2  = 'INVE' p2p6 ;
p2p3  = p2 'DROIT' 'DINI' d1 'DFIN' d1 p3 ;
p3p4  = p3 'DROIT' 'DINI' d1 'DFIN' d1 p4 ;
p4p8  = p4 'DROIT' 'DINI' d1 'DFIN' d1 p8 ;
p8p5  = p8 'DROIT' 'DINI' d1 'DFIN' d1 p5 ;
p5p6  = p5 'DROIT' 'DINI' d1 'DFIN' d1 p6 ;
mesh2 = 'DALL' (p6p2 'ET' p2p3) p3p4 (p4p8 'ET' p8p5) p5p6 'PLAN' ;
*
*
*=======================
* Modèles NAVIER_STOKES
*=======================
*
*
* Définition des équations
* $DOMTOT : Modèle volumique défini sur le maillage complet
*
* Conditions aux limites en vitesse :
* $INLET  : Modèle surfacique défini à l'entrée fluide (Poiseuille)
* $WALL   : Modèle surfacique défini sur les murs (adhérence en paroi)
*
* Conditions aux limites en pression :
* $OUTLET : Modèle surfacique défini à la sortie fluide (sert à imposer
*           la pression de sortie en pression continue MSOMMET)
*
* Post-traitement
* $BOTTOM : Modèle surfacique défini sur le plancher après la marche
*           (sert à évaluer la position du point de réattachement)
*
DOMTOT  = 'CHAN' 'QUAF' (mesh1 ET mesh2) ;
$DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' DISCR ;
$INLET  = 'MODE' p7p1   'NAVIER_STOKES' DISCR ;
$OUTLET = 'MODE' (p4p8 'ET' p8p5) 'NAVIER_STOKES' DISCR ;
$BOTTOM = 'MODE' p3p4 'NAVIER_STOKES' DISCR ;
$WALL   = 'MODE' (p1p2 ET p2p3 ET p3p4 ET p5p6 ET p6p7)
                 'NAVIER_STOKES' DISCR ;
*
* Elimination ad hoc
* (En 2D, il faut éliminer les points centres des modèles surfaciques
* avec les points faces des modèles volumiques à cause des MACROs)
FDOMTOT = 'DOMA' $DOMTOT 'FACE'   ;
CINLET  = 'DOMA' $INLET  'CENTRE' ;
COUTLET = 'DOMA' $OUTLET 'CENTRE' ;
CBOTTOM = 'DOMA' $BOTTOM 'CENTRE' ;
CWALL   = 'DOMA' $WALL   'CENTRE' ;
'ELIM' (FDOMTOT 'ET' CINLET 'ET' COUTLET 'ET' CWALL 'ET' CBOTTOM) EPS0 ;
*
* On écrase les anciens maillages afin d'éviter toute ambiguitée
DOMTOT = 'DOMA' $DOMTOT 'MAILLAGE' ;
INLET  = 'DOMA' $INLET  'MAILLAGE' ;
OUTLET = 'DOMA' $OUTLET 'MAILLAGE' ;
BOTTOM = 'DOMA' $BOTTOM 'MAILLAGE' ;
WALL   = 'DOMA' $WALL   'MAILLAGE' ;
*
* Maillage pour d'éventuelles conditions aux limites en pression
'SI' ('EGA' kpres 'MSOMMET') ;
   OUTLETP = 'DOMA' $outlet KPRES ;
'SINON' ;
   OUTLETP = 'DOMA' $domtot KPRES ;
'FINSI' ;
*
*===========================
* Description des équations
*===========================
*
* Grandeurs adimensionnées
Umax = 1.0  ;
Re   = 150. ;
*
* Profil de vitesse parabolique à l'entrée
YINLET = 'COOR' 2 INLET ;
YMAX   = 'MAXI' YINLET ;
YMIN   = 'MINI' YINLET ;
UIN = (YINLET '-' YMAX) '*' (YINLET '-' YMIN) ;
UIN = UIN '*' (-4.0*Umax/((YMAX-YMIN)*(YMAX-YMIN))) ;
UIN = 'NOMC' 'UX' UIN 'NATU' 'DISCRET' ;
VIN = 'KCHT' $INLET 'SCAL' 'SOMMET' 'COMP' 'UY' 0.  ;
*
* Description du système en vitesse-pression
RV = 'EQEX' $DOMTOT 'ITMA' 50000 'ALFA' 1. 'FIDT' 100
           'OPTI' 'EF' 'IMPL'  KSUPG KPRES
           'ZONE' $DOMTOT 'OPER' residu FREQ0 EPS0 GRAPH
           'ZONE' $DOMTOT 'OPER' 'NS' 1. 'UN' (1./Re) 'INCO' 'UN'
           'ZONE' $DOMTOT 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PRES'
*                                          Conditions aux limites par défaut
*                           'ZONE' $OUTLET 'OPER' 'TOIM' (0. 0.) 'INCO' 'UN'
      'OPTI' KMASS KSUPG
      'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'UN' DT 'UN' (1./Re) 'INCO' 'UN'
;
*
* Description des conditions aux limites
RV = 'EQEX' RV
           'CLIM' 'UN' 'UIMP' WALL  0.  'UN' 'VIMP' WALL  0.
                  'UN' 'UIMP' INLET UIN 'UN' 'VIMP' INLET VIN
;
*
* Déclaration des inconnues et initialisations (table INCO)
RV . 'INCO'          = 'TABLE' 'INCO' ;
RV . 'INCO' . 'UN'   = 'KCHT' $DOMTOT 'VECT' 'SOMMET' (0. 0.) ;
* BUG QUAF/CENTRE rv . 'INCO' . 'PRES' = 'KCHT' $DOMTOT 'SCAL' 'CENTREP0' 0. ;
RV . 'INCO' . 'PRES' = 'KCHT' $DOMTOT 'SCAL' KPRES 0. ;
*
* Champs supplémentaires pour la procédure residu
RV . 'INCO' . 'IT'   = 'PROG' ;
RV . 'INCO' . 'TI'   = 'PROG' ;
RV . 'INCO' . 'ER'   = 'PROG' ;
RV . 'INCO' . 'POSI' = 'PROG' ;
RV . 'INCO' . 'UN2'  = 'KCHT' $DOMTOT 'VECT' 'SOMMET' (1.D-3 1.D-3) ;
*
   rv . 'GPROJ' = 'TABLE' ;
   rv . 'GPROJ' . 'NOMVIT'  = 'UN' ;
   rv . 'GPROJ' . 'NOMPRES' = 'CHAINE' 'PRES' ;
*
EXEC RV ;
*
*=================
* Post-traitement
*=================
*
* Localisation du point de réattachement
* Tracés de la vitesse et de la pression
*
CNT1 = WALL ;
NLI0 = 14   ;
*
* Point de réattachement (point où dUx/dy=0)
V1    = attac ;
'SAUT' 10 'LIGNE' ;
'MESS' 'ABSCISSE DU POINT DE REATTACHEMENT :' V1 ;
'SI' GRAPH ;
  EV1 = 'EVOL' 'MANU' RV . 'INCO' . 'TI' RV . 'INCO' . 'POSI' ;
  'DESS' EV1 'MIMA' 'GRIL' 'TITR' 'Localisation reattachement = f(t)' ;
*
* Vitesse
  un  = RV . 'INCO' . 'UN' ;
  vun = 'VECT' UN 0.5 'UX' 'UY' 'JAUNE' ;
  trace vun DOMTOT CNT1 'TITR' 'Vitesse' ;
*
* Pression
  pe = RV . 'INCO' . 'PRES' ;
  'SI' ('EGA' kpres 'MSOMMET') ;
       mp = 'DOMA' $DOMTOT 'MMAIL' ;
       trace pe mp ('CONTOUR' mp) ;
  'SINON' ;
       pn = 'ELNO' $DOMTOT ('KCHT' $DOMTOT 'SCAL' kpres 0. pe) kpres ;
  'FINSI' ;
'FINSI' ;
*
*========================
* Test de non régression
*========================
*
'SI' COMPLET ;
   TEST = ('ABS' (V1 '-' 6.4842)) < 0.0005 ;
'SINON' ;
*   TEST = ('ABS' (V1 '-' 6.3556)) < 0.0005 ;
   TEST = ('ABS' (V1 '-' 6.4821)) < 0.0005 ;
'FINSI' ;
*'MESSAGE' 'Méthode de projection algebrique incrementale' ;
'MESSAGE' 'Méthode de double projection algebrique incrementale' ;
'MESSAGE' 'Element Vitesse  : ' DISCR ;
'MESSAGE' 'Element Pression : ' KPRES ;
'MESSAGE' 'Décentrement     : ' KSUPG ;
'MESSAGE' 'Matrice masse    : ' KMASS ;
'SI' TEST ;
   'ERREUR' 0 ;
'SINON' ;
   'ERREUR' 5 ;
'FINSI' ;
*
'FIN' ;
 
 
 
 
 
 
 
 
 
 
 
 

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