* fichier : BINGHAMp.dgibi ************************************************************************ * Section : Fluides Convection ************************************************************************ ** Fluide de BINGHAM : Ecoulement de POISEUILLE ** Comparaison à une solution analytique ** CANAL LONGUEUR 10. LARGEUR 4. ** Algorithme implicite (dans ce cas beaucoup plus rapide) ** Auteur : Isabelle Claudel Décembre 1997 GRAPH='N' ; COMPLET= FAUX ; ERR1=5.E-2 ; option dime 2 elem qua4 ; opti ISOV suli ; *************************************** ***PROCEDURE CALCUL DE LA VISCOSITE *** *************************************** DEBP CALCUL ; ARGU RX*TABLE ; iarg=rx.'IARG' ; si( non ( ega iarg 2)) ; mess 'Procedure CALCUL : nombre d arguments incorrect ' iarg ; quitter CALCUL ; finsi ; si ( ega ('TYPE' rx.'ARG1') 'MOT ') ; UN=rv.'INCO'.(rx.'ARG1') ; sinon ; mess 'Procedure CALCUL : type argument invalide ' ; quitter CALCUL ; finsi ; si ( ega ('TYPE' rx.'ARG2') 'MOT ') ; NU=rv.'INCO'.(rx.'ARG2') ; sinon ; mess 'Procedure CALCUL : type argument invalide ' ; quitter CALCUL ; finsi ; *mess ' Proc CALCUL ' ; UN1= exco ux UN ; UN2= exco uy UN ; UNN1=NOMC UN1 'SCAL' ; UNN2=NOMC UN2 'SCAL' ; unn1= kcht $mt scal sommet UNN1; unn2= kcht $mt scal sommet UNN2; GU1 = kops UNN1 'GRAD' $mt ; GU2 = kops UNN2 'GRAD' $mt ; ddudx= exco ux GU1 ; ddudy= exco uy GU1 ; ddvdx= exco ux GU2 ; ddvdy= exco uy GU2 ; dudx= NOMC ddudx 'SCAL' ; dudy= NOMC ddudy 'SCAL' ; dvdx= NOMC ddvdx 'SCAL' ; dvdy= NOMC ddvdy 'SCAL' ; dudx= kcht $mt scal centre dudx ; dudy= kcht $mt scal centre dudy ; dvdx= kcht $mt scal centre dvdx ; dvdy= kcht $mt scal centre dvdy ; D11= dudx ; D12= (kops (kops dudy '+' dvdx) '/' 2) ; D21= (D12) ; D22= (dvdy) ; *** Proten = Produit tensoriel D:D *** ProTen = kops (kops D11 '*' D11) '+' (kops D12 '*' D12) ; ProTen = kops ProTen '+' (kops D22 '*' D22) ; ProTen = kops ProTen '+' (kops D21 '*' D21) ; rv.'INCO'.'ProTen' = ProTen ; *list ProTen ; *** Norme de D = racine carree (Proten / 2) *** NormD = kops ProTen '/' 2. ; NormD = kops NormD '**' 0.5 ; **Comportement rhéologique de Bingham ** Seuil = kcht $mt scal centre 10. ; MUB = kcht $mt scal centre 10. ; MU0 = kcht $mt scal centre 1000. ; Dseuil = kops Seuil '/' (kops 2. '*' (kops MU0 '-' MUB) ) ; TESTi = kcht $mt scal centre (NormD MASQUE 'EGINFE' Dseuil) ; TESTs = kcht $mt scal centre (NormD MASQUE 'SUPERIEUR' Dseuil) ; NormD = kops ( kops NormD '*' TESTs ) '+' ( kops Dseuil '*' TESTi) ; MU = kops MUB '+' (kops Seuil '/' (kops 2. '*' NormD) ) ; ro = kcht $mt scal centre 1. ; NU = kops MU '/' ro ; rv.'INCO'.(rx.'ARG2')=NU ; rv.'INCO'.'GU1'=D12 ; as2 ama1 = 'KOPS' 'MATRIK' ; FINPROC as2 ama1 ; *********************** *******MAILLAGE******** *********************** nbe=16 ; nbv=16 ; nbee = nbe / 2 ; p1=0. 0.; p2=4. 0.; entree= p1 d nbe p2 ; c1= 5. 0. ; isa = 2. 0. ; entree2 = isa d nbee p2 ; 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 (0 10) ; q2=p2 plus (0 10) ; pp1=p1 d q1 nbv; pp2=p2 d q2 nbv; sortie=entree plus (0 10) ; finsi ; mientree = entree plus (0. 5.) ; 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 ; angle=0. ; mt=mt tour p1 angle ; entree=entree tour p1 angle ; sortie=sortie tour p1 angle ; pp1=pp1 tour p1 angle ; pp2=pp2 tour p1 angle ; elim (pp1 et pp2 et entree et sortie et mt et entree2) 0.0001 ; *ccc=2. 5. ; *option donn 5 ; mt = mt 'COUL' bleu ; *trace mt ; ************************ *** RESOLUTION *** ************************ mt0= mt ; mt=chan mt quaf ; $mt=mode mt 'NAVIER_STOKES' LINE ; doma $mt 'IMPR' ; entre1= chan entree quaf ; sorti1= chan sortie quaf ; elim (mt et entre1 et sorti1 et pp1 et pp2 ) 1.e-4 ; $entree = mode entre1 'NAVIER_STOKES' LINE ; $sortie = mode sorti1 'NAVIER_STOKES' LINE ; * doma $entree 'IMPR' ; * doma $sortie 'IMPR' ; ***dP/dy = - 20*** to = kcht $entree vect centre ( 0. 0.) ; tos = kcht $sortie vect centre ( 0. 200.) ; *to = kcht $entree vect centre ( 0. 6.) ; *tos = kcht $sortie vect centre ( 0. 0.) ; EPS=1.E-6 ; DT=1. ; RV = EQEX $mt OPTI EF IMPL 'CENTREE' ZONE $mt 'OPER' CALCUL 'UN' 'NU' ZONE $mt 'OPER' DUDW EPS INCO 'UN' OPTI EF IMPL 'CENTREE' ZONE $mt 'OPER' NS 1. 'UN' 'NU' INCO 'UN' ZONE $entree 'OPER' 'TOIMP' TO INCO 'UN' ZONE $sortie 'OPER' 'TOIMP' TOS INCO 'UN' ; RV=EQEX RV CLIM UN UIMP (pp1 et pp2) 0. UN VIMP (pp1 et pp2) 0. ; rv.'INCO'=table 'INCO' ; rv.'INCO'.'UN' = kcht $mt VECT SOMMET (0. 0.) ; rv.'INCO'.'NU' = kcht $mt SCAL CENTRE 10. ; rv.'INCO'.'ProTen'= kcht $mt SCAL SOMMET 0. ; DEBPROC POST ; ARGU RV*TABLE GRAPH*MOT ; *** POST TRAITEMENT *** evol2v10 = EVOL 'CHPO' (rv.'INCO'.'UN') UY (entree2) ; *list evol2v ; XX=extr evol2v10 'ABSC' ; vth10 = prog -2.25 -2.25 -2.25 -2.1875 -2. -1.6875 -1.25 -0.6875 0. ; vc=extr evol2v10 'ORDO' ; ER=SOMM( abs ((vth10 - vc)/2.5) ) ; mess ' Ecart sur profil de V : ' ER ; evol1v10 = EVOL 'MANU' Rayon xx Vitesse vth10 ; evolt = evol1v10 et evol2v10 ; TAB1=TABLE ; TAB1.'TITRE'=TABLE ; TAB1.1='MARQ ETOI REGU ' ; TAB1.'TITRE' . 1=MOT 'V_theorique ' ; TAB1.2='TIRR REGU ' ; TAB1.'TITRE' . 2=MOT 'V_calculee ' ; un=rv.INCO.'UN' ; ung1=vect un 0.5 ux uy jaune ; si (EGA GRAPH 'O' ); DESS evolt 'TITX' 'X (m)' 'TITY' 'V (m/s)' 'LEGE' TAB1 ; *list evolt ; mt0 = doma $mt maillage ; trace ung1 mt0 ; nu = elno (rv.'INCO'.'NU') $mt ; trace nu mt0 12 ; FINSI ; si ( er > err1 ) ; erreur 5 ; finsi ; FINPROC ; rv.'OMEGA'=0.9; rv.'NITER'=20 ; EXEC RV ; POST RV GRAPH ; FIN ;