Télécharger caculVF.dgibi

Retour à la liste

Numérotation des lignes :

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

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