Télécharger la plaquette Cast3M


fr en

Télécharger darcy3VF.dgibi

Retour à la liste

Numérotation des lignes : 

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

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