Exemples de jeux de données Cast3M

* fichier :  trajec.dgibi
'SAUTER' PAGE ;
'OPTION' 'ECHO' 0 ;
GRAPH = faux ;
*
'TITRE' ' HYDROCOIN 7b inclusion spherique' ;
*-----------------------------------------------------------------------
*=                       HYDROCOIN 7B
*= Cas isotrope avec une permeabilite p1 a l'interieur d'un cercle de
*= rayon RAY0, et une permeabilite p2 a l'exterieur.
*= Les dimensions sont (-50,50)*(0,30)
*= Le maillage est construit en symétrisant 5 zones.
*=
*=       Nombre d'élément = 1422   *       Nombre de faces  = 2901
*=
*= On va calculer les trajectoires pour les deux types de formulation.
*= Les vitesses sont données par une formule analytique.
*-----------------------------------------------------------------------
*
*
*
*    C12                                                       C11
*      C12--------------------------------------------------C11
*      |                                                      |
*      |Z5                                                    |
*      |                                                      |
*      |                             C8                       |
*      C9---------------------------*                         |
*      |Z4                          |                         |
*      |Z3                          |                         |
*      |                            |                         |
*      |Z2                          |                         |
*      |Z1                          |                         |
*      C2                           |                         |
*      |     C3                     |                         |
*      |                            |                         |
*      |                            |                         |
*      |  Z0    Z1   Z2  Z3    Z4   |           Z5            |
*      C0----C1----C4--C5--C6------C7-----------------------C10
*
*
*  Z0 : Délimité par les pts C0 a C3             |
*       Maillage IX0*IY0 avec DALLER             |   Zone ou regne une
*  Z1 : Délimité par les pts C1 , C4             |      permabilite
*       Maillage IZ1*SXY avec REGLER             |          p1
*  Z2 : Délimité par les pts C4 , C5             |
*       Maillage I45*SXY avec ROTA               |
*
*  Z3 : Délimité par les pts C5 , C6             |
*       Maillage I56*SXY avec ROTA               |
*  Z4 : Délimité par les pts C6R, C7 a C9        |   Zone ou regne une
*       Maillage IZ4*SXY avec REGLER             |      permabilite
*  Z5 : Délimité par les pts C7 a C12            |          p2
*       Maillage IZ5*SXY avec REGLER             |
*
*
*
'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
*
*- INITIALISATIONS pour la decoupe des droites support des zones
*
IX0 = 9         ;
IY0 = 15        ;
SXY = IX0 + IY0 ;
I45 = 8         ;
I56 = 4         ;
IZ1 = 1         ;
IZ4 = 7         ;
IZ5 = 4         ;
*
*- INITIALISATIONS des coordonnees des points support des droites
*
RAC2 = 2.**.5 ;
RAY0 = 10.0   ;
RS5  = RAY0 / 5.0  ;
PTS5 = RS5  / RAC2 ;
XC4  = RAY0 - RS5 / I45 / 2. + RS5 ;
XC6  = 12.0 ;
XC7  = 25.0 ;
YC8  = 20.0 ;
XC10 = 50.0 ;
YC11 = 30.0 ;
*
*- Definition des points de base
*
C0  = 0.0  0.0  ;
C1  = RS5  0.0  ;
C2  = 0.0  RS5  ;
C3  = PTS5 PTS5 ;
C4  = XC4  0.0  ;
C5  = RAY0 0.0  ;
C6  = XC6  0.0  ;
C7  = XC7  0.0  ;
C8  = XC7  YC8  ;
C9  = 0.0  YC8  ;
C10 = XC10 0.0  ;
C11 = XC10 YC11 ;
C12 = 0.0  YC11 ;
*
*- Definition de la partie carree pres de l'origine
*
C0C1 = C0 D IX0 C1                ;
C1C3 = C1 D IY0 C3                ;
C3C2 = C3 D IX0 C2                ;
C2C0 = C2 D IY0 C0                ;
Z0   = 'DALLER' C0C1 C1C3 C3C2 C2C0 ;
*
*- Definition de la zone jusqu'au changement de permeabilite
*
C4C5  = 'DROIT' C4 I45 C5         ;
Z2    = 'ROTATION' C4C5 SXY 90 C0 ;
*
*- Definition de la zone entre le carre et le cercle
*
C132 = C1C3 'ET' C3C2       ;
COT4 = 'COTE' 4 Z2          ;
COT4 = 'INVERSE' COT4          ;
Z1   = 'REGLER' IZ1 C132 COT4 ;
C1C4 = 'COTE' 4 Z1          ;
C1C4 = 'INVERSE' C1C4          ;
*
*- Definition du cercle au dela du changement de permeabilite
*
C5C6 = 'DROIT' C5 I56 C6         ;
Z3   = 'ROTATION' C5C6 SXY 90 C0 ;
*
*- Definition de la zone entre cercle et rectangle1
*
C7C8 = 'DROIT' C7 IX0 C8        ;
C8C9 = 'DROIT' C8 IY0 C9        ;
C789 = C7C8 'ET' C8C9       ;
COT2 = 'COTE' 2 Z3          ;
Z4   = 'REGLER' IZ4 C789 COT2 ;
C6C7 = 'COTE' 4 Z4          ;
*
*- Definition de la zone rectangle2
*
C10C11 = 'DROIT' C10 IX0 C11       ;
C11C12 = 'DROIT' C11 IY0 C12       ;
C1012  = C10C11 'ET' C11C12    ;
Z5     = 'REGLER' IZ5 C789 C1012 ;
*
* Définition des zones et symétrisation
*
ZONED1 = Z0 'ET' Z1 'ET' Z2                         ;
ZONED2 = Z3 'ET' Z4 'ET' Z5                         ;
ZONEG1 ZONEG2 = 'SYMETRIE' ZONED1 ZONED2 'DROIT' C0 C12 ;
ZONE1  = ZONED1 'ET' ZONEG1                       ;
ZONE2  = ZONED2 'ET' ZONEG2                       ;
*
*
*- Creation maillage GEOMETRIQUE
*
PTOT  = ZONE1 'ET' ZONE2 ;
'ELIMINATION' PTOT 0.001       ;
PTOTF = 'CHANGER' PTOT  'QUAF' ;
ZONF1 = 'CHANGER' ZONE1 'QUAF' ;
ZONF2 = 'CHANGER' ZONE2 'QUAF' ;
'ELIMINATION' 0.001 (PTOTF 'ET' ZONF1 'ET' ZONF2) ;
*
*- Initialisation des parametres physiques
*
AA    = RAY0      ;
GG    = 1.E+3     ;
KKB   = 1.E-15    ;
KKC   = 1.E-13    ;
MMU   = 1.E-3     ;
FIB   = 1.E-1     ;
FIC   = 1.E-1     ;
AA2   = AA * AA   ;
KCSM  = KKC / MMU ;
KBSM  = KKB / MMU ;
KCPKB = KKC + KKB ;
KCMKB = KKC - KKB ;
*
*----------
* VITESSE |
*----------
*
*      |                       2    2   2
*      |  kb           kc-kb  a  ( y - x )
*      |  -- GG ( 1 - ------- ------------  pour r>a
*      |  MU           kc+kb      r4
*   V  |
*    x |
*      |  kc      2kb
*      |  -- GG --------                    pour r<a
*      |  MU     kc+kb
*      |
*
*      |                       2    2   2
*      |  kb           kc-kb  a  ( y - x )
*      |  -- GG ( 1 - ------- ------------  pour r>a
*      |  MU           kc+kb      r4
*   V  |
*    y |
*      |  0.                                pour r<a
*      |
*      |
*
*                                              --------------------
*                                              définition du LACHER
*                                              --------------------
LACHER = 'TABLE' ;
LACHER.TEMPS_LACHER = 'PROG' 0.   ;
LACHER.TEMPS_LIMITE =  1.E+15 ;
LACHER.CFL =  0.05 ;
LACHER.DELTAT_SAUVE =  0. ;
LACHER.1 = (-49.99999  10.) 'ET' (-49.99999 12.) 'ET' (-49.99999 14.)
           'ET' (-49.99999 16.) 'ET' (-49.99999  18.)  ;
 
* ----------------------------------------------------------------
*       lorsque la trajectoire passera d'une zone à l'autre
*  on réinitialisera le lacher avec la derniere position calculée
* ----------------------------------------------------------------
*
*                           -------------------
*                            formulation DARCY
*                           -------------------
*
*- Creation MODELE HYBRIDE
*
MHYTOT = MODE PTOTF DARCY ISOTROPE                 ;
MHYZ1  = MODE ZONF1 DARCY ISOTROPE  ;
MHYZ2  = MODE ZONF2 DARCY ISOTROPE  ;
ZHYB1  = DOMA MHYZ1 SURFACE       ;
ZHYB2  = DOMA MHYZ1 NORMALE       ;
YHYB1  = DOMA MHYZ2  SURFACE      ;
YHYB2  = DOMA MHYZ2  NORMALE      ;
*
*                                               --------------------
*                                               = Charge aux FACES =
*                                               --------------------
*
XX YY = 'COOR' (DOMA MHYTOT FACE )  ;
XX1   = REDU XX (DOMA MHYZ1 FACE) ;
XX2   = REDU XX (DOMA MHYZ2 FACE) ;
YY2   = REDU YY (DOMA MHYZ2 FACE) ;
RX2   = XX2 * XX2         ;
RY2   = YY2 * YY2         ;
RR2   = RX2 + RY2         ;
RR4   = RR2 * RR2         ;
*
*
*                                     ------------------------------
*                                     = Vitesse 'ET' debit aux FACES =
*                                     ------------------------------
*
CSTZ1 = KCSM * KKB * GG * 2.0 / KCPKB    ;
UZ1   = MANU CHPO (DOMA MHYZ1 FACE) 1 SCAL CSTZ1 ;
UZ2   = (KBSM * KCMKB / KCPKB * GG * ((AA2 * (RX2 - RY2))/(RR4)))
      + (KBSM * GG)                      ;
VZ1   = MANU CHPO (DOMA MHYZ1 FACE) 1 SCAL 0.                    ;
VZ2   = KBSM * KCMKB/KCPKB * GG * 2. * ( XX2*YY2*AA2 / (RR4) ) ;
UZ1   = NOMC VX UZ1 ;
UZ2   = NOMC VX UZ2 ;
VZ1   = NOMC VY VZ1 ;
VZ2   = NOMC VY VZ2 ;
*
UZ1   = CHANGER 'ATTRIBUT' UZ1 'NATURE' 'DISCRET' ;
UZ2   = CHANGER 'ATTRIBUT' UZ2 'NATURE' 'DISCRET' ;
VZ1   = CHANGER 'ATTRIBUT' VZ1 'NATURE' 'DISCRET' ;
VZ2   = CHANGER 'ATTRIBUT' VZ2 'NATURE' 'DISCRET' ;
*
VN1   = UZ1 'ET' VZ1 ;
VN2   = UZ2 'ET' VZ2 ;
MOT1  = MOTS VX VY ;
MOT2  = MOTS UX UY ;
VAVN1 = PSCA VN1 ZHYB2 MOT1 MOT2 ;
VAVN2 = PSCA VN2 YHYB2 MOT1 MOT2 ;
VAVN1 = NOMC SCAL VAVN1 ;
VAVN2 = NOMC SCAL VAVN2 ;
QFACE1= VAVN1* ZHYB1 ;
QFACE1= NOMC FLUX QFACE1 ;
QFACE2= VAVN2* YHYB1 ;
QFACE2= NOMC FLUX QFACE2 ;
*                           --------------------------------------
*                           trajectoires formulation mixte hybride
*                           --------------------------------------
MOD2 UUU2 = 'TRAJ' 'CONVECTION_EXPLICITE' MHYZ2 QFACE2  LACHER ;
LATAB2 = 'EXTR' MOD2 'ZONE';
LACHER.TEMPS_LACHER = 'PROG' 5*0. ;
'REPETER' BLOC11 5 ;
  I    = &BLOC11 ;
  NBPT = 'NBEL' LATAB2. (2*I) ;
  PT1  = 'POINT' LATAB2. (2*I) FINAL ;
  TP1  = 'EXTR' UUU2 'TMPS' I NBPT 2 ;
  LACHER.I = 'MANU' 'POI1' PT1 ;
  'REMPLACER' LACHER.TEMPS_LACHER I TP1 ;
'FIN' BLOC11 ;
 
MOD1 UUU1 = 'TRAJ' 'CONVECTION_EXPLICITE' MHYZ1 QFACE1  LACHER ;
LATAB1 = 'EXTR' MOD1 'ZONE';
'REPETER' BLOC12 5 ;
  I    = &BLOC12 ;
  NBPT = 'NBEL' LATAB1. (2*I) ;
  PT1  = 'POINT' LATAB1. (2*I) FINAL ;
  TP1  = 'EXTR' UUU1 'TMPS' I NBPT 2 ;
  LACHER.I = 'MANU' 'POI1' PT1 ;
  'REMPLACER' LACHER.TEMPS_LACHER I TP1 ;
'FIN' BLOC12 ;
 
MOD3 UUU3 = 'TRAJ' 'CONVECTION_EXPLICITE' MHYZ2 QFACE2  LACHER ;
LATAB3 = 'EXTR' MOD3 'ZONE';
*                                       ------------------------
*                                        controle des résultats
*                                       ------------------------
BTTP  = 'PROG' 8.3324E+10 8.5795E+10 8.9149E+10 9.3984E+10 10.2586E+10 ;
VERRE = 0 ;
'REPETER' BLOC13 5 ;
  I       = &BLOC13 ;
  NBPT    = 'NBEL' LATAB3. (2*I) ;
  XI2  I2 = 'COOR' ('POINT' LATAB3. (2*I) 'FINAL');
  XI1 I1  = 'COOR' ('POINT' LATAB2. (2*I) 'INITIAL')   ;
  EPSS    = 'ABS' ( I1 - I2) ;
  'SI' ( EPSS > 5.E-4) ;
    VERRE=VERRE+1 ;
  'FINSI' ;
  TP1= EXTR UUU3 'TMPS' I NBPT 2 ;
  TP2= EXTR BTTP I ;
  EPST= ABS( TP1 - TP2) ;
  'SI' ( EPST > 1.E+8 ) ;
    VERRE=VERRE+1 ;
  'FINSI' ;
'FIN' BLOC13 ;
*                             ----------------------
*                             tracé des trajectoires
*                             ----------------------
'SI' GRAPH ;
   CONT1 = CONTOUR (DOMA MHYZ1 'MAILLAGE') ;
   CONT2 = CONTOUR (DOMA MHYZ2 'MAILLAGE') ;
   CROB1 = EXTR UUU1 'MAIL' ;
   CROB2 = EXTR UUU2 'MAIL' ;
   CROB3 = EXTR UUU3 'MAIL' ;
  'TITRE' ' Trajectoires   HYDROCOIN 7b inclusion spherique' ;
  'TRACER' (CROB1 'ET' CROB2 'ET' CROB3 'ET' CONT1 'ET' CONT2) ;
'FINSI' ;
*                           -------------------------
*                           formulation NAVIER_STOKES
*                           -------------------------
*
*
*- Creation MODELE NAVIER_STOKES
*
MNSTOT = MODE PTOTF NAVIER_STOKES LINE                 ;
MNSZ1  = MODE ZONF1 NAVIER_STOKES LINE  ;
MNSZ2  = MODE ZONF2 NAVIER_STOKES LINE  ;
*
*                                               --------------------
*                                               = Charge aux SOMMETS
*                                               --------------------
*
XXS YYS = COOR (DOMA MNSTOT 'MAILLAGE')   ;
XXS1    = REDU XXS (DOMA MNSZ1 'MAILLAGE') ;
XXS2    = REDU XXS (DOMA MNSZ2 'MAILLAGE') ;
YYS2    = REDU YYS (DOMA MNSZ2 'MAILLAGE') ;
RXS2    = XXS2 * XXS2         ;
RYS2    = YYS2 * YYS2         ;
RRS2    = RXS2 + RYS2         ;
RRS4    = RRS2 * RRS2         ;
*
*
*                                     -----------------------
*                                     = Vitesse AUX SOMMETS =
*                                     -----------------------
*
CSTZ1 = KCSM * KKB * GG * 2.0 / KCPKB    ;
UZS1  = 'MANU' 'CHPO' ('DOMA' MNSZ1 'MAILLAGE') 1 'SCAL' CSTZ1 ;
UZS2  = (KBSM * KCMKB / KCPKB * GG * ((AA2 * (RXS2 - RYS2))/(RRS4)))
      + (KBSM * GG)                      ;
VZS1  = 'MANU' 'CHPO' (DOMA MNSZ1 'MAILLAGE') 1 'SCAL' 0. ;
VZS2  = KBSM * KCMKB/KCPKB * GG * 2. * ( XXS2*YYS2*AA2 / (RRS4) ) ;
UZS1  = NOMC VX UZS1 ;
UZS2  = NOMC VX UZS2 ;
VZS1  = NOMC VY VZS1 ;
VZS2  = NOMC VY VZS2 ;
*
UZS1  = CHANGER 'ATTRIBUT' UZS1 'NATURE' 'DISCRET' ;
UZS2  = CHANGER 'ATTRIBUT' UZS2 'NATURE' 'DISCRET' ;
VZS1  = CHANGER 'ATTRIBUT' VZS1 'NATURE' 'DISCRET' ;
VZS2  = CHANGER 'ATTRIBUT' VZS2 'NATURE' 'DISCRET' ;
*
VNS1  = UZS1 'ET' VZS1 ;
VNS2  = UZS2 'ET' VZS2 ;
*
LACHER.1 = (-49.99999 10.) 'ET' (-49.99999 12.) 'ET' (-49.99999 14.)
      'ET' (-49.99999 16.) 'ET' (-49.99999 18.)  ;
LACHER.TEMPS_LACHER = 'PROG' 5*0. ;
*           
MODS2 UUUS2 = 'TRAJ' MNSZ2 VNS2 LACHER ;
 
LATABS2 = 'EXTR' MODS2 'ZONE';
'REPETER' BLOC21 5 ;
  I    = &BLOC21 ;
  NBPT = 'NBEL' LATABS2. (2*I) ;
  PT1  = 'POINT' LATABS2. (2*I) 'FINAL' ;
  TP1  = 'EXTR' UUUS2 'TMPS' I NBPT 2 ;
  LACHER.I = 'MANU' 'POI1' PT1 ;
  'REMPLACER' LACHER.TEMPS_LACHER I TP1 ;
'FIN' BLOC21 ;
*
MODS1 UUUS1= 'TRAJ' MNSZ1 VNS1 LACHER ;
 
LATABS1 = 'EXTR' MODS1 'ZONE';
LACHER.TEMPS_LACHER= 'PROG' 5*0. ;
'REPETER' BLOC22 5 ;
  I    = &BLOC22 ;
  NBPT = 'NBEL' LATABS1. (2*I) ;
  PT1  = 'POINT' LATABS1. (2*I) 'FINAL' ;
  TP1  = 'EXTR' UUUS1 'TMPS' I NBPT 2 ;
  LACHER.I = 'MANU' 'POI1' PT1 ;
  'REMPLACER' LACHER.TEMPS_LACHER I TP1 ;
'FIN' BLOC22 ;
 
MODS3 UUUS3 = 'TRAJ' MNSZ2 VNS2 LACHER ;
LATABS3 = 'EXTR' MODS3 'ZONE';
*                                       ------------------------
*                                        contrôle des résultats
*                                       ------------------------
BTTP = 'PROG' 8.3324E+10 8.5795E+10 8.9149E+10 9.3984E+10 10.2586E+10 ;
REPETER BLOC23 5 ;
  I      = &BLOC23 ;
  NBPT   = 'NBEL' LATABS3. (2*I) ;
  XI2 I2 = 'COOR' ('POINT' LATABS3. (2*I) 'FINAL');
  XI1 I1 = 'COOR' ('POINT' LATABS2. (2*I) 'INITIAL')   ;
  EPSS   = 'ABS' (I1 - I2) ;
  'SI' (EPSS > 5.E-4) ;
     VERRE = VERRE + 1 ;
  'FINSI' ;
  TP1  = 'EXTR' UUUS3 'TMPS' I NBPT 2;
  TP2  = 'EXTR' BTTP I ;
  EPST = 'ABS' ( TP1 - TP2)  ;
  'SI' ( EPST > 5.E+8 ) ;
    VERRE = VERRE + 1 ;
  'FINSI' ;
'FIN' BLOC23 ;
*                             ----------------------
*                             tracé des trajectoires
*                             ----------------------
'SI' GRAPH ;
  CONT1 = CONTOUR (DOMA MNSZ1 'MAILLAGE') ;
  CONT2 = CONTOUR (DOMA MNSZ2 'MAILLAGE') ;
  CROB1 = EXTR UUUS1 'MAIL' ;
  CROB2 = EXTR UUUS2 'MAIL' ;
  CROB3 = EXTR UUUS3 'MAIL' ;
  TITRE ' Trajectoires   HYDROCOIN 7b inclusion spherique' ;
  'TRACER' (CROB1 'ET' CROB2 'ET' CROB3 'ET' CONT1 'ET' CONT2) ;
'FINSI'  ;
 
'SI' (VERRE 'EGA' 0) ;
  'ERREUR' 0 ;
'SINON' ;
  'ERREUR' 5 ;
'FINSI' ;
 
'FIN'  ;
*