Télécharger cc3d3.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : cc3d3.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 P2 + bulle - P1 et Q2 - Q1 algorithme de projection
  12.  
  13. GRAPH = 'N' ;
  14. err1=1.e-6;
  15. DISCR='QUAF';
  16.  
  17. ds1=0.02 ; ds2=0.3 ;
  18. ITMAX=20 ;
  19.  
  20. option dime 2 elem qua4 ;
  21.  
  22. opti isov suli ;
  23. p1= 0 0 ; p12=0.5 0. ; p2= 1 0 ;
  24.  
  25. ab=p1 d dini ds1 dfin ds2 p12 d dini ds2 dfin ds1 p2 ;
  26. ab12= p12 d dini ds1 dfin ds2 (0.5 0.5) d dini ds2
  27. dfin ds1 (0.5 1.) ;
  28. mt0 = ab trans dini ds1 dfin ds2 (0 0.5) trans dini ds2
  29. dfin ds1 (0 0.5) ;
  30.  
  31. mt1= ab trans dini ds1 dfin ds2 (0 0.5);
  32.  
  33. opti elem tri3 ;
  34.  
  35. mt2= (cote 3 mt1) trans dini ds2
  36. dfin ds1 (0 0.5) ;
  37.  
  38. *cnt2=cont mt2 ;
  39. *mt2 = cnt2 surf ;
  40. mt = mt1 et mt2 ;
  41. elim (mt0 et mt) 1.e-5 ;
  42.  
  43. si (ega graph 'O');
  44. trace mt ;
  45. finsi ;
  46.  
  47. bc=cote 2 mt0 ;
  48. cd=cote 3 mt0 ;
  49. da=cote 4 mt0 ;
  50. ct=ab et bc et cd et da ;
  51. elim ct 1.e-3 ;
  52. cnt= ab et bc et cd et da ;
  53.  
  54.  
  55. opti dime 3 elem cub8;
  56.  
  57. prof=0.1;
  58. mth=mt plus (0 0 prof) ;
  59. ab= ab plus (0 0 prof) ;
  60. oeil = 10 10 100 ;
  61. nph=1 ;
  62. cav=mt volu nph mth ;
  63.  
  64. f1=face 1 cav ;
  65. f2=face 2 cav ;
  66. f3=face 3 cav ;
  67.  
  68. psup=(ab plus (0 1 0)) trans nph (0 0 ((-1.)*prof)) ;
  69. psup= CHAN QUAF psup ;
  70. ff3= CHAN QUAF f3 ;
  71. cav= chan cav quaf ;
  72. ab=chan ab quaf ;
  73. f1= chan QUAF f1 ;
  74. f2= chan QUAF f2 ;
  75.  
  76. srti=da plus (0.5 0. 0.) ;
  77. srti = inve srti ;
  78. elim (srti et cav) 1.e-5 ;
  79. srti =chan srti quaf ;
  80.  
  81. elim (cav et psup et ab et ff3 et f1 et f2 et srti) 1.e-5 ;
  82.  
  83. $mt=mode cav 'NAVIER_STOKES' DISCR ;
  84. $AB=mode AB 'NAVIER_STOKES' DISCR ;
  85. DOMA $mt IMPR ;
  86. $srti=mode srti 'NAVIER_STOKES' DISCR ;
  87. DOMA $mt IMPR ;
  88. srti = doma $srti maillage ;
  89.  
  90. DT=1. ;
  91. MU=1. ;
  92. RO= 400.;
  93. kpress='MSOMMET';
  94. prep1=doma $mt kpress;
  95. bcp=elem prep1 POI1 (lect 1 ) ;
  96.  
  97. rv= eqex 'OMEGA' 1. 'NITER' 1 ITMA ITMAX
  98. 'OPTI' 'EF' 'IMPL' KPRESS 'SUPG'
  99. ZONE $mt OPER NS 1. 'UN' (MU/RO) INCO 'UN'
  100. 'OPTI' 'EF' 'BDF2'
  101. ZONE $mt OPER DFDT 1. 'UNM' 'UNMM' DT 'UN' MU INCO 'UN'
  102. ;
  103.  
  104. rv=eqex rv
  105. CLIM
  106. UN UIMP ff3 0. UN VIMP ff3 0. UN WIMP ff3 0.
  107. UN UIMP PSUP 1. UN WIMP (f1 et f2) 0. ;
  108.  
  109. rv=eqex rv CLIM UN WIMP CAV 0. ;
  110.  
  111.  
  112. rv.inco= table inco ;
  113. rv.inco.tn = kcht $mt scal sommet 5. ;
  114. rv.inco.un = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
  115. rv.inco.'UNM' = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
  116. rv.inco.'UNMM' = kcht $mt vect sommet (1.e-5 1.e-5 1.e-5) ;
  117.  
  118. rv.inco.pres = kcht $mt scal KPRESS 1.e-5 ;
  119.  
  120. rv.'METHINV'.TYPINV=3 ;
  121. rv.'METHINV'.IMPINV=0 ;
  122. rv.'METHINV'.NITMAX=100;
  123. rv.'METHINV'.PRECOND=3 ;
  124. rv.'METHINV'.RESID =1.e-8 ;
  125.  
  126. RVP= EQEX
  127. 'OPTI' 'EF' KPRESS
  128. ZONE $mt OPER KBBT (-1.) INCO 'UN' 'PRES'
  129. CLIM PRES TIMP bcp 0.
  130. ;
  131.  
  132. rvp.'METHINV'.TYPINV=2 ;
  133. rvp.'METHINV'.IMPINV=0 ;
  134. rvp.'METHINV'.NITMAX=100;
  135. rvp.'METHINV'.PRECOND=3 ;
  136. rvp.'METHINV'.RESID =1.e-10;
  137. rvp.'METHINV' . 'FCPRECT'=100 ;
  138. rvp.'METHINV' . 'FCPRECI'=100 ;
  139.  
  140. RV.'PROJ'= RVP ;
  141.  
  142. exec rv ;
  143. da = chan da quaf ;
  144. elim ( da et cav) 1.e-5 ;
  145. srti=da plus (0.5 0. 0.) ;
  146. srti = inve srti ;
  147. elim (srti et cav) 1.e-5 ;
  148. evolV = EVOL 'CHPO' (rv.'INCO'.'UN') UX (srti ) ;
  149. evx=extr evolV 'ORDO' ;
  150. list evx ;
  151. evy=extr evolV 'ABSC' ;
  152. evolV= evol 'MANU' 'Vitesse' evx 'Hauteur' evy ;
  153.  
  154. si ('EGA' graph 'O' );
  155. TAB1=TABLE;
  156. TAB1.'TITRE'=TABLE ;
  157. TAB1 . 1 ='MARQ REGU ' ;
  158.  
  159. TAB1.'TITRE' . 1 = mot 'Composante_UY ' ;
  160. DESS evolV 'TITX' 'R (m)' 'TITY' 'V (m/s)' LEGE TAB1 ;
  161.  
  162. c1=vect (rv.inco.'UN') 0.3 ux uy uz jaune ;
  163. *trace c1 mt ;
  164. trace c1 cav ;
  165. trace (rv.inco.'PRESSION') mt cnt ;
  166. finsi ;
  167.  
  168. lrr=prog
  169. -9.45667E-40 -1.18439E-02 -2.33401E-02 -4.26135E-02 -6.08201E-02
  170. -9.03463E-02 -.11898 -.16633 -.18921 -.16629
  171. -8.99933E-02 -1.74638E-02 2.50559E-02 7.02040E-02 .12375
  172. .18613 .31337 .46512 .68187 .83580
  173. 1.0000;
  174.  
  175. lr = evx ;
  176.  
  177. Dgsrti = doma $srti 'XXDIAGSI' ;
  178. li= 'EXTR' (evol chpo Dgsrti srti) 'ORDO' ;
  179.  
  180. ER=( SOMM( ((lr - lrr)*(lr - lrr))*li) )** 0.5 ;
  181. mess ' Ecart sur CENTREP1 TRI7 QUADR ' er ;
  182. Si ( er > err1 ) ; erreur 5 ; finsi ;
  183. FIN ;
  184.  
  185.  
  186.  
  187.  
  188.  
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  

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