Télécharger la plaquette Cast3M


fr en

Télécharger darcy3EFMH.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  darcy3EFMH.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. ************************** CAS TEST : darcy3.dgibi ******************
  6. *
  7.  GRAPH = 'N'   ;
  8. 'SAUT' 'PAGE' ;
  9. NITER = 40 ;
  10. *
  11. *-------------------------------------------------------------------
  12. * TEST DARCY3
  13. *                  CALCUL DARCY ORTHOTROPE 3D
  14. *
  15. * TEST TRANGEOL en EFMH
  16. * staionnaire comme limite d'un transitoire
  17. *
  18. * On effectue trois calculs sur un cube, maillé par des cubes
  19. * réguliers.
  20. * Les conditions aux limites varient suivant le cas considéré :
  21. * On impose le flux ou la charge sur les cotés du domaine.
  22. *
  23. * La solution analytique en charge est un polynome de degré un,
  24. * la vitesse est constante et la conductivité hydraulique orthotrope.
  25. *             __           __
  26. *            |               |
  27. *            |  1    0    0  |
  28. *        K = |  0   3/4   0  |
  29. *            |  0    0   1/2 |
  30. *            |__           __|
  31. *
  32. *        H(x,y,z) = -45 x -80 y -60z + 200.
  33. *        V(x,y,z) = ( 45 ; 60 ; 30 )
  34. *
  35. * On s'attend à une précision de l'ordre de la précision machine.
  36. *
  37. *-------------------------------------------------------------------
  38. *
  39. 'SAUT' 'PAGE' ;
  40. *
  41. *- Options générales de calcul.
  42. *
  43. 'TITR' 'EFMH DARCY ORTHOTROPE 3D Lineaire : darcy3.dgibi' ;
  44. 'OPTI' 'DIME' 3 'ELEM' 'TET4'                             ; 
  45.  'OPTI' 'ECHO' 1                                           ;
  46. *
  47. *                                                      ------------
  48. *                                                      = MAILLAGE =
  49. *                                                      ------------
  50. OEIL  = 5.D0  6.D0  7.D0 ;
  51. VECX  = 1.D0  0.D0  0.D0 ;
  52. VECY  = 0.D0  1.D0  0.D0 ;
  53. VECZ  = 0.D0  0.D0  1.D0 ;
  54. *
  55. ENX = 6          ;
  56. ENY = 6          ;
  57. ENZ = 6          ;
  58. DX  = 1.D0 / ENX ;
  59. DY  = 1.D0 / ENY ;
  60. DZ  = 1.D0 / ENZ ;
  61. *
  62. *- Création des points
  63. *
  64. A0 = 0.D0  0.D0  0.D0 ;
  65. B0 = 1.D0  0.D0  0.D0 ;
  66. C0 = 1.D0  1.D0  0.D0 ;
  67. D0 = 0.D0  1.D0  0.D0 ;
  68. E0 = 0.D0  0.D0  1.D0 ;
  69. F0 = 1.D0  0.D0  1.D0 ;
  70. G0 = 1.D0  1.D0  1.D0 ;
  71. H0 = 0.D0  1.D0  1.D0 ;
  72. *
  73. *- Création des droites
  74. *
  75. AB = 'DROI' ENX A0 B0 ;
  76. AD = 'DROI' ENY A0 D0 ;
  77. AE = 'DROI' ENZ A0 E0 ;
  78. BC = 'DROI' ENY B0 C0 ;
  79. BF = 'DROI' ENZ B0 F0 ;
  80. CD = 'DROI' ENX C0 D0 ;
  81. CG = 'DROI' ENZ C0 G0 ;
  82. DH = 'DROI' ENZ D0 H0 ;
  83. EF = 'DROI' ENX E0 F0 ;
  84. EH = 'DROI' ENY E0 H0 ;
  85. FG = 'DROI' ENY F0 G0 ;
  86. GH = 'DROI' ENX G0 H0 ;
  87. *
  88. FE = 'INVE' EF ;
  89. EA = 'INVE' AE ;
  90. DC = 'INVE' CD ;
  91. HD = 'INVE' DH ;
  92. DA = 'INVE' AD ;
  93. HE = 'INVE' EH ;
  94. GF = 'INVE' FG ;
  95. FB = 'INVE' BF ;
  96. *
  97. *- Creation des faces du cube
  98. *
  99. SDRO = 'DALL' AB BF FE EA 'PLAN' ;
  100. SGAU = 'DALL' DC CG GH HD 'PLAN' ;
  101. SBAS = 'DALL' AB BC CD DA 'PLAN' ;
  102. SHAU = 'DALL' EF FG GH HE 'PLAN' ;
  103. SDEV = 'DALL' BC CG GF FB 'PLAN' ;
  104. SDER = 'DALL' AD DH HE EA 'PLAN' ;
  105. *
  106. *- Création maillage géométrique
  107. *
  108. ENXM  = ENX + ENY + ENZ              ;
  109. ELI0  = 1.D0 / ENXM / 10.D0          ;
  110. 'SI' ('EGA' ('VALEUR' 'ELEM') 'CUB8') ;
  111.    CUBE1 = 'PAVE' SDER SDEV SBAS SHAU SGAU SDRO ;
  112. 'SINON' ;
  113.    CUBE1 = 'VOLU' (SDER 'ET' SDEV 'ET' SBAS 'ET' SHAU
  114.                    'ET' SGAU 'ET' SDRO) ;
  115. 'FINSI' ;
  116. QFTOT = CHANGE CUBE1 QUAF ;
  117. QFGAU = CHANGE SGAU QUAF ;
  118. QFDRO = CHANGE SDRO QUAF ;
  119. QFHAU = CHANGE SHAU QUAF ;
  120. QFBAS = CHANGE SBAS QUAF ;
  121. QFDEV = CHANGE SDEV QUAF ;
  122. QFDER = CHANGE SDER QUAF ;
  123.  ELIM ELI0 (QFTOT ET QFGAU ET QFDRO ET QFHAU ET QFBAS ET QFDEV ET
  124.               QFDER )  ;
  125. *
  126. *- Création maillage HYBRIDE et sous-objets (conditions aux limites)
  127. *
  128. MODHYB = MODE QFTOT 'DARCY' 'ANISOTROPE'     ;
  129. MODGAU = MODE QFGAU 'DARCY' 'ANISOTROPE'     ;
  130. MODDRO = MODE QFDRO 'DARCY' 'ANISOTROPE'     ;
  131. MODHAU = MODE QFHAU 'DARCY' 'ANISOTROPE'     ;
  132. MODBAS = MODE QFBAS 'DARCY' 'ANISOTROPE'     ;
  133. MODDEV = MODE QFDEV 'DARCY' 'ANISOTROPE'     ;
  134. MODDER = MODE QFDER 'DARCY' 'ANISOTROPE'     ;
  135. C11   = 'DOMA' MODHYB 'VOLUME'       ;
  136. CHYB1 = 'DOMA' MODHYB 'SURFACE'       ;
  137. CHYB2 = 'DOMA' MODHYB 'NORMALE'       ;
  138. CEGAU  = 'DOMA' MODGAU 'CENTRE' ;
  139. CEDRO  = 'DOMA' MODDRO 'CENTRE' ;
  140. CEHAU  = 'DOMA' MODHAU 'CENTRE' ;
  141. CEBAS  = 'DOMA' MODBAS 'CENTRE' ;
  142. CEDEV  = 'DOMA' MODDEV 'CENTRE' ;
  143. CEDER  = 'DOMA' MODDER 'CENTRE' ;
  144. *
  145. *- Solution analytique
  146. *
  147. XX  YY  ZZ  = 'COOR' (DOMA MODHYB 'FACE' )  ;
  148. XXC YYC ZZC = 'COOR' (DOMA MODHYB 'CENTRE') ;
  149. *
  150. VKX = 1.D0             ;
  151. VKY = 0.75D0           ;
  152. VKZ = 0.5D0            ;
  153. VVKX = MANU 'CHPO'  (DOMA MODHYB 'CENTRE') 'K11' 1.D0 ;
  154. VVKX = CHANGER ATTRIBUT VVKX  'NATU' DISCRET  ;
  155. VVKY = MANU 'CHPO'  (DOMA MODHYB 'CENTRE') 'K22' 0.75D0 ;
  156. VVKY = CHANGER ATTRIBUT VVKY  'NATU' DISCRET  ;
  157. VVKZ = MANU 'CHPO'  (DOMA MODHYB 'CENTRE') 'K33' 0.5D0 ;
  158. VVKZ = CHANGER ATTRIBUT VVKZ  'NATU' DISCRET  ;
  159. AA  = -45.D0           ;
  160. BB  = -80.D0           ;
  161. CC  = -60.D0           ;
  162. DD  = 200.D0           ;
  163. AAA = -1.D0 * VKX * AA ;
  164. BBB = -1.D0 * VKY * BB ;
  165. CCC = -1.D0 * VKZ * CC ;
  166. *
  167. PANAF = (AA * XX)  + (BB * YY)  + (CC * ZZ)  + DD                 ;
  168. PANAC = (AA * XXC) + (BB * YYC) + (CC * ZZC) + DD                 ;
  169. VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 3 'VX' AAA 
  170.          'VY' BBB 'VZ' CCC ;
  171. VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE')   3 'VX' AAA 
  172.           'VY' BBB 'VZ' CCC ;
  173. *
  174. *                                                  --------------
  175. *                                                  = RESOLUTION =
  176. *                                                  --------------
  177. *
  178. MATI3  = (NOMC 'K11' VVKX) et
  179.          (NOMC 'K22' VVKY) et
  180.          (NOMC 'K33' VVKZ) et
  181.          (NOMC 'K21' (0.D0 * VVKX)) et
  182.          (NOMC 'K31' (0.D0 * VVKX)) et
  183.          (NOMC 'K32' (0.D0 * VVKX))           ;
  184. *                                                        ;
  185. *- Conditions aux limites
  186. *
  187. BBGAU = 'BLOQ' CEGAU 'TH' ;
  188. BBDRO = 'BLOQ' CEDRO 'TH' ;
  189. BBHAU = 'BLOQ' CEHAU 'TH' ;
  190. BBBAS = 'BLOQ' CEBAS 'TH' ;
  191. BBDEV = 'BLOQ' CEDEV 'TH' ;
  192. BBDER = 'BLOQ' CEDER 'TH' ;
  193. *
  194. *- TH imposée
  195. *
  196. TTIMP = 'REDU' PANAF CEGAU ;
  197. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  198. EEGAU = 'DEPI' BBGAU TTIM2          ;
  199. TTIMP = 'REDU' PANAF CEDRO ;
  200. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  201. EEDRO = 'DEPI' BBDRO TTIM2          ;
  202. TTIMP = 'REDU' PANAF CEBAS ;
  203. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  204. EEBAS = 'DEPI' BBBAS TTIM2          ;
  205. TTIMP = 'REDU' PANAF CEHAU ;
  206. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  207. EEHAU = 'DEPI' BBHAU TTIM2          ;
  208. TTIMP = 'REDU' PANAF CEDEV ;
  209. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  210. EEDEV = 'DEPI' BBDEV TTIM2          ;
  211. TTIMP = 'REDU' PANAF CEDER ;
  212. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  213. EEDER = 'DEPI' BBDER TTIM2          ;
  214. *
  215. *- Flux imposé
  216. *
  217. FLDRO = -60.D0 * (nomc 'FLUX' (doma moddro VOLUME));
  218. FLGAU  = 60.D0  * (nomc 'FLUX' (doma modgau VOLUME));
  219. FLHAU  =  30.D0 * (nomc 'FLUX' (doma modhau VOLUME));
  220. FLBAS  = -30.D0  * (nomc 'FLUX' (doma modbas VOLUME));
  221. FLDEV  =  45.D0 * (nomc 'FLUX' (doma moddev VOLUME)); 
  222. FLDER  = -45.D0  * (nomc 'FLUX' (doma modder VOLUME));
  223. *
  224. *
  225. *- Assemblage et résolution en TH
  226. *
  227.  
  228.  
  229. h_lim = 'RESOUD' (BBGAU 'ET' BBDEV)
  230.                  (EEGAU 'ET' EEDEV) ;
  231. h_lim = 'EXCO' 'TH' h_lim;
  232. CHCLIM = TABLE;
  233. CHCLIM . 'NEUMANN' = ('NOMC' 'I35' ( FLHAU 'ET' FLBAS
  234.                                  'ET' FLDRO 'ET' FLDER));
  235. CHCLIM . 'DIRICHLET' = 'NOMC' 'I35' h_lim;
  236.  
  237.  
  238. GEOL1 = TABLE;
  239. GEOL1 . 'CONCENTRATION'     = 'NOMC' 'I35' (0.D0 * PANAC)  ;
  240. GEOL1 . 'LUMP'              = FAUX                  ;
  241. GEOL1 . 'TYPDISCRETISATION' = 'EFMH'                ;
  242. GEOL1 . 'THETA_DIFFUSION'   = 1.0D0                 ;
  243. GEOL1 . 'THETA_CONVECTION'  = 1.0D0                 ;
  244. GEOL1 . 'DECENTREMENT'      = FAUX                  ;
  245. GEOL1 . 'DELTAT'            = 10.D0 / NITER         ;
  246. GEOL1 . 'DIFFUSIVITE'       = mati3                 ;
  247. GEOL1 . 'SOLVEUR'           = 3                     ;
  248. GEOL1 . 'PRECONDITIONNEUR'  = 3                     ;
  249. GEOL1 . 'POROSITE'          = 'MANU' 'CHPO' 
  250.                 (doma MODHYB CENTRE) 'CK' 1.D0 ;
  251. GEOL1 . 'CLIMITES'          = CHCLIM                ;
  252. GEOL1 . 'RECALCUL'          = VRAI                  ;
  253. *
  254. GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1;
  255.  
  256.  
  257.  
  258. * boucle transitoire
  259.  
  260.  
  261.  
  262. REPE blocc NITER   ;
  263.    mess 'boucle' &blocc;
  264.    GEOL1 . 'RECALCUL'          = FAUX                  ;
  265.    GEOL1 GEOL2 = TRANGEOL Modhyb GEOL1 GEOL2;
  266. FIN blocc;
  267.  
  268.  
  269.  
  270.  
  271. CHTER3 = 'NOMC' 'TH' GEOL2 . 'TRACE_CONC';
  272. PCEN3 = 'NOMC' 'H' GEOL1 . 'CONCENTRATION' ;
  273. QFACE3 = 'NOMC' 'FLUX' GEOL1 . 'FLUXDIFF' ;
  274.  
  275. *
  276. *- Calcul de V
  277. *
  278. VCENT3 = 'HVIT' MODHYB  QFACE3         ;
  279. QFACE3 = 'EXCO' QFACE3 'FLUX' 'SCAL'             ;
  280. VFACE3 = QFACE3 * CHYB2 / CHYB1                  ;
  281. *
  282. *                                                   -----------------
  283. *                                                   = Calcul ERREUR =
  284. *                                                   -----------------
  285. * ERReur relative en Trace de charge TH aux faces des éléments
  286. * ERReur relative en charge H au centre des éléments
  287. * Erreur relative sur la vitesse au centre des éléments
  288. *
  289. ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
  290. ERRTP3 = ERRTP3 - PANAF / PANAF    ;
  291. ERRTP3 = 'ABS' ERRTP3              ;
  292. ERRP3  = 'EXCO' PCEN3 'H' 'SCAL'   ;
  293. ERRP3  = ERRP3 - PANAC / PANAC     ;
  294. ERRP3  = 'ABS' ERRP3               ;
  295. *
  296. MOT1 = 'MOTS' 'VX' 'VY' 'VZ'        ;
  297. VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
  298. *
  299. VD3  = VANAC - VCENT3            ;
  300. VC3  = 'PSCA' VD3 VD3 MOT1 MOT1  ;
  301. SDC3 = 'ABS' ( VC3 / VDVD )      ;
  302. SDC3 = SDC3 '**' 0.5             ;
  303. *
  304. *                                                 -------------------
  305. *                                                 = Tracé resultats =
  306. *                                                 -------------------
  307. 'SI' ('NEG' GRAPH 'N') ;
  308. *
  309. *- Transformation des quantités aux centres en MCHAML constant.
  310. *
  311. ERRP3 = 'KCHA' MODHYB 'CHAM' ERRP3 ;
  312. SDS3  = 'KCHA' MODHYB 'CHAM' SDC3  ;
  313. *
  314. *   L'erreur relative sur la charge au centre
  315. *   L'erreur relative sur la Vitesse au centre
  316. *
  317. 'TITR' 'darcy3/3 : Erreur relative sur la charge'  ;
  318. 'TRAC' MODHYB ERRP3                                ;
  319. 'TITR' 'darcy3/3 : Erreur relative sur la vitesse' ;
  320. 'TRAC' MODHYB SDS3                                 ;
  321. *
  322. 'FINSI' ;
  323. *
  324. *                                                 -------------------
  325. *                                                 = Gestion ERREURS =
  326. *                                                 -------------------
  327. MAXTP3 = 'MAXI' ERRTP3 ;
  328. MAXP3  = 'MAXI' ERRP3  ;
  329. MAXV3  = 'MAXI' SDC3   ;
  330. *
  331. 'SAUT' 'PAGE'                                                       ;
  332. 'SAUT' 2 'LIGNE'                                                    ;
  333. 'MESS' '                             ERREURS RELATIVES            ' ;
  334. 'SAUT' 1 'LIGNE'                                                    ;
  335. 'MESS' ' cas test          TH                  H                 V' ;
  336. 'SAUT' 1 'LIGNE'                                                    ;
  337. 'MESS' ' numero3     ' maxtp3  '   ' maxp3 '   ' maxv3              ;
  338. 'SAUT' 2 'LIGNE'                                                    ;
  339. *
  340. EPS0 = 1.D-7 ;
  341. LOG3 = MAXTP3 > EPS0 ;
  342. LOG6 = MAXP3  > EPS0 ;
  343. LOG9 = MAXV3  > EPS0 ;
  344. LTP0 = LOG3 ;
  345. LP0  = LOG6 ;
  346. LV0  = LOG9 ;
  347. L0   = LTP0 'OU' LP0  'OU' LV0  ;
  348. 'SI' ( L0 ) ;
  349.    'ERRE' 5 ;
  350. 'SINO'      ;
  351.    'ERRE' 0 ;
  352. 'FINSI'     ;
  353. *
  354. 'FIN' ;
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  

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