Numérotation des lignes :

* 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 100OPTI EFM1 SUPGZONE$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 $mtzone$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=SOMM( abs (vcr - vc) ) ;mess ' Ecart sur profil de V : ' ER ; si ( er > err1 ) ; erreur 5 ; finsi ; FINSI ; FIN ;

© Cast3M 2003 - Tous droits réservés.
Mentions légales