Télécharger coudep.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : coudep.dgibi
  2. * Ecoulement dans un Coude
  3. * Test 3D pression continue methode de projection
  4. * teste en 3D VNIMP FPU NS
  5. * Ce test pose des problemes en QUAF ??
  6. * Il faut un minimum de mailles pour converger a fort Reynolds
  7. *
  8. opti dime 2 elem QUA8 ;
  9.  
  10. GRAPH=VRAI ;
  11. GRAPH=FAUX ;
  12. COMPLET = VRAI ;
  13. COMPLET = FAUX ;
  14.  
  15. KPRES='MSOMMET' ;
  16.  
  17. Si COMPLET ;
  18. R=0.5 ; nl = 4 ;n1= 5 ; lgd1=1.; lgd2=2.; n2=5 ;
  19. NU=1.e-5 ;
  20. NUT=5.e-3;
  21. YP=5.e-3 ;
  22. DT=1. ;
  23. Sinon ;
  24. R=0.5 ; nl = 2 ;n1= 2 ; lgd1=1.; lgd2=2.; n2=2 ;
  25. NU=1.e-5 ;
  26. NUT=5.e-2;
  27. YP=5.e-2 ;
  28. DT=1. ;
  29. Finsi ;
  30.  
  31. P0=0 0 ; p1= R 0 ;
  32. p2 = p1 tour 90 p0 ;
  33. p3 = p2 tour 90 p0 ;
  34. p4 = p3 tour 90 p0 ;
  35. mct= p1 c p0 n1 p2 c p0 n1 p3 c p0 n1 p4 c p0 n1 p1;
  36. ct1= p1 c p0 n1 p2;
  37. ct2= p2 c p0 n1 p3;
  38. ct3= p3 c p0 n1 p4;
  39. ct4= p4 c p0 n1 p1;
  40.  
  41. *ct2= p2 d n1 p3;
  42. *ct4= p4 d n1 p1;
  43.  
  44. *mct= p1 c p0 n1 p2 d n1 p0 d n1 p1;
  45. *mc= surf mct ;
  46. mc= daller ct1 ct2 ct3 ct4 ;
  47. *trace mc ;
  48.  
  49. opti dime 3 elem seg2 ;
  50.  
  51. q1= p0 plus (0 0 lgd1);
  52. ax1 = q1 plus (1 0 0) ;
  53. ax2 = q1 plus (1 1 0) ;
  54. q2= q1 tour 90. ax1 ax2 ;
  55. q3= q2 plus (lgd2 0 0);
  56.  
  57. lg= p0 d nl q1 c ax1 n2 q2 d nl q3 ;
  58. *lg= p0 d nl q1 ;
  59.  
  60. opti elem CU20 ;
  61.  
  62. *trace lg ;
  63. lg1= chan lg 'POI1' ;
  64. nbp= (nbno lg1) - 2 ;
  65. pi1= p0 ;
  66. si1= mc ;
  67. pi2= poin 2 lg1 ;
  68. vh1=pi2 moins pi1 ;
  69. si2= si1 plus vh1 ;
  70. Entree = si1 ;
  71. vc=si1 volu 2 si2;
  72. *trace cache vc ;
  73.  
  74. repeter bloc1 nbp ;
  75. pi1=pi2;
  76. si1=si2 ;
  77. pi2= poin (&bloc1 + 2) lg1 ;
  78. vh2=pi2 moins pi1 ;
  79. vhr=vh2 ;
  80.  
  81. axe =vh1 pvec vhr ;
  82. alfa= atg (norm(axe / (vh1 psca vhr) ));
  83. Si ( alfa >EG 1.) ;
  84. si2= si1 tour alfa pi1 (pi1 plus axe) ;
  85. Sinon ;
  86. si2= si1 ;
  87. Finsi ;
  88. si2=si2 plus vh2 ;
  89. vh1=vh2 ;
  90. vc = vc et (si1 volu 2 si2);
  91. fin bloc1 ;
  92. Sortie = si2 ;
  93. oeil = 0. -1.e5 -100. ;
  94. *trace cache vc oeil ;
  95. Sortie= Sortie coul rouge ;
  96. Entree= Entree coul rouge ;
  97.  
  98. Mvc = chan vc QUAF ;
  99. Mentree = chan entree QUAF ;
  100. Msortie = chan sortie QUAF ;
  101.  
  102. ELIM (Mvc et Mentree et Msortie) 1.e-5 ;
  103.  
  104. DEBPROC TEST DISCR*MOT ;
  105.  
  106. $vc = mode Mvc 'NAVIER_STOKES' DISCR ;
  107. $sortie = mode Msortie 'NAVIER_STOKES' DISCR ;
  108. sortip=doma $sortie KPRES;
  109. $entree = mode Mentree 'NAVIER_STOKES' DISCR ;
  110. entree = doma $entree maillage ;
  111. Menvlp = doma $vc 'ENVELOPPE' ;
  112. Mparoi = diff Menvlp (Mentree et Msortie) ;
  113.  
  114. doma $vc 'IMPR' ;
  115.  
  116. * On extrait le champ des normales de la surface d'entree
  117. * en prenant garde a l'orientation des normales
  118. nentr= doma $entree 'NORMALEV';
  119.  
  120. Mcte = doma $entree 'ENVELOPPE';
  121. $cte= mode Mcte 'NAVIER_STOKES' DISCR ;
  122. cte=doma $cte maillage ;
  123.  
  124. * On extrait les composantes nx ny etc
  125. nx= exco 'UX' nentr;
  126. ny= exco 'UY' nentr;
  127. nz= exco 'UZ' nentr;
  128. * on cree les componsantes cartesiennes du champ de vitesse
  129. * qui seront imposees en condition limite type Dirichlet
  130. U0=1. ;
  131. usx=u0*nx ;
  132. usy=u0*ny ;
  133. usz=u0*nz ;
  134.  
  135. $envlp = mode Menvlp 'NAVIER_STOKES' DISCR ;
  136. envlp=doma $envlp maillage ;
  137. $paroi = mode Mparoi 'NAVIER_STOKES' DISCR ;
  138. paroi=doma $paroi maillage ;
  139.  
  140.  
  141. rv= eqex 'ITMA' 5 'NITER' 1 'OMEGA' 1.
  142. OPTI EF 'SUPG' 'IMPL' KPRES
  143. 'ZONE' $paroi 'OPER' 'FPU' 1. 'UN' NU 'UET' YP 'INCO' 'UN'
  144. 'ZONE' $vc 'OPER' 'NS' 1. 'UN' NUT 'INCO' 'UN'
  145. OPTI EF 'CENTREE'
  146. 'ZONE' $vc 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  147. 'CLIM'
  148. 'UN' UIMP entree usx 'UN' VIMP entree usy 'UN' WIMP entree usz
  149. * 'UN' UIMP paroi 0. 'UN' VIMP paroi 0. 'UN' WIMP paroi 0.
  150. ;
  151.  
  152. rv.'METHINV'.TYPINV=3 ;
  153. rv.'METHINV'.IMPINV=0 ;
  154. rv.'METHINV'.NITMAX=100;
  155. rv.'METHINV'.PRECOND=3 ;
  156. rv.'METHINV'.RESID =1.e-10 ;
  157.  
  158. rv.'INCO'= TABLE 'INCO' ;
  159. rv.'INCO'.'UN'=kcht $vc vect sommet (0. 0. 0.) ;
  160. rv.'INCO'.'PRES'=kcht $vc scal KPRES 0. ;
  161.  
  162. RVP= EQEX 'OPTI' 'EF' KPRES
  163. 'ZONE' $vc 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
  164. 'ZONE' $paroi 'OPER' 'VNIMP' $vc 0. 'INCO' 'UN' 'PRES'
  165. 'CLIM' 'PRES' TIMP sortip 0.
  166. ;
  167.  
  168. rvp.'METHINV'.TYPINV=3 ;
  169. rvp.'METHINV'.IMPINV=0 ;
  170. rvp.'METHINV'.NITMAX=100;
  171. rvp.'METHINV'.PRECOND=3 ;
  172. rvp.'METHINV'.RESID =1.e-10 ;
  173.  
  174.  
  175. RV.'PROJ'= RVP ;
  176.  
  177. exec rv ;
  178.  
  179. un=rv.inco.'UN' ;
  180. un1=redu un envlp ;
  181. ung=vect un1 0.1 ux uy uz jaune ;
  182.  
  183. Si GRAPH ;
  184. trace ung envlp oeil ;
  185. Finsi;
  186.  
  187. q= dbit un $envlp ;
  188. list q ;
  189. qe=dbit un $entree ;
  190. qs=dbit un $sortie ;
  191. qp=dbit un $paroi ;
  192. mess 'qe=' qe 'qs=' qs 'qp=' qp 'qs+qp=' (qs+qp) ;
  193.  
  194.  
  195. yplus=yp*(rv.inco.'UET')/nu ;
  196. mess 'Mini yplus ' (mini yplus) ' Maxi yplus ' (maxi yplus) ;
  197.  
  198. pn=exco rv.inco.'PRESSION' 'PRES' ;
  199. mtp= doma $vc 'MMAIL' ;
  200.  
  201. Si GRAPH ;
  202. trace pn mtp ;
  203. Finsi;
  204.  
  205. mess 'Mini Pres ' (mini pn) ' Maxi Pres ' (maxi pn) ;
  206.  
  207. FINPROC ;
  208.  
  209. *test 'QUAF' ;
  210.  
  211. test 'MACRO';
  212.  
  213. test 'LINE' ;
  214.  
  215. FIN ;
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  
  224.  

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