Accueil | Cast3M

* fichier : cyltest6.dgibi
************************************************************************
************************************************************************
'OPTI' 'ECHO' 0 ;
************************************************************************
*                                                                      *
*             Ecoulement autour d'un cylindre circulaire               *
* Résolution des équations de Navier Stokes laminaires instationnaires *
*                                                                      *
* Pompé sur cyltest.dgibi                                              *
* Ajouts Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)                     *
*        mél : gounand@semt2.smts.cea.fr                               *
*        en date du  20/12/2007                                        *
*        teste une méthode couplée et une méthode de projection        *
*        algébrique incrémentale : vérification du bilan des forces    *
*        locales et de la force s'exercant sur le cylindre             *
* SG 08/10/2024 : passage aux quadratiques pour la discrétisation      *
*                 du cercle en TRIA RAFT                               *
*                 Ajout double projection                              *
************************************************************************
'SAUTER' 2 'LIGNE' ;
'MESSAGE' ' Execution de cyltest6.dgibi' ;
'SAUTER' 2 'LIGNE' ;
*
interact = FAUX ;
graph    = FAUX ;
*
*****************************************
*                                       *
*           Procedure tracvit           *
*                                       *
*****************************************
'DEBPROC' TRACVIT;
'ARGUMENT' RVX*'TABLE';
rv = rvx . 'EQEX' ;
umax = 2. ;
CHUN = 'VECTEUR' (rv.inco.un)  (1. '/' (1. * umax)) UX UY JAUN;
'TRACER' CHUN mt cmt 'TITRE' 'Champ de vitesse'
         'NCLK' ;
as2 ama1 = 'KOPS' 'MATRIK';
'RESPRO' as2 ama1;
'FINP';
*****************************************
*                                       *
*   Proc    Bloc de pression nul        *
*                                       *
*****************************************
 'DEBPROC' BLOCPNUL ;
 'ARGUMENT' rvx*'TABLE' ;
 RV = RVX . 'EQEX';
*
 'SI' ('NON' ('EXISTE' rv 'MATPNULL')) ;
*
* Creation d'un bloc de pression nulle
*
   'MESSAGE' ('CHAINE' 'Calcul du bloc de pression') ;
   $mt = rvx . 'DOMZ' ;
   rt = 'EQEX' 'OPTI' 'EF' 'IMPL' 'CENTREE' 'CENTREP1'
               'ZONE' $mt
               'OPER' 'KMAB' 0. 'INCO' 'UN' 'PN' ;
   Rt.'INCO' = RV . 'INCO' ;
   smbb matb = kmab rt . '1KMAB' ;
   chpv = 'KCHT' $mt 'VECT' 'SOMMET' 'COMP' '1UN' '2UN' (0. 0.) ;
   matpres = 'KOPS' 'CMCT' matb chpv matb ;
   rv . 'MATPNULL' = matpres ;
 'SINON' ;
   matpres =rv . 'MATPNULL'  ;
 'FINSI' ;
 smbvide matvide = 'KOPS' 'MATRIK' ;
'FINPROC' smbvide matpres ;
*
* Options globales
*
'OPTION' 'DIME' 2 'ELEM' 'QUA8' 'ECHO' 0 ;
'SI' interact ;
   'OPTION' 'TRAC' 'X'  ;
'SINON' ;
   'OPTION' 'TRAC' 'PSC';
'FINSI' ;
'OPTION' isov suli;
disv  = 'QUAF' ;
DISP  = 'CENTREP1';
* iraff   : raffinement du maillage
* mproj   : méthode de projection ou couplée
* miter   : méthode itérative ou directe
iraff   = 1 ;
impkres = 0 ;
*
lok = VRAI ;
*
'REPETER' imet 6 ;
   'SAUTER' 1 'LIGN' ;
   iimet = &imet ;
   'SI' ('EGA' iimet 1) ;
      mproj = FAUX ;
      miter = FAUX ;
   'FINSI' ;
   'SI' ('EGA' iimet 2) ;
      mproj = FAUX ;
      miter = VRAI ;
   'FINSI' ;
   'SI' ('EGA' iimet 3) ;
      mproj = VRAI ; lprec = FAUX ;
      miter = FAUX ;
   'FINSI' ;
   'SI' ('EGA' iimet 4) ;
      mproj = VRAI ; lprec = FAUX ;
      miter = VRAI ;
   'FINSI' ;
   'SI' ('EGA' iimet 5) ;
      mproj = VRAI ; lprec = VRAI ;
      miter = FAUX ;
   'FINSI' ;
   'SI' ('EGA' iimet 6) ;
      mproj = VRAI ; lprec = VRAI ;
      miter = VRAI ;
   'FINSI' ;
   'SI' mproj ;
      'SI' lprec ;
         'MESSAGE' 'Methode projection preconditionnee' ;
      'SINON' ;
         'MESSAGE' 'Methode projection' ;
      'FINSI' ;
   'SINON' ;
      'MESSAGE' 'Methode directe' ;
   'FINSI' ;
   'SI' miter ;
      'MESSAGE' 'Solveur itératif' ;
   'SINON' ;
      'MESSAGE' 'Solveur direct' ;
   'FINSI' ;
   'SAUTER' 1 'LIGN' ;
*============================
*  Construction du maillage:
*============================
*  Construction des points:
nc = 4 ;
nx = 11 ; ny = 6 ;
nc = '*' nc iraff ;
nx = '*' nx iraff ;
ny = '*' ny iraff ;
*
p0 = 0. 0. ;
p1 = 0.5 0. ;
p2 = 0. 0.5 ;
p3 = -0.5 0. ;
p4 =  0. -0.5 ;
*
c1 = 'CERCLE' nc p1 p0 p2 ;
c2 = 'CERCLE' nc p2 p0 p3 ;
c3 = 'CERCLE' nc p3 p0 p4 ;
c4 = 'CERCLE' nc p4 p0 p1 ;
*
ligc = c1 'ET' c2 'ET' c3 'ET' c4 ;
ligc = 'INVERSE' ligc ;
*
*
pA = -6. -6. ;
pB = 15. -6. ;
pC = 15. 6. ;
pD = -6. 6. ;
bas = pA 'DROIT' nx pB ;
dro = pB 'DROIT' ny pC ;
hau = pC 'DROIT' nx pD ;
gau = pD 'DROIT' ny pA ;
*
lige = bas 'ET' dro 'ET' hau 'ET' gau ;
*
cntt = ligc 'ET' lige ;
*
mt = tria cntt ;
xmt ymt = 'COORDONNEE' mt ;
rmt = '**' ('+' ('**' xmt 2) ('**' ymt 2)) 0.5D0 ;
*rmt = '*' rmt 3. ;
rmt = '/' rmt 3. ;
rmt = '/' rmt iraff ;
mtr = 'RAFT' rmt mt ;
'SI' graph ;
   'TRACER' mtr 'TITR' ('CHAINE' 'nbl=' ('NBEL' mtr)) ;
'FINSI' ;
*'OPTION' 'DONN' 5 ;
cercle = ligc ;
* Création du domaine total:
mt = mtr;
cmt = 'CONTOUR' mt;
*==============================================
*  Changement des éléments du maillage en QUAF:
*==============================================
 
_mt     = 'CHANGER' mt 'QUAF';
_cmt    = 'CHANGER' cmt 'QUAF';
_cercle = 'CHANGER' cercle 'QUAF';
_entree = 'CHANGER' gau 'QUAF';
_sortie = 'CHANGER' dro 'QUAF';
_hau  = 'CHANGER' hau           'QUAF';
_bas  = 'CHANGER' bas           'QUAF';
'ELIM' 1.E-5 (_mt et _cercle et _entree 'ET' _sortie
              'ET' _cmt 'ET' _hau 'ET' _bas);
 
NN   = 'NBEL' mt;
'MESSAGE' 'Nombre d éléments du maillage' NN;
 
*=======================================
*  Formulation du domaine Navier Stokes:
*=======================================
 
$mt     = 'MODE' _mt 'NAVIER_STOKES' disv ;
$cmt    = 'MODE' _cmt 'NAVIER_STOKES' disv ;
$cercle = 'MODE' _cercle 'NAVIER_STOKES' disv ;
$entree = 'MODE' _entree 'NAVIER_STOKES' disv ;
$sortie = 'MODE' _sortie 'NAVIER_STOKES' disv ;
$hau    = 'MODE' _hau    'NAVIER_STOKES' disv ;
$bas    = 'MODE' _bas 'NAVIER_STOKES' disv ;
*
mt     = 'DOMA' $mt     'MAILLAGE' ;
cmt    = 'DOMA' $cmt    'MAILLAGE' ;
cercle = 'DOMA' $cercle 'MAILLAGE' ;
entree = 'DOMA' $entree 'MAILLAGE' ;
sortie = 'DOMA' $sortie 'MAILLAGE' ;
hau    = 'DOMA' $hau    'MAILLAGE' ;
bas    = 'DOMA' $bas    'MAILLAGE' ;
 
*====================================
*  Création des tables de résolution:
*====================================
 
dt     = .5D0 ;
*rey    = 30. ;
rey    = 100.;
nu     = 1./rey;
arelax = 1.;
*
niter = 6 ;
nitime = 3 ;
*
RV = 'EQEX' 'OMEGA' 1.0 'NITER' niter  'ITMA' 1  'FIDT' 1 ;
*
RV = 'EQEX' RV
       'OPTI' 'EF' 'CENTREE' 'IMPL'
*        'ZONE' $mt 'OPER' 'NS' 'INV_RE' 'INCO' 'UN'
       'ZONE' $mt 'OPER' 'KONV' 1. 'UN' 'INV_RE' 'INCO' 'UN'
       'OPTI' 'EF' 'CENTREE' 'IMPL'
       'ZONE' $mt 'OPER' 'LAPN' 'INV_RE' 'INCO' 'UN'
       'OPTI' 'EF' 'CENTREE' 'IMPL'
       'ZONE' $mt 'OPER' 'DFDT' 1. 'UNM' 'DT' 'INCO' 'UN'
      'OPTI' 'EF' 'IMPL' 'CENTREE' DISP
      'ZONE' $mt
      'OPER' 'KBBT' 1.    'INCO' 'UN' 'PN' ;
 
'SI' graph ;
   rv = 'EQEX' rv
     'ZONE' $mt 'OPER' TRACVIT ;
'FINSI' ;
 
*  Implantation des conditions aux limites:
 
RV = 'EQEX' RV
      'CLIM' 'UN' 'UIMP' cercle 0. 'UN' 'VIMP' cercle 0.
             'UN' 'UIMP' entree 1. 'UN' 'VIMP' entree 0.
             'UN'                       'VIMP' hau  0.
             'UN'                       'VIMP' bas  0. ;
*
* Choix de la méthode de résolution (couplée, projection)
* et des solveurs (direct, itératif)
*
'SI' ('NON' mproj) ;
*  Solveur pour le système vitesse-pression
   rvm = rv . 'METHINV' ;
   rvm . 'SCALING' = 1 ;
   rvm . 'IMPINV'  = impkres ;
   'SI' ('NON' miter) ;
      rvm . 'TYPINV' = 1 ;
   'SINON' ;
*  On rajoute le bloc de pression nul
      rv = 'EQEX' rv  'ZONE' $mt 'OPER' BLOCPNUL ;
*
      rvm . 'TYPINV' = 4 ;
      rvm . 'LBCG'   = 2 ;
      rvm . 'PRECOND' = 5 ;
      rvm . 'ILUTLFIL' = 3.D0 ;
   'FINSI' ;
'SINON' ;
   rv . 'GPROJ' = 'TABLE' ;
   rv . 'GPROJ' . 'NOMVIT'  = 'UN' ;
   rv . 'GPROJ' . 'NOMPRES' = 'PN' ;
   'SI' ('NON' lprec) ;
      rv . 'GPROJ' . 'NOPREC' = VRAI ;
   'FINSI' ;
*!   rv . 'GPROJ' . 'dblproj' = faux ;
*  Solveur pour les composantes de la vitesse
   rvm = rv . 'METHINV' ;
   rvm . 'SCALING' = 1 ;
   rvm . 'IMPINV'  = impkres ;
   'SI' ('NON' miter) ;
      rvm . 'TYPINV' = 1 ;
   'SINON' ;
      rvm . 'TYPINV' = 4 ;
      rvm . 'LBCG'   = 2 ;
      rvm . 'PRECOND' = 3 ;
   'FINSI' ;
*  Solveur pour l'équation de pression
   rv . 'GPROJ'. 'METHINV' = 'TABLE' 'METHINV' ;
   rvgm = rv . 'GPROJ' . 'METHINV' ;
   rvgm . 'SCALING' = 1 ;
   rvgm . 'IMPINV'  = impkres ;
   'SI' ('NON' miter) ;
      rvgm . 'TYPINV' = 1 ;
   'SINON' ;
      rvgm . 'TYPINV' = 4 ;
      rvgm . 'LBCG'   = 2 ;
      rvgm . 'PRECOND' = 3 ;
*      rvgm . 'PRECOND' = 5 ;
*      rvgm . 'ILUTLFIL' = 1.D0 ;
   'FINSI' ;
'FINSI' ;
*  Création de la table des inconnues et initialisation:
RV.INCO         = 'TABLE' INCO;
*
RV.INCO.'UN'     = 'KCHT' $mt 'VECT' 'SOMMET' (1.E-3 1.E-3);
RV.INCO.'UNM'    = 'KCHT' $mt 'VECT' 'SOMMET' (1.E-3 1.E-3);
rv . 'INCO' . 'G' = (0. -9.8) ;
RV.INCO.'DT'    = dt  ;
RV.INCO.'INV_RE'= nu;
RV.INCO.'PN'  = 'KCHT' $mt 'SCAL' DISP 0.;
*
* Table avec les résultats en vitesse
*
res = 'TABLE' ;
itres = 0 ;
res . itres = 'COPIER' (RV.INCO.'UN') ;
***************************************************
*
* Boucle en temps faite main
*
'TEMPS' 'ZERO' ;
'REPETER' iitime nitime ;
   rv . 'INCO' . 'UNM' = 'COPIER' (res .  itres) ;
   EXEC RV;
   itres = '+' itres 1 ;
   res . itres = 'COPIER' (RV.INCO.'UN') ;
'FIN' iitime ;
TABTPS = TEMP 'NOEC';
tcpu = TABTPS.'TEMPS_CPU'.'INITIAL' ;
*
***************************************************
lastime = '-' ('DIME' res) 1 ;
vit   = res . lastime ;
vitm1 = res . ('-' lastime 1) ;
umax = 2. ;
*
*
* Créer les matrices
*
bidon  matkon  = 'KONV' (rv . '1KONV') ;
bidon  matlap  = 'LAPN' (rv . '2LAPN') ;
chdfdt matdfdt = 'DFDT' (rv . '3DFDT') ;
*
rt = 'EQEX' 'ZONE' $mt
      'OPTI' 'EF' 'IMPL' 'CENTREE' DISP
      'OPER' 'KMBT' -1.    'INCO' 'PN' 'UN'
        ;
rt.'INCO' = rv . 'INCO' ;
bidon matgrp = 'KMBT' (rt . '1KMBT') ;
*
pre = rv . 'INCO' . 'PN' ;
nomavv = 'MOTS' 'UX' 'UY' ;
nomapv = 'MOTS' '1UN' '2UN' ;
nomavp = 'MOTS' 'SCAL' ;
nomapp = 'MOTS' 'PN' ;
vit   = 'NOMC' nomavv nomapv vit ;
vitm1 = 'NOMC' nomavv nomapv vitm1 ;
pre = 'NOMC' nomavp nomapp pre ;
forlap = 'KOPS' matlap '*' vit ;
forkon = 'KOPS' matkon '*' vit ;
forgrp  = 'KOPS' matgrp '*' pre ;
fordfdt = 'KOPS' matdfdt '*' ('-' vit vitm1) ;
*
forlap = ('NOMC' nomapv nomavv forlap) ;
forkon = ('NOMC' nomapv nomavv forkon) ;
forgrp = ('NOMC' nomapv nomavv forgrp) ;
fordfdt = ('NOMC' nomapv nomavv fordfdt) ;
*
fort = forlap '+' forkon '+' forgrp '+' fordfdt ;
*
fmax = 'MAXIMUM' fort 'ABS' ;
echf = '/' 0.5 fmax ;
*
vlap  = 'VECTEUR' forlap echf UX UY VERT ;
vkon  = 'VECTEUR' forkon echf UX UY ROUG ;
vgrp  = 'VECTEUR' forgrp echf UX UY JAUN ;
vdfdt = 'VECTEUR' fordfdt echf UX UY TURQ ;
vf    = 'VECTEUR' fort   echf UX UY BLAN ;
vtot  =  vlap 'ET' vkon 'ET' vgrp 'ET' vdfdt 'ET' vf  ;
*
'SI' graph ;
tit = 'CHAINE' 'vert=lap;roug=kon;jaun=grp;turq=dfdt;blan=reac;fmax='
                fmax ;
'TRACER' vtot _mt ('CONTOUR' _mt) 'TITRE' tit ;
'FINSI' ;
*
pmt = 'CHANGER' mt 'POI1' ;
pcmt = 'CHANGER' cmt 'POI1' ;
pint = 'DIFF' pmt pcmt ;
fortint = 'REDU' fort pint ;
fmaxint = 'MAXIMUM' fortint 'ABS' ;
'MESSAGE' ('CHAINE' 'fmax=' fmax) ;
'MESSAGE' ('CHAINE' 'fmaxint=' fmaxint) ;
errbf   = '*' ('/' fmaxint fmax) 100. ;
'MESSAGE' ('CHAINE' 'erreur bilan forces =' errbf ' %') ;
forcer = 'REDU' fort cercle ;
rfor   = 'RESULT' forcer ;
'MESSAGE' 'Resultante des forces sur le cylindre' ;
rfx = 'MAXIMUM' ('EXCO' 'UX' rfor) ;
rfy = 'MAXIMUM' ('EXCO' 'UY' rfor) ;
'MESSAGE' ('CHAINE' '  en x =' rfx) ;
'MESSAGE' ('CHAINE' '  en y =' rfy) ;
'MESSAGE' ('CHAINE' 'tcpu=' tcpu) ;
*
* Tests
*   Erreur sur le bilan des forces au 08/10/2024
*     2.e-3 % pour la méthode directe
*     3.5   % pour la méthode projection algébrique
*     5.3e-2  % pour la méthode double projection algébrique
*  Valeur de la force de trainée (fx)
*     -7.261-1 pour la méthode directe
*     -7.029e-1 pour la méthode projection algébrique
*     -7.262e-1 pour la méthode double projection algébrique
*
'SI' mproj ;
*   vrbf = 8. ;
* vrbf = 2. ;
   vrbf = 0.1 ;
'SINON' ;
*   vrbf = 56.D-3 ;
   vrbf = 3.D-3 ;
'FINSI' ;
*vrft = -6.9e-1 ;
*tol  =  0.3e-1 ;
*vrft = -6.94e-1 ;
*tol  =  0.02e-1 ;
vrft = -7.261e-1 ;
tol  =  0.002e-1 ;
*
tst = '<' errbf vrbf ;
'SI' ('NON' tst) ;
   cherr = 'CHAINE' '!!! Erreur bilan force > ' vrbf ' %' ;
   'MESSAGE' cherr ;
'FINSI' ;
lok = lok 'ET' tst ;
*
tst = '<' ('ABS' ('-' rfx vrft)) tol ;
'SI' ('NON' tst) ;
   cherr = 'CHAINE' '!!! Erreur force trainee <> ' vrft ' a ' tol
    ' pres' ;
   'MESSAGE' cherr ;
'FINSI' ;
lok = lok 'ET' tst ;
'FIN' imet ;
'SAUTER' 2 'LIGNE' ;
'SI' lok ;
   'MESSAGE' 'Tout sest bien passe' ;
'SINON' ;
   'MESSAGE' 'Il y a eu des erreurs' ;
'FINSI' ;
'SAUTER' 2 'LIGNE' ;
*
'SI' interact ;
   'OPTION' 'DONN' 5 'ECHO' 1 ;
'FINSI' ;
'SI' ('NON' lok) ;
   'ERREUR' 5 ;
'FINSI' ;
* 'FIN' du jeu de données.
fin;