* fichier : dvisp.dgibi * CAVITE CARREE ­ VAHL­DAVIS Methode de projection implicite * * Malvina Renesson Aout 1999 * *************************************************************** GRAPH=VRAI ; GRAPH=FAUX ; COMPLET=VRAI; COMPLET=FAUX ; TCPT = FAUX; TKPR = VRAI; TRESOU = VRAI; IMPARA = FAUX; Si (NON Complet); itmax= 80 ; * itmax= 2 ; *itmax= 20 ; *itmax= 8 ; DISCR='LINE' ; DISCR='MACRO'; *DISCR='QUAF'; d1 = 0.04; d2 = 0.12; d1 = 0.06; d2 = 0.18; *d1 = 0.01; d2 = 0.03; ** d1 = 0.2 ; d2 = 0.3 ; d1 = 0.01 ; d2 = 0.02 ; *d1 = 0.005; d2 = 0.01 ; DT=.25; Ra = 1.e6 ; Sinon ; *Ra = 1.e5 ; *itmax= 600 ; Ra = 1.e6 ; itmax= 600 ; DISCR='LINE' ; d1 = 0.01; d2 = 0.07; DT=.25; Finsi ; KSUPG ='CENTREE' ; KPRESS='MSOMMET' ; *KPRESS='CENTREP1' ; *KPRESS='CENTRE' ; ****************************** *ESTIMATION DE LA CONVERGENCE* ****************************** DEBPROC CALCUL; ARGU RVX*'TABLE' MDZ/'MMODEL' MDZ/'TABLE'; RV = RVX.'EQEX'; DD = RV.PASDETPS.'NUPASDT'; NN = DD/5; LO = (DD-(5*NN)) EGA 0; SI (LO); UN = RV.INCO.'UN'; UNM1 = RV.INCO.'UNM1'; unx = kcht $mt scal sommet (exco 'UX' un); unm1x = kcht $mt scal sommet (exco 'UX' unm1); uny = kcht $mt scal sommet (exco 'UY' un); unm1y = kcht $mt scal sommet (exco 'UY' unm1); ERRX = KOPS unx '-' unm1x; ERRY = KOPS uny '-' unm1y; ELIX = MAXI ERRX 'ABS'; ELIY = MAXI ERRY 'ABS'; ELIX = (LOG (ELIX + 1.0E-10))/(LOG 10.); ELIY = (LOG (ELIY + 1.0E-10))/(LOG 10.); MESSAGE 'ITER' RV.PASDETPS.'NUPASDT' 'ERREUR LINF' ELIX ELIY; RV.INCO.'UNM1' = KCHT $MT 'VECT' 'SOMMET' (RV.INCO.'UN'); IT = PROG RV.PASDETPS.'NUPASDT'; ER = PROG ELIY; RV.INCO.'IT' = (RV.INCO.'IT') ET IT; RV.INCO.'ER' = (RV.INCO.'ER') ET ER; FINSI; as2 ama1 = 'KOPS' 'MATRIK' ; 'FINPROC' as2 ama1 ; ************ * MAILLAGE * ************ OPTI DIME 2; OPTI ELEM QUA8; crit = 0.0000001; xxll = 1.; p1 = 0. 0.; p15 = (xxll/2.) 0.; p2 = xxll 0.; p25 = xxll 0.5; p3 = xxll 1.; p35 = (xxll/2.) 1.; p4 = 0. 1.; p45 = 0. 0.5; bas1 = p1 d dini d1 dfin d2 p15; bas2 = p15 d dini d2 dfin d1 p2; bas = (bas1 et bas2);elim crit bas; cdro = p2 d dini d1 dfin d2 p25 d dini d2 dfin d1 p3; haut = p3 d dini d1 dfin d2 p35 d dini d2 dfin d1 p4; cgau = p4 d dini d1 dfin d2 p45 d dini d2 dfin d1 p1; coupx = bas plus (0. 0.5) ; coupy=cgau plus (0.5 0.) ; cnt = bas et cdro et haut et cgau; mt = bas cdro haut cgau daller; bord = cont mt; *mt = surf bord; elim (mt et coupx) 1.e-5; Si GRAPH ; trace mt ; Finsi ; *option donn 5; orienter mt; ******** * DOMA * ******** Mmt = chan QUAF mt ; $mt = MODE Mmt 'NAVIER_STOKES' DISCR ; mt = doma $mt maillage; Mcgau=chan cgau QUAF ; $cgau=mode Mcgau 'NAVIER_STOKES' DISCR ; Mcdro=chan cdro QUAF ; $cdro=mode Mcdro 'NAVIER_STOKES' DISCR ; Mcnt=chan cnt QUAF ; $cnt=mode Mcnt 'NAVIER_STOKES' DISCR ; Mcoupx = chan coupx 'QUAF' ; $coupx=mode Mcoupx 'NAVIER_STOKES' DISCR ; Mcoupy = chan coupy 'QUAF' ; $coupy=mode Mcoupy 'NAVIER_STOKES' DISCR ; Elim (Mmt et Mcgau et Mcdro et Mcnt et Mcoupx et Mcoupy) 1.e-5 ; cgau=doma $cgau maillage ; cdro=doma $cdro maillage ; cnt=doma $cnt maillage ; coupx=doma $coupx maillage ; coupy=doma $coupy maillage ; *option donn 5; ************** * PARAMETRES * ************** Pr = 0.71; Gr = Ra/Pr; NU = 1/(Gr**0.5); ALF= NU/Pr; gb = kcht $mt vect centre (0. -1.); uref = 0.2 ; *********************** * CREATION DES TABLES * *********************** RV = EQEX 'OMEGA' 1. 'NITER' 1 ITMA itmax 'FIDT' 1 'OPTI' 'EF' 'IMPL' KSUPG 'ZONE' $MT 'OPER' CALCUL 'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'GB' 'TN' 0.5 'INCO' 'UN' 'ZONE' $MT 'OPER' 'TSCAL' 1. 'UN' ALF 0. 'INCO' 'TN' 'OPTI' 'CENTREE' ZONE $mt OPER DFDT 1. 'UN' DT INCO 'UN' ZONE $mt OPER DFDT 1. 'TN' DT INCO 'TN' ; **** RV = EQEX RV 'CLIM' 'UN' 'UIMP' cnt 0. 'UN' 'VIMP' cnt 0. 'TN' 'TIMP' cgau 1. 'TN' 'TIMP' cdro 0.; rv.'METHINV'.TYPINV=1 ; rv.'METHINV'.IMPINV=0 ; rv.'METHINV'.NITMAX=400; rv.'METHINV'.PRECOND=3 ; rv.'METHINV'.RESID =1.e-10; rv. 'METHINV' . 'FCPRECT'=1 ; rv. 'METHINV' . 'FCPRECI'=1 ; **** * La cavité est fermée il faut imposer la pression en 1 point ! prep1=doma $mt kpress; bcp=elem prep1 POI1 (lect 1 ) ; RVP = EQEX 'OPTI' 'EF' KPRESS ZONE $mt OPER KBBT -1. INCO 'UN' 'PRES' CLIM PRES TIMP bcp 0. ; rvp.'METHINV'.TYPINV=1 ; rvp.'METHINV'.IMPINV=0 ; rvp.'METHINV'.NITMAX=300; rvp.'METHINV'.PRECOND=3 ; rvp.'METHINV'.RESID =1.e-8 ; rvp.'METHINV' . 'FCPRECT'=100 ; rvp.'METHINV' . 'FCPRECI'=100 ; RV.INCO = TABLE INCO; RV.INCO.'UN' = kcht $MT VECT SOMMET (0. 0.); mt=doma $mt 'MAILLAGE'; cx=coor 1 mt; cy=coor 2 mt; un1=(nomc cx 'UY') - (nomc cy 'UX'); RV.INCO.'UN1'= un1; ung= vect un1 0.5 ux uy jaune ; *trace ung mt; *RV.INCO.'TN' = kcht $MT SCAL SOMMET (1.-corsx); RV.INCO.'TN' = kcht $MT SCAL SOMMET 0.5 ; RV.INCO.'PRES' = kcht $MT SCAL KPRESS 0. ; RV.INCO.'UNM1' = kcht $MT VECT SOMMET (1.E-5 1.E-5); RV.INCO.'IT' = PROG 1; RV.INCO.'ER' = PROG 0.; RV.INCO.'GB' = gb ; RV.'PROJ' = RVP; RV.'TCPT' = TCPT ; RV.'TKPR' = TKPR ; RV.'TRESOU'= TRESOU; RV.'IMPARA'= IMPARA; rv.'ITMA'=1; temps; REPETER BLOC 8 ; EXEC RV ; FIN BLOC; temps; *opti donn 5; ************* * RESULTATS * ************* *option sauv 'SAV'; *sauv Pr Ra Gr NU ALF gb uref RV; ************* OPTI ISOV SULI; un=(RV.INCO.'UN') * (Pr*(Gr**0.5)) ; unch = vect un (0.3/(Pr*(Gr**0.5))) UX UY JAUNE; u1y = exco un 'UX' ; u2x = exco un 'UY' ; evu=evol 'CHPO' u2x coupx ; evv=evol 'CHPO' u1y coupy ; EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG---E---inf' (RV.INCO.'ER'); GRADT = KOPS RV.INCO.'TN' 'GRAD' $MT; DTDX = KCHT $MT 'SCAL' 'CENTRE' (EXCO 'UX' GRADT); DTDX = (-1)*(ELNO $MT DTDX); EVOLd = EVOL 'CHPO' DTDX SCAL cdro ; EVOLg = EVOL 'CHPO' DTDX SCAL cgau ; evol1=evold et evolg ; LISTE EVOL1; TAB1 = TABLE; TAB1.1 = 'MARQ LOSA'; TAB1.2 = 'MARQ CROI'; Si GRAPH ; DESS EVOL1 'XBOR' 0. 1. 'GRIL' TAB1 'TITRE' ' Nusselt' ; trac RV.INCO.'TN' mt (cont mt) 14 'TITRE' ' Isothermes' ; trace unch mt (cont mt) 'TITRE' ' Vitesses' ; dess evu 'TITRE' 'Coupe ox à y=1/2 (uy)' ; dess evv 'TITRE' 'Coupe oy à x=1/2 (ux)' ; dess EVOL4 'XBOR' 0. itmax 'YBOR' -10.0 0.0; Finsi; Nusd = (redu DTDX cdro); Nusg = (redu DTDX cgau); Numg = somt (Nusg*(doma $cgau 'XXDIAGSI')) ; Numd = somt (Nusd*(doma $cdro 'XXDIAGSI')) ; un2=KOPS un '*' (Pr*(Gr**0.5)); sw = KOPS un 'ROT' $mt ; rk = EQEX $mt 'OPTI' 'EF' 'IMPL' ZONE $mt OPER LAPN -1. INCO 'PSI' *ZONE $mt OPER KONV 1. 'UN' 1. INCO 'PSI' ZONE $mt OPER FIMP sw INCO 'PSI' 'CLIM' 'PSI' 'TIMP' cnt 0.; rk.'INCO'=table 'INCO' ; rk.'INCO'.'PSI'=kcht $mt scal sommet 0. ; rk.'METHINV'.TYPINV=1 ; RK.'TCPT'=TCPT; EXec rk ; psi=rk.'INCO'.'PSI'; psi2=kops psi '*' (Pr*(Gr**0.5)); Si GRAPH ; trac psi mt CNT 14 'TITRE' 'Fonction de courant '; Finsi; opti 'ECHO' 0 ; MESS '----------- Ra = ' Ra '------------------------------------'; MESS ' NUSSELT MOYEN NUSSELT MAX NUSSELT MIN' ; MESS 'CAST3M : Gauche ' NUMG ' ' (MAXI Nusg) ' ' (MINI Nusg) ; MESS 'CAST3M : Droite ' NUMD ' ' (MAXI Nusd) ' ' (MINI Nusd) ; Si (EGA Ra 1.e5); MESS 'Vahl Davis : Ra 1.e6 4.519 7.717 0.729'; Finsi ; Si (EGA Ra 1.e6); MESS 'Vahl Davis : Ra 1.e6 8.8 17.925 0.989'; Finsi ; MESS '---------------------------------------------------------------'; MESS ' Umax Vmax ' ; MESS 'CAST3M : ' (MAXI (exco un 'UX')) (MAXI (exco un 'UY')); Si (EGA Ra 1.e5); MESS 'Vahl Davis : Ra 1.e6 34.73 68.59 ' ; Finsi ; Si (EGA Ra 1.e6); MESS 'Vahl Davis : Ra 1.e6 64.63 219.36 ' ; Finsi ; MESS '---------------------------------------------------------------'; MESS ' Psi max ' ; MESS 'CAST3M : ' (MAXI PSI) ; Si (EGA Ra 1.e5); MESS 'Vahl Davis : Ra 1.e6 9.612 ' ; Finsi; Si (EGA Ra 1.e6); MESS 'Vahl Davis : Ra 1.e6 16.75 ' ; Finsi; MESS '---------------------------------------------------------------'; opti 'ECHO' 1 ; *opti donn 5; Si (NON COMPLET) ; ern=abs( NUMD - 6.32)/6.32 ; ernm=abs((MAXI Nusg) - 9.31)/9.31 ; eru=abs((MAXI (exco un 'UY')) - 154.)/154.; erp=abs((MAXI PSI) - 11.1)/11.1; mess 'ern = ' ern ' ernm = ' ernm ' eru = ' eru ' erp = ' erp ; Si (ern > 0.005); erreur 5 ; Finsi ; Si (ernm > 0.005); erreur 5 ; Finsi ; Si (eru > 0.005); erreur 5 ; Finsi ; Si (erp > 0.005); erreur 5 ; Finsi ; Finsi ; FIN ;