Datasets examples

* fichier :  tran4.dgibi
complet = faux;
*
*************************************************
*  pour calcul complet mettre complet à : vrai;
*************************************************
************************************************************************
*       Test de l'option changement de phase de l'opérateur capacité  *
*       et de l'algorithme de transit3                                *
************************************************************************
*                       PROBLEME DE STEPHAN UNIDIMENSIONNEL
*
*    Transport de la chaleur par diffusion uniquement.
*
*    T(0,t) = 0 si t<0
*    T(0,t) = 2 si t>0
*        --------------------------------------
*   x=0  |       BARRE                          |   T=0° à l'infini
*        --------------------------------------
*
*    Changement de phase à 1° avec chaleur latente massique L.
*    Tous les autres caractéristiques sont egales à 1. 
*
*    On cherche l'avancée du front de changement de phase 
*    en fonction du temps (entre 0 et 1 seconde)
*
* REFERENCE pour la solution analytique
* "Conduction of Heat in Solids"  H S Carslaw and J C Jaeger 
*  Oxford Clarendon Press 1959
*  ch XI p291
************************************************************************
* 
*==============
* paramètres pour l'integration numérique
* largeur du changement de phase
EPSILON = 0.1 ;
*  chaleur latente
L = 4.7427 ;
*  constante de la solution analytique
BETA = 0.25 ;
*  pas de temps
DT = 0.05 ;
* pas d'espace
DX = 0.05 ;
*==============
*
*   MAILLAGE avec des qua4
opti dime 2 elem seg2;
opti trac ps ;
dens DX ;
p1 =   0.  0. ;
p4 =   0.  DX ;
dens (4.*DX) ;
p2 =   10.  0. ;
p3 =   10.  DX ;
*
li1 = d    p1 p2 ;
li2 = d 1  p2 p3 ;
li3 = d    p3 p4 ;
li4 = d 1  p4 p1 ;
opti elem qua4 ;
su1 = dall  li1  li2  li3  li4 'PLAN' ;
*trace su1 ;
*
*  MODELE DE CALCUL
*
mod1 = modele su1 thermique phase qua4  ;
*
*  CARACTERISTIQUES
*
mat = mate mod1 'C' 1. 'K' 1. 'RHO' 1.   tpha 1. QLAT L;
*
*   température imposée en x=0
bl1 = bloq 'T' li4 ;
timp1 = depimp bl1 2. ;
ev1 = evol manu abs (prog 0. 1.) ord (prog 1. 1.);
cha = char TIMP timp1 ev1;

*
* table pour thermic
xtab = table ;
xtab.temperatures = table;
xtab.temperatures . 0 = ( manu chpo su1 1 'T' 0. ) + 
              ( manu chpo li4 1 'T' 2. )  ;
xtab.'MODELE' = mod1 ;
xtab.'CARACTERISTIQUES' = mat;
xtab.'CHARGEMENT' = cha;
xtab.'BLOCAGES_THERMIQUES' = bl1 ;
si complet ;
   tfin = 1.;
sinon;
   tfin = 0.25;
finsi;
xtab.'TEMPS_CALCULES' = prog 0. pas 0.05 tfin;
xtab.'PROCEDURE_THERMIQUE' = DUPONT;
si complet;
xtab.'CRITERE' = 1E-2;
sinon;
xtab.'CRITERE' = 5E-2;
finsi;
xtab.'RELAXATION_THETA' = 1.;
*
*   CALCUL
*
pasapas xtab ;
*
*   ANALYSE DES RESULTATS
*
*     le front avance comme 2*beta*srqt(t)
*     beta et Chaleur latente sont liés
*
ltps = prog 0.;
lfr = prog 0. ;
lfrtheo = prog 0.;
lfrsans = prog 0. ;
*
*      boucle sur les pas de temps
*
npas = (dime (xtab.temps)) - 1;
repeter analyse npas;
ite = &analyse;
temp1 = xtab. temps . ite ;
ev1 = evol chpo (xtab. temperatures . ite ) ( li1) 'T';
lt = extr ev1 'T' ;
*
*      boucle sur les points de la droite li1
*      on cherche le point dont la température vaut 1
repeter trak ( dime lt) ;
tt2 = extrai lt (&trak) ;
si ( tt2 < 1. ) ;
      xx1 = extr ( extr ev1 'ABSC' ) (&trak - 1) ;
      tt1 =  extrai lt (&trak - 1) ;
      xx2 = extr ( extr ev1 'ABSC' ) &trak ;
*       on interpole entre deux noeuds
*       xx est la position du front
      xx = xx1 + ((1. - tt1 ) * ( xx2 - xx1 ) / ( tt2 - tt1 ));
      quitter trak ;
finsi ;
fin trak ;
*
*
xxtheo = 2. * beta * ( temp1 ** 0.5 ) ;
xxsans = 0.95 * ( temp1 ** 0.5 ) ;
*titre 'Instant' (xtab. temps . ite) ;
*dessin ev1 ;
mess 'TEMPS' temp1 ' Front =' xx 'Valeur Theorique' xxtheo;
ltps = ltps et ( prog temp1 ) ;
lfr = lfr et ( prog xx) ;
lfrtheo = lfrtheo et ( prog xxtheo ) ;
lfrsans = lfrsans et ( prog xxsans ) ;
fin analyse ;
*
*   message d'erreur
*   pour le test on conseille dt=0.025 dx=0.01 epsilon=0.1
erreu = abs(xx - xxtheo) / xx ;
list erreu;
si ( erreu > 0.05 ) ;
  mess 'Résultat incorrect' ;
  erreur 5 ;
finsi ;
*
*   VISUALISATION
*
evfr = evol manu 'TEMPS' ltps 'ABSCISSE' lfr ;
evfrtheo = evol manu 'TEMPS' ltps 'ABSCISSE' lfrtheo ;
evfrsans = evol manu 'TEMPS' ltps 'ABSCISSE' lfrsans ;
titre 'Avance du front.        Epsilon=' epsilon 'Dt=' dt 'Dx='
dx ;
*dess ( evfr et evfrtheo et evfrsans ) ;
*
*opti sauv 'test.res' ;
*sauver evfr evfrtheo evfrsans ;
fin ;