* fichier :  phase_02.dgibi
************************************************************************
* Section : Changement_De_Phase Changement_De_Phase
* Section : Thermique Conduction
************************************************************************
graph= faux;
*
*  Cas test :  phase2d.dgibi
*
*  barreau ayant une temperature variant de 0 à 250 jusqu'a son milieu
*  puis de 250 à 0°C à son extremité. On suppose une température de
* fusion à 200°C, on laisse évoluer la température
* tout devient solide.
*
*


OPTI DIME 2 ELEM qua4 MODE plan;
*
*************** Temperature de changement de phase ************************
Tfma = 200. ;
*
*************** Geometrie *************************************************
*p1 = 0. 0. 0. ; p11 = (0. 0. 10.e-3); p12=(0. 0. -10.e-3);
*p2 = 10.e-3 0. 0. ;
*p3 = 10.e-3 10.e-3 0. ;
*p4 = 0. 10.e-3 0. ;
p1 = 0. 0. ; p11 = (0.  10.e-3); p12=(0.  -10.e-3);
p2 = 10.e-3  0. ;
*p3 = 10.e-3 10.e-3 0. ;
*p4 = 0. 10.e-3 0. ;
n=1;
l1 = p12 droi 100  p1 droi   100 p11;
si graph;trac l1;finsi;
su1= l1 trans n p2;
*vol1= su1 volu trans n p4;
*trac cach  vol1;
**************** Modeles et caracteristiques *********************************
modth1 = MODE su1 'THERMIQUE' 'CONDUCTION' ;
modph1 = MODE su1 'CHANGEMENT_PHASE' 'INCO' 'T' 'Q';

latmin= 16.25e19;  latint = 16.85;latmax= 25.63e19;

LAT= 16.5e19;

LAT=1.e5;
CAPCALO= 900.;
matth1 = MATE modth1 'C' capcalo 'K' 100. 'RHO' 4500. ;
matph1 = MATE modph1  'PRIM' Tfma 'DUAL' (4500. * LAT);
*
*************** Pas de temps "grossier" et fin ********************************
* Durée totale
TT1 = 1.;
* pas de temps "grossier"
dt1 = 5.e-2 ;
* pas de temps "fin"
dt2 = 0.5e-6 ;
* pas de temps "très fin" juste pour avant la disparition totale de la zone fondue
dt3 = 0.5e-8 ;
*
* pas de temps du début "grossier" jusqu'à Temps1 et ensuite fin jusqu'à Temps2 ;
*
Temps1 = 0.350e-3 ;
Temps2 = 0.378e-3 ;
Temps3 = 0.385e-3 ;
listT = PROG 0. PAS dT1 Temps1 PAS dt2 Temps2 PAS dt3 Temps3
        PAS dt1 TT1 ;
listt= prog 0. pas dt1 tt1;
n1 = DIME listT ;
listS = PROG 0. pas 1. (n1-1) ;
listS = (1./(n1-1.)) * listS ;
EvS = EVOL MANU 'TEMP' listT 'SOUR' listS ;
si graph; dess EvS TITR 'evolution chargement'; finsi;
*
**************** Définition du champ de temperatures initial ****************
*xco = COOR 3 vol1 ;
xco = COOR 2 su1 ;
xco = ABS xco ;
max1 = MAXI xco ;
xco = max1 - xco ;
xco = (1./max1) * xco ;
*
* eventuellement, pour un pic de temperature plus localisé...
*xco = xco * xco ;
*
Txco = Tfma * 1.25 * xco ;
Txco = CHAN COMP 'T' Txco;
Txco = CHAN 'ATTRIBUT' Txco NATURE 'DIFFUS' ;
*TRAC vol1 Txco 'TITR' ' Temperature initiale' ;
si graph ;TRAC su1 Txco 'TITR' ' Temperature initiale' ; finsi;
evch = evol chpo l1 Txco;
 si graph ; dess evch; finsi;
********************************************************************
* chargement nul !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
********************************************************************
*sour0 = MANU CHPO vol1 1 'Q' 0. 'NATURE' 'DISC';
sour0 = MANU CHPO su1 1 'Q' 0. 'NATURE' 'DISC';
chthQ0 = CHAR 'Q' sour0 EvS ;
********************************************************************
xtab = table ;
*xtab.'CELSIUS'= VRAI ;
xtab.'TEMPERATURES' = table ;
xtab.'TEMPERATURES' . 0 = Txco ;
xtab.'MODELE' = modth1 et modph1 ;
xtab.'CARACTERISTIQUES' = matth1 et matph1 ;
xtab.'CHARGEMENT' = chthQ0 ;
xtab.'PROCEDURE_PERSO2' = VRAI ;
DEBP PERSO2 ;
  chT = PECHE xtab.'TEMPERATURES' ;
cht = exco cht 'T' 'T';
*  TRAC vol1 chT ;
  Tmax1 = MAXI chT ;
  titre ' profil de temperature le long du barreau';
  evch=  evol chpo l1 cht ;
  si graph ;    dess evch   nclk; finsi;
FINP ;
xtab.'TEMPS_CALCULES' = listT ;
*xtab.'TEMPS_SAUVES' = (PROG 0. TT1) ;
*xtab.'PAS_AJUSTE' = vrai ;
*xtab.'RELAXATION_THETA' = 1.;
*xtab.'CRITERE' = 1.e-3 ;
*
PASAPAS xtab ;



* on compare energie fournie/dépensée par chaleur latente dql
* à la différence d'énergie calorique % à la solution homogène 125°C
*  obtenue sans chaleur latente (dqt).

tder = peche xtab temperatures ;
tder= tder exco 'T';
evt = evol chpo tder l1;
so =( ( INTG evt) ) / 2.e-2;
list so;
evini= evol chpo txco l1;
sini= ( ( INTG evini) ) / 2.e-2;

dqt =  capcalo * ( abs( so - sini));

*proportions de phase, initiale, prise dans le cacul
pr0= xtab.proportions_phase . 0;
zo= intg pr0 modph1;
zz = manu chml modph1 'SCAL' 1.;

zotot= intg zz modph1;
na = dime xtab.proportions_phase ;
prfi= xtab . proportions_phase . (na - 1);
zfi= intg prfi modph1;
zmod = abs ( zfi - zo ) ;

propchzo=abs( zmod / zotot);
dql = LAT  * propchzo  ;

mess ' dqt'  dqt  ' dql ' dql;

err= abs (  ( dql- dqt ) / dql);
si ( err > 1.e-4) ;
mess ' erreur dans phase2d.dgibi' ;
erreur 5 ; finsi;

fin;