TITLE 'Solidificazione di una lastra di acqua'
COORDINATES cartesian1  { monodimensionale}
VARIABLES
temp(0.01)
solid(0.01)
! SELECT         { method controls }
DEFINITIONS    { parametri }
Tw=253    { Temperatura della parete (solido) }
Tinf=293  { Temperatura iniziale del fuso }
Tm=273  { Temperatura di transizione }
T0=1   { Ampiezza della Gaussiana (il materiale e' solido a Tm-2*T0, e' fuso a Tm+2*T0 }
S0 =( 1-erf((Temp-Tm)/T0))/2    { gaussiana }
lambda=2.2*s0+.6*(1-s0)   { conducibilita' (valori per ghiaccio e acqua) W/mK}
rho=900*s0+1000*(1-s0)   { densita' (valori per ghiaccio e acqua) Kg/m3}
Cp=2100*s0+4200*(1-s0)   { calore specifico (valori per ghiaccio e acqua) J/KgK}
alfas=2.2/900/2100   { diffusivita' termica del solido  m2/s}
alfal=0.6/1000/4200   { diffusivita' termica del fuso  m2/s}
Qm=335000      { calore latente J/Kg  }
L=.1   { Spessore della lastra in m  }
 
{ Soluzione analitica  }
r=(alfas/alfal)^.5
gam=0.23936
delta=gam*(4*alfas*t)^.5
delta1=delta*.82
Tsol=Tw+(Tm-Tw)*ERF(x/(4*alfas*t)^.5)/ERF(gam)
Tmel=Tinf-(Tinf-Tm)*ERFC(x/(4*alfal*t)^.5)/ERFC(gam*r)
Tsoluz=if x<delta then Tsol else Tmel
solsol=if x<delta then 1 else 0
solsol1=if x<delta1 then 1 else 0
 
{ valori per soluzione numerica  }
response = 0.01   {tempo di risposta per la soluzione numerica}
eps=1e-6    {coeff. for numerical diffusion term}
 
INITIAL VALUES
temp=tinf
solid =(1-erf((Tinf-Tm)/T0))/2
 
EQUATIONS
 
temp:   rho*cp*dt(temp) - div(lambda*grad(temp))  = rho*Qm*dt(solid)
 
! l'equazione seguente costringe solid ad avvicinarsi il piu' possibile a S0,
! rendendola pero' funzione del tempo secondo la discretizzazione data dal programma
!se "solid" si allontana da S0 (il valore che dovrebbe avere), il secondo membro lo forza
!a ritornare sul valore corretto
 
solid:   dt(solid) - eps*div(grad(solid)) =  (S0-solid)/response
 
 
BOUNDARIES       { The domain definition }
  REGION 1
START(0) POINT VALUE(temp)=Tw  {parete di sinistra a Tw}
LINE TO (L) natural(temp)=0   {parete di destra adiabatica}
 
TIME 0 TO 3600
 
MONITORS         { show progress }
For cycle=10
elevation(temp,Tsoluz) from (0) to (L)
elevation(solid, solsol)  from (0) to (L)
elevation(Cp)  from (0) to (L)
elevation(lambda)  from (0) to (L)
elevation(rho)  from (0) to (L)
PLOTS
 
END