Accueil | Cast3M

* 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 ;