Download fcourant.dgibi

Back to the list

* fichier : fcourant.dgibi
'OPTI' 'ECHO' 0 ;
************************************************************************
* NOM         : FCOURANT
* DESCRIPTION : Teste la procédure calculant la fonction de courant
*               sur un poiseuille en 2D et 2D axi
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 21/12/2007, version initiale
* HISTORIQUE : v1, 21/12/2007, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
'SAUTER' 2 'LIGNE' ;
'MESSAGE' ' Execution de fcourant.dgibi' ;
'SAUTER' 2 'LIGNE' ;
*
interact = FAUX ;
graph    = FAUX ;
*
'OPTION' 'DIME' 2 'ELEM' 'TRI6' ;
*
'SI' ('NON' interact) ;
  'OPTION' 'TRAC' 'PS' ;
'SINON' ;
  'OPTION' 'TRAC' 'X' ;
'FINSI' ;  
*
lok = VRAI ;
*
'REPETER' iaxi 2 ;
   iiaxi = &iaxi ;
   'SAUTER' 1 'LIGNE' ;
   'SI' ('EGA' iiaxi 1) ;
      'OPTI' 'MODE' 'PLAN' ;
      'MESSAGE' 'Mode plan ' ;
   'SINON' ;
      'OPTI' 'MODE' 'AXIS' ;      
      'MESSAGE' 'Mode axi ' ;
   'FINSI' ;
   'SAUTER' 1 'LIGNE' ;
*
laxi = ('EGA' ('VALEUR' 'MODE') 'AXIS') ;
*
* Maillage
*
a = 1.2 ; b = 2.5 ;
nx = 5 ; ny = 7 ;
*nx = 10 ; ny = 14 ;
dpdz = -1.5 ; mu = 0.9 ;
vtoim = 0. ('*' dpdz b) ;
disv = 'QUAF' ;
disp = 'CENTREP1' ;
*
pA = 0. 0. ; pB = a 0. ; pC = a b ; pD = 0. b ;
bas = 'DROIT' nx pA pB ; dro = 'DROIT' ny pB pC ;
hau = 'DROIT' nx pC pD ; gau = 'DROIT' ny pD pA ;
*
cmt = bas 'ET' dro 'ET' hau 'ET' gau ;
mt  = 'SURFACE' cmt ;
_mt = 'CHANGER' mt 'QUAF' ;
*
$mt  = 'MODE' _mt  'NAVIER_STOKES' disv ;
$hau = 'MODE' hau 'NAVIER_STOKES' disv ;
*
mt = 'DOMA' $mt 'MAILLAGE' ;
*
* Solution exacte 2D 
*    u_x = 0
*    u_y = vmax (1 - x^2/a^2)
*          vmax = -dp/dz (a^2 / 2 mu)
*    p   = dp/dz y
*    psi = vmax x (1 - x^2/3a^2)
*
'SI' ('NON' laxi) ;
   vmax = '*' ('*' dpdz -1.) ('/' ('*' a a) ('*' mu 2)) ;
   x = 'COORDONNEE' 1 mt ;
   y = 'COORDONNEE' 2 ('DOMA' $mt disp) ;
   cy = '+' ('*' ('/' ('*' x x) ('*' a a)) -1.) 1. ;
   unx = 'NOMC' 'UY' ('*' vmax cy) ;
   pnx = '*' dpdz y ;
   cpsi = '+' ('*' ('/' ('*' x x) ('*' ('*' a a) 3)) -1.) 1. ;
   psix = '*' vmax ('*' x cpsi) ;
'SINON' ;   
*
* Solution exacte 2D axi
*    u_r = 0
*    u_z = vmax (1 - r^2/a^2)
*          vmax = -dp/dz (a^2 / 4 mu)
*    p   = dp/dz y
*    psi = vmax x (1 - x^2/3a^2)
*
   vmax = '*' ('*' dpdz -1.) ('/' ('*' a a) ('*' mu 4)) ;
   r = 'COORDONNEE' 1 mt ;
   z = 'COORDONNEE' 2 ('DOMA' $mt disp) ;
   cz = '+' ('*' ('/' ('*' r r) ('*' a a)) -1.) 1. ;
   unx = 'NOMC' 'UY' ('*' vmax cz) ;
   pnx = '*' dpdz z ;
   cpsi = '+' ('*' ('/' ('*' r r) ('*' ('*' a a) 2)) -1.) 1. ;
   psix = '*' ('*' vmax PI) ('*' ('*' r r) cpsi) ;
'FINSI' ;   
*
* Poiseuille
*
rv = 'EQEX' 'ITMA' 1 'ALFA' 1. 'NITER' 1
     'OPTI' 'EF' 'IMPL' 'CENTREE' disp
     'ZONE' $mt  'OPER' 'NS' 1. 'UN' mu 'INCO' 'UN'
     'ZONE' $mt  'OPER' 'KBBT' 1. 'INCO' 'UN' 'PN'
     'ZONE' $hau 'OPER' 'TOIM' vtoim 'INCO' 'UN'
     'CLIM' 'UN' 'UIMP' dro 0. 'UN' 'VIMP' dro 0.
     'CLIM' 'UN' 'UIMP' gau 0. 
   ;
rv . 'INCO' = 'TABLE' 'INCO' ;
RV . 'INCO' . 'UN'   = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ;
RV . 'INCO' . 'PN' = 'KCHT' $mt 'SCAL' disp 0. ;
*
EXEC rv ;
*
un  = rv . 'INCO' . 'UN' ;
pn  = rv . 'INCO' . 'PN' ;
clpsi = 'MANUEL' 'CHPO' pA 1 'SCAL' 0. ;
psi = FCOURANT _mt disv un clpsi ;
* 
'SI' graph ;
   titg = 'CHAINE' ' mod = ' ('VALEUR' 'MODE') ;
   vun  = 'VECT' un 'UX' 'UY' 'JAUN' ;
   tit = 'CHAINE' 'Vitesse' titg ;
   'TRACER' vun mt cmt 'TITR' tit ;
   ppn = 'ELNO' $mt pn disp ;
   tit = 'CHAINE' 'Pression' titg ;
   'TRACER' ppn mt  'TITR' tit ;
   tit = 'CHAINE' 'Fonction de courant' titg ;
   'TRACER' psi mt  'TITR' tit ;
'FINSI' ;
*
erru = 'MAXIMUM' ('-' un unx) 'ABS' ;
errp = 'MAXIMUM' ('-' pn pnx) 'ABS' ;
errf = 'MAXIMUM' ('-' psi psix) 'ABS' ;
*
'MESSAGE' ('CHAINE' 'erru = ' erru) ;
'MESSAGE' ('CHAINE' 'errp = ' errp) ;
'MESSAGE' ('CHAINE' 'errf = ' errf) ;
*
* le 21/12/2007 sur ma machine :
* Mode plan
*   erru =  4.21885E-15
*   errp =  3.15303E-14
*   errf =  6.36767E-04
* Mode axi
*   erru =  2.66454E-15
*   errp =  2.62013E-14
*   errf =  5.45703E-03
*
tolup = 1.D-12 ;
'SI' laxi ;
   tolf = 6.D-3 ;
'SINON' ;
   tolf = 7.D-4 ;
'FINSI' ;
*
tst = ('<' erru tolup) ;
'SI' ('NON' tst) ;
   chmes = 'CHAINE' '!!! Erreur, on na pas erru < ' tolup ;
   'MESSAGE' chmes ;
'FINSI' ;
lok = lok 'ET' tst ;   
*
tst = ('<' errp tolup) ;
'SI' ('NON' tst) ;
   chmes = 'CHAINE' '!!! Erreur, on na pas errp < ' tolup ;
   'MESSAGE' chmes ;
'FINSI' ;
lok = lok 'ET' tst ;   
*
tst = ('<' errf tolf) ;
'SI' ('NON' tst) ;
   chmes = 'CHAINE' '!!! Erreur, on na pas errf < ' tolf ;
   'MESSAGE' chmes ;
'FINSI' ;
lok = lok 'ET' tst ;
*
'FIN' iaxi ;
*
'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' ;   
*
* End of dgibi file FCOURANT
*
'FIN' ;

 

© Cast3M 2003 - All rights reserved.
Disclaimer