Télécharger la plaquette Cast3M


fr en

Télécharger darcy5.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  darcy5.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. ********************  CAS TEST : darcy5.dgibi  **********************
  6. *
  7.  GRAPH = FAUX   ;
  8. *GRAPH = VRAI   ;
  9. 'SAUT' 'PAGE' ;
  10. *
  11. *-------------------------------------------------------------------
  12. * TEST DARCY5
  13. *
  14. *    Validation de la résolution des équations de DARCY formulées en
  15. *             (VITESSE - PRESSION) par comparaison  avec 
  16. * la résolution des équations de DARCY formulées en (VITESSE - CHARGE)
  17. *
  18. *
  19. *    (Résolution par une méthode d'éléments finis mixtes hybrides)
  20. *
  21. * Le maillage est composés de triangles et de quadrangles et contient
  22. * des lignes de maillages courbes 
  23. *
  24. * la massif est supposé isotrope, indéformable et saturé
  25. *
  26. * Pour permettre la comparaison, la densité de l'eau reste constante 
  27. *
  28. *-------------------------------------------------------------------
  29. *
  30.  OPTI ECHO 0                                              ;
  31. *
  32. 'SAUT' 'PAGE' ;
  33. *
  34. *- Options générales de calcul
  35. *
  36. TITRE 
  37. 'EFMY DARCY FORMULATION (VITESSE - PRESSION) 2D : darcy5.dgibi' ;
  38.  
  39. OPTI DIME 2 ELEM QUA4                                    ;
  41.  
  42. *
  43. *- Données physiques
  44. *
  45. *- Densité de l'eau
  46. RHO0 = 1000. ;
  47. *
  48. *- Pression athmosphérique
  49. P0 = 1.013E5 ;
  50. *
  51. *- Composante du vecteur gravité
  52. Gx0 = 0. ;
  53. Gy0 = -9.81 ;
  54. *
  55. *- Perméabilité intrinsèque
  56. LA0 = 1.e-9 ;
  57. *
  58. *======================================================================
  59. *
  60. *- Construction du maillage
  61. *
  62. *======================================================================
  63. *
  64. *- Discrétisation spatiale
  65. *
  66. ENX = 15          ;
  67. ENY0 = 5         ;
  68. ENY1 = 5        ;
  69. ENY2 = 5         ;
  70. *
  71. *- Création des points
  72. *
  73. A0 = 0.d0 -2.d3 ; B0 = 5.d3 -2.d3 ;
  74. *
  75. A1 = 0.d0 -1.d3 ; B1 = 5.d3 -1.d3 ; C1 = 2.d3 -1.25d3 ;
  76. *
  77. A2 = 0.d0 0.0d3 ; B2 = 5.d3 0.0d3 ; C2 = 2.d3 0.25d3 ;
  78. *
  79. A3 = 0.d0 1.5d3 ; B3 = 5.d3 1.d3 ; C3 = 2.5d3 1.5d3 ;
  80. *
  81. *- Création des droites
  82. *
  83. AB0 = DROI ENX A0 B0 ;
  84. AB1 = A1 cer3 C1 ENX B1 ;
  85. AB2 = A2 cer3 C2 ENX B2 ;
  86. AB3 = A3 cer3 C3 ENX B3 ;
  87. *
  88. *- Creation des SURFACES
  89. *
  90. NIVEAU0 = AB0 regle ENY0  AB1 ;
  91. OPTI ELEM TRI3 ;
  92. NIVEAU1 = AB1 regle ENY1  AB2 ;
  93. OPTI ELEM QUA4 ;
  94. NIVEAU2 = AB2 regle ENY2  AB3 ;  
  95. MASSIF0 = NIVEAU0 ET NIVEAU1 ET NIVEAU2  ;
  96. HAUT  = INVE AB3 ;
  97. MASSIF0 = ORIENT MASSIF0        ;
  98. QMASSIF = CHANGE MASSIF0 QUAF ;
  99. QFHAUT  = CHANGE HAUT QUAF ;
  100.  ELIM 0.01 (QMASSIF ET QFHAUT ) ;
  101. *
  102. *- Domaines
  103. *
  104. *HYTOT = DOMA MASSIF0 0.01   ;
  105. *CHYB1 = DOMA HYTOT 'SURFACE' ;
  106. *CHYB2 = DOMA HYTOT 'NORMALE' ;
  107. *MCHYB = DOMA HYTOT 'ORIENTAT' ;
  108. *HYHAUT = DOMA HAUT INCL HYTOT 0.01 ;
  109. *
  110. *- Modélisation
  111. *
  112. *MODHYB = MODE HYTOT DARCY ISOTROPE HYQ4  HYT3                 ;
  113.  MODHYB = MODE QMASSIF DARCY ISOTROPE                 ;
  114.  MODHAU = MODE QFHAUT DARCY ISOTROPE                 ;
  115.  CHYB1 = DOMA MODHYB 'SURFACE' ;
  116.  CHYB2 = DOMA  MODHYB 'NORMALE' ;
  117.  CEHAUT= DOMA MODHAU 'CENTRE' ;
  118.  CETOT= DOMA MODHYB 'CENTRE' ;
  119. *
  120. *- perméabilité intrinsèque
  121. *
  122. LMME =  MANU CHPO CETOT 1 'K' la0 'NATURE' 'DIFFU' ;
  123. *
  124. *======================================================================
  125. *
  126. *- Résolution en  en (VITESSE ; PRESSION)
  127. *
  128. *======================================================================
  129. *
  130. LMMA = KCHA MODHYB 'CHAM' LMME                       ;
  131. MATI_P  = MATERIAU  MODHYB 'K' LMMA                 ;
  132. *
  133. ****  Masse HYBride
  134. CND1A_P  = MHYBR MODHYB  MATI_P                ;
  135. *
  136. ****  Masse FORce
  137. M_P  = MHYBR MODHYB 'MASSE'                ;
  138. *
  139. ****  MAtrice en Trace de pression TP
  140. HND1A_P  = MATP MODHYB  CND1A_P                ;
  141. *
  142. *- Conditions aux limites
  143. *
  144. BHAUT_P = BLOQUE CEHAUT 'TH' ;
  145. *
  146. *- TP imposee
  147. *
  148. PIMP = MANU CHPO CEHAUT 1 'TH' P0 'NATURE' 'DISCRET' ;
  149. EHAUT_P = DEPI BHAUT_P PIMP  ; 
  150. *
  151. *-  contribution des forces de volume au second membre 
  152. *
  153. RHO = MANU CHPO CETOT 1 'SCAL' RHO0 'NATURE' 'DISCRET' ;
  154. RCH =  MANU CHPO CETOT 2 'FX' gx0 'FY' gy0  'NATURE' 'DISCRET' ;
  155. RCH =  RHO * RCH ;
  156. GRAV2TP = SQTP MODHYB  HND1A_P M_P RCH ;
  157. *
  158. *- Assemblage matrice et second membre
  159. *
  160. CCC1_P = HND1A_P ET BHAUT_P ;
  161. * second membre du système matriciel en trace de charge
  162. FFF1_P = GRAV2TP et EHAUT_P ;
  163. *
  164. *- Resolution en trace de pression
  165. *
  166. CHTER1_P = RESOUDRE CCC1_P FFF1_P ;
  167. CHTER1_P = exco CHTER1_P 'TH' 'TH';
  168. *
  169. *- Calcul des pressions élémentaires
  170. *
  171. PCEN1 = HYBP MODHYB  CND1A_P CHTER1_P M_P RCH  ;
  172. *
  173. *- Calcul des débits
  174. *
  175. QFACE1_P = HDEBI MODHYB  CND1A_P PCEN1 CHTER1_P  M_P RCH;
  176. VCENT1_P = HVIT  MODHYB  QFACE1_P  ;
  177. *VV = VECT VCENT1_P 0.5e9 'VX' 'VY'  VERT ;
  178. *titre 'VITESSES AUX CENTRES DES ELEMENTS : ECHL=1e9:2' ;
  179. *TRAC VV CETOT ;
  180. *
  181. *=======================================================================
  182. *
  183. *- Résolution en (VITESSE ; CHARGE)
  184. *
  185. *=======================================================================
  186. *
  187. *
  188. *- conductivité hydraulique
  189. *
  190. LMME = LMME * (ABS Gy0) * RHO0  ;
  191. LMMA = KCHA MODHYB 'CHAM' LMME                       ;
  192. MATI_H  = MATERIAU  MODHYB 'K' LMMA                 ;
  193. *
  194. ****  Masse HYBride
  195. CND1A_H  = MHYBR MODHYB  MATI_H                ;
  196. *
  197. ****  MAtrice en Trace de charge TH
  198. HND1A_H  = MATP MODHYB  CND1A_H                ;
  199. *
  200. *- Conditions aux limites
  201. *
  202. BHAUT_H = BLOQUE CEHAUT 'TH' ;
  203. *
  204. *- TH imposee
  205. *
  206. PIMP = NOMC 'TH' (CEHAUT COOR 2) ;
  207. EHAUT_H = DEPI BHAUT_H PIMP  ;   
  208. *
  209. *- Assemblage matrice et second membre
  210. *
  211. CCC1_H = HND1A_H ET BHAUT_H ;
  212. * second membre du système matriciel en trace de charge
  213. FFF1_H = EHAUT_H ;
  214. *
  215. *- Resolution en trace de charge
  216. *
  217. CHTER1_H = RESOUDRE CCC1_H FFF1_H ;
  218. CHTER1_H = exco CHTER1_H 'TH' 'TH';
  219. *
  220. *- Calcul de la charge
  221. *
  222. TCEN1 = HYBP MODHYB  CND1A_H CHTER1_H   ;
  223. *
  224. *- Calcul de la vitesse
  225. *
  226. QFACE1_H = HDEBI MODHYB  CND1A_H TCEN1 CHTER1_H  ;
  227. VCENT1_H = HVIT  MODHYB  QFACE1_H  ;
  228. *VV = VECT VCENT1_H 0.5e9 'VX' 'VY'  VERT ;
  229. *titre 'CALCUL EN CHARGE : VITESSES AUX CENTRES DES ELTS : ECHL=1e9:2';
  230. *TRAC VV CETOT ;
  231. *
  232. *=======================================================================
  233. *
  234. *- Comparaison des deux résolutions
  235. *
  236. *=======================================================================
  237. *
  238. *
  239. *- comparaison des traces de charges
  240. *
  241. TP2 =  (NOMC CHTER1_P SCAL) - P0 / (rho0 * gy0)  ;
  242. TP2 = ((DOMA MODHYB 'FACE') COOR 2) - TP2 ;
  243. ER2 = (NOMC CHTER1_H SCAL)/TP2 - 1. ;
  244. maxtp1 = ER2 ABS MAXI ;
  245. *
  246. *
  247. *- comparaison des charges
  248. *
  249. PCEN2 = (NOMC PCEN1 SCAL) - P0 / (rho0 * gy0)  ;
  250. PCEN2 = CETOT COOR 2 - PCEN2 ;
  251. ER2 = (NOMC TCEN1 SCAL)/PCEN2 - 1. ;
  252. maxp1 = ER2 ABS MAXI ;
  253.  
  254. 'SI' GRAPH  ;
  255. TITRE 'VARIATIONS RELATIVES DES CHARGES DH/H' ; 
  256. TRAC (KCHA MODHYB 'CHAM' er2) MODHYB (CONTOUR
  257.  (DOMA MODHYB 'MAILLAGE')) ;
  258. 'FINSI' ;
  259. *
  260. *- comparaison des vitesse
  261. *
  262. QFACE1_P = NOMC SCAL QFACE1_P                ;
  263. VFACE1_P = QFACE1_P * CHYB2 / CHYB1          ;
  264. QFACE1_H = NOMC SCAL QFACE1_H                ;
  265. VFACE1_H = QFACE1_H * CHYB2 / CHYB1          ;
  266. maxvf1 = (kops QFACE1_P '-' QFACE1_H) / (maxi QFACE1_H) ABS MAXI ;
  267.  
  268. MOT1 = 'MOTS' 'VX' 'VY'             ;
  269. VDVD = 'PSCA' VCENT1_H VCENT1_H MOT1 MOT1 ;
  270. VD1  = VCENT1_P - VCENT1_H           ;
  271. VC1  = 'PSCA' VD1 VD1 MOT1 MOT1 ;
  272. SDC1 = 'ABS' ( VC1 / VDVD )     ;
  273. SDC1 = SDC1 '**' 0.5D0          ;
  274. maxvc1 = MAXI SDC1 ;
  275.  
  276. 'SI' GRAPH  ;
  277. TITRE ' VARIATIONS RELATIVES DES VITESSES |DV|/|V|' ; 
  278. TRAC ('KCHA' MODHYB 'CHAM' SDC1) MODHYB (CONTOUR
  279.  (DOMA MODHYB 'MAILLAGE')) ;
  280. 'FINSI' ;
  281.  
  282. *
  283. 'SAUT' 'PAGE'                                                        ;
  284. 'SAUT' 2 'LIGNE'                                                     ;
  285. 'MESS' '                       VARIATIONS RELATIVES                '   ;
  286. 'SAUT' 1 'LIGNE'                                                     ;
  287. 'MESS' '      TH                H                Vf              Vc' ;
  288. 'SAUT' 1 'LIGNE'                                                     ;
  289. 'MESS' ' ' maxtp1  ' ' maxp1 ' ' maxvf1 ' ' maxvc1                   ;
  290. 'SAUT' 2 'LIGNE'                                                     ;
  291. *
  292. EPS0 = 1.E-7 ; 
  293. LOG1 = MAXTP1 > EPS0                      ;
  294. LOG2 = MAXP1  > EPS0                      ;
  295. LOG3 = MAXVF1 > EPS0                      ;
  296. LOG4 = MAXVC1 > EPS0                      ;
  297. L0   = LOG1 'OU' LOG2 'OU' LOG3 'OU' LOG4 ;
  298. *
  299. 'SI' ( L0 ) ;
  300.    'ERRE' 5 ;
  301. 'SINO'      ;
  302.    'ERRE' 0 ;
  303. 'FINS'      ;
  304.  
  305. FIN ;
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  

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