* fichier : bc30.dgibi
*$$$ BC30
*X BC30 (Jeux de donnees)
**
** BC30
**
** CANAL CHAUFFE INCLINE ANGL = 30. LONGUEUR H0=10 LARGEUR 1.
** K-E PLAN --- RE=2.E+5
**
** teste NSKE(SUPG) TSCAL(PSI) LAPN FPU ECHI et VNIMP
** Formulation EFM1 algorithme semi implicite
** teste aussi les méthodes d'inversion de la pression :
** Choleski et gradient conjugué.
**
**
** --- 26 JANVIER 1987 --- modifie --- 26 MAI 1993 ---
** modifie --- 17 DECEMBRE 1997 ---
**
** A noter :
**
** 1/ la manière dont sont donnés les différents coefficients
** du fait que df/dt a comme coefficient 1. pour toutes les
** équations.
** 2/ La manière dont est réalisé le couplage : construction
** du maillage de raccord (uniquement valable en 2D)
** et procédure CALCUL .
** 3/ Le bilan thermique dans ce cas est assez mal vérifié .
** Les causes : maillage grossier, flux thermique tres faible
** par rapport au debit enthalpique, terme d'accumulation
** d teta / dt non pris en compte dans le bilan.
**
**
GRAPH = 'N' ;
COMPLET = FAUX ;
NITMA=200;
SI (NON COMPLET) ;
err1=5.E-4;
NITMA = 40 ;
FINSI ;
option dime 2 elem tri6 ;
option dime 2 elem qua8 ;
p1=0 0 ;
p2=1. 0.;
nbpe=7 ; nbpp=10 ;
*nbpe=40 ; nbpp=80 ;
*nbpe=50 ; nbpp=150 ;
tole=1.e-5 ;
H0=10.;
entree= p1 d nbpe p2 ;
*entree= p1 d dini 0.05 dfin 0.2 p2 ;
c1= 5. 0. ;
q1=p1 tour c1 -90. ;
q2=p2 tour c1 -90. ;
pp1=p1 c c1 q1 nbpp ;
pp2=p2 c c1 q2 nbpp ;
q1=p1 plus (0 H0) ;
q2=p2 plus (0 H0) ;
pp1=p1 d q1 nbpp ;
eps=0.1 ;
pri=pp1 plus ((-1.*eps) 0) ;
pp2=p2 d q2 nbpp ;
sortie=entree tour c1 -90 ;
sortie=entree plus (0 H0) ;
pp1=inve pp1 ;
sortie=inve sortie;
elim tole (sortie et pp1 et pp2 );
cont=entree et pp2 et sortie et pp1 ;
*bell=surf cont ;
bell=daller entree pp2 sortie pp1 ;
paroi=pri trans 4 (-0.1 0) ;
pr= cote 3 paroi ;
entref=entree elem seg3 (lect 2 pas 1 6) ;
bell=bell tour c1 -20. ;
entree=entree tour c1 -20. ;
entref=entref tour c1 -20. ;
sortie=sortie tour c1 -20. ;
pp1=pp1 tour c1 -20. ;
pp2=pp2 tour c1 -20. ;
paroi=paroi tour c1 -20. ;
paroi=paroi coul rouge ;
pr= pr tour c1 -20. ;
pri= pri tour c1 -20. ;
MT=bell et paroi ;
pp1=inve pp1 ;
elim (pp1 et pp2 et entree et entref et sortie et bell) tole ;
elim ( pri et pr et paroi ) tole ;
titred= 'PLAQUES // COUDEES : TEST TRIO-EF NSKE FPU VNIMP' ;
titre= 'PLAQUES // COUDEES : TEST TRIO-EF NSKE FPU VNIMP' ;
Mmt = CHAN mt QUAF ;
Mbell = CHAN bell QUAF ;
Mparoi = CHAN paroi QUAF ;
Mentree= CHAN entree QUAF ;
Mentref= CHAN entref QUAF ;
Msortie= CHAN sortie QUAF ;
Mpri = CHAN pri QUAF ;
Mpp1 = CHAN pp1 QUAF ;
Mpp2 = CHAN pp2 QUAF ;
elim (Mmt et Mbell et Mparoi et Mentree et Mentref et Msortie
et Mpri et Mpp1 et Mpp2) tole ;
$mt = MODE MMT 'NAVIER_STOKES' 'MACRO' ;
$bell = MODE Mbell 'NAVIER_STOKES' 'MACRO' ;
$paroi= MODE Mparoi 'NAVIER_STOKES' 'MACRO' ;
paroi = doma $paroi 'MAILLAGE' ;
DOMA $mt 'IMPR' ; DOMA $bell 'IMPR' ; DOMA $paroi 'IMPR' ;
$pri = MODE Mpri 'NAVIER_STOKES' 'MACRO' ;
$pp1 = MODE Mpp1 'NAVIER_STOKES' 'MACRO' ;
$pp2 = MODE Mpp2 'NAVIER_STOKES' 'MACRO' ;
$entree = MODE Mentree 'NAVIER_STOKES' 'MACRO' ;
$entref = MODE Mentref 'NAVIER_STOKES' 'MACRO' ;
$sortie = MODE Msortie 'NAVIER_STOKES' 'MACRO' ;
pp1c=pp1 elem (lect 2 pas 1 (nbel pp1));
pric=pri elem (lect 2 pas 1 (nbel pri));
Mpp1c = CHAN pp1c QUAF ;
Mpric = CHAN pric QUAF ;
elim (Mpric et Mpp1c et Mmt ) tole ;
$pp1c=MODE Mpp1c 'NAVIER_STOKES' 'MACRO' ;
$pric=MODE Mpric 'NAVIER_STOKES' 'MACRO' ;
pp1l=doma $pp1c MAILLAGE ;
pril=doma $pric MAILLAGE ;
rfp=racc eps pp1l pril ;
rfpl= rfp chan ligne ;
rfp= diff rfpl (pp1l et pril) ;
rfpi= inve rfp ;
************************* PROCEDURE CALCUL *****************
DEBP CALCUL ;
ARGU RX*TABLE ;
rv=rx.'EQEX' ;
tn=rv.'INCO'.'TN' ;
un=rv.'INCO'.'UN' ;
TF=redu tn pp1c ;
rv.'INCO'.'TF'= kcht $pric scal sommet (tf kpro rfp) ;
TP=redu tn pric ;
rv.'INCO'.'TP'= kcht $pp1c scal sommet (tp kpro rfpi) ;
rv.'INCO'.'UF'= kcht $bell vect sommet (redu un bell) ;
as2 ama1 = 'KOPS' 'MATRIK' ;
finproc as2 ama1 ;
************************* PROCEDURE CALCUL *****************
*option donn 5 ;
NU =7.3E-7 ;
LAMBDAP=15.2 ;
ROCPP=3.7E6 ;
ALFP=LAMBDAP/ROCPP ;
*ALFP=4.06E-6 ;
LAMBDAE=0.62 ;
ROCPE=4.15E6 ;
ALFE=LAMBDAE/ROCPE ;
*ALFE=1.5E-7 ;
mess ' ALFE=' ALFE ' ALFP=' ALFP ;
YP =5.E-2 ;
R0 =0.5 ;
TPSC =1500.;
TREF=0. ;
ZERO =1.E-20 ;
gb = 0. 0. ;
H=1.E4 ;
HF=H/ROCPE ;
HP=H/ROCPP ;
HP=H/ROCPE ;
CMU = 0.09 ;
UREF = 4. ;
KE = (0.05 * UREF)*(0.05 * UREF) ;
EE = (KE**1.5) / R0 ;
nut = 0.09 * KE * KE / EE ;
mess ' NUT entree ' nut ;
UET1 = KCHT $PP1 SCAL CENTRE 5.E-2 ;
UET2 = KCHT $PP2 SCAL CENTRE 5.E-2 ;
NUT = KCHT $BELL SCAL CENTRE nut ;
RV= EQEX 'DUMP' 'ITMA' NITMA 'ALFA' 0.9 OPTI 'SUPG'
'ZONE' $MT 'OPER' 'CALCUL'
'ZONE' $BELL 'OPER' 'NSKE' NU NUT 'INCO' 'UN' 'KN' 'EN'
'ZONE' $PP1 'OPER' 'FPU' NU UET1 YP 'INCO' 'UN' 'KN' 'EN'
'ZONE' $PP2 'OPER' 'FPU' NU UET2 YP 'INCO' 'UN' 'KN' 'EN'
;
RV=EQEX RV OPTI EFM1 'PSI'
'ZONE' $BELL 'OPER' 'TSCA' ALFE 'UF' 0. NUT 0.7 'INCO' 'TN'
'ZONE' $PAROI 'OPER' 'LAPN' ALFP 'INCO' 'TN'
'ZONE' $PRIC 'OPER' 'ECHI' HP 'TF' 'INCO' 'TN'
'ZONE' $PP1C 'OPER' 'ECHI' HF 'TP' 'INCO' 'TN'
;
RV=EQEX RV OPTI EFM1 'CENTREE'
'ZONE' $BELL 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
'ZONE' $MT 'OPER' 'DFDT' 1. 'TN' 'DELTAT' 'INCO' 'TN'
'ZONE' $BELL 'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
'ZONE' $BELL 'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN'
;
RV=EQEX RV CLIM
'KN' TIMP ENTREE KE 'EN' TIMP ENTREE EE
'TN' TIMP PR 100. 'TN' TIMP ENTREE 20.
;
RV.INCO=TABLE 'INCO' ;
UF = kcht $BELL VECT SOMMET (UREF*(0.5 0.9));
RV.INCO.'UN'=kcht $BELL VECT SOMMET (0. 0.) UF ;
RV.INCO.'KN'=kcht $BELL SCAL SOMMET KE ;
RV.INCO.'EN'=kcht $BELL SCAL SOMMET EE ;
tp=kcht $paroi scal sommet 100. ;
RV.INCO.'TN'=kcht $MT SCAL SOMMET 20. tp ;
* Suivi temporel de quelques points
lh= elem (doma $bell sommet) 'POI1' (lect 1 pas 20 101) ;
lc1= elem (doma $pp1 sommet) 'POI1' (lect 4 pas 4 20) ;
lc2= elem (doma $pri sommet) 'POI1' (lect 4 pas 4 20) ;
lc=lc1 et lc2 ;
his = khis 'UN' 1 lh
'UN' 2 lc1
'KN' lh
'EN' lh
'TN' lc1 ;
rv.'HIST'=his ;
DEBPROC TEST ityp*entier ksto*entier rv*table ;
* 1 : Choleski (valeur par defaut)
**** Le reste est non maintenu et ne fonctionne pas sur machine 64bytes**
* 2 : gradient conjugue sans preconditionnement
* 3 : gradient conjugue avec preconditionnement diagonal
* 4 : gradient conjugue avec preconditionnement par Choleski incomplet
* 7 : gradient conjugue avec preconditionnement diagonal (variante)
*
*Pour 2,3,4, la matrice de pression est calculee et stockee une fois
*pour toutes. Pour 7, elle n'est jamais stockee explicitement et permet
*de ce fait d'atteindre les plus grosses tailles de problemes.
*
*ityp=1 ;
* obtenir ityp*entier ;
*
*Pour 2,3,4, on propose deux types de stockage : morse et compresse. Le
*second type n'est interessant que pour les machines vectorielles.
* ksto=0 morse ksto=1 compresse
* ksto=0 Valeur par defaut
*
*ksto=0 ;
*obtenir ksto*entier ;
RVP = EQPR $BELL KTYPI ITYP 'ZONE' $BELL 'OPER' 'PRESSION' 0.
'ZONE' $PP1 'OPER' 'VNIMP' 0. 'ZONE' $PP2 'OPER' 'VNIMP' 0.
'ZONE' $ENTREE 'OPER' 'VNIMP' UREF 'ZONE' $ENTREF 'OPER' 'VTIMP' 0.;
si (non (ega ityp 1));
rvp.METHODE.KSTOCK=ksto ;
rvp.METHODE.JTEPS=5.E-4;
rvp.METHODE.'IPAT'=1;
finsi ;
RV.'PRESSION' =RVP ;
EXEC RV ;
mess ' Fin PROC ' ;
FINPROC RV ;
DEBPROC POST RV*TABLE GRAPH*MOT ;
un=rv.'INCO'.'UN' ;
tn=rv.'INCO'.'TN' ;
qe=dbit un $entree ;
qs=dbit un $sortie ;
dq=(abs qs) - (abs qe) ;
dqr=dq / (abs qe) ;
mess ' Erreur sur la masse ' dqr ;
ER1= abs (abs(dqr) - 1.E-10) ;
srti=inve (doma $sortie MAILLAGE) ;
mu2=kops un 'PSCA' un ;
mu= mu2 ** 0.5 ;
evolV = EVOL 'CHPO' mu SCAL srti ;
evy=extr evolV 'ORDO' ;
list evy ;
ley=prog 3.4573 3.8455 4.0470 4.0806 4.0848 4.0855 4.0853 4.0852
4.0853 4.0855 4.0848 4.0806 4.0470 3.8455 3.4573 ;
ER2=resu( (abs (ley - evy)) / 4. ) ;
mess ' Ecart sur profil de vitesse en sortie ' er2 ;
evolT = EVOL 'CHPO' tn SCAL srti ;
evt=extr evolT 'ORDO' ;
list evt ;
let=prog 20.985 20.501 20.179 20.078 20.045 20.026 20.015 20.008
20.005 20.002 20.001 20.001 20.000 20.000 20.000 ;
ER3=resu( (abs (let - evt)) / 20. ) ;
mess ' Ecart sur profil de température en sortie ' er3 ;
D=doma $pp1c 'XXDIAGSI' ;
TF= kcht $pp1c scal sommet (redu tn pp1c) ;
TP=rv.inco.'TP' ;
DTFP=kops tp '-' tf ;
FP=somt ( hf* (kops d '*' (kops tp '-' tf) ) ) ;
mess 'FP=' fp ;
uteta= kops ('KCHT' $bell 'SCAL' 'SOMMET' tn) '*' un ;
ute=dbit uteta $entree ;
uts=dbit uteta $sortie ;
dut=(abs uts) - (abs ute) - fp ;
dutr=dut/FP ;
mess (' Bilan : erreur relative = ') dutr ;
ER4= abs (abs(dutr) - .58735) ;
si (non COMPLET ) ;
si ( ER1 > 1.E-9 ) ; erreur 5 ; finsi ;
si ( ER2 > err1 ) ; erreur 5 ; finsi ;
si ( ER3 > err1 ) ; erreur 5 ; finsi ;
si ( ER4 > err1 ) ; erreur 5 ; finsi ;
finsi ;
SI ('EGA' GRAPH 'O' );
TAB1=TABLE;
TAB1.'TITRE'=TABLE ;
TAB1 . 1 ='MARQ REGU ' ;
TAB1.'TITRE' . 1 = mot 'Module de U ' ;
DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
TAB1.'TITRE' . 1 = mot 'Temperature sortie ' ;
DESS evolT 'TITX' 'R (m)' 'TITY' 'T °C' LEGE TAB1 ;
ung1=vect un 0.1 ux uy jaune;
pn=elno (kcht $bell scal centre (rv.'PRESSION'.pression)) $bell ;
his=rv.'HIST' ;
dessin his.'TABD' his.'1UN' ;
dessin his.'TABD' his.'2UN' ;
dessin his.'TABD' his.'TN' ;
dessin his.'TABD' his.'KN' ;
dessin his.'TABD' his.'EN' ;
trace ung1 bell ;
trace pn bell ;
trace tn bell ;
trace tn paroi;
kn=rv.inco.'KN' ;
en=rv.inco.'EN' ;
nut=0.09*kn*kn/en ;
trace nut bell ;
trace kn bell ;
trace en bell ;
FINSI ;
FINPROC ;
******************** TESTS **************************
ityp=1 ;
ksto = 0 ;
rv = test ityp ksto rv ;
* ityp=4 ;non maintenu, ne marche pas sur machine 64bytes
ityp=1 ;
ksto = 0 ;
rv = test ityp ksto rv ;
* ityp=4 ;non maintenu, ne marche pas sur machine 64bytes
ityp=1 ;
ksto = 1 ;
rv = test ityp ksto rv ;
* ityp=3 ;non maintenu, ne marche pas sur machine 64bytes
ityp=1 ;
ksto = 1 ;
rv = test ityp ksto rv ;
* ityp=2 ;non maintenu, ne marche pas sur machine 64bytes
ityp=1 ;
ksto = 1 ;
rv = test ityp ksto rv ;
post rv graph ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer