Melting

Previous  Top  Next

{ MELTING.PDE }
{ ************************************************************
This problem shows the application of FlexPDE to the melting of metal.

We choose as our system variables the temperature, "temp", and the
fraction of material which is solid at any point, "solid".

The temperature is given by the heat equation,

rho*cp*dt(temp) - div(lambda*grad(temp)) = Source

where cp is the heat capacity, rho the density and lambda the conductivity.

The latent heat, Qm, is an amount of heat released as "Solid" changes from
zero to one. We have Qm = integral[0,1]((dH/dSolid)*dSolid), or assuming
dH/dSolid is constant, dH/dSolid = Qm.
Then heat source from freezing is

dH/dt = (dH/dSolid)*(dSolid/dt) = Qm*dt(Solid).

We assume that the solid fraction can be represented by a function

solid = 0.5*erfc((temp-Tm)/T0)

where Tm is the melting temperature, and T0 is a temperature range over
which the melting transition occurs. Since there are no spatial derivatives
in this equation, we introduce a diffusion term with small coefficient to act
as a noise filter.

The particular problem addressed here is a disk of cold solid material
immersed in a bath of liquid. The initial temperatures are such that material
first freezes onto the disk, but after equilibrium is reached all the material
is liquid. The outer boundary is insulated.

Since the initial condition is a discontinuous distribution, we use a separate
REGION to define the cold initial disk, so that the grid lines will follow the
shape. We also add a FEATURE bounding the disk to help the gridder
define the abrupt transition. SELECT SMOOTHINIT helps minimize
interpolator overshoots.

****************************************************** }

TITLE
'Melting Metal'

COORDINATES
ycylinder('r','z')

SELECT
cubic
smoothinit

VARIABLES
temp(threshold=10)
solid(threshold=0.1)

DEFINITIONS

Qm= 225000 { latent heat }
Tm=1850 { Melting temperature }
T0= 20 { Melting interval +- T0 }
temp_liq=2000 { initial liquid temperature }
temp_sol=400 { initial solid temperature }
Tinit
sinit

R_inf = 0.7 { Domain Radius m}

{ plate }
d=0.05
dd = d/5 { a defining layer around discontinuity }
R_Plate=0.15


lambda = 30+4.5e-5*(temp-1350)^2 { Conductivity }
rho=2500 { Density kg/m3 }
cp = 700 { heat capacity }

INITIAL VALUES
temp=Tinit
solid = 0.5*erfc((tinit-Tm)/T0)

EQUATIONS
temp:   rho*cp*dt(temp) - div(lambda*grad(temp)) = Qm*dt(solid)
solid:   solid - 1e-6*div(grad(solid)) = 0.5*erfc((temp-Tm)/T0)

BOUNDARIES

region 'Outer'
Tinit = temp_liq
sinit = 0

start 'outer' (0,-R_inf)
value(temp)= temp_liq   arc(center=0,0) angle 180
natural(temp)=0      line to close

region 'Plate'
Tinit = temp_sol
sinit = 1

start(0,0)
mesh_spacing=dd
line to (R_Plate,0) to (R_Plate,d) to (0,d) to close

TIME 0 by 1e-5 to 600

MONITORS
for cycle=10
grid(r,z) zoom (0,-0.1,0.25,0.25)
elevation(temp) from(0.1,-0.1) to (0.1,0.15) range=(0,2000)
elevation(solid) from(0.1,-0.1) to (0.1,0.15) range=(0,1)

PLOTS
for t= 0 1e-4 1e-3 1e-2 0.1 1 10 by 10 to 100 by 100 to 300 by 300 to endtime

contour( temp) range=(0,2000)
contour( temp) zoom (0,-0.2,0.45,0.45) range=(0,2000)
elevation(temp) from(0.1,-0.1) to (0.1,0.15) range=(0,2000)
contour(solid) range=(0,1)
contour(solid) zoom (0,-0.2,0.45,0.45) range=(0,1)
elevation(solid) from(0.1,-0.1) to (0.1,0.15) range=(0,1)

HISTORIES

history(temp) at (0.051,d/2) at (0.075,d/2) at (R_plate,d/2)
history(temp) at (0.051,d) at (0.075,d) at (R_plate,d)
history(solid) at (0.051,d/2) at (0.075,d/2) at (R_plate,d/2)
history(solid) at (0.051,d) at (0.075,d) at (R_plate,d)
history(integral(cp*temp+Qm*(1-solid))) as "Total Energy"


END