Télécharger la plaquette Cast3M


fr en

Télécharger ODWp.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  ODWp.dgibi
  2. **  Fluide PSEUDOPLASTIQUE : Ecoulement de POISEUILLE
  3. **  loi d'OSTWALD de WAELE  nn= 0.8, 0.6 , 0.4
  4. **  CANAL LONGUEUR 10. LARGEUR 4.
  5. **  Algorithme semi - implicite
  6. **  Auteur : Isabelle Claudel  Décembre 1997
  7.  
  8. GRAPH='N' ;
  9. COMPLET= FAUX ;
  10. ERR1=5.E-5 ;
  11.  
  12. OPTION DIME 2 ELEM QUA4 ;
  13.  
  14.  
  15. ***************************************
  16. ***PROCEDURE CALCUL DE LA VISCOSITE ***
  17. ***************************************
  18.  
  19. DEBP CALCUL ;
  20. *** loi d'OSTWALD de WAELE ***
  21. ARGU RX*TABLE ;
  22. iarg=rx.'IARG' ;
  23.  
  24. si( non ( ega iarg 3)) ;
  25. mess 'Procedure CALCUL : nombre d arguments incorrect ' iarg ;
  26. quitter CALCUL ;
  27. finsi ;
  28.  
  29. si ( ega ('TYPE' rx.'ARG1') 'MOT     ') ;
  30. UN=rv.'INCO'.(rx.'ARG1') ;
  31. sinon ;
  32. mess 'Procedure CALCUL : type argument invalide ' ;
  33. quitter CALCUL ;
  34. finsi ;
  35.  
  36. si ( ega ('TYPE' rx.'ARG2') 'MOT     ') ;
  37. NU=rv.'INCO'.(rx.'ARG2') ;
  38. sinon ;
  39. mess 'Procedure CALCUL : type argument invalide ' ;
  40. quitter CALCUL ;
  41. finsi ;
  42.  
  43. si ( ega ('TYPE' rx.'ARG3') 'MOT     ') ;
  44. NN=rv.'INCO'.(rx.'ARG3') ;
  45. sinon ;
  46. mess 'Procedure CALCUL : type argument invalide ' ;
  47. quitter CALCUL ;
  48. finsi ;
  49.  
  50. UN1=  exco 'UX' UN ;
  51. UN2=  exco 'UY' UN ;
  52.  
  53. UNN1=NOMC UN1 'SCAL' ;
  54. UNN2=NOMC UN2 'SCAL' ;
  55.  
  56. unn1= kcht $mt scal sommet UNN1;
  57. unn2= kcht $mt scal sommet UNN2;
  58.  
  59. GU1 = kops  UNN1 'GRAD'  $mt ;
  60. GU2 = kops  UNN2 'GRAD'  $mt ;
  61.  
  62. ddudx= exco 'UX' GU1 ;
  63. ddudy= exco 'UY' GU1 ;
  64. ddvdx= exco 'UX' GU2 ;
  65. ddvdy= exco 'UY' GU2 ;
  66.  
  67. dudx= NOMC ddudx 'SCAL' ;
  68. dudy= NOMC ddudy 'SCAL' ;
  69. dvdx= NOMC ddvdx 'SCAL' ;
  70. dvdy= NOMC ddvdy 'SCAL' ;
  71.  
  72. dudx= kcht $mt scal centre dudx ;
  73. dudy= kcht $mt scal centre dudy ;
  74. dvdx= kcht $mt scal centre dvdx ;
  75. dvdy= kcht $mt scal centre dvdy ;
  76.  
  77.  
  78. D11= dudx ;
  79. D12= (kops (kops dudy '+' dvdx) '/' 2) ;
  80. D21= (D12)  ;
  81. D22= (dvdy) ;
  82.  
  83. ProTen = kops (kops D11 '*' D11) '+' (kops D12 '*' D12)  ;
  84. ProTen = kops ProTen '+' (kops D22 '*' D22)  ; 
  85. ProTen = kops ProTen '+' (kops D21 '*' D21)  ;  
  86. ProTen = kops ProTen '/' 2. ;
  87. ProTen = kops ProTen '**' 0.5 ;
  88.  
  89. *list ProTen ;
  90.  
  91. **Comportement rhéologique pseudoplastique **
  92.  
  93. MUi   =  10. ;
  94.  
  95. MU0 =  100. ;
  96. nexp =  nn - 1. ;
  97. inv =  1. ;
  98. inexp = inv / nexp  ;
  99. Dseuil =  ( MU0 / MUi) ** inexp ;   
  100.  
  101. TESTi = kcht $mt scal centre (ProTen MASQUE 'EGINFE' Dseuil) ;
  102. TESTs = kcht $mt scal centre (ProTen MASQUE 'SUPERIEUR' Dseuil) ;
  103. ProTen=kops ( kops ProTen '*' TESTs ) '+' ( kops Dseuil '*' TESTi);
  104.  
  105. MU =  kops  MUi '*' (kops (kops ProTen '*' 2) '**' nexp) ;
  106.  
  107. ro = kcht $mt scal centre 1. ;
  108. NU = kops MU '/' ro ;
  109.  
  110.  
  111. rv.'INCO'.(rx.'ARG2')=NU ;
  112.  
  113. rv.'INCO'.'GU1'=D12 ;
  114. as2 ama1  = 'KOPS' 'MATRIK' ;
  115. FINPROC as2 ama1 ;
  116.  
  117.  
  118.  
  119. ***MAILLAGE***
  120.  
  121. nbe=16 ; nbv=16 ;
  122. nbee = nbe / 2 ;
  123.  
  124. p1=0. 0.;
  125. p2=4. 0.;
  126. entree= p1 d nbe p2 ;
  127. c1= 5. 0. ;
  128. isa = 2. 0. ;
  129.  
  130. entree2 = isa d nbee p2 ;
  131.  
  132. ikas=0 ;
  133. *mess 'ikas= 0 --> DROIT (option par defaut)  ikas=1 --> COURBE  ? ';
  134. *obtenir ikas*entier ;
  135.  
  136. si (EGA ikas 1) ;
  137.  
  138. q1=p1 tour c1 -90. ;
  139. q2=p2 tour c1 -90. ;
  140. pp1=p1 c c1 q1 nbv;
  141. pp2=p2 c c1 q2 nbv;
  142. sortie=entree tour c1 -90 ;
  143.  
  144. sinon ;
  145.  
  146. q1=p1 plus (0 10)  ;
  147. q2=p2 plus (0 10)  ;
  148. pp1=p1 d q1 nbv;
  149. pp2=p2 d q2 nbv;
  150. sortie=entree plus (0 10) ;
  151. finsi ;
  152. mientree = entree plus (0. 5.) ;
  153.  
  154. pp1=inve pp1 ;
  155. sortie=inve sortie;
  156. elim 0.0001 (sortie et pp1 et pp2 et entree );
  157. cnt=entree et pp2 et sortie et pp1 ;
  158.  
  159. *mt=surf cnt ;
  160. mt=daller entree pp2 sortie pp1 ;
  161.  
  162. angle=0. ;
  163.  
  164. mt=mt tour p1 angle ;
  165. entree=entree tour p1 angle ;
  166. sortie=sortie tour p1 angle ;
  167. pp1=pp1 tour p1 angle ;
  168. pp2=pp2 tour p1 angle ;
  169.  
  170. elim (pp1 et pp2 et entree et sortie et mt et entree2) 0.0001 ;
  171. *ccc=2. 5. ;
  172.  
  173. *option donn 5 ;
  174. mt = mt 'COUL' bleu ;
  175.  
  176. *trace mt ;
  177.  
  178.  
  179. *** RESOLUTION ***
  180.  
  181.  mtQ=chan mt QUAF ;
  182. $mt=mode mtQ 'NAVIER_STOKES' LINE ;
  183.  
  184.  entreeQ=chan entree QUAF ;
  185.  sortieQ=chan sortie QUAF ;
  186.  pp1Q   =chan pp1    QUAF ;
  187.  pp2Q   =chan pp2    QUAF ;
  188.  elim ( mtQ et entreeQ et sortieQ et pp1Q et pp2Q) 1.e-4 ;
  189.  
  190. $entree = mode entreeQ 'NAVIER_STOKES' LINE ;
  191. $sortie = mode sortieQ 'NAVIER_STOKES' LINE ;
  192. $pp1    = mode pp1Q 'NAVIER_STOKES' LINE ;
  193. $pp2    = mode pp2Q 'NAVIER_STOKES' LINE ;
  194.  
  195. ***dP/dy = - 20***
  196. to  = kcht $entree vect centre ( 0.  0.) ;
  197. tos = kcht $sortie vect centre ( 0.  200.) ;
  198.  
  199.  
  200. rv=eqex $mt 'DUMP' ALFA 0.9 ITMA 100
  201. OPTI EFM1 SUPG
  202. ZONE  $mt    'OPER'   CALCUL 'UN' 'NU' 'NN'
  203. ZONE  $mt    'OPER'    NS    'NU'            INCO 'UN'
  204. OPTI EFM1 'CENTREE'
  205. ZONE  $mt    'OPER'    DFDT 1. 'UN' 'DELTAT' INCO 'UN'
  206. ZONE  $entree  'OPER'   'TOIMP' TO                INCO 'UN'
  207. ZONE  $sortie  'OPER'   'TOIMP' TOS               INCO 'UN'
  208. ;
  209.  
  210.  
  211. rvp = eqpr $mt
  212. zone $mt   oper PRESSION  0.
  213. zone $pp1    oper VNIMP     0.
  214. zone $pp1    oper VTIMP     0. 
  215. zone $pp2    oper VNIMP     0.
  216. zone $pp2    oper VTIMP     0.
  217. ;
  218.  
  219. rv.'INCO'=table 'INCO' ;
  220. rv.'INCO'.'UN' = kcht $mt VECT SOMMET  (0. 0.) ;
  221. rv.'INCO'.'NU' = kcht $mt SCAL CENTRE  10. ; 
  222. rv.'INCO'.'GU1'= kcht $mt SCAL SOMMET  0.  ;
  223. rv.pression=rvp ;
  224.  
  225. lh= (manu poi1 ((doma $mt maillage) poin proc(2. 2.) ) ) ;
  226. his = khis 'UN' 2 lh ; 
  227. rv.'HIST' = his ;
  228. rv.'TFINAL' = 10.e30 ;
  229.  
  230.  
  231. DEBPROC POST RV*TABLE GRAPH*MOT ;
  232.  
  233. evol2v08 = EVOL 'CHPO' (rv.'INCO'.'UN') UY entree2 ;
  234.  
  235. *** Solution analytique ***
  236. x = prog 0. PAS 0.25 2. ;
  237. unite = prog 9.* 1. ;
  238. nn = rv.'INCO'.'NN' ;
  239. L = 2. ;
  240.  
  241. invn = 1. / nn ;
  242. list invn ;
  243.  
  244. nn1 = (nn + 1.) / nn ;
  245. list nn1 ;
  246. nn2 = 1. / nn1 ;
  247. list nn2 ;
  248.  
  249. MUi = 10. ;
  250.  
  251. ddp = 20. / (MUi) ;
  252. list ddp ;
  253.  
  254. ddpn = EXP (invn * (LOG ddp)) ;
  255. list ddpn ;
  256.  
  257. LLn = L ** nn1 ;
  258. list LLn ; 
  259.  
  260. Vx = nn2 * ddpn * LLn * (unite - ((x / L) ** nn1)) ;
  261. Vx = Vx * -1. ;
  262.  
  263.  
  264. ****** Graphe représentant la vitesse calculée avec la vitesse analytique *****
  265.  
  266. evol1v08 = EVOL 'MANU' largeur x Vitesse Vx ;
  267. evolt = evol1v08 et evol2v08;
  268. rv.'RET'=evolt ;
  269. rv.'Vcal'=evol2v08 ;
  270. TAB1=TABLE;
  271. TAB1.'TITRE'=TABLE ;
  272. TAB1.1='MARQ ETOI REGU '     ;
  273. TAB1.'TITRE' . 1=mot ' Vitesse_théorique '     ;
  274. TAB1.2='TIRR      REGU ';
  275. TAB1.'TITRE' . 2=mot ' Vcalculée ';
  276.  
  277. si (EGA GRAPH 'O' );
  278. dessin (his.'TABD') (his.'2UN') ;
  279. ung1 = vect 0.5 (rv.'INCO'.'UN') ux uy jaune ;
  280. trace ung1 mt ;
  281. DESS evolt 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1 ;
  282. FINSI ;
  283.  
  284. FINPROC ;
  285.  
  286.  
  287. SI COMPLET ;
  288.  
  289. RV.INCO.'NN'=0.8 ;
  290. RV.ITMA= 1500 ;
  291. EXEC RV ;
  292.  
  293. post rv GRAPH ;
  294. evolt=rv.'RET' ;
  295.  
  296. RV.INCO.'NN'=0.6 ;
  297. RV.ITMA= 3500 ;
  298. EXEC RV ;
  299.  
  300. post rv GRAPH ;
  301. evolt=evolt et rv.'RET' ;
  302.  
  303. TAB1=TABLE;
  304. TAB1.'TITRE'=TABLE ;
  305. TAB1.1='MARQ ETOI REGU '     ;
  306. TAB1.'TITRE' . 1=mot ' V_théorique N=0.8'     ;
  307. TAB1.2='TIRR      REGU ';
  308. TAB1.'TITRE' . 2=mot ' Vcalculée N=0.8';
  309. TAB1.3='MARQ PLUS REGU '     ;
  310. TAB1.'TITRE' . 3=mot ' V_théorique N=0.6'     ;
  311. TAB1.4='TIRR      REGU ';
  312. TAB1.'TITRE' . 4=mot ' Vcalculée N=0.6';
  313. DESS evolt 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE  TAB1 ;
  314.  
  315. SINON ;
  316. RV.INCO.'NN'=0.8 ;
  317. RV.ITMA= 40 ;
  318. EXEC RV ;
  319.  
  320. post rv GRAPH ;
  321. ev=rv.'Vcal' ;
  322. vc=extr ev 'ORDO' ;
  323. list vc ;
  324. vcr=prog -.30694 -.30664 -.30531 -.30203 -.29464
  325.          -.27818 -.24186 -.16326 1.88015E-17  ;
  326.  
  327. ER=SOMM( abs (vcr - vc) ) ;
  328. mess ' Ecart sur profil de V : ' ER ;
  329.  si ( er > err1 ) ; erreur 5 ; finsi ;
  330.  
  331. FINSI ;
  332.  
  333. FIN ;
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  

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