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.  

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