Accueil | Cast3M

* 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=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 ;