Télécharger la plaquette Cast3M


fr en

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

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