Melting

{  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