* fichier : Comte-Bellot.dgibi ************************************************************************************ * * * Ecoulement entre deux plaques planes infinies (Symétrie) * * * * Penser a augmenter YP qd Rey diminue * * Rey = 4.E6 -> YP = 1.e-3 Rey = 4.e4 -> YP=1.E2 * * Nu = 1.e-7 Rey=4.e6 * ************************************************************************************ OPTI TRAC 'PSC' ; GRAPH = 'N' ; OPTION DIME 2 ELEM QUA8; COMPLET=VRAI; COMPLET=FAUX; sqrt2 = 2.**0.5 ; tole = 1e-7; ****************************************************************** * Paramètres géométriques du sous-canal ****************************************************************** * Demi-espacement entre les deux plaques h = 0.25; * Longueur de l'espace simulé (20 fois le diamètre hydraulique) Si COMPLET; L = 20; Sinon; L = 10; Finsi; ****************************************************************** * Paramètres du maillage ****************************************************************** * nombre de points sur la verticale nv0 = 15 ; * nombre de points sur la longueur nh0 = 35; nh0 = 70; ****************************************************************** * Paramètres lié à l'écoulement ****************************************************************** * Constante du modèle de turbulence * CMU = 0.09 ; * Viscosité cinématique du fluide NU = 1.e-7; * Diamètre hydraulique du sous-canal (DH = 1 car adim) DH = 4.*h ; * Vitesse initiale UI = 0.4626; * Constante cnu cnu = 0.09 ; ****************************************************************** * Paramètres liés à la résolution numérique ****************************************************************** KPRESS = 'CENTREP1' ; DISCR = 'MACRO'; * Nombre d'itérations * Pas de temps Si Complet; NITMA =800; DT = 0.1 ; Sinon; NITMA =10; DT = 2. ; Finsi ; ****************************************************************** * Construction du maillage ****************************************************************** * Points pour construire le maillage élémentaire p0 = 0. 0. ; p1 = L 0. ; p2 = L h ; p3 = 0. h ; * Construction des cotés élémentaires bsym = DROI nh0 p0 p1 ; sortie = DROI (-1*nv0) p1 p2 'DINI' 0.15 'DFIN' 0.1 ; bwall = DROI nh0 p2 p3 ; cot34 = DROI (-1*nv0) p0 p3 'DINI' 0.15 'DFIN' 0.1 ; entree = INVE cot34 ; dom = daller bsym sortie bwall entree ; *tracer dom; * Description du domaine de calcul et des frontières NN = 'NBEL' dom; 'MESSAGE' 'Nombre d éléments du maillage' NN; ******************************************************************** * Changement des éléments du maillage en DISCR: ******************************************************************** £dom = 'CHANGER' dom QUAF ; £entree = 'CHANGER' entree QUAF ; £sortie = 'CHANGER' sortie QUAF ; £bsym = 'CHANGER' bsym QUAF ; £bwall = 'CHANGER' bwall QUAF ; elim tole (£dom et £entree et £sortie et £bsym et £bwall) ; ************************************************************** * Formulation du domaine Navier Stokes: ************************************************************** $dom = 'MODE' £dom 'NAVIER_STOKES' DISCR ; $entree = 'MODE' £entree 'NAVIER_STOKES' DISCR ; $sortie = 'MODE' £sortie 'NAVIER_STOKES' DISCR ; $bsym = 'MODE' £bsym 'NAVIER_STOKES' DISCR ; $bwall = 'MODE' £bwall 'NAVIER_STOKES' DISCR ; DOMA $dom 'IMPR' ; entree = 'DOMA' $entree 'MAILLAGE' ; sortie = 'DOMA' $sortie 'MAILLAGE' ; bwall = 'DOMA' $bwall 'MAILLAGE' ; bsym = 'DOMA' $bsym 'MAILLAGE' ; ***************************************************************** * Création de la table des inconnues et initialisation: ***************************************************************** * Rey = 4.E6 -> YP = 1.e-3 Rey = 4.e4 -> YP=1.E2 YP = 1.e-2 ; YP = 1.e-3 ; KE = (0.02 * UI)*(0.02 * UI)/2.; EE = (KE**1.5)/h ; * Vitesse de frottement : * Viscosité turbulente : NUT = 0.09 * KE * KE / EE ; Rey=UI*DH/NU; 'MESSAGE' 'NUT =' NUT ' Rey=' Rey; RV = EQEX 'ITMA' NITMA 'OPTI' 'EF' 'IMPL' 'SUPG' 'ZONE' $bwall 'OPER' 'FPU' 1. 'UN' NU 'UET' YP 'INCO' 'UN' 'KN' 'EN' 'ZONE' $dom 'OPER' 'KEPSILON' 1. 'UN' NU 'DT' 'INCO' 'KN' 'EN' 'ZONE' $dom 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN' ; * Condition limite sur le mur *RV = EQEX RV 'CLIM' 'UN' 'UIMP' bwall 0. 'UN' 'VIMP' bwall 0.; * Condition limite de symétrie RV = EQEX RV 'CLIM' 'UN' 'VIMP' bsym 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' entree UI 'UN' 'VIMP' entree 0.; RV = EQEX RV 'CLIM' 'KN' 'TIMP' entree KE; RV = EQEX RV 'CLIM' 'EN' 'TIMP' entree EE; RV = EQEX RV 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $dom 'OPER' DFDT 1. 'UN' DT 'INCO' 'UN'; RVP = EQEX 'OPTI' 'EF' KPRESS 'ZONE' $dom 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES' 'ZONE' $bwall 'OPER' 'VNIMP' $dom 0. 'INCO' 'UN' 'PRES'; RV.'PROJ'= RVP ; RV.INCO= 'TABLE' INCO; RV.INCO.'UN' = 'KCHT' $dom 'VECT' 'SOMMET' (UI 0.); RV.INCO.'UNM1' = 'KCHT' $dom 'VECT' 'SOMMET' (0.D0 0.D0); RV.INCO.'KN' = KCHT $dom SCAL SOMMET KE; RV.INCO.'EN' = KCHT $dom SCAL SOMMET EE; RV.INCO.'PRES' = 'KCHT' $dom 'SCAL' KPRESS 0.; RV.INCO.'DT' = DT ; *RV.INCO.'UET' = KCHT $bwall 'SCAL' 'SOMMET' 2.E-2; rv.inco.'MUF' = 200.*NU; temp place; exec rv ; temp place; * RESULTATS un= rv.inco.'UN' ; kn = rv.inco.'KN' ; en = rv.inco.'EN' ; prodt = rv.inco.'PRODT' ; nut =(rv.inco.'MUF')*(1./NU) ; pres = ELNO $dom (exco 'PRES' rv.inco.'PRESSION') KPRESS; nun = (pscal un (mots 'UX' 'UY') un (mots 'UX' 'UY'))**0.5 ; USORTIE = 'EVOL' 'CHPO' un 'UX' sortie; KSORTIE = EVOL 'CHPO' RV.INCO.'KN' sortie; ESORTIE = EVOL 'CHPO' RV.INCO.'EN' sortie; NSORTIE = EVOL 'CHPO' nut sortie; PTSORTIE = EVOL 'CHPO' prodt sortie; YPLUSI = YP * (rv.inco.'UET') / NU; evyplus=evol chpo yplusi (inve bwall); YSORTIE = COOR 2 sortie; * Vitesse de frottement à la sortie du domaine (au point p2 plus précisément) : UFROT = 'EXTR' (rv.inco.'UET') 'SCAL' p2; * Coordonnées yplus de la sortie YPLUS = ((0.25+YP)-YSORTIE) * UFROT / NU; * Champ de vitesse normé UPLUS = USORTIE / UFROT; Ksi= rv.inco.'MUF' * prodt *(inve rv.inco.'EN'); *--------- Référence calcul fin ---------------------------------------- *----------------------------------------------------------------------- si faux; nx = extr nsortie 'ABSC' ; ny = extr nsortie 'ORDO' ; ky = extr ksortie 'ORDO' ; uy = extr usortie 'ORDO' ; xypl=extr evyplus 'ABSC' ; yypl=extr evyplus 'ORDO' ; list nx; list ny; list ky; list uy; list xypl; list yypl; finsi; xnuts=prog 0.0000 2.38164E-02 4.76327E-02 6.71384E-02 8.66440E-02 0.10262 0.11859 0.13168 0.14476 0.15548 0.16619 0.17497 0.18375 0.19093 0.19812 0.20401 0.20989 0.21472 0.21954 0.22348 0.22743 0.23067 0.23390 0.23655 0.23920 0.24137 0.24354 0.24531 0.24709 0.24854 0.25000 ; ynuts=prog 260.31 366.25 901.62 1383.4 1818.7 2121.0 2354.5 2485.6 2555.9 2565.4 2532.7 2464.9 2358.3 2239.4 2091.6 1947.3 1781.9 1629.6 1462.3 1313.2 1153.8 1015.0 869.28 744.33 615.34 506.32 395.64 303.88 212.84 139.97 63.065 ; yky=prog 1.08988E-05 2.03752E-05 6.73979E-05 1.24382E-04 1.88152E-04 2.44708E-04 3.01914E-04 3.49488E-04 3.96502E-04 4.35483E-04 4.75142E-04 5.07189E-04 5.38433E-04 5.63879E-04 5.88796E-04 6.09357E-04 6.29581E-04 6.46532E-04 6.63290E-04 6.77661E-04 6.92061E-04 7.04813E-04 7.17758E-04 7.29813E-04 7.42301E-04 7.55025E-04 7.68504E-04 7.84668E-04 8.00892E-04 8.24540E-04 7.63839E-04 ; yuy=prog 0.50623 0.50569 0.50187 0.49693 0.49080 0.48504 0.47862 0.47285 0.46655 0.46091 0.45477 0.44928 0.44327 0.43787 0.43190 0.42647 0.42043 0.41488 0.40862 0.40279 0.39613 0.38983 0.38249 0.37538 0.36686 0.35831 0.34760 0.33619 0.32065 0.30182 0.26896 ; xypl=prog 0.0000 0.14286 0.28571 0.42857 0.57143 0.71429 0.85714 1.0000 1.1429 1.2857 1.4286 1.5714 1.7143 1.8571 2.0000 2.1429 2.2857 2.4286 2.5714 2.7143 2.8571 3.0000 3.1429 3.2857 3.4286 3.5714 3.7143 3.8571 4.0000 4.1429 4.2857 4.4286 4.5714 4.7143 4.8571; xypl=xypl et (prog 5.0000 5.1429 5.2857 5.4286 5.5714 5.7143 5.8571 6.0000 6.1429 6.2857 6.4286 6.5714 6.7143 6.8571 7.0000 7.1429 7.2857 7.4286 7.5714 7.7143 7.8571 8.0000 8.1429 8.2857 8.4286 8.5714 8.7143 8.8571 9.0000 9.1429 9.2857 9.4286 9.5714 9.7143 9.8571); xypl=xypl et (prog 10.0000 10.143 10.286 10.429 10.571 10.714 10.857 11.000 11.143 11.286 11.429 11.571 11.714 11.857 12.000 12.143 12.286 12.429 12.571 12.714 12.857 13.000 13.143 13.286 13.429 13.571 13.714 13.857 14.000 14.143 14.286 14.429 14.571 14.714 14.857); xypl=xypl et (prog 15.000 15.143 15.286 15.429 15.571 15.714 15.857 16.000 16.143 16.286 16.429 16.571 16.714 16.857 17.000 17.143 17.286 17.429 17.571 17.714 17.857 18.000 18.143 18.286 18.429 18.571 18.714 18.857 19.000 19.143 19.286 19.429 19.571 19.714 19.857 20.000) ; yypl=prog 227.14 200.86 186.99 180.11 175.86 172.94 170.42 169.16 167.07 166.45 164.82 164.58 163.12 162.93 161.80 161.67 160.69 160.50 159.77 159.58 158.94 158.71 158.22 158.01 157.57 157.33 156.99 156.79 156.46 156.25 155.99 155.81 155.56 155.37 155.17; yypl=yypl et (prog 155.01 154.81 154.65 154.49 154.35 154.19 154.05 153.92 153.80 153.67 153.55 153.44 153.34 153.22 153.13 153.03 152.95 152.85 152.77 152.68 152.62 152.53 152.47 152.39 152.34 152.26 152.21 152.14 152.10 152.03 151.99 151.93 151.90 151.84 151.80); yypl=yypl et (prog 151.75 151.73 151.67 151.65 151.60 151.58 151.54 151.52 151.48 151.46 151.43 151.41 151.38 151.37 151.33 151.33 151.30 151.29 151.26 151.26 151.23 151.23 151.20 151.21 151.18 151.19 151.16 151.17 151.14 151.16 151.13 151.15 151.12 151.14 151.12); yypl=yypl et (prog 151.13 151.11 151.13 151.12 151.13 151.12 151.13 151.13 151.14 151.14 151.14 151.15 151.15 151.17 151.15 151.19 151.16 151.21 151.17 151.23 151.18 151.26 151.19 151.29 151.21 151.33 151.23 151.36 151.26 151.39 151.30 151.41 151.35 151.42 151.42 151.38) ; evyplusr=evol manu xypl yypl ; evnutsr=evol manu xnuts ynuts ; evuyr=evol manu xnuts yuy ; evkyr=evol manu xnuts yky ; *----------------------------------------------------------------------- ERR=0; TAB1=TABLE; TAB1.'TITRE'= table ; TAB1.1 =' MARQ CROI REGU'; TAB1.2 ='TIRR MARQ TRIA REGU'; TAB1.'TITRE' . 1 ='Référence'; TAB1.'TITRE' . 2 ='Calcul grossier'; m=(somme evyplusr 'ABSO') resu ; dyplus=(resu (somme (evyplus - evyplusr) 'ABSO'))/m; mess ' YPLUS : m=' m ' dyplus=' dyplus; Si (dyplus > 0.07);ERR=ERR + 1;Finsi; m=(somme evuyr 'ABSO') resu ; duyr=(resu (somme (evuyr - usortie) 'ABSO'))/m; mess ' Usortie : m=' m ' duyr=' duyr; Si (duyr > 0.07);ERR=ERR + 1;Finsi; m=(somme evnutsr 'ABSO') resu ; dnut=(resu (somme (evnutsr - nsortie) 'ABSO'))/m; mess ' Nu sortie : m=' m ' dnut=' dnut; Si (dnut > 0.3);ERR=ERR + 1;Finsi; m=(somme evkyr 'ABSO') resu ; dks=(resu (somme (evkyr - ksortie) 'ABSO'))/m; mess ' K sortie : m=' m ' dks=' dks; Si (dks > 0.3);ERR=ERR + 1;Finsi; Si (EGA GRAPH 'O') ; DESS USORTIE TITR ' USORTIE ' ; DESS evyplus TITR 'Y+'; DESS KSORTIE TITR ' KSORTIE ' ; DESS ESORTIE TITR ' ESORTIE ' ; DESS NSORTIE TITR ' NUTSORTIE / NU ' ; DESS PTSORTIE TITR ' ProdtSORTIE ' ; trace kn dom (cont dom) TITR 'KN' ; trace en dom (cont dom) TITR 'EN' ; ung = vect un 0.1 ux uy jaune ; trace ung dom; trace rv.inco.'TKTE' dom (cont dom) TITR 'TKTE'; trace nut dom (cont dom) TITR 'nut/nu ' ; trace Ksi dom (cont dom) TITR ' Ksi'; trace prodt dom (cont dom) TITR ' P'; dess (evyplusr et evyplus) TITR ' Y+' LEGE TAB1 'GRIL'; dess (evuyr et usortie) TITR 'Usortie' LEGE TAB1 'GRIL'; dess (evnutsr et nsortie) TITR 'Nut/nu sortie' LEGE TAB1 'GRIL'; dess (evkyr et ksortie) TITR 'Ksortie' LEGE TAB1 'GRIL'; FINSI; Si (ERR > 0); mess 'NOMBRE d ERREURS ' ERR; ERREUR 5; FINSI; FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer