Télécharger transsatVF.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : transsat.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. GRAPH = vrai ;
  5. ******************** CAS TEST : transsat.dgibi ************************
  6. *
  7. *
  8. * Test de fonctionnement de DARCYSAT sur un probleme multizone.
  9. * ===================================================================
  10. * Infiltration d'eau dans une barrière ouvragée dans son site d'acceuil.
  11. *
  12. * - condition initiale :
  13. * > BO : désaturation uniforme correspondant a une succion de 7827,42 m ;
  14. * > SITE saturé sous 500m d'eau
  15. * - les paramètres des lois de perméabilité diffèrent entre le SITE et la BO
  16. *
  17. * ===================================================================
  18. * Les options de modélisation déclarées dans la table transmise à
  19. * la procedure DARCYSAT sont les suivantes :
  20. *
  21. * - les effets gravitationnels ne sont pas pris en compte : pression =
  22. * charge (indice GRAVITE absent ; valeur par defaut : FAUX);
  23. *
  24. * - Le pas de temps est automatique (indice TEMPS_CALCULES absent)
  25. * sont fournis :
  26. * > le pas de temps initial (indice 'DT_INITIAL'),
  27. * > le nombre d'itérations recherchées par pas de temps (indice 'NITER')
  28. * > le nombre de pas de temps maximum (indice 'NPAS')
  29. * > le nombre de CFL initial (indice 'CFL')
  30. *
  31. * - La méthode de pénalisation est désactivée (indice XI absent) ;
  32. *
  33. * - Homogénéisation : les propriétés physiques sont calculées à partir
  34. * des charges aux éléments (indice HOMOGENEISATION de valeur CENTRE)
  35. *
  36. *
  37. ************************************************************************
  38. 'OPTION' 'ECHO' 1 ;
  39. 'SAUTER' 'PAGE';
  40. *
  41.  
  42. TITRE 'Saturation d une barriere ouvragee : transsat.dgibi' ;
  43. 'OPTION' 'DIME' 2 'ELEM' 'QUA4';
  44. 'OPTION' 'ISOV' 'LIGN' 'TRACER' 'PSC' ;
  45.  
  46. *--------------------------------------------------------------------
  47. * Création maillage
  48. *
  49. *- Discrétisation :
  50. ENX = 1 ;
  51. ENY = 100 ;
  52. 'DENS' (3./eny) ;
  53.  
  54. *- Création des points et des droites
  55. *
  56. A0 = -.5 6.D0 ; B0 = .5 6.D0 ;
  57. A1 = -.5 4.D0 ; B1 = .5 4.D0 ;
  58. A2 = -.5 1.D0 ; B2 = .5 1.D0 ;
  59. A3 = -.5 0.D0 ; B3 = .5 0.D0 ;
  60.  
  61. *- Création des droites
  62. *
  63. AB0 = 'DROIT' ENX A0 B0 ;
  64. AB1 = 'DROIT' ENX A1 B1 ;
  65. AB2 = 'DROIT' ENX A2 B2 ;
  66. AB3 = 'DROIT' ENX A3 B3 ;
  67.  
  68. *- Création des surfaces
  69. *
  70. BO = 'REGLER' AB3 AB2 ;
  71. BO = 'COULEUR' BO 'VERT' ;
  72. BO = 'INVERSE' BO ;
  73. SITE = AB2 'REGLER' DINI (3./eny) DFIN 0.2 AB1 'REGLER' 10 AB0 ;
  74. SITE = 'COULEUR' SITE 'BLEU' ;
  75. SITE = 'INVERSE' SITE ;
  76. MASSIF0 = BO 'ET' SITE;
  77.  
  78. * entrée
  79. ENTR = 'COULEUR' ('INVERSE' AB0) 'ROUGE' ;
  80.  
  81. 'SI' GRAPH ;
  82. 'TRACER' (MASSIF0 'ET' ENTR) ;
  83. 'FINSI' ;
  84.  
  85. *- Création des maillages contenant tous les points
  86. *************************************************************************
  87. * IMPORTANT : ne pas concaténer de maillages QUAF *
  88. *************************************************************************
  89. *
  90. QMASSIF = 'CHANGER' MASSIF0 'QUAF' ;
  91. QSITE = 'CHANGER' SITE 'QUAF' ;
  92. QBO = 'CHANGER' BO 'QUAF' ;
  93. QFENTR = 'CHANGER' ENTR 'QUAF' ;
  94.  
  95. 'ELIMINATION' (QMASSIF 'ET' QFENTR 'ET' QBO 'ET' QSITE) 0.00001 ;
  96.  
  97. *- Modèles
  98. *************************************************************************
  99. * IMPORTANT : Il faut concaténer les modèles et ne surtout pas créer un *
  100. * modèle total à une seule zone *
  101. *************************************************************************
  102. MBO = 'MODELE' QBO 'DARCY' 'ISOTROPE' ;
  103. MSITE = 'MODELE' QSITE 'DARCY' 'ISOTROPE' ;
  104. MENTR = 'MODELE' QFENTR 'DARCY' 'ISOTROPE' ;
  105. MODHYB = MSITE 'ET' MBO;
  106.  
  107. CEENTR = 'DOMA' MENTR 'CENTRE' ;
  108. HYCEN = 'DOMA' MODHYB 'CENTRE' ;
  109. HYFAC = 'DOMA' MODHYB 'FACE';
  110.  
  111. *- Création ligne de suivi pour le post-traitement
  112. * ligne des points centres (cas 1D)
  113. *
  114. * on ordonnance de bas en haut les points centres
  115. orig = 0. 0. ;
  116. lcen1 = orig 'ET' hycen ;
  117. 'ORDO' lcen1 ;
  118. lcen2 = 'DIFF' lcen1 ('MANU' 'POI1' orig) ;
  119. * et on crée la ligne de seg2
  120. lcenc = 'QUELCONQUE' lcen2 'SEG2' ;
  121.  
  122. *- ordonnées
  123. *
  124. NDIME = 'VALEUR' 'DIME' ;
  125. ZCC = 'COOR' HYCEN NDIME ;
  126. ZFF = 'COOR' HYFAC NDIME ;
  127.  
  128. *--------------------------------------------------------------------
  129. *
  130. * Paramètres
  131. *
  132. unjour = 24 * 3600. ;
  133. unmois = 30 * unjour ;
  134. unan = unjour * 365.25 ;
  135. Tini = 0. ;
  136. Tfin = 1000 * unan ;
  137.  
  138. *--- profondeur du site
  139. Profsite = 500. ;
  140. T = 50 + 273.15 ;
  141. hc_bo = -7827.42 ;
  142. RHOg = 9810. ;
  143. poro_bo = 0.3 ;
  144.  
  145. *--- coefficient d'emmagasinement
  146. * C = poro * rhoeau * g *
  147. * (compressibilite_eau - compressibilite_solide + alpha/poro)
  148. * avec alpha = 3 (1 - 2 Nupoisson) / Eyoung
  149. * pour FOCA : nu = 0.25 ; E = 500 MPa
  150. * poro = 0.3 ; RHOg = 9810 ; compres_eau = 5E-10 Pa-1 ; compres_solide = 0
  151. emma_sit = 3.83E-6 ;
  152. nu_bo = 0.4 ;
  153. Young_bo = 150.E6 * 2 * (1 + nu_bo) ;
  154. alph_bo = 3 * (1 - (2 * nu_bo)) / Young_bo ;
  155. emma_bo = 1.D0 *poro_bo * RHOg * (5E-10 - 0 + (alph_bo/poro_bo)) ;
  156. emma0 = ('MANU' 'CHPO' ('DOMA' msite 'CENTRE') 1
  157. 'SCAL' emma_sit 'NATURE' 'DISCRET')
  158. 'ET' ('MANU' 'CHPO' ('DOMA' mbo 'CENTRE') 1
  159. 'SCAL' emma_bo 'NATURE' 'DISCRET') ;
  160. emma0 = 'KCHT' modhyb 'SCAL' 'CENTRE' 'COMP' 'SCAL' emma0 ;
  161.  
  162. *emma0 = 1.D-2 * emma0 ;
  163.  
  164.  
  165. *--------------------------------------------------------------------
  166. * Conditions aux limites
  167. * 1) flux nul sur la frontiere du colis => rien a faire
  168. * 2) flux nuls en latéral => rien a faire
  169. * 3.1) pression imposee en limite de site Pliq = h avec h =500m
  170. BENTR = 'BLOQUE' CEENTR 'TH' ;
  171. EENTR = 'MANUEL' 'CHPO' CEENTR 1 'TH' Profsite ;
  172.  
  173. * chargement des CLs
  174. LICALC = 'PROG' 0.D0 1.e10 ;
  175. LIST1 = 'PROG' 2 * 1.D0 ;
  176. VALI0 = 'CHAR' 'THIM' EENTR ('EVOL' 'MANU' LICALC LIST1) ;
  177.  
  178. *--------------------------------------------------------------------
  179. *- initialisation des inconnues
  180. * (doit être compatible avec les conditions aux limites)
  181. *
  182. Thi_site = 'MANU' 'CHPO' ('DOMA' Msite 'FACE') 1
  183. 'TH' Profsite 'NATURE' 'DISCRET' ;
  184. Thi_bo = 'MANU' 'CHPO' ('DOMA' Mbo 'FACE') 1
  185. 'TH' hc_bo 'NATURE' 'DISCRET' ;
  186. Th_ini = 'KCHT' modHYB 'SCAL' 'FACE' 'COMP' 'TH' Thi_site Thi_bo;
  187. *
  188. Hi_site = 'MANU' 'CHPO' ('DOMA' Msite 'CENTRE') 1
  189. 'SCAL' Profsite 'NATURE' 'DISCRET' ;
  190. Hi_bo = 'MANU' 'CHPO' ('DOMA' Mbo 'CENTRE') 1
  191. 'SCAL' hc_bo 'NATURE' 'DISCRET' ;
  192. H_ini = kcht modhyb SCAL CENTRE (Hi_site 'ET' Hi_bo) ;
  193. H_ini = nomc 'H' H_INI ;
  194. Flux_ini = 'MANU' 'CHPO' ('DOMA' modhyb 'FACE') 1
  195. 'FLUX' 0. 'NATURE' 'DISCRET' ;
  196.  
  197. * ---------------------------
  198. * = Table DARCY_TRANSITOIRE =
  199. * ---------------------------
  200. *- initialisation table
  201. SATUR = 'TABLE' ;
  202. SATUR. 'TEMPS' = 'TABLE' ;
  203. SATUR. 'CHARGE' = 'TABLE' ;
  204. SATUR. 'FLUX' = 'TABLE' ;
  205.  
  206. *- données géométriques
  207. SATUR. 'SOUSTYPE' = 'DARCY_TRANSATUR' ;
  208. SATUR. 'MODELE' = MODHYB ;
  209.  
  210. *- conditions initiales
  211. SATUR. 'TEMPS' . 0 = Tini ;
  212. SATUR. 'CHARGE' . 0 = H_ini ;
  213. SATUR. 'FLUX' . 0 = Flux_ini ;
  214.  
  215. *- conditions aux limites et chargements
  216. SATUR. 'TRACE_IMPOSE' = VALI0 ;
  217.  
  218.  
  219. *- définition de la loi de perméabilité
  220. SATUR.'LOI_PERMEABILITE' = 'TABLE' 'MULTIZONE' ;
  221.  
  222. * pour le site
  223. LoiP = 'TABLE' 'PUISSANCE';
  224. LoiP . 'ALPHA' = RHOg / 1.D7 ;
  225. LoiP . 'PERMSAT' = 1.D-12 ;
  226. LoiP . 'MODELE' = MSITE ;
  227. SATUR. 'LOI_PERMEABILITE'. 'SITE' = LoiP ;
  228.  
  229. * pour la BO
  230. LoiP = 'TABLE' 'PUISSANCE';
  231. LoiP .'ALPHA' = RHOg / (7.D6);
  232. LoiP .'PERMSAT' = 6.D-14 ;
  233. LoiP .'MODELE' = MBO ;
  234. SATUR.'LOI_PERMEABILITE'. 'BO' = LoiP ;
  235.  
  236. *- définition de la teneur en eau
  237. SATUR.'LOI_SATURATION' = 'TABLE' 'MULTIZONE' ;
  238.  
  239. * pour le site
  240. LoiS = 'TABLE' 'VAN_GENUCHTEN';
  241. LoiS . 'PORO' = 0.05 ;
  242. LoiS . 'TERESIDU' = 0. ;
  243. LoiS . 'NEXP' = 1.7 ;
  244. LoiS . 'MEXP' = 1 - (1 / 1.7) ;
  245. LoiS . 'BHETA' = RHOg / 10.D6 ;
  246. LoiS . 'MODELE' = MSITE ;
  247. SATUR.'LOI_SATURATION' . 'SITE' = LoiS ;
  248.  
  249. 'SI' GRAPH ;
  250. *- Vérification du choix du dp minimum pour le calcul de la capacité.
  251. * droite support des variables
  252. dx = 'DROIT' (0. -7827.4 ) 1000 (0. 0.) ;
  253. zc = 'COOR' dx 2 ;
  254. ev2 = 'EVOL' 'BLEU' 'CHPO' zc 'SCAL' dx ;
  255. px = 'EXTR' ev2 'ORDO' ;
  256. * calcul de la teneur en eau pour la pression zc
  257. s0 t0 dum = HT_PRO (SATUR.'LOI_SATURATION'.'SITE') ZC ;
  258. ev0 = 'EVOL' 'TURQ' 'CHPO' s0 'SCAL' dx ;
  259. evt = 'EVOL' 'VERT' 'MANU' px (100. * ('EXTR' ev0 'ORDO')) ;
  260. 'DESSIN' evt 'TITX' 'Pc(m)' 'TITY' 'S(%)'
  261. 'TITRE' 'Loi capillaire S(Pc)' ;
  262. * calcul de la teneur en eau pour la pression zc - dp
  263. dp = 1. ;
  264. s1 t1 dum = HT_PRO (SATUR.'LOI_SATURATION'.'SITE') ('KOPS' zc '-' dp);
  265. * représentation de la capacité
  266. c1 = (t0 - t1) / dp;
  267. ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
  268. evc = 'EVOL' 'TURQ' 'MANU' px ('EXTR' ev1 'ORDO') ;
  269. 'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
  270. 'TITRE' 'Capacite capillaire SITE' ;
  271. FINSI ;
  272.  
  273. * pour la BO
  274. LoiS = TABLE 'VAN_GENUCHTEN';
  275. LoiS . 'PORO' = poro_bo ;
  276. LoiS . 'TERESIDU' = 0. ;
  277. LoiS . 'NEXP' = 1. / (1. - 0.27) ;
  278. LoiS . 'MEXP' = 0.27 ;
  279. LoiS . 'BHETA' = RHOg / 17.E6 ;
  280. LoiS . 'MODELE' = MBO ;
  281. SATUR.'LOI_SATURATION' . 'BO' = LoiS ;
  282.  
  283. 'SI' GRAPH ;
  284. *- Vérification du choix du dp minimum pour le calcul de la capacité.
  285. *--- calcul de la teneur en eau pour p = zc
  286. * calcul de la teneur en eau pour la pression zc
  287. s0 t0 dum = HT_PRO (SATUR.'LOI_SATURATION'.'BO') ZC ;
  288. ev0 = 'EVOL' 'TURQ' 'CHPO' s0 'SCAL' dx ;
  289. evt = 'EVOL' 'VERT' 'MANU' px (100. * ('EXTR' ev0 'ORDO')) ;
  290. 'DESSIN' evt 'TITX' 'Pc(m)' 'TITY' 'S(%)'
  291. 'TITRE' 'Loi capillaire S(Pc)' ;
  292. * calcul de la teneur en eau pour la pression zc - dp
  293. dp = 1. ;
  294. s1 t1 dum = HT_PRO (SATUR.'LOI_SATURATION'.'BO') ('KOPS' zc '-' dp);
  295. * représentation de la capacité
  296. c1 = (t0 - t1) / dp;
  297. ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
  298. evc = 'EVOL' 'TURQ' 'MANU' px ('EXTR' ev1 'ORDO') ;
  299. 'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
  300. 'TITRE' 'Capacite capillaire B.O.' ;
  301. FINSI ;
  302. SATUR. 'COEF_EMMAGASINEMENT' = emma0 ;
  303.  
  304. * données numériques
  305. SATUR . 'ITMAX' = 40;
  306. SATUR. 'TEMPS_FINAL' = 1.e10 ;
  307. SATUR. 'HOMOGENEISATION' = 'CENTRE' ;
  308. SATUR. 'NPAS' = 100 ;
  309. SATUR . 'RESIDU_MAX' = 1.D-4;
  310. SATUR. 'DT_INITIAL' = 12 * 3600. ;
  311. SATUR. 'MESSAGES' = 1 ;
  312. SATUR . 'LUMP' = FAUX ;
  313.  
  314. SATUR . 'LUMP' = FAUX;
  315. SATUR . 'TYPDISCRETISATION' = 'VF' ;
  316. SATUR . 'DECENTR' = FAUX;
  317. TABRES = table METHINV;
  318. TABRES . 'TYPINV' = 3;
  319. TABRES . 'PRECOND' = 3;
  320.  
  321. SATUR . 'METHINV' = TABRES;
  322. SATUR . 'METHINV' . RESID = 1.D-20;
  323.  
  324. SATUR. 'TEMPS_SAUVES' = ('PROG' 1 'PAS' 1 41.)*36000. ;
  325. SATUR. 'TEMPS_CALCULES' = ('PROG' 1 'PAS' 1 41.)*36000. ;
  326.  
  327.  
  328. * ===============
  329. * Calcul
  330. * ===============
  331. DARCYSAT SATUR ;
  332.  
  333. * ===============
  334. * Post-traitement
  335. * ===============
  336. * Le test est effectué en vérifiant la solution
  337.  
  338. *- solution après 40 pas de temps
  339. lp = 'PROG' -7827.4 -7827.4 -7827.4 -7827.4 -7827.4 -7827.4 -7827.4
  340. -7827.4 -7827.4 -7827.4 -7827.3 -7827.2 -7826.9 -7826.3
  341. -7825.2 -7823.1 -7819.2 -7812.2 -7800.1 -7779.9 -7747.0
  342. -7695.1 -7615.8 -7498.0 -7328.3 -7090.4 -6765.7 -6334.1
  343. -5773.9 -5064.7 -4189.3 -3141.8 -1944.8 -1272.1 -1188.6
  344. -1101.0 -1009.7 -915.24 -818.60 -720.85 ;
  345. lp = lp 'ET' ('PROG'
  346. -623.36 -527.60 -435.04 -346.88 -263.88 -186.13 -113.58
  347. -47.363 14.082 73.448 130.79 185.41 236.60 283.74
  348. 326.27 363.79 396.07 423.06 444.92 462.03 474.90
  349. 484.16 490.51 494.64 497.16 498.60 499.36 499.72
  350. 499.88 499.95 499.98 499.99 500.00 500.00 500.00
  351. 500.00);
  352.  
  353.  
  354.  
  355. id = 40 ;
  356. ev = 'EVOL' 'CHPO' SATUR.'CHARGE'.id 'H' LCENC ;
  357. lp2 = 'EXTR' ev 'ORDO' ;
  358. err1 = SOMM ('ABS' (lp - lp2));
  359. err1 = err1 / (SOMM (ABS lp));
  360. mess 'régression : ' err1 ;
  361.  
  362. 'LISTE' lp2 ;
  363.  
  364.  
  365. 'SI' (err1 > 5.D-1);
  366. mess 'erreur anormale après ' (@ARR id 0) ' pas : ' err1 ;
  367. * 'ERREUR' 5;
  368. 'SINON' ;
  369. 'ERREUR' 0;
  370. 'FINSI' ;
  371.  
  372. * ===============
  373. * Post-traitement
  374. * ===============
  375.  
  376. 'SI' GRAPH ;
  377. *--- liste des indices des temps postraites
  378. NN = (('DIME' SATUR.TEMPS) - 1) / 3 ;
  379. LT = 'LECT' 0 'PAS' 3 (3 * NN) ;
  380. TDES = TRACHIS SATUR 'SATURATION' LCENC LT
  381. ('MOTS' 'SCAL') ('MOTS' ' ') 'PREF' ' ' ;
  382. DESTRA TDES 'MIMA' 'AXES' 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
  383. TDES = TRACHIS SATUR 'CHARGE' LCENC LT
  384. ('MOTS' 'H') ('MOTS' ' ') 'PREF' ' ' ;
  385. DESTRA TDES 'MIMA' 'AXES' 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
  386. 'FINSI' ;
  387.  
  388. 'FIN';
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  

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