Accueil | Cast3M

* fichier :  tran4.dgibi
************************************************************************
************************************************************************
 
*******************************************
*    Test tran4.dgibi: jeux de données    *
*    ---------------------------------    *
*******************************************
 
* SI GRAPH = N PAS DE GRAPHIQUE AFFICHE
* SINON SI GRAPH DIFFERENT DE N TOUS
* LES GRAPHIQUES SONT AFFICHES
 
GRAPH = 'N' ;
 
SAUT PAGE;
SI (NEG GRAPH 'N') ;
  OPTI ECHO 1   ;
  OPTI TRAC PSC   ;
SINO ;
  OPTI ECHO 0   ;
FINSI ;
 
SAUT PAGE;
 
*******************************************
*  pour calcul complet mettre complet=vrai;
*******************************************
complet = vrai;
 
*-----------------------------------------*
*                 TRAN4                   *
*                                         *
* 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                    *
*  -------------------------------------  *
* |       BARRE                         | *
*  -------------------------------------  *
*  x=0                               T=0° *
*                              à l'inFINi *
*                                         *
* Changement de phase à 1° avec chaleur   *
* latente massique L. Toutes les autres   *
* caractéristiques sont egales à 1.       *
*                                         *
* Recherche de 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
 
*   Pour le test on conseille
*   dt=0.025 dx=0.01 EPSILON=0.1
 
* 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.025 ;
* pas d'espace
DX = 0.05 ;
*==========================================
 
*   MAILLAGE avec des QUA4
OPTI DIME 2 ELEM SEG2;
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' ;
 
SI(NEG GRAPH 'N');
  TRAC 'QUAL' su1 ;
FINSI;
 
*  MODELE DE CALCUL
mod1 = MODE su1 THERMIQUE 'QUA4'  ;
mod2 = MODE su1 'CHANGEMENT_PHASE' 'INCO' 'T' 'Q' ;
 
*  CARACTERISTIQUES
mat  = MATE mod1 'C' 1. 'K' 1. 'RHO' 1. ;
mat2 = MATE mod2 'PRIM' 1. 'DUAL' L;
 
*   TEMPÉRATURE IMPOSÉE EN X=0
bl1   = BLOQ 'T' li4 ;
timp1 = DEPI bl1 2. ;
ev1   = EVOL MANU ABS (PROG 0. 1.) ORD (PROG 1. 1.);
cha   = CHAR TIMP timp1 ev1;
 
* TABLE POUR THERMIC
xtab              = TABL ;
xtab.temperatures = TABL;
xtab.temperatures . 0 = ( MANU CHPO su1 1 'T' 0. ) + ( MANU CHPO li4 1 'T' 2. )  ;
xtab.'MODELE'              = mod1 et mod2 ;
xtab.'CARACTERISTIQUES'    = mat  et mat2 ;
xtab.'CHARGEMENT'          = cha;
xtab.'BLOCAGES_THERMIQUES' = bl1 ;
 
SI complet ;
  tFIN = 1.;
SINON;
  tFIN = 0.25;
FINSI;
 
xtab.'TEMPS_CALCULES' = PROG 0. pas DT 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*(t**0.5)
*     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;
REPE 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 :
* recherche du point dont T vaut 1
  REPE trak ( DIME lt) ;
    tt2 = EXTR lt (&trak) ;
    SI ( tt2 < 1. ) ;
      xx1 = EXTR ( EXTR ev1 'ABSC' ) (&trak - 1) ;
      tt1 =  EXTR lt (&trak - 1) ;
      xx2 = EXTR ( EXTR ev1 'ABSC' ) &trak ;
 
*     Interpolation 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 ) ;
 
* TITR 'Instant' (xtab. temps . ite) ;
* dess 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 ;
 
*   VISUALISATION
evfr     = EVOL BLEU MANU 'TEMPS' ltps 'ABSCISSE' lfr ;
evfrtheo = EVOL ROUG MANU 'TEMPS' ltps 'ABSCISSE' lfrtheo ;
evfrsans = EVOL VERT MANU 'TEMPS' ltps 'ABSCISSE' lfrsans ;
 
SI(NEG GRAPH 'N');
  TITR 'Avance du front : EPSILON=' EPSILON 'Dt=' dt 'Dx=' dx ;
  DESS ( evfr ET evfrtheo ) ;
FINSI;
 
 
*Retrait du premier terme nul des LISTREEL
LENT1 = LECT;
REPE SURIND ((DIME lfr) - 1);
  Indice = &SURIND;
  MESS Indice;
  LENT1 = LENT1 ET (LECT (Indice+1));
 
FIN SURIND;
 
lfr2     = EXTR lfr LENT1;
lfrtheo2 = EXTR lfrtheo LENT1;
 
erreu=100. * (abs(lfr2 - lfrtheo2) / lfrtheo2);
 
MESS 'Erreur Maximale : ' (MAXI erreu) ;
 
*   Message d'erreur
SI ( (MAXI erreu) > 13. ) ;
  MESS 'Résultat incorrect' ;
  ERRE 5 ;
FINSI ;
 
*OPTI SAUV 'test.res'        ;
*SAUV evfr evfrtheo evfrsans ;
FIN ;