* fichier : ODWp.dgibi
** Fluide PSEUDOPLASTIQUE : Ecoulement de POISEUILLE
** loi d'OSTWALD de WAELE nn= 0.8, 0.6 , 0.4
** CANAL LONGUEUR 10. LARGEUR 4.
** Algorithme semi - implicite
** Auteur : Isabelle Claudel Décembre 1997
GRAPH='N' ;
COMPLET= FAUX ;
ERR1=5.E-5 ;
OPTION DIME 2 ELEM QUA4 ;
***************************************
***PROCEDURE CALCUL DE LA VISCOSITE ***
***************************************
DEBP CALCUL ;
*** loi d'OSTWALD de WAELE ***
ARGU RX*TABLE ;
iarg=rx.'IARG' ;
si( non ( ega iarg 3)) ;
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 ;
si ( ega ('TYPE' rx.'ARG3') 'MOT ') ;
NN=rv.'INCO'.(rx.'ARG3') ;
sinon ;
mess 'Procedure CALCUL : type argument invalide ' ;
quitter CALCUL ;
finsi ;
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 = kops (kops D11 '*' D11) '+' (kops D12 '*' D12) ;
ProTen = kops ProTen '+' (kops D22 '*' D22) ;
ProTen = kops ProTen '+' (kops D21 '*' D21) ;
ProTen = kops ProTen '/' 2. ;
ProTen = kops ProTen '**' 0.5 ;
*list ProTen ;
**Comportement rhéologique pseudoplastique **
MUi = 10. ;
MU0 = 100. ;
nexp = nn - 1. ;
inv = 1. ;
inexp = inv / nexp ;
Dseuil = ( MU0 / MUi) ** inexp ;
TESTi = kcht $mt scal centre (ProTen MASQUE 'EGINFE' Dseuil) ;
TESTs = kcht $mt scal centre (ProTen MASQUE 'SUPERIEUR' Dseuil) ;
ProTen=kops ( kops ProTen '*' TESTs ) '+' ( kops Dseuil '*' TESTi);
MU = kops MUi '*' (kops (kops ProTen '*' 2) '**' nexp) ;
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 ***
mtQ=chan mt QUAF ;
$mt=mode mtQ 'NAVIER_STOKES' LINE ;
entreeQ=chan entree QUAF ;
sortieQ=chan sortie QUAF ;
pp1Q =chan pp1 QUAF ;
pp2Q =chan pp2 QUAF ;
elim ( mtQ et entreeQ et sortieQ et pp1Q et pp2Q) 1.e-4 ;
$entree = mode entreeQ 'NAVIER_STOKES' LINE ;
$sortie = mode sortieQ 'NAVIER_STOKES' LINE ;
$pp1 = mode pp1Q 'NAVIER_STOKES' LINE ;
$pp2 = mode pp2Q 'NAVIER_STOKES' LINE ;
***dP/dy = - 20***
to = kcht $entree vect centre ( 0. 0.) ;
tos = kcht $sortie vect centre ( 0. 200.) ;
rv=eqex $mt 'DUMP' ALFA 0.9 ITMA 100
OPTI EFM1 SUPG
ZONE $mt 'OPER' CALCUL 'UN' 'NU' 'NN'
ZONE $mt 'OPER' NS 'NU' INCO 'UN'
OPTI EFM1 'CENTREE'
ZONE $mt 'OPER' DFDT 1. 'UN' 'DELTAT' INCO 'UN'
ZONE $entree 'OPER' 'TOIMP' TO INCO 'UN'
ZONE $sortie 'OPER' 'TOIMP' TOS INCO 'UN'
;
rvp = eqpr $mt
zone $mt oper PRESSION 0.
zone $pp1 oper VNIMP 0.
zone $pp1 oper VTIMP 0.
zone $pp2 oper VNIMP 0.
zone $pp2 oper VTIMP 0.
;
rv.'INCO'=table 'INCO' ;
rv.'INCO'.'UN' = kcht $mt VECT SOMMET (0. 0.) ;
rv.'INCO'.'NU' = kcht $mt SCAL CENTRE 10. ;
rv.'INCO'.'GU1'= kcht $mt SCAL SOMMET 0. ;
rv.pression=rvp ;
lh= (manu poi1 ((doma $mt maillage) poin proc(2. 2.) ) ) ;
his = khis 'UN' 2 lh ;
rv.'HIST' = his ;
rv.'TFINAL' = 10.e30 ;
DEBPROC POST RV*TABLE GRAPH*MOT ;
evol2v08 = EVOL 'CHPO' (rv.'INCO'.'UN') UY entree2 ;
*** Solution analytique ***
x = prog 0. PAS 0.25 2. ;
unite = prog 9.* 1. ;
nn = rv.'INCO'.'NN' ;
L = 2. ;
invn = 1. / nn ;
list invn ;
nn1 = (nn + 1.) / nn ;
list nn1 ;
nn2 = 1. / nn1 ;
list nn2 ;
MUi = 10. ;
ddp = 20. / (MUi) ;
list ddp ;
ddpn = EXP (invn * (LOG ddp)) ;
list ddpn ;
LLn = L ** nn1 ;
list LLn ;
Vx = nn2 * ddpn * LLn * (unite - ((x / L) ** nn1)) ;
Vx = Vx * -1. ;
****** Graphe représentant la vitesse calculée avec la vitesse analytique *****
evol1v08 = EVOL 'MANU' largeur x Vitesse Vx ;
evolt = evol1v08 et evol2v08;
rv.'RET'=evolt ;
rv.'Vcal'=evol2v08 ;
TAB1=TABLE;
TAB1.'TITRE'=TABLE ;
TAB1.1='MARQ ETOI REGU ' ;
TAB1.'TITRE' . 1=mot ' Vitesse_théorique ' ;
TAB1.2='TIRR REGU ';
TAB1.'TITRE' . 2=mot ' Vcalculée ';
si (EGA GRAPH 'O' );
dessin (his.'TABD') (his.'2UN') ;
ung1 = vect 0.5 (rv.'INCO'.'UN') ux uy jaune ;
trace ung1 mt ;
DESS evolt 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
FINSI ;
FINPROC ;
SI COMPLET ;
RV.INCO.'NN'=0.8 ;
RV.ITMA= 1500 ;
EXEC RV ;
post rv GRAPH ;
evolt=rv.'RET' ;
RV.INCO.'NN'=0.6 ;
RV.ITMA= 3500 ;
EXEC RV ;
post rv GRAPH ;
evolt=evolt et rv.'RET' ;
TAB1=TABLE;
TAB1.'TITRE'=TABLE ;
TAB1.1='MARQ ETOI REGU ' ;
TAB1.'TITRE' . 1=mot ' V_théorique N=0.8' ;
TAB1.2='TIRR REGU ';
TAB1.'TITRE' . 2=mot ' Vcalculée N=0.8';
TAB1.3='MARQ PLUS REGU ' ;
TAB1.'TITRE' . 3=mot ' V_théorique N=0.6' ;
TAB1.4='TIRR REGU ';
TAB1.'TITRE' . 4=mot ' Vcalculée N=0.6';
DESS evolt 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
SINON ;
RV.INCO.'NN'=0.8 ;
RV.ITMA= 40 ;
EXEC RV ;
post rv GRAPH ;
ev=rv.'Vcal' ;
vc=extr ev 'ORDO' ;
list vc ;
vcr=prog -.30694 -.30664 -.30531 -.30203 -.29464
-.27818 -.24186 -.16326 1.88015E-17 ;
ER=resu( abs (vcr - vc) ) ;
mess ' Ecart sur profil de V : ' ER ;
si ( er > err1 ) ; erreur 5 ; finsi ;
FINSI ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer