* fichier : canalKLbr.dgibi * Ce cas teste le modèle K-L bas Reynolds (Wolfshtein Yap) * sur l'écoulement turbulent dans un canal plan. opti trace 'PSC'; *mess ' Configuration Complete ' ; COMPLET = VRAI ; *mess ' Configuration POST ' ; GRAPH = VRAI ; *----------------------------------------------------------------------- *----------------------------------------------------------------------- * |Triangles | Quadrangles * --------------------------- * LINE | | OK * LINB | | * MACRO| ? | OK * QUAF | | OK *----------------------------------------------------------------------- DISCR=QUAF; DISCR=MACRO; *DISCR=LINE; KPRES=CENTREP1; *KPRES=MSOMMET ; *KPRES=CENTRE ; KSUPG=SUPG; KSUPG=CENTREE; KSUPG=SUPG ; ******=================================================********** ******=================================================********** ******=================================================********** ******=================================================********** ******=================================================********** Si COMPLET; NBIT =200; DT = 0.04; H=1. ; L=30.*H; dcl=H/1200; dm=H/40; dcl=H/400; dm=H/40; nm=20; dcx=H/15.; dmx=H/10.; dcx=H/20.; dmx=H/30.; Sinon; NBIT =7; DT = 0.1; H=1. ; L=10.*H; dcl=H/20; dm=H/5; nm=15; dcx=H/10.; dmx=H/5.; Finsi; P0 = 0. 0. ; P1 = L 0. ; P2 = L H ; P3 = 0. H ; PS32= (P3 MOINS PS31); entree= P0 d dini dcl dfin dm ps31 d nm ps32 d dini dm dfin dcl P3; *opti trace 'X'; *trace entree ; *entree= P0 d dini dcl dfin dm ps3 d dini dm dfin dcl P3; *entree= P0 d 12 P3; par1 = p0 d dini dcx dfin dmx p1 ; *par1 = p0 d dini dcl dfin dcl p1 ; pc= 1. 1.; ang= 43.; ang= 60.; ang= 0. ; mt = daller entree par2 sortie par1; v1=1. 0. ; v2=0. -1.; v3=-1. 0.; v4=0. 1. ; list vt1; list vt2; list vt3; list vt4; $mt = model Mmt 'NAVIER_STOKES' DISCR; $entree = model Mentree 'NAVIER_STOKES' DISCR; $sortie = model Msortie 'NAVIER_STOKES' DISCR; $par1 = model Mpar1 'NAVIER_STOKES' DISCR; $par2 = model Mpar2 'NAVIER_STOKES' DISCR; doma $mt 'IMPR'; doma $entree 'IMPR'; lm = (lm1*ilm1)+(lm2*ilm2) + lma; dparoi=dparoi*(H - dparoi); dparoi=dparoi+1.e-10; Si GRAPH ; trace mt TITR ' Maillage '; Finsi ; U0 = 1. ; *U0 = 4. ; Ufr = U0 ; *Ufr = ye*(H - ye)*U0*1.5*4. ; uux=(cos ang) * U0; uuy=(sin ang) * U0; MU = 3.14E-5; Ro = 1.; YP = 1.e-2 ; Nu = Mu /Ro ; Rey= Ufr *H * 2. / NU; 'OPTI' 'EF' 'IMPL' KSUPG KPRES ZONE $mt 'OPER' 'KEPSILON' 1. 'UN' 'NU' DT 'INCO' 'KN' 'EN' ZONE $mt 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN' 'OPTI' 'EFM1' 'IMPL' 'CENTREE' ; *rv.'ALGO_KEPSILON'=mots 'KL'; RV= EQEX RV CLIM * UN UIMP entre1 Ufr UN VIMP entree 1.e-4 UN UIMP entre1 uux UN VIMP entree (uuy + 1.e-4) UN UIMP (par1 et par2) 0. UN VIMP (par1 et par2) 0. KN TIMP entre1 ((0.02*Ufr)**2.) KN TIMP (par1 et par2) 0. EN TIMP (par1 et par2) 0. ; rv.'METHINV'.TYPINV=3 ; rv.'METHINV'.IMPINV=0 ; rv.'METHINV'.NITMAX=700; rv.'METHINV'.PRECOND=3 ; *rv.'METHINV'.PRECOND=2 ; rv.'METHINV'.ILUTLFIL=4 ; rv.'METHINV'.ILUTDTOL=0.; rv.'METHINV'.RESID =1.e-8; rv.'METHINV' . 'FCPRECT'=1 ; rv.'METHINV' . 'FCPRECI'=1 ; RVP= EQEX 'OPTI' 'EF' KPRES ; rvp.'METHINV'.TYPINV=1 ; rvp.'METHINV'.IMPINV=0 ; rvp.'METHINV'.NITMAX=5000; rvp.'METHINV'.PRECOND=3 ; rvp.'METHINV'.ILUTLFIL=10; rvp.'METHINV'.ILUTDTOL=0.; rvp.'METHINV'.RESID =1.e-10 ; rvp.'METHINV' . 'FCPRECT'=50000 ; rvp.'METHINV' . 'FCPRECI'=50000 ; rv.'INCO'= table inco ; rv.'INCO'.'NU'=NU; rv.'INCO'.'MU'=MU; rv.'Tps' = 1. ; rv.'DT' = DT; rv.'NUPAT' =1; Si COMPLET ; exec rv ; sauv ; exec rv ; sauv ; exec rv ; sauv ; exec rv ; sauv ; exec rv ; sauv ; sinon ; exec rv ; * opti rest 'QPC4/canalKLbr2_5'; * rest Finsi ; un = rv.inco.'UN' ; um=(1./H)*qs; unn=(1./um)*un; *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * POST-TRAITEMENT *----------------------------------------------------------------------- *--- Edition d'un profil de vitesse asymptotique de Référence ---------- *--- Modèle de Buleev Debproc ANABU H*'FLOTTANT' Rey*'FLOTTANT' UM*'FLOTTANT' ; H2=H*H; p0=0. 0. ; p1=0. H. ; Pm=0. (H/2.); dpar=H/Rey*0.05; dctr=H*0.05; de= p0 d dini dpar dfin dctr pm d dini dctr dfin dpar p1 ; Dh=2.*H; Rho=1.; mu=um/Rey*Rho*Dh; m=7.1; nu=mu/rho; uet=um; REPETER CCP 15 ; fp=(uet/Um)**2.*8. ; F=fp/2/Dh*rho*(um*um); xi = z*(1./H); xiet = (1. - ((1. - (nu*m/h/Uet*4))**0.5))/2.; u1=H2*F/2./mu*xi*(1. - xi); ax=xi*(1. - xi)/xiet/(1. - xiet) + 1.e-10; ax=m*m*F*H/Rho* (log ax); ax =(abs ax) + 1.e-10; ax0=(h*h*F/mu/2*xiet*(1. - xiet))**2.; u2=(ax0+ax)**0.5; iu1=iua+iub; uref = (iu1 * u1) + ((1. - iu1)*u2) ; w=0.9; uet=(w*(uet*um/umn)) + ((1. - w)*uet); dumn=abs ((umn - um)/um); ' dumn=' dumn; *dess evuref ; Si (dumn < 1.e-4) ; quitter CCP; Finsi ; Fin CCP; * Calcul des u+g et u+d evuplus=evuref*(1./uet); Resproc evuplug evuplud evuref uet; Finproc ; *--- Profil de Kader --------------------------------------------------- pas 500. 1.e4 ; ysr=yplusr/Ret; Gamma=((Pr*yplusr)**4.)/(unp+(5.*(Pr**3.)*yplusr)); Gamma=Gamma*(-0.01); tplusr=Pr*yplusr*(exp Gamma); beta=((3.85*(Pr**(1./3.)) - 1.3)**2.); beta=beta+(2.12*(log Pr)); aa1=1.5*(2.*unp - ysr); aa2=2.5*(2.*unp - ysr); bb1=(unp - ysr)*(unp - ysr)*2. + unp; bb2=(unp - ysr)*(unp - ysr)*4. + unp; aa1=aa1/bb1 ; aa2=aa2/bb2 ; cc1= ((2.12*(Log ((unp+yplusr)*aa1))) + (beta*unp)) *(exp(unp/Gamma)); cc2= ((2.12*(Log ((unp+yplusr)*aa2))) + (beta*unp)) *(exp(unp/Gamma)); tplusr1=tplusr + cc1; tplusr2=tplusr + cc2; *--- Profil de Kader --------------------------------------------------- FINPROC evt1 evt2; DEBPROC POST ; tabk = table ; evupd1 evupg1 evuref1 uet1 = anabu 1. 2000. 1. ; evupd2 evupg2 evuref2 uet2 = anabu 1. 1.e4 1. ; evupd3 evupg3 evuref3 uet3 = anabu 1. 1.e5 1. ; evupd4 evupg4 evuref4 uet4 = anabu 1. 1.e6 1. ; evupd5 evupg5 evuref5 uet5 = anabu 1. 1.e7 1. ; *--- Profil de Reichardt ----------------------------------------------- Ka = 0.41; pas 500. 1.e4 pas 1.e4 2.e5 ; E=9.5; uplusr=(1./Ka)*(Log (unp + (Ka*yplusr))); uplusr=uplusr + (7.8* (unp - (exp (yplusr*(-1./11)))) - ( (yplusr*(-1./11)) * (exp (yplusr*(-1./3.))) ) ); *--- Profil de Reichardt ----------------------------------------------- *--- Profil de Référence pour le KLbr ---------------------------------- xrr=prog 2.99401E-06 1.24325E-03 4.78060E-03 1.05194E-02 1.83637E-02 2.82181E-02 3.99867E-02 5.35741E-02 6.88838E-02 8.58208E-02 0.10429 0.12419 0.14544 0.16793 0.19156 0.21625 0.24190 0.26841 0.29568 0.32363 0.35214 0.38114 0.41052 0.44019 0.47005 0.50000 0.52995 0.55981 0.58948 0.61886 0.64786 0.67637 0.70432 0.73159 0.75810 ; 0.78375 0.80844 0.83207 0.85456 0.87581 0.89571 0.91418 0.93112 0.94643 0.96001 0.97178 0.98163 0.98948 0.99522 0.99875 1.00000); urr = prog 0.0000 9.36980E-02 0.35215 0.56346 0.68057 0.75493 0.80827 0.84978 0.88391 0.91302 0.93845 0.96103 0.98127 0.99955 1.0161 1.0311 1.0447 1.0569 1.0678 1.0774 1.0858 1.0927 1.0983 1.1024 1.1049 1.1058 1.1051 1.1027 1.0988 1.0934 1.0865 1.0783 1.0688 1.0580 1.0459 ; 1.0325 1.0177 1.0013 0.98319 0.96315 0.94083 0.91571 0.88702 0.85347 0.81285 0.76100 0.68940 0.57826 0.37914 0.11368 5.00378E-04 ); krr = prog 1.14242E-08 5.70695E-04 5.85036E-03 9.18812E-03 9.24352E-03 8.78015E-03 8.29986E-03 7.86419E-03 7.46515E-03 7.08857E-03 6.72369E-03 6.36405E-03 6.00640E-03 5.64984E-03 5.29542E-03 4.94455E-03 4.60029E-03 4.26604E-03 3.94641E-03 3.64662E-03 3.37300E-03 3.13278E-03 2.93380E-03 2.78421E-03 2.69096E-03 2.65941E-03 2.69096E-03 2.78421E-03 2.93380E-03 3.13278E-03 3.37300E-03 3.64662E-03 3.94641E-03 4.26604E-03 4.60029E-03 ; 4.94455E-03 5.29542E-03 5.64984E-03 6.00640E-03 6.36405E-03 6.72369E-03 7.08857E-03 7.46515E-03 7.86419E-03 8.29986E-03 8.78015E-03 9.24352E-03 9.18812E-03 5.85036E-03 5.70695E-04 1.14242E-08); mrr = prog 2.0000 2.0031 2.4276 4.7639 9.0737 15.067 22.284 30.207 38.403 46.558 54.454 61.929 68.847 75.089 80.550 85.148 88.826 91.561 93.372 94.328 94.549 94.214 93.548 92.808 92.241 92.030 92.241 92.808 93.548 94.214 94.549 94.328 93.372 91.561 88.826; 85.148 80.550 75.089 68.847 61.929 54.454 46.558 38.403 30.207 22.284 15.067 9.0737 4.7639 2.4276 2.0031 2.0000); *--- Profil de Référence pour le KLbr ---------------------------------- tabk.'Rey'=Rey; un = rv.inco.'UN' ; um=(1./H)*qs; unn=(1./um)*un; TAB1 = TABLE ; tab1.'TITRE' = TABLE ; 'GRIL' 'CARR' lege TAB1 titr tit1; *----------------------------------------------------------------------- uet=(mu*(abs (dudy + 1.e-5)))**0.5; tity 'm/s' titx 'm'; f1=8.*(((1./um)*uet1)**2.); f2=8.*(((1./um)*uet2)**2.); fiso = 0.202*(Rey**(-0.196)); TAB1 = TABLE ; tab1.'TITRE' = TABLE ; TITR ' Coefficient de frottement' 'GRIL' lege TAB1; tabk.'evfiso'=evfiso; tabk.'evf2' =evf2 ; tabk.'evf1' =evf1 ; tabk.'fmax' =fmax ; yplus1= yp1*uet1*(rv.inco.'MU'**(-1.)); yplus2= yp2*uet2*(rv.inco.'MU'**(-1.)); tity ' ' titx 'm'; xypls1=xypls1 enlev 1; yypls1=yypls1 enlev 1; xypls2=xypls2 enlev 1; yypls2=yypls2 enlev 1; TAB1 = TABLE ; tab1.'TITRE' = TABLE ; TITR ' U+ en sortie' 'GRIL' lege TAB1 tity ' ' titx 'm'; tabk.'evr1' =evr1 ; tabk.'evupl1'=evupl1; tabk.'evupl2'=evupl2; nutsnu=(rv.inco.'MUF')*(1./MU); npl1=npl1 enlev 1 ; npl2=npl2 enlev 1 ; TAB1 = TABLE ; tab1.'TITRE' = TABLE ; TITR ' Nut/nu en sortie' 'GRIL' lege TAB1 tity ' ' titx 'm'; tabk.'evnpl1'=evnpl1; tabk.'evnpl2'=evnpl2; trace ung mt 'TITR' ' Vitesses '; *----------------------------------------------------------------------- tity 'm/s' titx 'm'; *----------------------------------------------------------------------- *----------------------------------------------------------------------- tity ' ' titx 'm'; tabk.'evmuax'=evmuax; *----------------------------------------------------------------------- *$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ *profils de référence evupdr evupgr evusr uetr = anabu 1. Rey 1. ; evdu = (evusi - evusr)*(evusi - evusr); tab1 = table; tab1.'TITRE' = TABLE ; tity 'm/s' titx 'm'; tity 'm2/s2' titx 'm'; tity 'm2/s' titx 'm'; tabk.'evus'=evus; tabk.'evusr'=evusr; errn = (abs (du - 4.06931E-03)); Si (NON COMPLET); Si(errn > 0.0001); erreur 5 ; Finsi ; Finsi; FINPROC tabk; ******* Profils de T+ de Kader *********************************** Si FAUX; evt2p1 evt1p1=kader 0.7 1.e4 'VERTE'; evt2p2 evt1p2=kader 0.7 1.e5 'ROUGE'; evt2p3 evt1p3=kader 0.7 1.e6 'BLEUE'; tab1 = table; tab1.'TITRE' = TABLE ; LOGX XBOR 0.1 1.E5 tity 'T+' titx 'Y+' 'GRIL' LEGE TAB1 evt2p1 evt1p1=kader 0.1 1.e4 'VERTE'; evt2p2 evt1p2=kader 1. 1.e5 'ROUGE'; evt2p3 evt1p3=kader 10. 1.e6 'BLEUE'; LOGX XBOR 0.1 1.E5 tity 'T+' titx 'Y+' 'GRIL' LEGE TAB1 Finsi; ******* Profils de T+ de Kader FIN ****************************** *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *profils de référence evupdr evupgr evusr uetr = anabu 1. Rey 1. ; evdu = (evusi - evusr)*(evusi - evusr); Si GRAPH ; t1=POST; tab1 = table; tab1.'TITRE' = TABLE ; tity 'm/s' titx 'm'; Finsi; errn = (abs (du - 4.54462E-03)); Si (NON COMPLET); Si(errn > 0.0001); erreur 5 ; Finsi ; Finsi; FIN;
© Cast3M 2003 - Tous droits réservés.
Mentions légales