* fichier : BINGHAMp.dgibi
** 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=resu( 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 ;
© Cast3M 2003 - All rights reserved.
Disclaimer