Télécharger la plaquette Cast3M


fr en

Télécharger darcy3.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  darcy3.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. ************************** CAS TEST : darcy3.dgibi ******************
  6. *
  7.  GRAPH = 'N'   ;
  8. 'SAUT' 'PAGE' ;
  9. *
  10. *-------------------------------------------------------------------
  11. * TEST DARCY3
  12. *                  CALCUL DARCY ORTHOTROPE 3D
  13. *    Résolution par une méthode d'éléments finis mixtes hybrides.
  14. *
  15. * Ce test permet de vérifier le bon fonctionnement des opérateurs
  16. * utilisé afin de résoudre les équations de DARCY par une méthode
  17. * d'éléments finis mixtes hybrides dans CASTEM2000.
  18. *
  19. * On effectue trois calculs sur un cube, maillé par des cubes
  20. * réguliers.
  21. * Les conditions aux limites varient suivant le cas considéré :
  22. * On impose le flux ou la charge sur les cotés du domaine.
  23. *
  24. * La solution analytique en charge est un polynome de degré un,
  25. * la vitesse est constante et la conductivité hydraulique orthotrope.
  26. *             __           __
  27. *            |               |
  28. *            |  1    0    0  |
  29. *        K = |  0   3/4   0  |
  30. *            |  0    0   1/2 |
  31. *            |__           __|
  32. *
  33. *        H(x,y,z) = -45 x -80 y -60z + 200.
  34. *        V(x,y,z) = ( 45 ; 60 ; 30 )
  35. *
  36. * On s'attend à une précision de l'ordre de la précision machine.
  37. *
  38. *-------------------------------------------------------------------
  39. *
  40. 'SAUT' 'PAGE' ;
  41. *
  42. *- Options générales de calcul.
  43. *
  44. 'TITR' 'EFMH DARCY ORTHOTROPE 3D Lineaire : darcy3.dgibi' ;
  45. 'OPTI' 'DIME' 3 'ELEM' 'CUB8'                             ;
  46. *'OPTI' 'ECHO' 0                                           ;
  47. *
  48. *                                                      ------------
  49. *                                                      = MAILLAGE =
  50. *                                                      ------------
  51. OEIL  = 5.D0  6.D0  7.D0 ;
  52. VECX  = 1.D0  0.D0  0.D0 ;
  53. VECY  = 0.D0  1.D0  0.D0 ;
  54. VECZ  = 0.D0  0.D0  1.D0 ;
  55. *
  56. ENX = 5          ;
  57. ENY = 10         ;
  58. ENZ = 5          ;
  59. DX  = 1.D0 / ENX ;
  60. DY  = 1.D0 / ENY ;
  61. DZ  = 1.D0 / ENZ ;
  62. *
  63. *- Création des points
  64. *
  65. A0 = 0.D0  0.D0  0.D0 ;
  66. B0 = 1.D0  0.D0  0.D0 ;
  67. C0 = 1.D0  1.D0  0.D0 ;
  68. D0 = 0.D0  1.D0  0.D0 ;
  69. E0 = 0.D0  0.D0  1.D0 ;
  70. F0 = 1.D0  0.D0  1.D0 ;
  71. G0 = 1.D0  1.D0  1.D0 ;
  72. H0 = 0.D0  1.D0  1.D0 ;
  73. *
  74. *- Création des droites
  75. *
  76. AB = 'DROI' ENX A0 B0 ;
  77. AD = 'DROI' ENY A0 D0 ;
  78. AE = 'DROI' ENZ A0 E0 ;
  79. BC = 'DROI' ENY B0 C0 ;
  80. BF = 'DROI' ENZ B0 F0 ;
  81. CD = 'DROI' ENX C0 D0 ;
  82. CG = 'DROI' ENZ C0 G0 ;
  83. DH = 'DROI' ENZ D0 H0 ;
  84. EF = 'DROI' ENX E0 F0 ;
  85. EH = 'DROI' ENY E0 H0 ;
  86. FG = 'DROI' ENY F0 G0 ;
  87. GH = 'DROI' ENX G0 H0 ;
  88. *
  89. FE = 'INVE' EF ;
  90. EA = 'INVE' AE ;
  91. DC = 'INVE' CD ;
  92. HD = 'INVE' DH ;
  93. DA = 'INVE' AD ;
  94. HE = 'INVE' EH ;
  95. GF = 'INVE' FG ;
  96. FB = 'INVE' BF ;
  97. *
  98. *- Creation des faces du cube
  99. *
  100. SDRO = 'DALL' AB BF FE EA 'PLAN' ;
  101. SGAU = 'DALL' DC CG GH HD 'PLAN' ;
  102. SBAS = 'DALL' AB BC CD DA 'PLAN' ;
  103. SHAU = 'DALL' EF FG GH HE 'PLAN' ;
  104. SDEV = 'DALL' BC CG GF FB 'PLAN' ;
  105. SDER = 'DALL' AD DH HE EA 'PLAN' ;
  106. *
  107. *- Création maillage géométrique
  108. *
  109. ENXM  = ENX + ENY + ENZ              ;
  110. ELI0  = 1.D0 / ENXM / 10.D0          ;
  111. CUBE1 = 'PAVE' SDER SDEV SBAS SHAU SGAU SDRO ;
  112. QFTOT = CHANGE CUBE1 QUAF ;
  113. QFGAU = CHANGE SGAU QUAF ;
  114. QFDRO = CHANGE SDRO QUAF ;
  115. QFHAU = CHANGE SHAU QUAF ;
  116. QFBAS = CHANGE SBAS QUAF ;
  117. QFDEV = CHANGE SDEV QUAF ;
  118. QFDER = CHANGE SDER QUAF ;
  119.  ELIM ELI0 (QFTOT ET QFGAU ET QFDRO ET QFHAU ET QFBAS ET QFDEV ET
  120.               QFDER )  ;
  121. *
  122. *- Création maillage HYBRIDE et sous-objets (conditions aux limites)
  123. *
  124. *HYTOT = 'DOMA' CUBE1 ELI0            ;
  125. *C11   = 'DOMA' HYTOT 'VOLUME'       ;
  126. *CHYB1 = 'DOMA' HYTOT 'SURFACE'       ;
  127. *CHYB2 = 'DOMA' HYTOT 'NORMALE'       ;
  128. *MCHYB = 'DOMA' HYTOT 'ORIENTAT'      ;
  129. *
  130. *HYGAU = 'DOMA' SGAU 'INCL' HYTOT ELI0 ;
  131. *HYDRO = 'DOMA' SDRO 'INCL' HYTOT ELI0 ;
  132. *HYHAU = 'DOMA' SHAU 'INCL' HYTOT ELI0 ;
  133. *HYBAS = 'DOMA' SBAS 'INCL' HYTOT ELI0 ;
  134. *HYDEV = 'DOMA' SDEV 'INCL' HYTOT ELI0 ;
  135. *HYDER = 'DOMA' SDER 'INCL' HYTOT ELI0 ;
  136. MODHYB = MODE QFTOT 'DARCY' 'ORTHOTROPE'     ;
  137. MODGAU = MODE QFGAU 'DARCY' 'ORTHOTROPE'         ;
  138. MODDRO = MODE QFDRO 'DARCY' 'ORTHOTROPE'         ;
  139. MODHAU = MODE QFHAU 'DARCY' 'ORTHOTROPE'         ;
  140. MODBAS = MODE QFBAS 'DARCY' 'ORTHOTROPE'         ;
  141. MODDEV = MODE QFDEV 'DARCY' 'ORTHOTROPE'         ;
  142. MODDER = MODE QFDER 'DARCY' 'ORTHOTROPE'         ;
  143. C11   = 'DOMA' MODHYB 'VOLUME'       ;
  144. CHYB1 = 'DOMA' MODHYB 'SURFACE'       ;
  145. CHYB2 = 'DOMA' MODHYB 'NORMALE'       ;
  146. CEGAU  = 'DOMA' MODGAU 'CENTRE' ;
  147. CEDRO  = 'DOMA' MODDRO 'CENTRE' ;
  148. CEHAU  = 'DOMA' MODHAU 'CENTRE' ;
  149. CEBAS  = 'DOMA' MODBAS 'CENTRE' ;
  150. CEDEV  = 'DOMA' MODDEV 'CENTRE' ;
  151. CEDER  = 'DOMA' MODDER 'CENTRE' ;
  152. *
  153. *- Solution analytique
  154. *
  155. XX  YY  ZZ  = 'COOR' (DOMA MODHYB 'FACE' )  ;
  156. XXC YYC ZZC = 'COOR' (DOMA MODHYB 'CENTRE') ;
  157. *
  158. VKX = 1.D0             ;
  159. VKY = 0.75D0           ;
  160. VKZ = 0.5D0            ;
  161. AA  = -45.D0           ;
  162. BB  = -80.D0           ;
  163. CC  = -60.D0           ;
  164. DD  = 200.D0           ;
  165. AAA = -1.D0 * VKX * AA ;
  166. BBB = -1.D0 * VKY * BB ;
  167. CCC = -1.D0 * VKZ * CC ;
  168. *
  169. PANAF = (AA * XX)  + (BB * YY)  + (CC * ZZ)  + DD                 ;
  170. PANAC = (AA * XXC) + (BB * YYC) + (CC * ZZC) + DD                 ;
  171. VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 3 'VX' AAA 
  172.          'VY' BBB 'VZ' CCC ;
  173. VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE')   3 'VX' AAA 
  174.           'VY' BBB 'VZ' CCC ;
  175. *
  176. *                                                  --------------
  177. *                                                  = RESOLUTION =
  178. *                                                  --------------
  179. *
  180. *MODHYB = MODE HYTOT 'DARCY' 'ORTHOTROPE' 'HYC8'        ;
  181. MATI3  = 'MATE' MODHYB 'DIRECTION' VECX VECY 'PARALLELE'
  182.                 'K1' VKX 'K2' VKY 'K3' VKZ               ;
  183. CND1A  = 'MHYB' MODHYB  MATI3                       ;
  184. HND1A  = 'MATP' MODHYB  CND1A                       ;
  185. *                                                        ;
  186. *- Conditions aux limites
  187. *
  188. BBGAU = 'BLOQ' CEGAU 'TH' ;
  189. BBDRO = 'BLOQ' CEDRO 'TH' ;
  190. BBHAU = 'BLOQ' CEHAU 'TH' ;
  191. BBBAS = 'BLOQ' CEBAS 'TH' ;
  192. BBDEV = 'BLOQ' CEDEV 'TH' ;
  193. BBDER = 'BLOQ' CEDER 'TH' ;
  194. *
  195. *- TH imposée
  196. *
  197. TTIMP = 'REDU' PANAF CEGAU ;
  198. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  199. EEGAU = 'DEPI' BBGAU TTIM2          ;
  200. TTIMP = 'REDU' PANAF CEDRO ;
  201. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  202. EEDRO = 'DEPI' BBDRO TTIM2          ;
  203. TTIMP = 'REDU' PANAF CEBAS ;
  204. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  205. EEBAS = 'DEPI' BBBAS TTIM2          ;
  206. TTIMP = 'REDU' PANAF CEHAU ;
  207. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  208. EEHAU = 'DEPI' BBHAU TTIM2          ;
  209. TTIMP = 'REDU' PANAF CEDEV ;
  210. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  211. EEDEV = 'DEPI' BBDEV TTIM2          ;
  212. TTIMP = 'REDU' PANAF CEDER ;
  213. TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
  214. EEDER = 'DEPI' BBDER TTIM2          ;
  215. *
  216. *- Flux imposé
  217. *
  218. FLUD  = -60.D0 * DX * DZ ;
  219. FLUG  = -1.D0  * FLUD    ;
  220. FLUH  =  30.D0 * DX * DY ;
  221. FLUB  = -1.D0  * FLUH    ;
  222. FLUV  =  45.D0 * DY * DZ ;
  223. FLUR  = -1.D0  * FLUV    ;
  224. *
  225. FLGAU = 'MANU' 'CHPO' CEGAU 1 'FLUX' FLUG 'NATURE' 'DISCRET' ;
  226. FLDRO = 'MANU' 'CHPO' CEDRO 1 'FLUX' FLUD 'NATURE' 'DISCRET' ;
  227. FLBAS = 'MANU' 'CHPO' CEBAS 1 'FLUX' FLUB 'NATURE' 'DISCRET' ;
  228. FLHAU = 'MANU' 'CHPO' CEHAU 1 'FLUX' FLUH 'NATURE' 'DISCRET' ;
  229. FLDEV = 'MANU' 'CHPO' CEDEV 1 'FLUX' FLUV 'NATURE' 'DISCRET' ;
  230. FLDER = 'MANU' 'CHPO' CEDER 1 'FLUX' FLUR 'NATURE' 'DISCRET' ;
  231. *
  232. *- Assemblage et résolution en TH
  233. *
  234. CCC1 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBDRO 'ET'
  235.                   BBGAU 'ET' BBDEV 'ET' BBDER      ;
  236. FFF1 =            EEHAU 'ET' EEBAS 'ET' EEDRO 'ET'
  237.                   EEGAU 'ET' EEDEV 'ET' EEDER      ;
  238. *
  239. CCC2 = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO            ;
  240. FFF2 = FLDEV 'ET' FLDER 'ET' EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO ;
  241. *
  242. CCC3 = HND1A 'ET'                                  BBGAU 'ET' BBDEV ;
  243. FFF3 = FLHAU 'ET' FLBAS 'ET' FLDRO 'ET' FLDER 'ET' EEGAU 'ET' EEDEV ;
  244. *
  245. CHTER1 = 'RESO' CCC1 FFF1 ;
  246. CHTER2 = 'RESO' CCC2 FFF2 ;
  247. CHTER3 = 'RESO' CCC3 FFF3 ;
  248. *
  249. *- Calcul de H
  250. *
  251. PCEN1 = 'HYBP' MODHYB  CND1A CHTER1 ;
  252. PCEN2 = 'HYBP' MODHYB  CND1A CHTER2 ;
  253. PCEN3 = 'HYBP' MODHYB  CND1A CHTER3 ;
  254. *
  255. *- Calcul de V
  256. *
  257. QFACE1 = 'HDEB' MODHYB  CND1A PCEN1 CHTER1  ;
  258. VCENT1 = 'HVIT' MODHYB  QFACE1         ;
  259. QFACE1 = 'EXCO' QFACE1 'FLUX' 'SCAL'             ;
  260. VFACE1 = QFACE1 * CHYB2 / CHYB1                  ;
  261. *
  262. QFACE2 = 'HDEB' MODHYB  CND1A PCEN2 CHTER2  ;
  263. VCENT2 = 'HVIT' MODHYB  QFACE2         ;
  264. QFACE2 = 'EXCO' QFACE2 'FLUX' 'SCAL'             ;
  265. VFACE2 = QFACE2 * CHYB2 / CHYB1                  ;
  266. *
  267. QFACE3 = 'HDEB' MODHYB  CND1A PCEN3 CHTER3  ;
  268. VCENT3 = 'HVIT' MODHYB  QFACE3         ;
  269. QFACE3 = 'EXCO' QFACE3 'FLUX' 'SCAL'             ;
  270. VFACE3 = QFACE3 * CHYB2 / CHYB1                  ;
  271. *
  272. *                                                   -----------------
  273. *                                                   = Calcul ERREUR =
  274. *                                                   -----------------
  275. * ERReur relative en Trace de charge TH aux faces des éléments
  276. * ERReur relative en charge H au centre des éléments
  277. * Erreur relative sur la vitesse au centre des éléments
  278. *
  279. ERRTP1 = 'EXCO' CHTER1 'TH' 'SCAL' ;
  280. ERRTP1 = ERRTP1 - PANAF / PANAF    ;
  281. ERRTP1 = 'ABS' ERRTP1              ;
  282. ERRP1  = 'EXCO' PCEN1 'H' 'SCAL'   ;
  283. ERRP1  = ERRP1 - PANAC / PANAC     ;
  284. ERRP1  = 'ABS' ERRP1               ;
  285. *
  286. ERRTP2 = 'EXCO' CHTER2 'TH' 'SCAL' ;
  287. ERRTP2 = ERRTP2 - PANAF / PANAF    ;
  288. ERRTP2 = 'ABS' ERRTP2              ;
  289. ERRP2  = 'EXCO' PCEN2 'H' 'SCAL'   ;
  290. ERRP2  = ERRP2 - PANAC / PANAC     ;
  291. ERRP2  = 'ABS' ERRP2               ;
  292. *
  293. ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
  294. ERRTP3 = ERRTP3 - PANAF / PANAF    ;
  295. ERRTP3 = 'ABS' ERRTP3              ;
  296. ERRP3  = 'EXCO' PCEN3 'H' 'SCAL'   ;
  297. ERRP3  = ERRP3 - PANAC / PANAC     ;
  298. ERRP3  = 'ABS' ERRP3               ;
  299. *
  300. MOT1 = 'MOTS' 'VX' 'VY' 'VZ'        ;
  301. VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
  302. *
  303. VD1  = VANAC - VCENT1            ;
  304. VC1  = 'PSCA' VD1 VD1 MOT1 MOT1  ;
  305. SDC1 = 'ABS' ( VC1 / VDVD )      ;
  306. SDC1 = SDC1 '**' 0.5             ;
  307. VD2  = VANAC - VCENT2            ;
  308. VC2  = 'PSCA' VD2 VD2 MOT1 MOT1  ;
  309. SDC2 = 'ABS' ( VC2 / VDVD )      ;
  310. SDC2 = SDC2 '**' 0.5             ;
  311. VD3  = VANAC - VCENT3            ;
  312. VC3  = 'PSCA' VD3 VD3 MOT1 MOT1  ;
  313. SDC3 = 'ABS' ( VC3 / VDVD )      ;
  314. SDC3 = SDC3 '**' 0.5             ;
  315. *
  316. *                                                 -------------------
  317. *                                                 = Tracé resultats =
  318. *                                                 -------------------
  319. 'SI' ('NEG' GRAPH 'N') ;
  320. *
  321. *- Transformation des quantités aux centres en MCHAML constant.
  322. *
  323. ERRP1 = 'KCHA' MODHYB 'CHAM' ERRP1 ;
  324. ERRP2 = 'KCHA' MODHYB 'CHAM' ERRP2 ;
  325. ERRP3 = 'KCHA' MODHYB 'CHAM' ERRP3 ;
  326. SDS1  = 'KCHA' MODHYB 'CHAM' SDC1  ;
  327. SDS2  = 'KCHA' MODHYB 'CHAM' SDC2  ;
  328. SDS3  = 'KCHA' MODHYB 'CHAM' SDC3  ;
  329. *
  330. * Dans chaque cas on trace
  331. *   L'erreur relative sur la charge au centre
  332. *   L'erreur relative sur la Vitesse au centre
  333. *
  334. 'TITR' 'darcy3/1 : Erreur relative sur la charge'  ;
  335. 'TRAC' MODHYB ERRP1                                ;
  336. 'TITR' 'darcy3/1 : Erreur relative sur la vitesse' ;
  337. 'TRAC' MODHYB SDS1                                 ;
  338. *
  339. 'TITR' 'darcy3/2 : Erreur relative sur la charge'  ;
  340. 'TRAC' MODHYB ERRP2                                ;
  341. 'TITR' 'darcy3/2 : Erreur relative sur la vitesse' ;
  342. 'TRAC' MODHYB SDS2                                 ;
  343. *
  344. 'TITR' 'darcy3/3 : Erreur relative sur la charge'  ;
  345. 'TRAC' MODHYB ERRP3                                ;
  346. 'TITR' 'darcy3/3 : Erreur relative sur la vitesse' ;
  347. 'TRAC' MODHYB SDS3                                 ;
  348. *
  349. 'FINSI' ;
  350. *
  351. *                                                 -------------------
  352. *                                                 = Gestion ERREURS =
  353. *                                                 -------------------
  354. MAXTP1 = 'MAXI' ERRTP1 ;
  355. MAXTP2 = 'MAXI' ERRTP2 ;
  356. MAXTP3 = 'MAXI' ERRTP3 ;
  357. MAXP1  = 'MAXI' ERRP1  ;
  358. MAXP2  = 'MAXI' ERRP2  ;
  359. MAXP3  = 'MAXI' ERRP3  ;
  360. MAXV1  = 'MAXI' SDC1   ;
  361. MAXV2  = 'MAXI' SDC2   ;
  362. MAXV3  = 'MAXI' SDC3   ;
  363. *
  364. 'SAUT' 'PAGE'                                                       ;
  365. 'SAUT' 2 'LIGNE'                                                    ;
  366. 'MESS' '                             ERREURS RELATIVES            ' ;
  367. 'SAUT' 1 'LIGNE'                                                    ;
  368. 'MESS' ' cas test          TH                  H                 V' ;
  369. 'SAUT' 1 'LIGNE'                                                    ;
  370. 'MESS' ' numero1     ' maxtp1  '   ' maxp1 '   ' maxv1              ;
  371. 'SAUT' 1 'LIGNE'                                                    ;
  372. 'MESS' ' numero2     ' maxtp2  '   ' maxp2 '   ' maxv2              ;
  373. 'SAUT' 1 'LIGNE'                                                    ;
  374. 'MESS' ' numero3     ' maxtp3  '   ' maxp3 '   ' maxv3              ;
  375. 'SAUT' 2 'LIGNE'                                                    ;
  376. *
  377. EPS0 = 1.D-7 ;
  378. LOG1 = MAXTP1 > EPS0 ; LOG2 = MAXTP2 > EPS0 ; LOG3 = MAXTP3 > EPS0 ;
  379. LOG4 = MAXP1  > EPS0 ; LOG5 = MAXP2  > EPS0 ; LOG6 = MAXP3  > EPS0 ;
  380. LOG7 = MAXV1  > EPS0 ; LOG8 = MAXV2  > EPS0 ; LOG9 = MAXV3  > EPS0 ;
  381. LTP0 = LOG1 'OU' LOG2 'OU' LOG3 ;
  382. LP0  = LOG4 'OU' LOG5 'OU' LOG6 ;
  383. LV0  = LOG7 'OU' LOG8 'OU' LOG9 ;
  384. L0   = LTP0 'OU' LP0  'OU' LV0  ;
  385. 'SI' ( L0 ) ;
  386.    'ERRE' 5 ;
  387. 'SINO'      ;
  388.    'ERRE' 0 ;
  389. 'FINSI'     ;
  390. *
  391. 'FIN' ;
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  

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