* fichier : hy1.dgibi ************************************************************************ * Section : Fluides Transitoire ************************************************************************ *$$$ HY1 ** Exemple HY1 ** ___________ ** ** --- 2 OCTOBRE 1990 --- ** ** CANAL LONGUEUR 10. LARGEUR 1. ** test cas isotherme NS et TOIMP ** ** On considère l'ecoulement de poiseuille dans un canal plan vertical ** compris entre les plans x=0 et x=1. la hauteur est de 10. ** On maintient entre l'entrée et la sortie une différence de pression. ** En regime établi le profil de vitesse ne depend que de x et la ** pression que de y : V(x)=4Vm X(1-X) Vm vitesse max (x=1/2) et ** Vmoy=2/3Vm Vm=-1/8/Mu DP/DY Mu viscosité dynamique ** T0=Mu DV/DX =4Mu Vm et le Reynolds Re=1/12/Mu/Mu DP/DY ** (RO=1 Diam=1) ** On a choisi Mu=5.e-2 et DP/DY=0.6 => Re=20 Vm=1.5 Vmoy=1 ** Avec ces valeurs le temps d'établissement du régime est dominé par ** la diffusion. Il faut 12s pour arriver au régime établi (Vmoy=0.99) ** soit approximativement 800 pas de temps avec des QUA8 et le double avec ** des TRI6. ** Si on réduit la viscosité (d'un facteur 10 par exemple) la vitesse ** augmente pour une même perte de charge et le temps d'établissement ** du régime est dominé‚par la convection. ** Dans le cas ou le canal est coudé (IKAS=1) on peut noter une ** différence sur les isobares lorsque la viscosité est faible ** (Mu=5.e-3). Ceci est du à l'effet d'inertie dans le coude. ** ** On peut tester aussi le cas ou on impose la tension en paroi. ** GRAPH = 'N' ; COMPLET = FAUX ; SI ( COMPLET ) ; nbit3=2150 ; err1 =1.e-5 ; qe13= 0.9922268 ; qe14= 0.9922544 ; dq1 =1.e-14 ; dpes1 = 5.6997 ; erpe1 = 1.e-3 ; nbit4=1000 ; tfin=12. ; SINON ; nbit3=180 ; err1 =1.e-5 ; err1 =1.e-5 ; qe14 =0.3989106 ; qe13 =0.3956187 ; dq1 =1.e-14 ; dpes1 = 5.6997 ; erpe1 = 1.e-3 ; nbit4=80 ; tfin=1. ; FINSI ; typelt=mot qua8 ; nbe=7 ; nbv=10 ; angle=0. ; nitma=nbit4 ; nitma=100 ; tfinal=tfin ; DEBPROC HY1 ; ARGU TYPELT*MOT NBE*ENTIER NBV*ENTIER ANGLE*FLOTTANT NITMA*ENTIER TFINAL*FLOTTANT GRAPH*MOT ; OPTION DIME 2 ELEM TYPELT ; p1=0 0.; p2=1. 0.; VH=0 10; entree= p1 d nbe p2 ; c1= 5. 0. ; ikas=0 ; *mess 'ikas= 0 --> DROIT (option par defaut) ikas=1 --> COURBE ? '; *obtenir ikas*entier ; si (EGA ikas 1) ; q1=p1 tour c1 -90. ; q2=p2 tour c1 -90. ; pp1=p1 c c1 q1 nbv; pp2=p2 c c1 q2 nbv; sortie=entree tour c1 -90 ; sinon ; q1=p1 plus vh ; q2=p2 plus vh ; pp1=p1 d q1 nbv; pp2=p2 d q2 nbv; sortie=entree plus vh ; finsi ; pp1=inve pp1 ; sortie=inve sortie; elim 0.0001 (sortie et pp1 et pp2 et entree ); cnt=entree et pp2 et sortie et pp1 ; *mt=surf cnt ; mt=daller entree pp2 sortie pp1 ; *pp1=inve pp1 ; entree=inve entree ; sortie=inve sortie ; *pp2=inve pp2 ; *mt=daller entree pp1 sortie pp2 ; *entree=inve entree ; mt=mt tour c1 angle ; entree=entree tour c1 angle ; sortie=sortie tour c1 angle ; pp1=pp1 tour c1 angle ; pp2=pp2 tour c1 angle ; elim (pp1 et pp2 et entree et sortie et mt) 0.0001 ; bell = chan mt quaf ; entree = chan entree quaf ; sortie = chan sortie quaf ; pp1 = chan pp1 quaf ; pp2 = chan pp2 quaf ; elim (mt et entree et sortie et pp1 et pp2) 1.e-4 ; $bell = mode bell 'NAVIER_STOKES' MACRO ; doma $bell 'IMPR' ; $entree = mode entree 'NAVIER_STOKES' MACRO ; doma $entree 'IMPR' ; $sortie = mode sortie 'NAVIER_STOKES' MACRO ; doma $sortie 'IMPR' ; $pp1 = mode pp1 'NAVIER_STOKES' MACRO ; doma $pp1 'IMPR' ; $pp2 = mode pp2 'NAVIER_STOKES' MACRO ; doma $pp2 'IMPR' ; nble=(nbel (doma $entree 'MAILLAGE') ) - 1 ; entref=elem (doma $entree 'MAILLAGE') (lect 2 pas 1 nble) ; *?$entref=mode entref 'NAVIER_STOKES' MACRO ;doma $entref 'IMPR' ; mu=5.E-2 ; ro=1 ; nu=mu/ro ; to = kcht $entree vect centre ( 0. 6.) ; tos = kcht $sortie vect centre ( 0. 0.) ; tt1 = kcht $pp1 vect centre (-0.03 0.) ; tt2 = kcht $pp2 vect centre (-0.03 0.) ; dt=1.e-2 ; rv=eqex $bell 'DUMP' ALFA 0.3 ITMA nitma 'TFINAL' tfinal ZONE $BELL OPER NS NU INCO 'UN' OPTI 'CENTREE' ZONE $BELL OPER DFDT 1. 'UN' 'DELTAT ' INCO 'UN' ZONE $ENTREE OPER TOIMP TO INCO 'UN' ZONE $SORTIE OPER TOIMP TOS INCO 'UN' **ZONE $ENTREE OPER TOIMP TO 'UN' nu INCO 'UN' **ZONE $SORTIE OPER TOIMP TOS 'UN' nu INCO 'UN' * ZONE $PP1 OPER TOIMP TT1 INCO 'UN' * ZONE $PP2 OPER TOIMP TT2 INCO 'UN' ; ; rvp= eqpr $bell zone $PP1 oper VNIMP 0. zone $bell oper PRESSION 0. zone $PP1 oper VTIMP 0. zone $PP2 oper VNIMP 0. zone $PP2 oper VTIMP 0. * zone $ENTREF oper VNIMP 1. * zone $ENTREF oper VTIMP 0. ; rv.'INCO'=table 'INCO' ; rv.'INCO'.'UN' = kcht $bell VECT SOMMET (0 0) ; rv.pression=rvp ; lh= (noeu 10) et (noeu 20) et (noeu 30) et (noeu 40) et (noeu 50) et (noeu 60) ; lj= (manu poi1 ((doma $bell 'MAILLAGE') poin proc( 0.5 0.5) ) ) ; his = khis 'UN' 1 lh 'UN' 2 (lh et lj) ; rv.'HIST'=his ; exec rv ; pn=elno (kcht $bell scal centre (rvp.pression)) $bell ; ung1 = vect 0.5 (rv.'INCO'.'UN') ux uy jaune ; qe=dbit (rv.'INCO'.'UN') $entree ; qs=dbit (rv.'INCO'.'UN') $sortie ; dq=(abs qs )-(abs qe) ; mess (' Bilan : dq=') dq ; pet=somt (redu pn entree) / (nbno entree) ; pst=somt (redu pn sortie) / (nbno sortie) ; dpes=pet - pst ; mess ' pet pst dpes = ' pet pst dpes ; si ('EGA' graph 'O' ); dessin (his.'TABD') (his.'1UN') ; dessin (his.'TABD') (his.'2UN') ; trace ung1 bell ; trace pn bell (cont bell) ; finsi ; FINPROC RV dq qe dpes ; type=mot tri6 ; * mess 'Type d élément TRI6 ou QUA8 ? ' ; * obtenir TYPE*mot ; angle=0. ; * obtenir angle*flottant ; nbe=7 ; nbv=10 ; type=mot qua8 ; rv dq qe dpes = hy1 TYPE nbe nbv angle nbit4 tfin graph ; erq = abs (qe - qe14) ; mess ' Erreur sur le debit ' erq ; si ( erq > err1 ) ; erreur 5 ; finsi ; dq=abs dq ; si ( dq > dq1 ) ; erreur 5 ; finsi ; * On teste quand meme le delta P qui viendrait d'une erreur de discretisation ddp=abs ( dpes - dpes1) ; mess ' ddp erpe1 ' ddp erpe1 ; si ( ddp > erpe1 ) ; erreur 5 ; finsi ; *opti donn 5 ; type=mot tri6 ; rv dq qe dpes = hy1 TYPE nbe nbv angle nbit3 tfin graph ; erq = abs (qe - qe13) ; mess ' Erreur sur le debit ' erq ; si ( erq > err1 ) ; erreur 5 ; finsi ; dq=abs dq ; si ( dq > dq1 ) ; erreur 5 ; finsi ; FIN ;