* fichier : fsckei.dgibi ************************************************************************ * Section : Fluides Transitoire * Section : Fluides Transitoire ************************************************************************ ************************************************************************************ * * * 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 ;