Télécharger satnsathoriz.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : satnsathoriz.dgibi
  2. 'OPTI' 'TRAC' PSC ;
  3. GRAPH = FAUX ;
  4.  
  5.  
  6. ********************************************
  7. DEBPROC HHT_pro H1*FLOTTANT ;
  8. T1 = H1 * BET ** N + 1. ** m * TSR + TR ;
  9. FINPROC T1 ;
  10. ********************************************
  11.  
  12. ********************************************
  13. DEBPROC KKR_pro H0*FLOTTANT ;
  14. *K0 = H0 * ALF EXP * KS ;
  15. tt = HHT_PRO H0;
  16. K0 = (tt '-' tr) '/' TSR ** ALF * KS ;
  17. FINPROC K0 ;
  18. ********************************************
  19.  
  20. ********************************************
  21. DEBPROC TH_pro T1*FLOTTANT ;
  22. H1 = T1 - TR / TSR ** M1 - 1. ** N1 / BET ;
  23. FINPROC H1 ;
  24. ********************************************
  25.  
  26.  
  27.  
  28. ********************************************
  29. DEBPROC DHDT_pro H1*FLOTTANT T1*FLOTTANT ;
  30.  
  31. MESS 'Pour Frederic DEBUT';
  32. MESS '(BET*H1)' (BET*H1);
  33. MESS 'MOINSN' MOINSN;
  34. MESS 'Valeur =' (BET*H1 ** MOINSN);
  35. MESS ' ';
  36. MESS 'Il y a la une puissance negative FLOTTANT d un ';
  37. MESS 'FLOTTANT... au debut ==> 3.36 ** -2.0304';
  38. MESS 'Pour Frederic FIN';
  39. FIN; COMM 'A retirer';
  40.  
  41. DHDT0 = BET*H1 ** MOINSN + 1. * H1 / (T1 - TR) * N1 * M1 ;
  42. *MESS 'DHDT0' DHDT0;
  43. DHDT0 = 1./BET *N1*M1*1./(TSR) *
  44. (((T1 - TR / TSR) ** M1 - 1.)
  45. ** (N1 - 1.)) *((T1 - TR / TSR) ** (M1 - 1.)) ;
  46. *MESS 'DHDT0' DHDT0;
  47. FINPROC DHDT0 ;
  48. ********************************************
  49.  
  50.  
  51.  
  52. ********************************************
  53. DEBPROC KINPHI int00*FLOTTANT ;
  54. phi0 = 0. ;
  55. mess 'int00' int00;
  56. phit = (2* KS * PSI0) / (DT * int00) ;
  57. 'MESS' 'phit' phit ;
  58. * phi0 = 1.15;
  59. phi0 = phit;
  60. int0 = int00 - (phi0/2.); ;
  61. Int_p0 = prog int0 ;
  62. phi_p0 = prog phi0 ;
  63. repeter b1 (nb-1 );
  64. d0 = extr D2_prog &b1 ;
  65. phi0 = (2. * d0 / int0) + phi0 ;
  66. phi_p0 = phi_p0 et (prog phi0) ;
  67. int0 = int0 - phi0 ;
  68. Int_p0 = Int_p0 et (Prog int0) ;
  69. fin b1 ;
  70. d0 = extr D2_prog NB;
  71. phi0 = 2. * d0 / iNT0 + phi0 ;
  72. phi_p0 = phi_p0 et (prog phi0) ;
  73. * list phi_p0 ;
  74. * list int_p0 ;
  75. PNM1 = EXTR phi_p0 NB ;
  76. MESS 'PNM' PNM1 ;
  77. int = EXTR int_p0 NB ;
  78. MESS 'int' int ;
  79. FINPROC phi_p0 Int_p0 ;
  80. ********************************************
  81.  
  82.  
  83.  
  84. ********************************************
  85. DEBPROC KPROG T00*FLOTTANT TN0*FLOTTANT NB0*ENTIER MO1*MOT ;
  86. DT = T00 - TN0 / NB0 ;
  87.  
  88. SI (EGA MO1 DEMI ) ;
  89. MESS 'DEMI = ' MO1 ;
  90. T = T00 - (dt*0.5) ;
  91. h_pr0 = prog ;
  92. K_pr0 = prog ;
  93. D_pr0 = prog ;
  94. C_pr0 = prog ;
  95. T_pr0 = prog ;
  96. H = TH_pro T ;
  97. K = KKR_pro H ;
  98. DHDT = DHDT_pro H T ;
  99. D = K*DHDT ;
  100. C = 1./DHDT ;
  101. h_pr0 = prog H;
  102. K_pr0 = prog K;
  103. D_pr0 = prog D;
  104. C_pr0 = prog C;
  105. T_pr0 = prog T;
  106. T = T - DT ;
  107.  
  108. SINON ;
  109. MESS 'TOTA = ' MO1 ;
  110. T = T00 ;
  111. H = TH_pro T ;
  112. K = KKR_pro H ;
  113. DHDT = DHDT_pro H T ;
  114. D = K*DHDT ;
  115. C = 1./DHDT ;
  116. h_pr0 = prog H;
  117. K_pr0 = prog K;
  118. D_pr0 = prog D;
  119. C_pr0 = prog C;
  120. T_pr0 = prog T;
  121. T = T - DT ;
  122. FINSI;
  123.  
  124. repe b0 nb0 ;
  125. H = TH_pro T ;
  126. K = KKR_pro H ;
  127. DHDT = DHDT_pro H T ;
  128. D = K*DHDT ;
  129. C = 1./DHDT ;
  130. h_pr0 = h_pr0 ET (prog H) ;
  131. K_pr0 = K_pr0 ET (prog K) ;
  132. D_pr0 = D_pr0 ET (prog D) ;
  133. C_pr0 = C_pr0 ET (prog C) ;
  134. T_pr0 = T_pr0 ET (prog T);
  135. T = T - DT ;
  136. FIN B0 ;
  137. FINPROC h_pr0 K_pr0 D_pr0 C_pr0 T_pr0 ;
  138. ********************************************
  139.  
  140.  
  141.  
  142.  
  143. ********************************************
  144. DEBPROC A x*FLOTTANT ;
  145. si (x < 2.9 ) ;
  146. X2 = (-1.)*X*X ;
  147. A0 = X2 EXP * X * (1. - (ERF X) ** (-1.)) * RA2PI + (2.*X2) ;
  148. sinon ;
  149. X2 = X ** (-2.) * -0.5 ;
  150. A0 = 108830.*X2+8162.*X2+706.*X2+74.*X2+10.*X2+2.*X2+1.;
  151. *A0 = 8162.*X2+706.*X2+74.*X2+10.*X2+2.*X2+1.;
  152. *A0 = 706.*X2+74.*X2+10.*X2+2.*X2+1.;
  153. *A0 = 74.*X2+10.*X2+2.*X2+1.;
  154. FINSI ;
  155. FINPROC A0 ;
  156. ********************************************
  157.  
  158.  
  159.  
  160. ********************************************
  161. DEBPROC INTN phi*LISTREEL ;
  162. PNM1 = EXTR phi NB ;
  163. DD0 = EXTR D2_PROG NB ;
  164. XX = PNM1 * 0.5D0 * (DD0 ** -0.5D0) ;
  165. IC0 = PNM1*0.5D0 + (2.D0*DD0/PNM1*(A XX)) ;
  166. FINPROC IC0 ;
  167. ********************************************
  168.  
  169.  
  170.  
  171. 'DEBPROC' CALPHILI ;
  172.  
  173. h_prog K_prog D_prog C_prog T_prog = KPROG T0 TN NB 'TOTA';
  174.  
  175. dinf = extr D_prog 1 ;
  176. mess 'dinf' dinf;
  177. dinf = extr D_prog 2 ;
  178. mess 'dinf2' dinf;
  179. T1 = extr T_prog 1 ;
  180. mess 'T1' T1;
  181. TN1 = extr T_prog (DIME T_prog) ;
  182. T01 = extr T_prog 1 ;
  183. db0 = 2. / ((T01 - TN1) ** 2.) ;
  184.  
  185. Y = T_prog - (prog (DIME T_prog)*TN1) * K_prog ;
  186. *-------------------------- Dmoy par K
  187. Db = (-1.) * ((INTG KfH_evol) extr 1) * db0 ;
  188. *-------------------------- calcul de I 1/2 (j = 1)
  189. *????????????????????????????????????
  190. Int = Db/pi**0.5*2.*nb ; 'MESS' 'Int' Int ;
  191. *Int = 1.e-3/(ts - tn)*nb ;list Int ;
  192. *Int = int * 2. ;
  193. *opti donn 5 ;
  194.  
  195. *-------------------------- calcul de D (I + 1/2)
  196. h2_prog K2_prog D2_prog C2_prog T2_prog = KPROG T0 TN nb DEMI ;
  197.  
  198.  
  199. *-------------------------- calcul de PHI (I + 1/2)
  200. phi_prog Int_prog = KINPHI Int ;
  201.  
  202.  
  203. I2_evol = evol manu T_prog phi_prog ;
  204. I2mi = I2_evol INTG * (-1.) / dt ; LIST I2mi ;
  205.  
  206. *I2_evol = evol manu phi_prog T_prog ;
  207. *DESS I2_evol ;
  208. *opti donn 5 ;
  209.  
  210. *-------------------------- calcul de ICHA (N - 1/2)
  211. ICHA = INTN phi_prog ;
  212.  
  213. *-------------------------- DELTA (J=1)
  214. DELTA = (EXTR Int_prog NB) - ICHA ;
  215.  
  216. *-------------------------- calcul de I 1/2 (j = 2)
  217. int1 = int ;
  218. Int = int1 - (DELTA * 0.5d0) ;
  219. deltI = Int - Int1 ;
  220.  
  221. repeter b2 20 ;
  222. *-------------------------- calcul de PHI (I + 1/2)
  223. phi_prog Int_prog = KINPHI Int ;
  224.  
  225. 'MESSAGE' 'valeur de xphi';
  226. 'LISTE' ('EXTRAIRE' 1 phi_prog);
  227.  
  228. I2_evol = evol manu T_prog phi_prog ;
  229. I2mi = I2_evol INTG * (-1.) / dt ;
  230. *-------------------------- calcul de ICHA (N - 1/2)
  231. ICHA = INTN phi_prog ;
  232. 'iter :' &b2 'icha' ICHA 'DELTA' DELTA 'Itot' (maxi I2mi) 'int' int;
  233.  
  234. *-------------------------- DELTA (J=1)
  235. delta1 = DELTA ;
  236. DELTA = (EXTR Int_prog NB) - ICHA ;
  237.  
  238. *-------------------------- calcul de I 1/2
  239. int1 = int ;
  240. *int = int1 - (DELTA * delta1 / (DELTA + delta1)) ;
  241. *int = int1 - (DELTA * deltI / (DELTA + delta1)) ;
  242. Int = int1 - (DELTA * 0.5d0 ) ;
  243. deltI = Int - Int1 ;
  244. si ( (DELTA abs) < (ICHA*1.e-4))
  245. quitter b2 ;
  246. finsi;
  247. menage ;
  248. fin b2 ;
  249.  
  250.  
  251. 'FINPROC' phi_prog h_prog T_prog ;
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262. 'MESS' 'Calcul de la solution de Philip en milieu saturé non saturé';
  263. xf = 0.0;
  264. 'MESS' 'xf (frontiére saturé / non saturé) ' xf;
  265. PSI0 = 1.D0;
  266.  
  267. 'MESS' 'Quelle est la pression d injection dans le milieu saturé?' ;
  268. mess 'PSI0' PSI0;
  269.  
  270. ALF = 0.08 * 100. ;
  271. KS = 5.85E-2 / 3600. ;
  272. TS = 0.3 ;
  273. TR = 0.055 ;
  274. TSR = TS - TR ;
  275. N = 2.0304 ; N1 = 1./N ; MOINSN = (-1.)*N ;
  276. M = -0.5075 ; M1 = 1./M ;
  277. BET = -0.029227 * 100. ;
  278.  
  279. nb = 500 ;
  280. H00 = -1.E-20 ;
  281. *H00 = 0;
  282. Hinf = -1.15D0 ;
  283.  
  284. RA2PI = PI ** (-0.5) * 2. ;
  285.  
  286. phi0 = 0. ;
  287.  
  288. HN = HINF ;
  289.  
  290. TN = HHT_pro HN ; LIST TN ;
  291.  
  292. HN = TH_pro TN ; LIST HN ;
  293.  
  294. KN = KKR_pro HN ; LIST KN ;
  295.  
  296. DHDTN = DHDT_pro HN TN ; LIST DHDTN ;
  297.  
  298. DN = KN*DHDTN ;
  299.  
  300. CN = 1./DHDTN ;
  301.  
  302. T0 = HHT_pro H00 ;
  303.  
  304. LIST (HN - (TH_pro (TN - 1.E-6)) / 1.E-6) ;
  305.  
  306. DT = T0 - TN / NB ;
  307.  
  308.  
  309.  
  310.  
  311. LTEMPS = 'PROG' 1000. 'PAS' 1000. (1000. * 3.) ;
  312. ntemps = 'DIME' LTEMPS ;
  313. lastemps = 0. ;
  314. TETSAT = TS ;
  315. TETNSAT = HHT_pro HINF ;
  316. OLDXDTET = (TETSAT - TETNSAT) * xf ;
  317.  
  318. *
  319. *
  320. * CALCUL DE PHI, H et TETA
  321. *
  322. phi_prog h_prog t_prog = CALPHILI ;
  323.  
  324.  
  325.  
  326. xfini = xf;
  327. 'REPETER' BTEMPS ntemps ;
  328. ttemps = 'EXTR' LTEMPS &BTEMPS ;
  329. DELTEMP = ttemps - lastemps ;
  330. * profile de la charge dans la zone saturée
  331. LXSAT = 'PROG' 0. ;
  332. LHSAT = 'PROG' PSI0 ;
  333. LTETSAT = 'PROG' TETSAT ;
  334.  
  335.  
  336. LXNSAT = phi_prog * ((TTEMPS)** 0.5) ;
  337.  
  338. ID = 'PROG' ('DIME' LXNSAT)*1. ;
  339. phif = EXTR phi_prog 1;
  340. mess 'phif' phif;
  341. xfr = xfini + (phif*(ttemps**0.5));
  342. mess 'xfr= 'xfr;
  343. LXNSAT = LXNSAT + (ID*xfini) ;
  344. LHNSAT = h_prog ;
  345. LTETNSAT = t_prog ;
  346. *
  347. LX = LXSAT ET LXNSAT ;
  348. LTET = LTETSAT ET LTETNSAT ;
  349. LH = LHSAT ET LHNSAT ;
  350.  
  351. 'SI' (&BTEMPS 'EGA' 1) ;
  352. * (LH = 'ORDO' LH CROISSANT);
  353. uun = 'PROG' ('DIME' LX) * 1.D0;
  354. courb = ('EVOL' 'MANU'(uun '-' LX) LH);
  355. 'SINON' ;
  356. * (LH = 'ORDO' LH CROISSANT) ;
  357. courb = courb 'ET' ('EVOL' 'MANU' (uun '-' LX) LH);
  358. 'FINSI' ;
  359. *
  360. 'FIN' BTEMPS ;
  361.  
  362.  
  363. TAB1 = 'TABLE' ;
  364. TAB1 . 'TITRE' = 'TABLE' ;
  365. TAB1 . 'TITRE' . 1 = 'exacte 1000' ;
  366. TAB1 . 'TITRE' . 2 = 'exacte 2000' ;
  367. TAB1 . 'TITRE' . 3 = 'exacte 3000' ;
  368. TAB1 . 1 = 'TIRR ';
  369. TAB1 . 2 = 'TIRR ';
  370. TAB1 . 3 = 'TIRR ';
  371.  
  372.  
  373. 'SI' (GRAPH);
  374. dess (courb) LEGE TAB1;
  375. 'FINSI' ;
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384. ******************** CAS TEST : cacul.dgibi ************************
  385. *
  386. *
  387. * Test de fonctionnement de DARCYSAT en 2D sans effet de gravité.
  388. * ===================================================================
  389. * Infiltration d'eau dans une colonne horizontale de sable uniformément
  390. * désaturé.
  391. *
  392. * - condition initiale : désaturation uniforme correspondant à une
  393. * succion de 1,15 m ;
  394. * - a l'instant initial, une extrémité de la colonne est noyée.
  395. * La frontière est supposée rester a pression nulle ;
  396. *
  397. * ===================================================================
  398. * Les options de modélisation declarées dans la table transmise à
  399. * la procédure DARCYSAT sont les suivantes :
  400. *
  401. *
  402. * - les effets gravitationnels ne sont pas pris en compte : pression =
  403. * charge (indice GRAVITE absent ; valeur par défaut : FAUX);
  404. *
  405. * - Une liste de temps de sauvegarde est fournie en valeur de l'indice
  406. * TEMPS_SAUVES ;
  407. *
  408. * - Le pas de temps est d'abord automatique (indice TEMPS_CALCULES absent)
  409. * L'utilisateur fournit
  410. * > le pas de temps initial (indice 'DT_INITIAL'),
  411. * > le nombre d'itérations recherché par pas de temps (indice 'NITER')
  412. * > le nombre de pas de temps (indice 'NPAS')
  413. *
  414. * - Le calcul est ensuite refais avec une liste de temps de calcul
  415. * donnée à l'indice TEMPS_CALCULES ;
  416. *
  417. * - L'homogenéisation spatiale des propriétes physiques s'effectue
  418. * par moyenne arithmétique des valeurs obtenues aux faces (indice
  419. * HOMOGENEISATION de valeur DECENTRE)
  420. *
  421. *!!!!!!!!!!!!!!!!!!!!!!!!
  422. * - La précision de convergence demandée est de 0,5 mm (indice RESIDU_MAX)
  423. * cas test tiré du rapport DMT 97/25 :
  424. * "Implémentation dans CASTEM 2000 d'un modèle de transfert hydrique
  425. * en milieu poreux non saturé"
  426. *
  427. *
  428. ********************* CAS TEST : cacul.dgibi ************************
  429.  
  430. 'OPTION' 'ECHO' 1 ;
  431. 'SAUTER' 'PAGE';
  432. *
  433.  
  434. 'TITRE' 'infiltration horizontale dans le sable : cacul.dgibi' ;
  435. 'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
  436. *'OPTION' 'ISOV' 'LIGN' ;
  437. *'TRACER' 'PSC' ;
  438.  
  439. ************************************************************************
  440. * Fonction qui calcule la perméabilité en fonction de la
  441. * saturation réduite.
  442. *
  443. * Loi puissance :
  444. * Perméabilité : K = Ks S^B
  445. *
  446. * Paramètres physiques à définir dans la table LOI :
  447. * --------------------------------------------------
  448. * ALPHA : coef. B (s.d.)
  449. * PERMSAT : coef. Ks, perméabilité à saturation (m/s)
  450. *
  451. * Entrée :
  452. * --------
  453. * LOI : table de données décrite ci-dessus
  454. * SAT : saturation réduite
  455. * TEST : mot facultatif 'NOTEST' qui permet de shunter les tests.
  456. *
  457. * Sortie :
  458. * --------
  459. * K1 : perméabilité totale en eau (m/s)
  460. *
  461. ************************************************************************
  462.  
  463.  
  464. *--------------------------------------------------------------------
  465. * Création maillage
  466. *
  467. *- Discrétisation :
  468. ENX = 1 ;
  469. ENY = 50 ;
  470. 'DENS' (1./eny) ;
  471.  
  472. *- Création des points et des droites
  473. *
  474. A0 = -0.5 1.0D0; B0 = 0.5 1.0D0;
  475. 'DENS' (1./eny) ;
  476. A1 = -0.5 0.5D0; B1 = 0.5 0.5D0;
  477. A2 = -0.5 0.25D0 ; B2 = 0.5 0.25D0 ;
  478. A3 = -0.5 0.0D0; B3 = 0.5 0.0D0 ;
  479.  
  480. *- Création des droites
  481. *
  482. AB0 = 'DROIT' ENX A0 B0 ;
  483. AB1 = 'DROIT' ENX A1 B1 ;
  484. AB2 = 'DROIT' ENX A2 B2 ;
  485. AB3 = 'DROIT' ENX A3 B3 ;
  486.  
  487. *- Creation des surfaces
  488. *
  489. MASSIF0 = AB3 'REGLER' AB2 'REGLER' AB1 'REGLER' AB0 ;
  490. ENY = 'NBEL' MASSIF0 ;
  491. ENTR = 'COULEUR' ('INVERSE' AB0) 'ROUGE' ;
  492. SORT = 'COULEUR' (AB3) 'ROUGE' ;
  493.  
  494. 'SI' GRAPH ;
  495. 'TRACER' (MASSIF0 'ET' ENTR 'ET' SORT) ;
  496. 'FINSI' ;
  497.  
  498. *- Creation des maillages contenant tous les points
  499. *
  500. QFTOT = 'CHANGER' MASSIF0 'QUAF' ;
  501. QFSORT = 'CHANGER' SORT 'QUAF' ;
  502. QFENTR = 'CHANGER' ENTR 'QUAF' ;
  503.  
  504.  
  505. 'ELIMINATION' 0.00001 (QFTOT 'ET' QFSORT 'ET' QFENTR) ;
  506.  
  507. *- Modèles
  508. *
  509. MODHYB = 'MODELE' QFTOT 'DARCY' 'ISOTROPE' ;
  510. MODENTR = 'MODELE' QFENTR 'DARCY' 'ISOTROPE' ;
  511. MODSORT = 'MODELE' QFSORT 'DARCY' 'ISOTROPE' ;
  512. CEENTR = 'DOMA' MODENTR 'CENTRE' ;
  513. CESORT = 'DOMA' MODSORT 'CENTRE' ;
  514. HYCEN = 'DOMA' MODHYB 'CENTRE' ;
  515. HYFAC = 'DOMA' MODHYB 'FACE';
  516.  
  517. *- Création ligne de suivi pour le post-traitement et le test
  518. * ligne des points centres (cas 1D)
  519. *
  520. 'REPETER' BCL (ENY - 1) ;
  521. IP = &BCL ;
  522. JP = IP + 1 ;
  523. PI = 'POINT' HYCEN IP ;
  524. PJ = 'POINT' HYCEN JP ;
  525. 'SI' (IP 'EGA' 1);
  526. LCENC = ('MANU' 'SEG2' PI PJ) ;
  527. 'SINON' ;
  528. LCENC = LCENC 'ET' ('MANU' 'SEG2' PI PJ) ;
  529. 'FINSI' ;
  530. FIN BCL ;
  531.  
  532. *--------------------------------------------------------------------
  533. *- pression initiale (metre d'eau) dans le sable
  534. HN = -1.15 ;
  535.  
  536. *--------------------------------------------------------------------
  537. *- Conditions aux limites
  538.  
  539.  
  540. * GBM modifie supprime bloque
  541. *- frontière en limite du domaine de calcul (milieu désaturé)
  542. ESORT = 'MANU' 'CHPO' CESORT 1 'TH' HN 'NATURE' 'DISCRET';
  543.  
  544. *- frontière mouillée
  545. EENTR = 'MANU' 'CHPO' CEENTR 1 'TH' 1. 'NATURE' 'DISCRET';
  546.  
  547.  
  548. *- chargement des CLs
  549. LICALC = 'PROG' 0.D0 1.e20 ;
  550. LIST1 = 'PROG' 2 * 1.D0 ;
  551. VALI0 = 'CHAR' 'THIM' (ESORT et EENTR) ('EVOL' 'MANU' LICALC LIST1) ;
  552.  
  553. *--------------------------------------------------------------------
  554. *- initialisation des inconnues
  555. * (doit être compatible avec les conditions aux limites)
  556. *
  557.  
  558. *- trace de charge d'eau
  559. *TH0 = 'MANU' 'CHPO' ('DIFF' CEENTR HYFAC) 1 'TH' HN 'NATURE' 'DISCRET';
  560. *TH0 = TH0 + SIMPE ;
  561.  
  562. *- charge d'eau
  563. H0 = 'MANU' 'CHPO' HYCEN 1 'H' HN 'NATURE' 'DIFFUS';
  564. H0 = 'MASQUE' ('COORDONNEE' 2 ('DOMA' modhyb centre))
  565. SUPERIEUR 0.98D0;
  566. H0 = 1 '-' H0;
  567. H0 = (HN * H0) + ((1.D0 - H0) * 1.D0);
  568. H0 = 'NOMC' 'H' H0 NATURE DISCRET;
  569. 'TRACER' ('KCHA' MODHYB H0 'CHAM') modhyb;
  570.  
  571.  
  572. *- flux
  573. QFACE0 = 'MANU' 'CHPO' HYFAC 1 'FLUX' 0.D0 'NATURE' 'DISCRET' ;
  574.  
  575.  
  576. * ---------------------------
  577. * = Table DARCY_TRANSITOIRE =
  578. * ---------------------------
  579. *- initialisation table
  580. SATUR = 'TABLE' ;
  581. SATUR. 'TEMPS' = 'TABLE' ;
  582. SATUR. 'CHARGE' = 'TABLE' ;
  583. SATUR. 'FLUX' = 'TABLE' ;
  584.  
  585. SATUR . 'ITMAX' = 45;
  586.  
  587.  
  588. *- données géommétriques
  589. SATUR. 'SOUSTYPE' = 'DARCY_TRANSATUR' ;
  590. SATUR. 'MODELE' = MODHYB ;
  591.  
  592. *- instant initial
  593. SATUR. 'TEMPS' . 0 = 0. ;
  594. SATUR. 'CHARGE' . 0 = 'COPIER' H0 ;
  595. SATUR. 'FLUX' . 0 = 'COPIER' QFACE0 ;
  596.  
  597. *********************************************
  598. * GBM
  599. *- conditions aux limites et chargements
  600. *SATUR. 'BLOCAGES_DARCY' = BSORT 'ET' BENTR ;
  601. *SATUR. 'CHARGEMENT' = VALI0 ;
  602. SATUR . 'TRACE_IMPOSE' = VALI0 ;
  603.  
  604. SATUR . 'LUMP' = FAUX;
  605. SATUR . 'TYPDISCRETISATION' = 'VF' ;
  606.  
  607.  
  608.  
  609. * GBM MODIFIE
  610.  
  611. TABRES = table METHINV;
  612. TABRES . 'TYPINV' = 1;
  613. TABRES . 'PRECOND' = 5;
  614.  
  615. SATUR . 'METHINV' = TABRES;
  616. SATUR . 'METHINV' . RESID = 1.D-20;
  617.  
  618. *******************************************
  619.  
  620.  
  621.  
  622. *******************************************
  623.  
  624.  
  625.  
  626. *- données physiques
  627. * loi de perméabilité
  628. LoiP = 'TABLE' 'PUISSANCE';
  629. LoiP. 'ALPHA' = 8. ;
  630. LoiP. 'PERMSAT' = 5.85E-2 / 3600. ;
  631. SATUR.'LOI_PERMEABILITE' = 'TABLE' LoiP ;
  632.  
  633. * loi de succion
  634. LoiS = 'TABLE' 'VAN_GENUCHTEN';
  635. LoiS. 'PORO' = 0.3 ;
  636. LoiS. 'TERESIDU' = 0.055 ;
  637. LoiS. 'NEXP' = 2.0304 ;
  638. LoiS. 'MEXP' = 0.5075 ;
  639. LoiS. 'BHETA' = 0.029227 * 100. ;
  640. SATUR.'LOI_SATURATION' = 'TABLE' LoiS ;
  641.  
  642. SATUR . 'COEF_EMMAGASINEMENT' = 0.D-6;
  643. *SATUR . 'XI' = 1.D-5;
  644.  
  645. *- données numériques
  646. SATUR. 'TEMPS_FINAL' = 2000.D0 ;
  647. SATUR. 'HOMOGENEISATION' = 'CHAINE' 'DECENTRE' ;
  648.  
  649.  
  650. SATUR.'SOUS_RELAXATION' = 1.;
  651. SATUR. 'NPAS' = 10000 ;
  652. SATUR. 'RESIDU_MAX' = 1.D-4 ;
  653. SATUR. 'NITER' = 10 ;
  654. SATUR. 'DT_INITIAL' = 0.1D0 ;
  655. SATUR. 'TEMPS_SAUVES' = ('PROG' 1 'PAS' 1 2.)*1000. ;
  656. *SATUR. 'TEMPS_CALCULES' = ('PROG' 1 'PAS' 1 3.)*1000. ;
  657. SATUR. 'MESSAGES' = 2 ;
  658.  
  659.  
  660.  
  661. *--------------------------------------------------------------------
  662. *- Vérification du choix du dp minimum pour le calcul de la capacité.
  663. * droite support des variables
  664. dx = 'DROIT' (0. -1.15 ) 1000 (0. 0.) ;
  665. zc = 'COOR' dx 2 ;
  666. ev2 = 'EVOL' 'BLEU' 'CHPO' zc 'SCAL' dx ;
  667. px = 'EXTR' ev2 'ORDO' ;
  668. * calcul de la teneur en eau pour la pression zc
  669. s0 t0 cap = HT_PRO (SATUR.'LOI_SATURATION') ZC ;
  670. ev0 = 'EVOL' 'TURQ' 'CHPO' s0 'SCAL' dx ;
  671. evt = 'EVOL' 'VERT' 'MANU' px (100. * ('EXTR' ev0 'ORDO')) ;
  672.  
  673. 'SI' GRAPH ;
  674. 'DESSIN' evt 'TITX' 'Pc(m)' 'TITY' 'S(%)'
  675. 'TITRE' 'Loi capillaire S(Pc)' ;
  676. 'FINSI' ;
  677.  
  678. * calcul de la teneur en eau pour la pression zc - dp
  679. dp = 1.e-4 ;
  680. s1 t1 cap = HT_PRO (SATUR.'LOI_SATURATION')
  681. ('KOPS' zc '-' dp) ;
  682. * représentation de la capacité
  683. c1 = (t0 - t1) / dp;
  684. ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
  685. evc = 'EVOL' 'JAUNE' 'MANU' px ('EXTR' ev1 'ORDO') ;
  686.  
  687. 'SI' GRAPH ;
  688. 'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
  689. 'TITRE' 'Capacite capillaire' ;
  690. 'FINSI' ;
  691.  
  692. * ===========================
  693. * | 1er CALCUL |
  694. * ===========================
  695. *- fonctionnement avec temps automatiques
  696.  
  697.  
  698.  
  699. DARCYSAT SATUR ;
  700.  
  701.  
  702. * Post-traitement
  703. * ===============
  704.  
  705. LT = 'LECT' 0 PAS 1 ((dime SATUR.TEMPS) - 1) ;
  706. liopt = 'MOTS' 'MIMA' 'AXES';
  707. TDES = TRACHIS SATUR 'CHARGE' LCENC LT 'PREF' ' ' 'UNIT' 's' ;
  708.  
  709.  
  710.  
  711.  
  712. * le fichier lu est produit par resanalyphilhorsatnsat.dgibi
  713.  
  714.  
  715. LT = 'LECT' 1 'PAS' 1 (('DIME' SATUR.'TEMPS') - 1) ;
  716. liopt = 'MOTS' 'MIMA' 'AXES';
  717. TDES = TRACHIS SATUR 'CHARGE' LCENC LT 'PREF' ' ' 'UNIT' 's' ;
  718. TDES . 3 = 'TABLE' ;
  719. TDES . 3 . VALEUR = 'EXTRAIRE' courb 'COUR' 1;
  720. TDES . 3 . LEGEND2 = TAB1 . TITRE . 1;
  721. TDES . 3 . LEGEND1 = ' ' ;
  722. TDES . 4 = 'TABLE' ;
  723. TDES . 4 . VALEUR = 'EXTRAIRE' courb 'COUR' 2;
  724. TDES . 4 . LEGEND2 = TAB1 . TITRE . 2;
  725. TDES . 4 . LEGEND1 = ' ' ;
  726. * TDES . 5 = 'TABLE' ;
  727. * TDES . 5 . VALEUR = 'EXTRAIRE' courb 'COUR' 3;
  728. * TDES . 5 . LEGEND2 = TAB1 . TITRE . 3;
  729. * TDES . 5 . LEGEND1 = ' ' ;
  730.  
  731. 'SI' GRAPH ;
  732. DESTRA TDES liopt 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
  733. 'FINSI' ;
  734.  
  735.  
  736. toto = 'PRIM' TDES . 2 . VALEUR;
  737. inta1 = 'EXTRAIRE' toto 'ORDO';
  738. inta1 = 'EXTRAIRE' ('DIME' inta1) inta1;
  739.  
  740. 'LISTE' inta1;
  741.  
  742. err1 = 'ABS' (inta1);
  743. 'LISTE' err1;
  744.  
  745. 'SI' ((err1 '-' 2.95D-2) > 1.D-3) ;
  746. 'ERREUR' 5;
  747. 'FINSI' ;
  748.  
  749.  
  750. *opti donn 5;
  751.  
  752. * ===========================
  753. * | 2nd CALCUL |
  754. * ===========================
  755. *- fonctionnement avec une liste de temps calculés
  756. SATUR. 'TEMPS' = 'TABLE' ;
  757. SATUR. 'TRACE_CHARGE' = 'TABLE' ;
  758. SATUR. 'CHARGE' = 'TABLE' ;
  759. SATUR. 'FLUX' = 'TABLE' ;
  760. SATUR. 'TEMPS' . 0 = 0. ;
  761. SATUR. 'TRACE_CHARGE' = TABLE ;
  762. SATUR. 'CHARGE' . 0 = H0 ;
  763. SATUR. 'FLUX' . 0 = QFACE0 ;
  764. SATUR . 'TYPDISCRETISATION' = 'EFMH' ;
  765.  
  766.  
  767.  
  768.  
  769. DARCYSAT SATUR ;
  770.  
  771. * Post-traitement
  772. * ===============
  773. *-- Test de non régression
  774. * Le test est effectué en vérifiant la solution
  775.  
  776. *- solution après 10 pas de temps
  777.  
  778.  
  779. *-- Tracé (tous les temps)
  780.  
  781. LT = 'LECT' 0 'PAS' 1 (('DIME' SATUR.'TEMPS') - 1) ;
  782. liopt = 'MOTS' 'MIMA' 'AXES';
  783. TDES = TRACHIS SATUR 'CHARGE' LCENC LT 'PREF' ' ' 'UNIT' 's' ;
  784. LT = 'LECT' 1 'PAS' 1 (('DIME' SATUR.'TEMPS') - 1) ;
  785. liopt = 'MOTS' 'MIMA' 'AXES';
  786. TDES = TRACHIS SATUR LCENC 'CHARGE' LT 'PREF' ' ' 'UNIT' 's' ;
  787. TDES . 3 = 'TABLE' ;
  788. TDES . 3 . VALEUR = 'EXTRAIRE' courb 'COUR' 1;
  789. TDES . 3 . LEGEND2 = TAB1 . TITRE . 1;
  790. TDES . 3 . LEGEND1 = ' ' ;
  791. TDES . 4 = 'TABLE' ;
  792. TDES . 4 . VALEUR = 'EXTRAIRE' courb 'COUR' 2;
  793. TDES . 4 . LEGEND2 = TAB1 . TITRE . 2;
  794. TDES . 4 . LEGEND1 = ' ' ;
  795. 'SI' GRAPH ;
  796. DESTRA TDES liopt 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
  797. 'FINSI' ;
  798.  
  799.  
  800.  
  801. toto = 'PRIM' TDES . 2 . VALEUR;
  802. inta2 = 'EXTRAIRE' toto 'ORDO';
  803. inta2 = 'EXTRAIRE' ('DIME' inta2) inta2;
  804.  
  805. 'MESSAGE' inta1 inta2;
  806.  
  807. err2 = 'ABS' (inta2);
  808. 'LISTE' err2;
  809.  
  810. 'SI' ((err2 '-' 1.22D-2) > 1.D-3) ;
  811. 'ERREUR' 5;
  812. 'FINSI' ;
  813.  
  814. 'MESSAGE' err1 err2;
  815.  
  816. 'FIN';
  817.  
  818.  
  819.  
  820.  
  821.  
  822.  
  823.  
  824.  
  825.  
  826.  
  827.  
  828.  
  829.  
  830.  

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