Télécharger cc3d2.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : cc3d2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. **
  5. ** --- 08 Novembre 1999 ---
  6. **
  7. ** TEST CAVITE CUBIQUE
  8. **
  9. ** teste KBBT NS en 3D + le Bi CG
  10. ** teste les elements a pression continue sur des hexahedres
  11. ** et prisme P1-P1 et Q1-Q1 algorithme de projection
  12. ** ces elements ne verifient pas LBB mais si DT pas trop petit
  13. ** ca marche quand meme (pression en damier)
  14.  
  15. GRAPH = 'N' ;
  16. DISCR='LINE' ;
  17.  
  18. ds1=0.02 ; ds2=0.1 ;
  19. ITMAX=20 ;
  20. err1=1.e-6;
  21.  
  22. option dime 2 elem qua4 ;
  23. *option dime 2 elem tri3 ;
  24.  
  25. opti isov suli ;
  26. p1= 0 0 ; p12=0.5 0. ; p2= 1 0 ;
  27.  
  28. ab=p1 d dini ds1 dfin ds2 p12 d dini ds2 dfin ds1 p2 ;
  29. ab12= p12 d dini ds1 dfin ds2 (0.5 0.5) d dini ds2
  30. dfin ds1 (0.5 1.) ;
  31. mt0 = ab trans dini ds1 dfin ds2 (0 0.5) trans dini ds2
  32. dfin ds1 (0 0.5) ;
  33.  
  34. mt1= ab trans dini ds1 dfin ds2 (0 0.5);
  35.  
  36. opti elem tri3 ;
  37.  
  38. mt2= (cote 3 mt1) trans dini ds2
  39. dfin ds1 (0 0.5) ;
  40.  
  41. *cnt2=cont mt2 ;
  42. *mt2 = cnt2 surf ;
  43. mt = mt1 et mt2 ;
  44.  
  45. elim (mt0 et mt) 1.e-5 ;
  46.  
  47. si ('EGA' graph 'O' );
  48. trace mt ;
  49. finsi ;
  50.  
  51. bc=cote 2 mt0 ;
  52. cd=cote 3 mt0 ;
  53. da=cote 4 mt0 ;
  54. ct=ab et bc et cd et da ;
  55. elim ct 1.e-3 ;
  56. cnt= ab et bc et cd et da ;
  57. * mt=ab bc cd da daller ;
  58. * mt=cnt surf ;
  59.  
  60. mt=orie mt ;
  61.  
  62. opti dime 3 elem cub8;
  63.  
  64. prof=0.1 ;
  65. mth=mt plus (0 0 prof) ;
  66. ab= ab plus (0 0 prof) ;
  67. oeil = 10 10 100 ;
  68. nph=1 ;
  69. cav=mt volu nph mth ;
  70.  
  71. f1=face 1 cav ;
  72. f2=face 2 cav ;
  73. f3=face 3 cav ;
  74.  
  75. psup=(ab plus (0 1 0)) trans nph (0 0 ((-1.)*prof)) ;
  76. *psup=(CHAN TRI3 psup) ;
  77. *ff3= (CHAN TRI3 f3) ;
  78. ff3 = f3 ;
  79.  
  80. ELIM (PSUP ET FF3) 1.e-4 ;
  81.  
  82.  
  83. cav= chan cav quaf ;
  84. ab=chan ab quaf ;
  85. srti=da plus (0.5 0. 0.) ;
  86. srti = inve srti ;
  87. elim (srti et cav) 1.e-5 ;
  88. srti =chan srti quaf ;
  89.  
  90.  
  91. elim (cav et psup et ab et ff3 et srti) 1.e-5 ;
  92.  
  93. $mt=mode cav 'NAVIER_STOKES' DISCR ;
  94. $AB=mode AB 'NAVIER_STOKES' DISCR ;
  95. $srti=mode srti 'NAVIER_STOKES' DISCR ;
  96. DOMA $mt IMPR ;
  97. srti = doma $srti maillage ;
  98.  
  99. DT=1.;
  100. MU=1. ;
  101. RO= 400. ;
  102. kpress='MSOMMET';
  103. prep1=doma $mt kpress;
  104. bcp=elem prep1 POI1 (lect 1 ) ;
  105.  
  106. rv= eqex 'OMEGA' 1. 'NITER' 1 ITMA ITMAX
  107. 'OPTI' 'EF' 'IMPL' KPRESS 'SUPG' 'NODIV'
  108. ZONE $mt OPER NS 1. 'UN' (MU/RO) INCO 'UN'
  109. * 'OPTI' 'EF' 'BDF2'
  110. * ZONE $mt OPER DFDT 1. 'UNM' 'UNMM' DT 'UN' MU INCO 'UN'
  111. ZONE $mt OPER DFDT 1. 'UNM' DT 'UN' MU INCO 'UN'
  112. ;
  113.  
  114. rv=eqex rv
  115. CLIM
  116. UN UIMP ff3 0. UN VIMP ff3 0. UN WIMP ff3 0.
  117. UN UIMP PSUP 1. UN WIMP (f2 et f1) 0.
  118. ;
  119.  
  120.  
  121. rv.inco= table inco ;
  122. rv.inco.tn = kcht $mt scal sommet 5. ;
  123. rv.inco.un = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
  124. rv.inco.'UNM' = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
  125. rv.inco.'UNMM' = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
  126.  
  127. rv.inco.pres = kcht $mt scal KPRESS 1.e-5 ;
  128.  
  129. rv.'METHINV'.TYPINV=3 ;
  130. rv.'METHINV'.IMPINV=0 ;
  131. rv.'METHINV'.NITMAX=100;
  132. rv.'METHINV'.PRECOND=3 ;
  133. rv.'METHINV'.RESID =1.e-8 ;
  134.  
  135. RVP= EQEX
  136. 'OPTI' 'EF' KPRESS
  137. ZONE $mt OPER KBBT (-1.) INCO 'UN' 'PRES'
  138. CLIM PRES TIMP bcp 0.
  139. ;
  140.  
  141. rvp.'METHINV'.TYPINV=2 ;
  142. rvp.'METHINV'.IMPINV=0 ;
  143. rvp.'METHINV'.NITMAX=100;
  144. rvp.'METHINV'.PRECOND=3 ;
  145. rvp.'METHINV'.RESID =1.e-10;
  146. rvp.'METHINV' . 'FCPRECT'=100 ;
  147. rvp.'METHINV' . 'FCPRECI'=100 ;
  148.  
  149. RV.'PROJ'= RVP ;
  150.  
  151. exec rv ;
  152.  
  153. evolV = EVOL 'CHPO' (rv.'INCO'.'UN') UX (srti ) ;
  154. evx=extr evolV 'ORDO' ;
  155. list evx ;
  156. evy=extr evolV 'ABSC' ;
  157. evolV= evol 'MANU' 'Vitesse' evx 'Hauteur' evy ;
  158.  
  159. si ('EGA' graph 'O' );
  160. TAB1=TABLE;
  161. TAB1.'TITRE'=TABLE ;
  162. TAB1 . 1 ='MARQ REGU ' ;
  163.  
  164. TAB1.'TITRE' . 1 = mot 'Composante_UY ' ;
  165. DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
  166.  
  167. c1=vect (rv.inco.'UN') 0.3 ux uy uz jaune ;
  168. menvf= doma $mt 'ENVELOPPE' ;
  169. trace c1 menvf ;
  170. trace (rv.inco.'PRESSION') menvf cnt ;
  171. finsi ;
  172.  
  173. list evx ;
  174. lrr=prog
  175. -2.85788E-38 -1.61082E-02 -3.22813E-02 -4.99973E-02 -7.04952E-02
  176. -9.56495E-02 -.12703 -.16245 -.18934 -.18725
  177. -.12259 -4.13788E-02 3.10240E-02 8.66209E-02 .12808
  178. .16266 .21239 .31615 .50353 .75483
  179. 1.0000;
  180.  
  181.  
  182. lr = evx ;
  183.  
  184. Dgsrti = doma $srti 'XXDIAGSI' ;
  185. li= 'EXTR' (evol chpo Dgsrti srti) 'ORDO' ;
  186.  
  187. ER=( SOMM( ((lr - lrr)*(lr - lrr))*li) )** 0.5 ;
  188. mess ' Ecart sur MSOMMET PRI6 CUB8 ' er ;
  189. Si ( er > err1 ) ; erreur 5 ; finsi ;
  190. FIN ;
  191.  
  192.  
  193.  
  194.  
  195.  
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202.  

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