Télécharger la plaquette Cast3M


fr en

Télécharger BINGHAMp.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  BINGHAMp.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. **  Fluide de BINGHAM : Ecoulement de POISEUILLE
  5. **  Comparaison à une solution analytique
  6. **  CANAL LONGUEUR 10. LARGEUR 4.
  7. **  Algorithme implicite (dans ce cas beaucoup plus rapide)
  8. **  Auteur : Isabelle Claudel  Décembre 1997
  9.  
  10. GRAPH='N' ;
  11. COMPLET= FAUX ;
  12. ERR1=5.E-2 ;
  13.  
  14.  
  15. option dime  2 elem qua4 ;
  16. opti ISOV suli ;
  17.  
  18.  
  19. *************************************** 
  20. ***PROCEDURE CALCUL DE LA VISCOSITE ***
  21. ***************************************
  22.  
  23. DEBP CALCUL ;
  24. ARGU RX*TABLE ;
  25. iarg=rx.'IARG' ;
  26.  
  27. si( non ( ega iarg 2)) ;
  28. mess 'Procedure CALCUL : nombre d arguments incorrect ' iarg ;
  29. quitter CALCUL ;
  30. finsi ;
  31.  
  32. si ( ega ('TYPE' rx.'ARG1') 'MOT     ') ;
  33. UN=rv.'INCO'.(rx.'ARG1') ;
  34. sinon ;
  35. mess 'Procedure CALCUL : type argument invalide ' ;
  36. quitter CALCUL ;
  37. finsi ;
  38.  
  39. si ( ega ('TYPE' rx.'ARG2') 'MOT     ') ;
  40. NU=rv.'INCO'.(rx.'ARG2') ;
  41. sinon ;
  42. mess 'Procedure CALCUL : type argument invalide ' ;
  43. quitter CALCUL ;
  44. finsi ;
  45.  
  46. *mess ' Proc CALCUL ' ;
  47.  
  48. UN1=  exco ux UN ;
  49. UN2=  exco uy UN ;
  50.  
  51. UNN1=NOMC UN1 'SCAL' ;
  52. UNN2=NOMC UN2 'SCAL' ;
  53.  
  54. unn1= kcht $mt scal sommet UNN1;
  55. unn2= kcht $mt scal sommet UNN2;
  56.  
  57. GU1 = kops  UNN1 'GRAD'  $mt ;
  58. GU2 = kops  UNN2 'GRAD'  $mt ;
  59.  
  60. ddudx= exco ux GU1 ;
  61. ddudy= exco uy GU1 ;
  62. ddvdx= exco ux GU2 ;
  63. ddvdy= exco uy GU2 ;
  64.  
  65. dudx= NOMC ddudx 'SCAL' ;
  66. dudy= NOMC ddudy 'SCAL' ;
  67. dvdx= NOMC ddvdx 'SCAL' ;
  68. dvdy= NOMC ddvdy 'SCAL' ;
  69.  
  70. dudx= kcht $mt scal centre dudx ;
  71. dudy= kcht $mt scal centre dudy ;
  72. dvdx= kcht $mt scal centre dvdx ;
  73. dvdy= kcht $mt scal centre dvdy ;
  74.  
  75.  
  76. D11= dudx ;
  77. D12= (kops (kops dudy '+' dvdx) '/' 2) ;
  78. D21= (D12)  ;
  79. D22= (dvdy) ;
  80.  
  81. *** Proten = Produit tensoriel D:D ***
  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. rv.'INCO'.'ProTen' = ProTen ; 
  87. *list ProTen ;
  88.  
  89.  
  90. *** Norme de D = racine carree (Proten / 2) ***
  91.  
  92. NormD = kops ProTen '/' 2. ;
  93. NormD = kops NormD '**' 0.5 ; 
  94.  
  95.  
  96. **Comportement rhéologique de Bingham **
  97.  
  98. Seuil = kcht $mt scal centre 10. ;
  99. MUB   = kcht $mt scal centre 10. ;
  100.  
  101. MU0 = kcht $mt scal centre 1000. ;
  102. Dseuil = kops Seuil '/' (kops 2. '*' (kops MU0 '-' MUB) ) ;   
  103.  
  104. TESTi = kcht $mt scal centre (NormD MASQUE 'EGINFE' Dseuil)    ;
  105. TESTs = kcht $mt scal centre (NormD MASQUE 'SUPERIEUR' Dseuil) ;
  106. NormD = kops ( kops NormD '*' TESTs ) '+' ( kops Dseuil '*' TESTi) ;
  107.  
  108.  
  109. MU = kops MUB '+' (kops Seuil '/' (kops 2. '*' NormD) ) ;
  110. ro = kcht $mt scal centre 1. ;
  111. NU = kops MU '/' ro ;
  112.  
  113. rv.'INCO'.(rx.'ARG2')=NU ;
  114.  
  115. rv.'INCO'.'GU1'=D12 ;
  116. as2 ama1  = 'KOPS' 'MATRIK' ;
  117. FINPROC as2 ama1 ;
  118.  
  119.  
  120.  
  121.  
  122. ***********************
  123. *******MAILLAGE********
  124. ***********************
  125.  
  126. nbe=16 ; nbv=16 ;
  127. nbee = nbe / 2 ;
  128.  
  129. p1=0. 0.;
  130. p2=4. 0.;
  131. entree= p1 d nbe p2 ;
  132. c1= 5. 0. ;
  133. isa = 2. 0. ;
  134.  
  135. entree2 = isa d nbee p2 ;
  136. ikas=0 ;
  137. *mess 'ikas= 0 --> DROIT (option par defaut)  ikas=1 --> COURBE  ? ';
  138. *obtenir ikas*entier ;
  139.  
  140. si (EGA ikas 1) ;
  141.  
  142. q1=p1 tour c1 -90. ;
  143. q2=p2 tour c1 -90. ;
  144. pp1=p1 c c1 q1 nbv;
  145. pp2=p2 c c1 q2 nbv;
  146. sortie=entree tour c1 -90 ;
  147.  
  148. sinon ;
  149.  
  150. q1=p1 plus (0 10)  ;
  151. q2=p2 plus (0 10)  ;
  152. pp1=p1 d q1 nbv;
  153. pp2=p2 d q2 nbv;
  154. sortie=entree plus (0 10) ;
  155. finsi ;
  156. mientree = entree plus (0. 5.) ;
  157.  
  158. pp1=inve pp1 ;
  159. sortie=inve sortie;
  160. elim 0.0001 (sortie et pp1 et pp2 et entree );
  161. cnt=entree et pp2 et sortie et pp1 ;
  162.  
  163. *mt=surf cnt ;
  164. mt=daller entree pp2 sortie pp1 ;
  165.  
  166. angle=0. ;
  167.  
  168. mt=mt tour p1 angle ;
  169. entree=entree tour p1 angle ;
  170. sortie=sortie tour p1 angle ;
  171. pp1=pp1 tour p1 angle ;
  172. pp2=pp2 tour p1 angle ;
  173.  
  174. elim (pp1 et pp2 et entree et sortie et mt et entree2) 0.0001 ;
  175. *ccc=2. 5. ;
  176.  
  177. *option donn 5 ;
  178. mt = mt 'COUL' bleu ;
  179.  
  180. *trace mt ;
  181.  
  182.  
  183. ************************
  184. *** RESOLUTION ***
  185. ************************
  186. mt0= mt ;
  187. mt=chan mt quaf ;
  188. $mt=mode mt 'NAVIER_STOKES' LINE ;
  189. doma $mt 'IMPR' ;
  190. entre1= chan entree quaf ;
  191. sorti1= chan sortie quaf ;
  192. elim (mt et entre1 et sorti1 et pp1 et pp2 ) 1.e-4 ;
  193. $entree = mode entre1 'NAVIER_STOKES' LINE ;
  194. $sortie = mode sorti1 'NAVIER_STOKES' LINE ;
  195. * doma $entree 'IMPR' ;
  196. * doma $sortie 'IMPR' ;
  197.  
  198. ***dP/dy = - 20***
  199.  to  = kcht $entree vect centre ( 0.  0.) ;
  200.  tos = kcht $sortie vect centre ( 0.  200.) ;
  201. *to  = kcht $entree vect centre ( 0.  6.) ;
  202. *tos = kcht $sortie vect centre ( 0.  0.) ;
  203.  
  204. EPS=1.E-6  ;
  205. DT=1. ;
  206. RV = EQEX $mt  OPTI EF IMPL 'CENTREE'
  207.  ZONE  $mt    'OPER'   CALCUL 'UN' 'NU'
  208.  ZONE  $mt    'OPER'    DUDW   EPS           INCO 'UN'
  209.  OPTI EF IMPL 'CENTREE'
  210. ZONE  $mt    'OPER'    NS 1. 'UN' 'NU'        INCO 'UN'
  211.  ZONE  $entree  'OPER'   'TOIMP' TO           INCO 'UN'
  212.  ZONE  $sortie  'OPER'   'TOIMP' TOS          INCO 'UN'
  213. ;
  214. RV=EQEX RV CLIM
  215.  UN UIMP (pp1 et pp2) 0.
  216.  UN VIMP (pp1 et pp2) 0.
  217. ;
  218. rv.'INCO'=table 'INCO' ;
  219. rv.'INCO'.'UN' = kcht $mt VECT SOMMET  (0. 0.) ;
  220. rv.'INCO'.'NU' = kcht $mt SCAL CENTRE  10.  ;
  221. rv.'INCO'.'ProTen'= kcht $mt SCAL SOMMET  0.  ;
  222.  
  223.  
  224.  
  225.  
  226. DEBPROC POST ;
  227. ARGU RV*TABLE GRAPH*MOT ;
  228. *** POST TRAITEMENT ***
  229.  
  230. evol2v10 = EVOL 'CHPO' (rv.'INCO'.'UN') UY (entree2) ;
  231. *list evol2v ;
  232. XX=extr evol2v10 'ABSC' ;
  233.  
  234. vth10 = prog  
  235.  -2.25
  236.  -2.25
  237.  -2.25
  238.  -2.1875
  239.  -2.
  240.  -1.6875
  241.  -1.25
  242.  -0.6875
  243.  0.
  244.  ;
  245.  
  246. vc=extr evol2v10 'ORDO' ;
  247. ER=SOMM( abs ((vth10 - vc)/2.5) ) ;
  248. mess ' Ecart sur profil de V : ' ER ;
  249.  
  250. evol1v10 = EVOL 'MANU' Rayon xx Vitesse vth10 ;
  251. evolt = evol1v10 et evol2v10 ;
  252. TAB1=TABLE ;
  253. TAB1.'TITRE'=TABLE ;
  254. TAB1.1='MARQ ETOI REGU '      ;
  255. TAB1.'TITRE' . 1=MOT 'V_theorique '      ;
  256. TAB1.2='TIRR      REGU ' ;
  257. TAB1.'TITRE' . 2=MOT 'V_calculee ' ;
  258.  
  259. un=rv.INCO.'UN' ;
  260. ung1=vect un 0.5  ux uy jaune ;
  261.  
  262. si (EGA GRAPH 'O' );
  263. DESS evolt 'TITX' 'X (m)' 'TITY' 'V (m/s)' 'LEGE' TAB1  ;
  264. *list evolt ;
  265. mt0 = doma $mt maillage ;
  266. trace ung1 mt0 ;
  267. nu = elno (rv.'INCO'.'NU') $mt ; 
  268. trace nu mt0 12 ;
  269. FINSI ;
  270.  
  271.  si ( er > err1 ) ; erreur 5 ; finsi ;
  272. FINPROC ;
  273.  
  274. rv.'OMEGA'=0.9;
  275. rv.'NITER'=20  ;
  276.  
  277. EXEC RV ;
  278.  
  279. POST RV GRAPH ;
  280.  
  281.  
  282.  FIN ;
  283.  
  284.  
  285.  
  286.  
  287.  
  288.  
  289.  
  290.  
  291.  
  292.  
  293.  

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