* fichier : fsckei.dgibi
************************************************************************************
* *
* Maillage d'un sous-canal d'un faisceau de tube *
* à pas carré *
* *
************************************************************************************
OPTION DIME 2 ELEM QUA8;
COMPLET=FAUX;
*COMPLET=VRAI;
NBIT = 10 ;
GRAPH = FAUX ;
DT = 1. ;
Si COMPLET;
GRAPH = VRAI ;
NBIT = 1000;
DT = 0.2;
Finsi ;
sqrt2 = 2.**0.5 ;
tole = 1e-5;
******************************************************************
* Paramètres géométriques du sous-canal
******************************************************************
* Espacement entre deux centres de cylindres
b0 = 8.;
* Rayon d'un cylindre
r0 = 4. ;
* demi-diagonale "fluide"
diag1 = (b0*sqrt2)+((-1.)*r0) ;
******************************************************************
* Paramètres du maillage
******************************************************************
* distance entre le point p0 (centre du sous canal) et le point pint01
* c'est une portion (entre 0 et 1) de diag1
alpha01 = 0.7 ;
* distance entre le point p3 et pint23 (portion entre 0 et 1 de la distance entre p2 et p3)
alpha23 = 0.5 ;
* nombre de points sur le cylindre
n2 = 8 ;
* nombre de points sur la diagonale proche du cylindre
n12 = 10 ;
* nombre de points sur la diagonale loin du cylindre
n11 = 5 ;
******************************************************************
* Paramètres lié à l'écoulement
******************************************************************
* Constante du modèle de turbulence
CMU = 0.09 ;
* Viscosité cinématique du fluide
NU = 1.e-6;
* Vitesse d'entrée
UI = 1.e-6;
UI = 0.1 ;
* Diamètre hydraulique du sous-canal
DH = ((4*b0**2)+((-1.*PI)*(r0**2)))/(2*PI*r0) ;
******************************************************************
* Paramètres liés à la résolution numérique
******************************************************************
KPRESS = 'CENTREP1' ;
DISCR = 'MACRO' ;
* Nombre d'itérations
NITMA = NBIT;
* Pas de temps
* Distance à la paroi
YP = 1e-3;
******************************************************************
* PROCEDURE
******************************************************************
DEBPROC CALCUL;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
Tps = rv.'Tps' ;
DT = rv.'DT' ;
NUPAT = rv.'NUPAT' + 1;
Tps = Tps + DT ;
rv.'Tps' = Tps ;
rv.'NUPAT'=NUPAT ;
un = rv.'INCO'.'UN' ;
rv.ltps=rv.ltps et (prog tps) ;
rc.'pts1_1'=rc.'pts1_1'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_1')));
rc.'pts1_2'=rc.'pts1_2'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_2')));
rc.'pts1_3'=rc.'pts1_3'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_3')));
rc.'pts1_4'=rc.'pts1_4'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_4')));
rc.'pts1_5'=rc.'pts1_5'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_5')));
rc.'pts1_6'=rc.'pts1_6'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_6')));
rc.'pts1_7'=rc.'pts1_7'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_7')));
rc.'pts1_8'=rc.'pts1_8'
et (PROG (extr tn 'SCAL' (TABTEC.'pts1_8')));
rc.'pts1_1u'=rc.'pts1_1u'
et (PROG (extr un 'UX' (TABTEC.'pts1_1')));
rc.'pts1_2u'=rc.'pts1_2u'
et (PROG (extr un 'UX' (TABTEC.'pts1_2')));
rc.'pts1_3u'=rc.'pts1_3u'
et (PROG (extr un 'UX' (TABTEC.'pts1_3')));
rc.'pts1_4u'=rc.'pts1_4u'
et (PROG (extr un 'UX' (TABTEC.'pts1_4')));
rc.'pts2_1'=rc.'pts2_1'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_1')));
rc.'pts2_2'=rc.'pts2_2'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_2')));
rc.'pts2_3'=rc.'pts2_3'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_3')));
rc.'pts2_4'=rc.'pts2_4'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_4')));
rc.'pts2_5'=rc.'pts2_5'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_5')));
rc.'pts2_6'=rc.'pts2_6'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_6')));
rc.'pts2_7'=rc.'pts2_7'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_7')));
rc.'pts2_8'=rc.'pts2_8'
et (PROG (extr tn 'SCAL' (TABTEC.'pts2_8')));
rc.'pts3_1'=rc.'pts3_1'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_1')));
rc.'pts3_2'=rc.'pts3_2'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_2')));
rc.'pts3_3'=rc.'pts3_3'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_3')));
rc.'pts3_4'=rc.'pts3_4'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_4')));
rc.'pts3_5'=rc.'pts3_5'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_5')));
rc.'pts3_6'=rc.'pts3_6'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_6')));
rc.'pts3_7'=rc.'pts3_7'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_7')));
rc.'pts3_8'=rc.'pts3_8'
et (PROG (extr tn 'SCAL' (TABTEC.'pts3_8')));
rc.'pts3_1v'=rc.'pts3_1v'
et (PROG (extr un 'UY' (TABTEC.'pts3_1')));
rc.'pts3_2v'=rc.'pts3_2v'
et (PROG (extr un 'UY' (TABTEC.'pts3_2')));
rc.'pts3_3v'=rc.'pts3_3v'
et (PROG (extr un 'UY' (TABTEC.'pts3_3')));
rc.'pts3_4v'=rc.'pts3_4v'
et (PROG (extr un 'UY' (TABTEC.'pts3_4')));
* evs1 = (evol manu rv.ltps 'Temps (s)' rc.'pts3_1v' 'UY')
* et (evol manu rv.ltps 'Temps (s)' rc.'pts3_2v' 'UY')
* et (evol manu rv.ltps 'Temps (s)' rc.'pts3_3v' 'UY')
* et (evol manu rv.ltps 'Temps (s)' rc.'pts3_4v' 'UY');
* dess evs1 'NCLK' TITR 'Composante UY de la vitesse section 3 ' ;
Mess ' Pas de temps ' NUPAT ' Temps ' tps ;
mess 'Ufr=' Ufr 'Uc=' Uc ;
ai = modulo NUPAT 100;
si (EGA ai 0);
Mess ' *******************************************' ;
Mess ' Pas de temps ' NUPAT ' Temps ' tps ;
mess 'Ufr=' Ufr 'Uc=' Uc ;
Mess ' *******************************************' ;
oeil =1.e4 -1.e4 1.e4 ;
oeil2=-1.e4 -1.e4 1.e4 ;
THIS = FAUX ;
TECT = FAUX ;
POST THIS TECT ;
Finsi ;
aii = modulo NUPAT 1000;
si (EGA aii 0);
THIS = VRAI ;
TECT = VRAI ;
POST THIS TECT ;
THIS = FAUX ;
TECT = FAUX ;
Finsi ;
qc=dbit un $schaud ;
qf=(dbit un $sfroid)*(-1.) ;
qs=dbit un $sortie ;
mess 'Tps = ' Tps ' QC =' qc ' QF ' qf ' QS ' qs ;
*mess ' FIN CALCUL ' ;
as2 ama1 = 'KOPS' 'MATRIK' ;
finproc as2 ama1 ;
******************************************************************
* Construction du maillage
******************************************************************
* Centre du cylindre
c0 = (-1.*b0) b0;
* Points pour construire le maillage élémentaire
p0 = 0. 0. ;
p1 = (-1.*diag1*0.5*sqrt2) (diag1*0.5*sqrt2) ;
p2 = (-1.*b0) (b0+((-1.)*r0));
p3 = (-1.*b0) 0 ;
pint01 = (-1.*alpha01*diag1*0.5*sqrt2) (alpha01*diag1*0.5*sqrt2);
pint23 = (-1.*b0) ((-1.)*(1.-alpha01)*diag1 + b0 + (-1.*r0));
* Points pour créer les maillages symétriques
p4 = 0. b0 ;
p5 = b0 0. ;
p6 = (r0 + (-1.*b0)) b0;
* Construction des cotés élémentaires
cot11 = DROI n11 p0 pint01 ;
cot12 = DROI (-1*n12) pint01 p1 'DINI' 0.1 'DFIN' 0.1 ;
cot2 = CERC n2 p1 c0 p2 ;
cot31bis = DROI (-1*n12) pint23 p2 'DINI' 0.1 'DFIN' 0.1 ;
cot32 = DROI n11 pint23 p3;
cot31 = INVE cot31bis;
cot4 = DROI n2 p3 p0 ;
elim tole (cot11 et cot12 et cot2 et cot31 et cot32 et cot4) ;
* Reconstitution de 4 cotés pour utiliser la fonction DALLER
cot1 = cot11 et cot12;
cot3 = cot31 et cot32;
dom1 = daller cot1 cot2 cot3 cot4 ;
dom2 = dom1 SYME 'DROIT' c0 p0 ;
dom3 = (dom1 et dom2) SYME 'DROIT' p0 p4 ;
dom4 = (dom1 et dom2 et dom3) SYME 'DROIT' p3 p5 ;
* Description du domaine de calcul et des frontières
domtot = (dom1 et dom2 et dom3 et dom4) ;
ouest = cot3 et (cot3 SYME 'DROIT' p3 p0) ;
est = ouest SYME 'DROIT' p0 p4 ;
nord = ouest SYME 'DROIT' p0 c0 ;
sud = nord SYME 'DROIT' p3 p5 ;
bordNO = CERC (2*n2) p6 c0 p2 ;
bordNE = bordNO SYME 'DROIT' p0 p4 ;
bordSO = bordNO SYME 'DROIT' p0 p3 ;
bordSE = bordSO SYME 'DROIT' p0 p4 ;
elim tole (domtot et bordNO et bordNE et bordSO et bordSE
et ouest et est et nord et sud) ;
*tracer (domtot et ouest et est et nord et sud et bordNO
* et bordNE et bordSO et bordSE) ;
NN = 'NBEL' domtot;
'MESSAGE' 'Nombre d éléments du maillage' NN;
WALL0 = bordNO et bordNE et bordSO et bordSE;
********************************************************************
* Changement des éléments du maillage en DISCR:
********************************************************************
£domtot = 'CHANGER' domtot QUAF ;
£bordNO = 'CHANGER' bordNO QUAF ;
£bordNE = 'CHANGER' bordNE QUAF ;
£bordSO = 'CHANGER' bordSO QUAF ;
£bordSE = 'CHANGER' bordSE QUAF ;
£ouest = 'CHANGER' ouest QUAF ;
£est = 'CHANGER' est QUAF ;
£nord = 'CHANGER' nord QUAF ;
£sud = 'CHANGER' sud QUAF ;
elim tole (£domtot et £bordNO et £bordNE et £bordSO et £bordSE
et £ouest et £est et £nord et £sud) ;
**************************************************************
* Formulation du domaine Navier Stokes:
**************************************************************
$domtot = 'MODE' £domtot 'NAVIER_STOKES' DISCR ;
$bordNO = 'MODE' £bordNO 'NAVIER_STOKES' DISCR ;
$bordNE = 'MODE' £bordNE 'NAVIER_STOKES' DISCR ;
$bordSO = 'MODE' £bordSO 'NAVIER_STOKES' DISCR ;
$bordSE = 'MODE' £bordSE 'NAVIER_STOKES' DISCR ;
$ouest = 'MODE' £ouest 'NAVIER_STOKES' DISCR ;
$est = 'MODE' £est 'NAVIER_STOKES' DISCR ;
$nord = 'MODE' £nord 'NAVIER_STOKES' DISCR ;
$sud = 'MODE' £sud 'NAVIER_STOKES' DISCR ;
$WALL0 = $bordNO et $bordNE et $bordSE et $bordSO;
*****************************************************************
* Création de la table des inconnues et initialisation:
*****************************************************************
KE = (0.05 * UI)*(0.05 * UI)/2.;
EE = (KE**1.5)/r0 ;
mess 'UI=' ui 'KE=' KE 'EE=' EE 'Teta=' (KE/EE) ;
RV = EQEX $domtot 'ITMA' NITMA
'OPTI' 'EF' 'IMPL' 'CENTREE'
'ZONE' $WALL0 'OPER' 'FPU' 1. 'UN' NU 'UF' YP 'INCO' 'UN' 'KN' 'EN'
'ZONE' $domtot 'OPER' 'KEPSILON' 1. 'UN' NU 'DT' 'INCO' 'KN' 'EN'
'ZONE' $domtot 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN'
'ZONE' $domtot 'OPER' 'TSCA' 1. 'UN' 'MUF' 'INCO' 'TN'
;
* Condition limite sur les surfaces des cylindres
*RV = EQEX RV 'CLIM' 'UN' 'UIMP' bordNO 0. 'UN' 'VIMP' bordNO 0.;
*RV = EQEX RV 'CLIM' 'UN' 'UIMP' bordSO 0. 'UN' 'VIMP' bordSO 0.;
*RV = EQEX RV 'CLIM' 'UN' 'UIMP' bordNE 0. 'UN' 'VIMP' bordNE 0.;
*RV = EQEX RV 'CLIM' 'UN' 'UIMP' bordSE 0. 'UN' 'VIMP' bordSE 0.;
* Condition limite sur les frontières libres
* Remarque : pas besoin d'imposer une condition à la sortie car en utilisant la fonction 'CENTREP1'
* pour la pression, si aucune valeur n'est précisé, son intégrale sur cette frontière est maintenue nulle.
RV = EQEX RV 'CLIM' 'UN' 'UIMP' ouest UI;
RV = EQEX RV 'CLIM' 'UN' 'VIMP' ouest 0.;
RV = EQEX RV 'CLIM' 'TN' 'TIMP' ouest 1.;
RV = EQEX RV 'CLIM' 'TN' 'TIMP' nord 0.;
RV = EQEX RV 'CLIM' 'KN' 'TIMP' ouest KE;
RV = EQEX RV 'CLIM' 'EN' 'TIMP' ouest EE;
RV = EQEX RV 'CLIM' 'UN' 'VIMP' nord 0.;
RV = EQEX RV 'CLIM' 'UN' 'VIMP' sud 0.;
RV = EQEX RV 'OPTI' 'EFM1' 'CENTREE'
'ZONE' $domtot 'OPER' DFDT 1. 'UN' 'DT' 'INCO' 'UN';
*opti donn 5;
RVP = EQEX 'OPTI' 'EF' KPRESS
'ZONE' $domtot 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
'ZONE' $bordNO 'OPER' 'VNIMP' $domtot 0. 'INCO' 'UN' 'PRES'
'ZONE' $bordNE 'OPER' 'VNIMP' $domtot 0. 'INCO' 'UN' 'PRES'
'ZONE' $bordSO 'OPER' 'VNIMP' $domtot 0. 'INCO' 'UN' 'PRES'
'ZONE' $bordSE 'OPER' 'VNIMP' $domtot 0. 'INCO' 'UN' 'PRES'
;
RV.'PROJ'= RVP ;
RV.INCO= 'TABLE' INCO;
RV.INCO.'UN' = 'KCHT' $domtot 'VECT' 'SOMMET' (0.D0 0.D0);
RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' (0.D0 0.D0);
RV.INCO.'KN' = KCHT $domtot 'SCAL' 'SOMMET' KE;
RV.INCO.'TN' = KCHT $domtot 'SCAL' 'SOMMET' 0.;
RV.INCO.'EN' = KCHT $domtot 'SCAL' 'SOMMET' EE;
RV.INCO.'TT' = KCHT $domtot 'SCAL' 'SOMMET' (KE/EE);
RV.INCO.'PRES' = 'KCHT' $domtot 'SCAL' KPRESS 0.;
* Vitesse de frottement
RV.INCO.'UF' = 'KCHT' $WALL0 'SCAL' 'SOMMET' 5.E-2;
RV.INCO.'DT' = DT ;
RV.INCO.'MUF' = NU;
exec rv ;
un= rv.inco.'UN' ;
pres = exco rv.inco.'PRESSION' 'PRES';
pn = elno $domtot pres 'CENTREP1';
u2= pscal un (mots '1UN' '2UN') un (mots '1UN' '2UN');
Si GRAPH ;
trace pn domtot (cont domtot) TITR ' Pression' ;
pp= pn - (0.5*u2) ;
trace pp domtot (cont domtot) TITR ' Pression - u2' ;
nun = (pscal un (mots 'UX' 'UY') un (mots 'UX' 'UY'))**0.5 ;
trace nun domtot (cont domtot) TITR 'Norme de la vitesse' ;
*trace pres domtot;
ung = vect un 5. ux uy jaune ;
trace ung domtot;
nut = rv.inco.'MUF' - NU;
trace (nut*(1./NU)) domtot (cont domtot) TITR 'Nut/NU' ;
Finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer