* fichier : bc30.dgibi ************************************************************************ * Section : Fluides Thermique ************************************************************************ *$$$ 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) ; lcu = extr un 'COMP'; mu2=un lcu 'PSCA' un lcu; 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=SOMM( (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=SOMM( (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 ;