Télécharger caculVFconservatif.dgibi

Retour à la liste

Numérotation des lignes :

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

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