* fichier : allee.dgibi
'OPTION' echo 0;
************************************************************************
* NOM : allee.dgibi
* DESCRIPTION :
*
* Ce jeu de données est dédié à l'écoulement d'un fluide
* incompressible autour d'un cylindre de section
* circulaire. L'écoulement est supposé laminaire. Le but est de
* retrouver les lâchers tourbillonnaires de von Karman et les
* fréquences de ces derniers.
*
* Références :
*
* U. Frisch. Turbulence. Cambridge Universuty Press, 1995.
* H. Schlichting. Boundary-layer theory. McGraw-Hill book company,
* 1968.
*
* LANGAGE : GIBIANE
* AUTEUR : F. PIGEONNEAU (CEA/DRN/DMT/SEMT/LTMF)
* e-mail : frankp@semt2.smts.cea.fr
************************************************************************
* APPELES :
************************************************************************
* SYNTAXE GIBIANE : Uniquement.
************************************************************************
* Description des paramètres et variables locaux:
************************************************************************
*
*
************************************************************************
* VERSION : v1, 15/12/98, version initiale.
* HISTORIQUE : v1, 15/12/98, création.
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de cette procédure afin de faciliter
* la maintenance !
************************************************************************
COMPLET = FAUX ;
'SI' (COMPLET) ;
* Temps final de calcul :
tfin=1. ;
* Définitions des nombres d'éléments:
nx1 = 10 ;
nx2 = 20 ;
nx3 = 40 ;
* Définitions des densites:
dbeg = 0.1 ;
dend = 0.5 ;
'SINON' ;
* Temps final de calcul :
tfin = 0.1 ;
* Définitions des nombres d'éléments:
nx1 = 2 ;
nx2 = 4 ;
nx3 = 5 ;
* Définitions des densites:
dbeg = 0.4 ;
dend = 0.5 ;
'FINSI' ;
*=====================
* Domaine de calcul 2D:
*=====================
'OPTION' 'DIME' 2 'ELEM' 'QUA8' ;
*****************************************
* *
* Procedure calcul *
* *
*****************************************
'DEBP' CALCUL ;
* Calcul de la convergence.
* Evolution de la vitesse en un point quelconque.
'ARGUMENT' RVX*'TABLE' ;
RV = RVX.'EQEX' ;
dd = RV.PASDETPS.'NUPASDT' ;
nn = dd/2 ;
'SI' ((dd-(2*nn))'EGA'0) ;
* Champs de vitesse aux instants n et n-1:
un = RV.INCO.'UN' ;
unm1 = RV.INCO.'UNM1' ;
* Extraction des composantes de la vitesse au pas de temps n et n-1:
unx = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' un) ;
unm1x = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' unm1) ;
uny = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' un) ;
unm1y = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' unm1) ;
* Sauvegarde de la vitesse en un point fixe du domaine à l'instant n:
z = 'REDU' unx RV.INCO.'PoiObs' ;
zz = 'PROG' ('EXTR' z 'SCAL' RV.INCO.'PoiObs') ;
RV.INCO.'Uvit'= RV.INCO.'Uvit' et zz ;
* Calcul de la convergence:
errx = 'KOPS' unx '-' unm1x ;
erry = 'KOPS' uny '-' unm1y ;
elix = 'MAXIMUM' errx 'ABS' ;
eliy = 'MAXIMUM' erry 'ABS' ;
elix = 'LOG'(elix + 1.0E-10)/('LOG' 10.) ;
eliy = 'LOG'(eliy + 1.0E-10)/('LOG' 10.) ;
erx = 'PROG' elix ;
ery = 'PROG' eliy ;
RV.INCO.'ERX' = RV.INCO.'ERX' et erx ;
RV.INCO.'ERY' = RV.INCO.'ERY' et ery ;
it = 'PROG' RV.PASDETPS.'NUPASDT' ;
temps = 'PROG' RV.PASDETPS.'TPS' ;
RV.INCO.'IT' = RV.INCO.'IT' et it ;
RV.INCO.'tps' = RV.INCO.'tps' et temps ;
RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET'
(RV.INCO.'UN') ;
'FINS' ;
as2 ama1 = 'KOPS' 'MATRIK' ;
'RESPRO' as2 ama1 ;
'FINP' ;
*============================
* Construction du maillage:
*============================
* Construction des points:
p0 = 0. 0. ;
p1 = -6. 0. ;
p2 = -0.5 0. ;
p3 = 0. 0.5 ;
p4 = 0.5 0. ;
p5 = 15. 0. ;
p6 = 15. 6. ;
p7 = 3. 6. ;
p8 = -3. 6. ;
p9 = -6. 6. ;
* construction des points de la partie supérieure:
p1p2 = p1 d nx1 p2 ;
p2p3 = 'CERCLE' (nx2/2) p2 p0 p3 ;
p3p4 = 'CERCLE' (nx2/2) p3 p0 p4 ;
p4p5 = p4 d nx3 p5 ;
p5p6 = p5 d dini dbeg dfin dend p6 ;
p6p7 = p6 d nx3 p7 ;
p7p8 = p7 d nx2 p8 ;
p8p9 = p8 d nx1 p9 ;
p9p1 = p9 d dini dend dfin dbeg p1 ;
* Création du contour:
peri1 = p1p2 et p2p3 et p3p4 et p4p5 ;
peri2 = p6p7 et p7p8 et p8p9 ;
s1 = 'DALL' peri1 p5p6 peri2 p9p1 'PLAN' ;
* construction des points de la partie inférieure:
p10 = -6. -6. ;
p11 = -3. -6. ;
p12 = 3. -6. ;
p13 = 15. -6. ;
p14 = 0. -0.5 ;
p1p10 = p1 d dini dbeg dfin dend p10 ;
p10p11 = p10 d nx1 p11 ;
p11p12 = p11 d nx2 p12 ;
p12p13 = p12 d nx3 p13 ;
p13p5 = p13 d dini dend dfin dbeg p5 ;
p5p4 = 'INVERSE' p4p5 ;
p4p14 = 'CERCLE' (nx2/2) p4 p0 p14 ;
p14p2 = 'CERCLE' (nx2/2) p14 p0 p2 ;
p2p1 = 'INVERSE' p1p2 ;
* Création du contour:
peri3 = p10p11 et p11p12 et p12p13 ;
peri4 = p5p4 et p4p14 et p14p2 et p2p1 ;
s2 = 'DALL' peri3 p13p5 peri4 p1p10 'PLAN' ;
* Création du domaine total:
'ELIMINER' 1.E-3 s1 s2 ;
domtot = s1 et s2 ;
*==============================================
* Changement des éléments du maillage en QUAF:
*==============================================
£domtot = 'CHANGER' domtot 'QUAF' ;
£cercle = 'CHANGER' (p2p3 et p3p4 et p4p14 et p14p2) 'QUAF' ;
£entree = 'CHANGER' (peri2 et p9p1 et p1p10 et peri3) 'QUAF' ;
'ELIM' 1.E-3 (£domtot et £cercle et £entree) ;
*=======================================
* Formulation du domaine Navier Stokes:
*=======================================
$domtot = 'MODE' £domtot 'NAVIER_STOKES' 'MACRO' ;
mdomtot = 'DOMA' $domtot maillage ;
*====================================
* Création des tables de résolution:
*====================================
maxit = 10000000 ;
rey = 20. ;
nu = 1./rey ;
arelax = 0.7 ;
RV = 'EQEX' $domtot 'ITMA' maxit 'ALFA' arelax
'ZONE' $domtot
'OPER' CALCUL ;
RV = 'EQEX' RV
'OPTI' 'EFM1' 'CENTREE' 'EXPL'
'ZONE' $domtot
'OPER' 'NS' 'INV_RE' 'INCO' 'UN' ;
RV = 'EQEX' RV
'OPTI' 'EF' 'CENTREE'
'ZONE' $domtot
'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN' ;
RVP = 'EQPR' $domtot
'ZONE' $domtot
'OPER' 'PRESSION' 0. ;
* 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. ;
* Création de la table des inconnues et initialisation:
RV.INCO = 'TABLE' INCO ;
RV.INCO.'UN' = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3) ;
RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3) ;
RV.INCO.'INV_RE'= nu ;
RV.'PRESSION'= RVP ;
RV.TFINAL = tfin ;
* Détermination d'un point d'observation:
obs = mdomtot 'POINT' 'PROC' (0.8 0.5) ;
RV.INCO.'PoiObs'= obs ;
* Initialisation des inconnues complémentaires:
RV.INCO.'Uvit'= 'PROG' 0. ;
RV.INCO.'ERX' = 'PROG' 0. ;
RV.INCO.'ERY' = 'PROG' 0. ;
RV.INCO.'IT' = 'PROG' 0 ;
RV.INCO.'tps' = 'PROG' 0. ;
RV.INCO.'entree'= £entree ;
RV.INCO.'cercle'= £cercle ;
*========================
* Résolution du problème:
*========================
EXEC RV ;
* Récupération de la vitesse au point PoiObs :
unx = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' RV.'INCO'.'UN') ;
z = 'REDU' unx RV.INCO.'PoiObs' ;
Uvit ='EXTR' z 'SCAL' RV.INCO.'PoiObs' ;
ecart = 'ABS' (Uvit - 1.0209) ;
'SI' ('NON' COMPLET) ;
'SI' (ecart > 1.E-4) ;
'ERREUR' 5 ;
'FINSI' ;
'FINSI' ;
* 'FIN' du jeu de données.
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer