Télécharger la plaquette Cast3M


fr en

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 ;
  282. 'iter :' &b2 'icha' ICHA 'DELTA' DELTA 'Itot' (maxi I2mi) 'int' int;
  283.  
  284. *-------------------------- DELTA (J=1)
  285. delta1 = DELTA ;
  286. DELTA = (EXTR Int_prog NB) - ICHA ;
  287.  
  288. *-------------------------- calcul de I 1/2    
  289. int1 = int ;
  290. *int = int1 - (DELTA * delta1 / (DELTA + delta1)) ;
  291. *int = int1 - (DELTA * deltI / (DELTA + delta1)) ;
  292. Int = int1 - (DELTA * 0.5d0 ) ;
  293. deltI = Int - Int1 ; 
  294. si ( (DELTA abs) < (ICHA*1.e-6))
  295. quitter  b2 ;
  296. finsi;
  297. menage ;
  298. fin b2 ;
  299.  
  300. *'LISTE' h_prog;
  301. *dess ( evol manu phi_prog T_prog);
  302. *dess ( evol manu phi_prog h_prog);
  303.  
  304.  
  305.  
  306. nlist = 'DIME' h_prog;
  307. nnlist = 2 * nlist;
  308. xabs = 'PROG' 0.D0 PAS (1.D0 '/' nnlist) 1.D0;
  309. * on inverse l'axe des x
  310. *'ORDONNER' CROISSANT  h_prog ;
  311. *'LISTE' h_prog;
  312. * il faut tronquer la derniere valeur de phi_prog et mettre
  313. * un grand nombre à la place.
  314.  
  315. roro = 'ENLEVER' phi_prog nlist ;
  316. roro = roro 'ET' ('PROG' 10.D0);
  317.  
  318. tabcool = 'TABLE' ;
  319. tabcool . 1 = 'TURQ';
  320. tabcool . 2 = 'TURQ';
  321. tabcool . 3 = 'TURQ';
  322. tabcool . 4 = 'TURQ';
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329. ***********************************************************************************
  330.  
  331.  
  332. *           LE CAS TEST
  333.  
  334.  
  335. **********************************************************************************
  336.  
  337.  
  338.  
  339.  
  340.  
  341. ******************** CAS TEST : cacul.dgibi ************************
  342. *
  343. *
  344. * Test de fonctionnement de DARCYSAT en 2D sans effet de gravité.
  345. * ===================================================================
  346. * Infiltration d'eau dans une colonne horizontale de sable uniformément
  347. * désaturé.
  348. *
  349. * -  condition initiale : désaturation uniforme correspondant à une
  350. *    succion de 1,15 m ;
  351. * -  a l'instant initial, une extrémité de la colonne est noyée.
  352. *    La frontière est supposée rester a pression nulle ;
  353. *
  354. * ===================================================================
  355. *    Les options de modélisation declarées dans la table transmise à
  356. *    la procédure DARCYSAT sont les suivantes :
  357. *
  358. *
  359. * -  les effets gravitationnels ne sont pas pris en compte : pression =
  360. *    charge (indice GRAVITE absent ; valeur par défaut : FAUX);
  361. *
  362. * -  Une liste de temps de sauvegarde est fournie en valeur de l'indice
  363. *    TEMPS_SAUVES ;
  364. *
  365. * -  Le pas de temps est d'abord automatique (indice TEMPS_CALCULES absent)
  366. *    L'utilisateur fournit
  367. *    > le pas de temps initial (indice 'DT_INITIAL'),
  368. *    > le nombre d'itérations recherché par pas de temps (indice 'NITER')
  369. *    > le nombre de pas de temps (indice 'NPAS')
  370. *
  371. * -  Le calcul est ensuite refais avec une liste de temps de calcul
  372. *    donnée à l'indice TEMPS_CALCULES ;
  373. *
  374. * -  L'homogenéisation spatiale des propriétes physiques s'effectue
  375. *    par moyenne arithmétique des valeurs obtenues aux faces (indice
  376. *    HOMOGENEISATION de valeur DECENTRE)
  377. *
  378.  
  379. * cas test tiré du rapport DMT 97/25 :
  380. * "Implémentation dans CASTEM 2000 d'un modèle de transfert hydrique
  381. *  en milieu poreux non saturé"
  382. *
  383. *
  384. ********************* CAS TEST : cacul.dgibi ************************
  385.  
  386. 'OPTION' 'ECHO' 1 ;
  387. 'SAUTER' 'PAGE';
  388. *
  389.  
  390. 'TITRE' 'infiltration horizontale dans le sable : cacul.dgibi'   ;
  391. 'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
  392. *'OPTION' 'ISOV' 'LIGN' ;
  393. *'TRACER' 'PSC' ;
  394.  
  395. ************************************************************************
  396. *  Fonction qui calcule la perméabilité en fonction de la
  397. *  saturation réduite.
  398. *
  399. *  Loi puissance :
  400. *  Perméabilité : K = Ks S^B
  401. *
  402. *  Paramètres physiques à définir dans la table LOI :
  403. *  --------------------------------------------------
  404. *    ALPHA   : coef. B (s.d.)
  405. *    PERMSAT : coef. Ks, perméabilité à saturation (m/s)
  406. *
  407. *  Entrée :
  408. *  --------
  409. *    LOI     : table de données décrite ci-dessus
  410. *    SAT     : saturation réduite
  411. *
  412. *  Sortie :
  413. *  --------
  414. *    K1      : perméabilité totale en eau (m/s)
  415. *
  416. ************************************************************************
  417.  
  418. ************************************************************************
  419.  
  420.  
  421. *--------------------------------------------------------------------
  422. * Création maillage
  423. *
  424. *- Discrétisation :
  425. ENX = 1      ;
  426. ENY = 200    ;
  427. 'DENS' (1./eny) ;
  428.  
  429. *- Création des points et des droites
  430. *
  431. A0 = -0.5   1.0D0;   B0 = 0.5   1.0D0;
  432.  
  433. dens (4./eny);
  434. A1 = -0.5   0.75D0;   B1 = 0.5   0.75D0;
  435. A2 = -0.5   0.5D0 ;  B2 = 0.5   0.5D0 ;
  436. A3 = -0.5   0.0D0;   B3 = 0.5   0.0D0 ;
  437.  
  438. *- Création des droites
  439. *
  440. AB0 = 'DROIT' ENX A0 B0 ;
  441. AB1 = 'DROIT' ENX A1 B1 ;
  442. AB2 = 'DROIT' ENX A2 B2 ;
  443. AB3 = 'DROIT' ENX A3 B3 ;
  444.  
  445. *- Creation des surfaces
  446. *
  447. MASSIF0 = AB3 'REGLER' AB2 'REGLER' AB1;
  448. MASSIF0 = MASSIF0 'REGLER' AB0 ;
  449. ENY  = 'NBEL' MASSIF0 ;
  450. ENTR = 'COULEUR' ('INVERSE' AB0) 'ROUGE' ;
  451. SORT = 'COULEUR' (AB3) 'ROUGE' ;
  452.  
  453. 'SI' GRAPH  ;
  454.   'TRACER' (MASSIF0 'ET' ENTR 'ET' SORT) ;
  455. 'FINSI' ;
  456.  
  457. *- Creation des maillages contenant tous les points
  458. *
  459. QFTOT  = 'CHANGER' MASSIF0 'QUAF' ;
  460. QFSORT = 'CHANGER' SORT    'QUAF' ;
  461. QFENTR = 'CHANGER' ENTR    'QUAF' ;
  462.  
  463.  
  464. 'ELIMINATION' 0.00001 (QFTOT 'ET' QFSORT 'ET' QFENTR) ;
  465.  
  466. *- Modèles
  467. *
  468. MODHYB  = 'MODELE' QFTOT  'DARCY' 'ISOTROPE'        ;
  469. MODENTR = 'MODELE' QFENTR 'DARCY' 'ISOTROPE'        ;
  470. MODSORT = 'MODELE' QFSORT 'DARCY' 'ISOTROPE'        ;
  471. CEENTR  = 'DOMA'  MODENTR 'CENTRE' ;
  472. CESORT  = 'DOMA'  MODSORT 'CENTRE' ;
  473. HYCEN   = 'DOMA'  MODHYB  'CENTRE' ;
  474. HYFAC   = 'DOMA'  MODHYB  'FACE';
  475.  
  476. *- Création ligne de suivi pour le post-traitement et le test
  477. * ligne des points centres (cas 1D)
  478. *
  479. 'REPETER' BCL (ENY - 1) ;
  480.     IP = &BCL ;
  481.     JP = IP + 1 ;
  482.     PI = 'POINT' HYCEN IP ;
  483.     PJ = 'POINT' HYCEN JP ;
  484.    'SI' (IP 'EGA' 1);
  485.       LCENC = ('MANU' 'SEG2' PI PJ) ;
  486.    'SINON' ;
  487.       LCENC = LCENC 'ET' ('MANU' 'SEG2' PI PJ) ;
  488.    'FINSI' ;
  489.  FIN BCL ;
  490.  
  491. *--------------------------------------------------------------------
  492. *- pression initiale (metre d'eau) dans le sable
  493. HN    = -1.15 ;
  494.  
  495. *--------------------------------------------------------------------
  496. *- Conditions aux limites
  497.  
  498.  
  499. * GBM modifie supprime bloque
  500. *- frontière en limite du domaine de calcul (milieu désaturé)
  501. ESORT = 'MANU' 'CHPO' CESORT 1 'TH' HN 'NATURE' 'DISCRET';
  502.  
  503. *- frontière mouillée
  504. EENTR =  'MANU' 'CHPO' CEENTR 1 'TH' 0. 'NATURE' 'DISCRET';
  505.  
  506.  
  507. *- chargement des CLs
  508. LICALC = 'PROG' 0.D0  1.e20 ;
  509. LIST1  = 'PROG' 2 * 1.D0 ;
  510. VALI0  = 'CHAR' 'THIM' (ESORT et EENTR)  ('EVOL' 'MANU' LICALC LIST1) ;
  511.  
  512. *--------------------------------------------------------------------
  513. *- initialisation des inconnues
  514. *
  515.  
  516. *- trace de charge d'eau
  517. *TH0  = 'MANU' 'CHPO' ('DIFF' CEENTR HYFAC) 1 'TH' HN 'NATURE' 'DISCRET';
  518. *TH0  = TH0 + SIMPE ;
  519.  
  520. *- charge d'eau
  521. H0   = 'MANU' 'CHPO' HYCEN 1 'H' HN 'NATURE' 'DIFFUS';
  522.  
  523. *- flux
  524. QFACE0 = 'MANU' 'CHPO' HYFAC 1 'FLUX' 0.D0 'NATURE' 'DISCRET' ;
  525.  
  526.  
  527. *                                         ---------------------------
  528. *                                         = Table DARCY_TRANSITOIRE =
  529. *                                         ---------------------------
  530. *- initialisation table
  531. SATUR                     = 'TABLE' ;
  532. SATUR. 'TEMPS'            = 'TABLE' ;
  533. SATUR. 'CHARGE'           = 'TABLE' ;
  534. SATUR. 'FLUX'             = 'TABLE' ;
  535.  
  536. SATUR . 'ITMAX' = 45;
  537.  
  538.  
  539. *- données géommétriques
  540. SATUR. 'SOUSTYPE'         = 'DARCY_TRANSATUR' ;
  541. SATUR. 'MODELE'           = MODHYB              ;
  542.  
  543.  
  544. SATUR . 'XI' = 1.D-10;
  545. SATUR . 'COEF_EMMAGASINEMENT' = 0.D0 ;
  546. *- instant initial
  547. SATUR. 'TEMPS' . 0        = 0.                  ;
  548. SATUR. 'CHARGE' . 0       = H0                  ;
  549. SATUR. 'FLUX' . 0         = QFACE0              ;
  550.  
  551. *********************************************
  552. *- conditions aux limites et chargements
  553. SATUR . 'TRACE_IMPOSE' = VALI0 ;
  554.  
  555. SATUR . 'LUMP'   = FAUX;
  556. SATUR . 'TYPDISCRETISATION'   = 'VF' ;
  557. SATUR . 'CONSERVATIF'   = VRAI  ;
  558.  
  559.  
  560.  
  561. TABRES = table METHINV;
  562. TABRES . 'TYPINV' = 2;
  563. TABRES . 'PRECOND' = 5;
  564.  
  565. SATUR . 'METHINV' = TABRES;
  566. SATUR . 'METHINV' . RESID = 1.D-20;
  567.  
  568. *******************************************
  569.  
  570.  
  571.  
  572. *******************************************
  573.  
  574.  
  575.  
  576. *- données physiques
  577. * loi de perméabilité
  578. LoiP                      = 'TABLE' 'PUISSANCE';
  579. LoiP. 'ALPHA'             = 8.D0 ;
  580. LoiP. 'PERMSAT'           = 5.85E-2 / 3600. ;
  581. SATUR.'LOI_PERMEABILITE'  = 'TABLE' LoiP ;
  582.  
  583. * loi de succion
  584. LoiS                      = 'TABLE' 'VAN_GENUCHTEN';
  585. LoiS. 'PORO'              = 0.3 ;
  586. *LoiS. 'PORO'              = 1. ;
  587.  
  588. LoiS. 'TERESIDU'          = 0.055 ;
  589. *LoiS. 'TERESIDU'          = 0.0 ;
  590.  
  591. LoiS. 'NEXP'              = 2.0304  ;
  592. *LoiS. 'NEXP'              = 1.D0 ;
  593.  
  594. LoiS. 'MEXP'              = 0.5075 ;
  595. *LoiS. 'MEXP'              = -1.D0 ;
  596.  
  597. LoiS. 'BHETA'             = 0.029227 * 100. ;
  598. SATUR.'LOI_SATURATION'    = 'TABLE' LoiS ;
  599.  
  600. *- données numériques
  601. 'SI' (COMPLET) ;
  602. SATUR. 'TEMPS_FINAL'     = 10000.D0 ;
  603. 'REPETER' bbb 4 ;
  604.   tfin = 2500.D0 * &bbb ;
  605. * var de Boltzman
  606.   znew = xabs '/' (tfin '**' 0.5D0);
  607.   ress = 'IPOL' znew roro h_prog;
  608.   'ORDONNER' CROISSANT ress;
  609.   'SI' ('EGA' &bbb 1) ;
  610.      courb = evol tabcool . 1 manu xabs ress ;
  611.   'SINON'     ;
  612.      courb = courb 'ET'
  613.      ('EVOL' tabcool . &bbb manu xabs ress);
  614.   'FINSI' ;
  615. 'FIN' bbb;
  616.  
  617. TAB1 = 'TABLE' ;
  618. TAB1 . 'TITRE' = 'TABLE' ;
  619. TAB1 . 'TITRE' . 1 = 'exacte 2500' ;
  620. TAB1 . 'TITRE' . 2 = 'exacte 5000' ;
  621. TAB1 . 'TITRE' . 3 = 'exacte 7500' ;
  622. TAB1 . 'TITRE' . 4 = 'exacte 10000' ;
  623. TAB1 . 1 = 'TIRR ';
  624. TAB1 . 2 = 'TIRR ';
  625. TAB1 . 3 = 'TIRR ';
  626. TAB1 . 4 = 'TIRR ';
  627. 'SINON' ;
  628. 'REPETER' bbb 4 ;
  629.   tfin = 500.D0 * &bbb ;
  630. * var de Boltzman
  631.   znew = xabs '/' (tfin '**' 0.5D0);
  632.   ress = 'IPOL' znew roro h_prog;
  633.   'ORDONNER' CROISSANT ress;
  634.   'SI' ('EGA' &bbb 1) ;
  635.      courb = evol tabcool . 1 manu xabs ress ;
  636.   'SINON'     ;
  637.      courb = courb 'ET'
  638.      ('EVOL' tabcool . &bbb manu xabs ress);
  639.   'FINSI' ;
  640. 'FIN' bbb;
  641.  
  642. TAB1 = 'TABLE' ;
  643. TAB1 . 'TITRE' = 'TABLE' ;
  644. TAB1 . 'TITRE' . 1 = 'exacte 500' ;
  645. TAB1 . 'TITRE' . 2 = 'exacte 1000' ;
  646. TAB1 . 'TITRE' . 3 = 'exacte 1500' ;
  647. TAB1 . 'TITRE' . 4 = 'exacte 2000' ;
  648. TAB1 . 1 = 'TIRR ';
  649. TAB1 . 2 = 'TIRR ';
  650. TAB1 . 3 = 'TIRR ';
  651. TAB1 . 4 = 'TIRR ';
  652. SATUR. 'TEMPS_FINAL'     = 2000.D0 ;
  653. 'FINSI' ;
  654.  
  655. SATUR. 'HOMOGENEISATION' = 'CHAINE' 'DECENTRE' ;
  656. SATUR. 'NPAS'            = 4000 ;
  657. SATUR. 'RESIDU_MAX'      = 1.D-4 ;
  658. SATUR. 'NITER'           = 10 ;
  659. SATUR. 'DT_INITIAL'      = 1.D0 ;
  660. 'SI' COMPLET ;
  661. SATUR. 'TEMPS_SAUVES'    = ('PROG' 250 'PAS' 250 1000)*10. ;
  662. 'SINON'  ;
  663. SATUR. 'TEMPS_SAUVES'    = ('PROG' 500 'PAS' 500 2000) ;
  664. 'FINSI' ;
  665. *SATUR. 'TEMPS_CALCULES'    = ('PROG' 1 'PAS' 1 10.)*360. ;
  666. SATUR. 'MESSAGES'        = 2 ;
  667.  
  668.  
  669. 'SI' GRAPH ;
  670. *--------------------------------------------------------------------
  671. *-  Vérification du choix du dp minimum pour le calcul de la capacité.
  672. * droite support des variables
  673.   dx  = 'DROIT' (0. -1.15 ) 1000  (0. 0.) ;
  674.   zc  = 'COOR' dx 2 ;
  675.   ev2 = 'EVOL' 'BLEU' 'CHPO' zc 'SCAL' dx ;
  676.   px  = 'EXTR' ev2 'ORDO' ;
  677. * calcul de la teneur en eau pour la pression zc
  678.   s0 t0 cap = HT_PRO (SATUR.'LOI_SATURATION') ZC  ;
  679.   ev0 = 'EVOL' 'TURQ' 'CHPO' s0 'SCAL' dx ;
  680.   evt = 'EVOL' 'VERT' 'MANU' px (100. * ('EXTR' ev0 'ORDO')) ;
  681.   'DESSIN' evt 'TITX' 'Pc(m)' 'TITY' 'S(%)'
  682.                'TITRE' 'Loi capillaire S(Pc)' ;
  683. * calcul de la teneur en eau pour la pression zc - dp
  684.   dp  = 1.e-4 ;
  685.   s1 t1 cap = HT_PRO (SATUR.'LOI_SATURATION') 
  686.               ('KOPS' zc '-'  dp) ;
  687. * représentation de la capacité
  688.   c1  = (t0 - t1) / dp;
  689.   ev1 = 'EVOL' 'ROUGE' 'CHPO' c1 'SCAL' dx ;
  690.   evc = 'EVOL' 'JAUNE' 'MANU' px ('EXTR' ev1 'ORDO') ;
  691.   'DESSIN' evc 'TITX' 'Pc(m)' 'TITY' 'Capa(1/m)'
  692.                'TITRE' 'Capacite capillaire' ;
  693. 'FINSI' ;
  694.  
  695. *                                           ===========================
  696. *                                           |       1er CALCUL        |
  697. *                                           ===========================
  698. *- fonctionnement avec temps automatiques
  699.  
  700. DARCYSAT SATUR ;
  701.  
  702. * Post-traitement
  703. * ===============
  704.  
  705.   LT        = 'LECT' 1  PAS 1 ((dime SATUR.TEMPS) - 1)  ;
  706.   liopt     = 'MOTS' 'MIMA' 'AXES';
  707.   TDES      = TRACHIS SATUR 'CHARGE' LCENC LT 'PREF' ' ' 'UNIT' 's' ;
  708.   TDES . 5 = 'TABLE' ;
  709.   TDES . 5 . VALEUR = 'EXTRAIRE'  courb 'COUR' 1;
  710.   TDES . 5 . LEGEND2 = TAB1 . TITRE . 1;
  711.   TDES . 5 . LEGEND1 = ' ' ;
  712.   TDES . 6 = 'TABLE' ;
  713.   TDES . 6 . VALEUR = 'EXTRAIRE'  courb 'COUR' 2;
  714.   TDES . 6 . LEGEND2 = TAB1 . TITRE . 2;
  715.   TDES . 6 . LEGEND1 = ' ' ;
  716.   TDES . 7 = 'TABLE' ;
  717.   TDES . 7 . VALEUR = 'EXTRAIRE'  courb 'COUR' 3;
  718.   TDES . 7 . LEGEND2 = TAB1 . TITRE . 3;
  719.   TDES . 7 . LEGEND1 = ' ' ;
  720.   TDES . 8 = 'TABLE' ;
  721.   TDES . 8 . VALEUR = 'EXTRAIRE'  courb 'COUR' 4;
  722.   TDES . 8 . LEGEND2 = TAB1 . TITRE . 4;
  723.   TDES . 8 . LEGEND1 = ' ' ;
  724.  
  725. 'SI' GRAPH ;  
  726.     DESTRA TDES liopt 'TITX' 'z (m)' 'TITY' 'Pw (m)' ;
  727. 'FINSI' ;
  728.  
  729. * On compare l'integrale de la saturation au temps final entre les
  730. * solutions analytiques et non analytiques.
  731.  
  732. toto = 'PRIM' TDES . 8 . VALEUR;
  733. inta = 'EXTRAIRE' toto 'ORDO';
  734. inta = 'EXTRAIRE'  ('DIME' inta) inta; 
  735.  
  736. toto = 'PRIM' TDES . 4 . VALEUR;
  737. intb = 'EXTRAIRE' toto 'ORDO';
  738. intb = 'EXTRAIRE'  ('DIME' intb) intb;
  739.  
  740. err = 'ABS' (inta '-' intb '/' inta);
  741. 'LISTE' err;
  742.  
  743. 'SI' (err > 7.D-2) ;
  744.    'ERREUR'  5;
  745. 'FINSI'  ;
  746.  
  747. 'FIN';
  748.  
  749.  
  750.  
  751.  
  752.  
  753.  
  754.  
  755.  
  756.  
  757.  
  758.  
  759.  

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