Phasechange.py
From GeoMod
To run the code, once you have VPython installed, you can either.
- Download this file: water_phasechange_06.py
- or, copy and paste the following code to a file (named "water_phasechange_06.py" for example) using your favorite text editor.
from visual import *
from visual.graph import * # import graphing features
class beaker_label:
"event details"
def __init__(self, x, y, z):
dx = 0.1
self.lhash = curve(color=(1,0,0))
self.lhash.append(pos=(x, y, z))
self.lhash.append(pos=(x+dx, y, z))
self.lab = label(pos=(x-dx,y,z), text=str(y),height=10, border=2, box = 0, opacity=0,color=(1,0,0))
def water_density(T,S,P):
"based on Temperature (Celcius) and Salinity (ppt) and Pressure (bars) from http://web.maths.unsw.edu.au/~bxs/DieCAST/MANUAL/density.html"
p1 = P + 5884.81703666 + T*(39.803732+T*(-0.3191477+T*0.0004291133))+ 2.6126277*S
p2 = 1747.4508988 + T*(11.51588 -0.046331033*T) - S*(3.85429655+0.01353985*T)
RHO = p1/(p2+0.7028423*p1)
return RHO
def ice_density(T):
"density of ice based on Temperature"
alpha = 50e-6
rho_0 = 0.917
return rho_0 + T * alpha
def slide_val(bx, val_min, val_max, xmin, xmax, step = 1.0):
"returns the value from a slider"
r = (val_max - val_min) * (bx-xmin) / (xmax-xmin)
tr = round(r/step,0) * step
bn = val_min + (tr)
return bn
def fill_beaker(T, i, w):
## print `i`
## print `w`
## print ' '
amp = 1.0
if T == 0.0:
if w <= 0.0: #water first since it is denser
water.axis = (0,-1.0,0)
water.visible = 0
elif w > 0.0:
water.axis = (0,w / water_density(T,0.0,1.0),0)
water.visible = 1
if i <= 0.0:
ice.axis = (0,-1.0,0)
ice.visible = 0
elif i > 0.0:
ice.axis = (0,i / ice_density(T),0)
ice.visible = 1
elif T < 0.0:
#print 'less'
water.axis = (0,-1.0,0)
water.visible = 0
ice.axis = (0,i / ice_density(T),0)
ice.visible = 1
elif T > 0.0 and T < 100.0:
ice.axis = (0,-1.0,0)
ice.visible = 0
water.axis = (0,w / water_density(T, 0.0, 1.0),0)
water.visible = 1
elif T == 100.0:
water.axis = (0,w / water_density(T, 0.0, 1.0),0)
#print `T` + " " + `ice.axis.y` + " " + `ice_density(T)`
#amplify
water.axis.y = water.axis.y * amp
ice.axis.y = ice.axis.y * amp
if water.axis.y <= 0.0:
ice.pos.y = 0.0
else:
ice.pos.y = water.axis.y
w = 640+4+4
h = 480+24+4
scene.x = 0
scene.y = 0
scene.width = w
scene.height = h
hotplate = cylinder(pos=(0,-0.1,0), axis=(0,-0.5,0), color=(1,0,0), radius = 1.5)
scene.center = (0,1.5,0)
scene.range = 3.75
Tctrl = curve(radius=0.05)
Tctrl.append(pos=(-hotplate.radius/4.0, hotplate.axis.y * 0.5, hotplate.radius), color=(0,0,1))
#Tctrl.append(color=(0,0,1))
Tctrl.append(pos=(0.0, hotplate.axis.y * 0.5, hotplate.radius), color=(1,1,1))
#Tctrl.append(color=(1,1,1))
Tctrl.append(pos=(hotplate.radius/4.0, hotplate.axis.y * 0.5, hotplate.radius), color=(1,0,0))
#Tctrl.append(color=(1,0,0))
Tbutt = sphere(pos=Tctrl.pos[1], radius = 0.1, color = (1,1,0))
Tzerobut = sphere(pos=Tctrl.pos[1]+(0,0.2,0), radius = 0.1, color = (0,1,0))
butstep = 10.0 # steps by which you can move the button
bh = 6 * hotplate.axis.y
brad = 1
bth = 0.01
beaker = frame()
beakbot = ring( radius = brad, pos=(0,0,0), thickness=bth, axis=hotplate.axis)
beaktop = ring( radius = brad, pos=bh*norm(hotplate.axis), thickness=bth, axis=hotplate.axis)
beakl = cylinder(pos = (-brad,0,0), radius=bth, axis = bh * norm(hotplate.axis))
beakr = cylinder(pos = (brad,0,0), radius=bth, axis = bh * norm(hotplate.axis))
tot_mass = 2.0
water = cylinder(pos=(0,0,0), axis=-hotplate.axis, radius = brad,color=(.5,.5,1), mass=0.0)
ice = cylinder(pos=(0,0,0), axis=-hotplate.axis, radius = brad, color=(1,1,1), mass = tot_mass)
nm = 3.0
dm = -bh / nm
print `dm`
measure = []
for n in arange(dm, -bh, dm):
print `n`
measure.append(beaker_label(0,n,beaktop.radius))
##oneline = curve()
##oneline.append(pos=(beakbot.radius*1.05, 0, 0))
##oneline.append(pos=(beakbot.radius*1.05, 1, 0))
##for i in arange(0,100,1):
## print `i` + " " + `water_density(i,0.0,1.0)`
##for i in arange(-100,0,1):
## print `i` + " " + `ice_density(i)`
#define specific heat capacities J/Kelvin.g
c_ice = 1.960
c_water = 4.22
l_fuse = 330 #latent heat of fusion/melting (J/g)
l_vape = 2500 #latent heat of vaporization (J/g)
Tinit = -20.0 #start temperature
framerate = 100.0
dt = 1.0/60.0 # timestep (seconds)
dH_min = -40.0
dH_max = 40.0
dH = slide_val(Tbutt.x, dH_min, dH_max, Tctrl.x[0], Tctrl.x[2])
#dH = dH_min + ((dH_max - dH_min) * (Tbutt.x-Tctrl.x[0]) / (Tctrl.x[2]-Tctrl.x[0])) #heat added (J/ second)
ice_phase = 1.0
water_phase = 0.0
T= Tinit
Hlabel = label(pos=Tctrl.pos[1] - (0,0.4,0),height=10, border=2, text="Heat = "+`round(dH)`+" J/s")
Tlabel = label(pos=(0,1.25 * beaktop.pos.y,0), text="Temperature = " + `round(T)`,height=20, border=2)
ilabel = label(pos=(1.4,0.5,0),height=10, border=2, text="Mass="+`round(tot_mass)`)
fill_beaker(T,ice.mass, water.mass)
#create graph display
gdisp = gdisplay(x=0, y=h, width=w, height=h/2,
ytitle="Temperature", xtitle="time")
# ymin=-50, ymax=110, xmax = 100)
Tcurve = gcurve(color=color.cyan) # a connected curve object
##print `T`
##print `dH`
##print `dt`
##print `c_ice`
pick = None
tm = 0
while 1:
rate(framerate)
tm = tm + dt
if scene.mouse.events:
m1 = scene.mouse.getevent() # obtain drag or drop event
if m1.drag and (m1.pick == Tbutt):
drag_pos = m1.pickpos
pick = m1.pick
scene.cursor.visible = 0 # make cursor invisible
elif m1.drop:
pick = None # end dragging
scene.cursor.visible = 1 # cursor visible
if m1.pick == Tzerobut and m1.release:
Tbutt.pos.x = round((Tctrl.x[0]+Tctrl.x[2])/2.0,0)
#dH = slide_val(Tbutt.x, dH_min, dH_max, Tctrl.x[0], Tctrl.x[2],10.0)
if pick:
new_pos = scene.mouse.project(normal=(0,0,1))
if new_pos != drag_pos:
pick.pos += new_pos - drag_pos
drag_pos = new_pos
if pick == Tbutt:
pick.pos.y = Tctrl.y[1]
pick.pos.z = Tctrl.z[1]
if pick.pos.x < Tctrl.x[0]:
pick.pos.x = Tctrl.x[0]
if pick.pos.x > Tctrl.x[2]:
pick.pos.x = Tctrl.x[2]
#print dH
dH = slide_val(Tbutt.pos.x, dH_min, dH_max, Tctrl.x[0], Tctrl.x[2],10.0)
#dH = dH * framerate/60.0 #adjustment for slow movie recording speed
Hlabel.text= "Heat = "+`dH`+" J/s"
if dH <= 0:
hotplate.color = (0.8*(1-dH/dH_min),1-dH/dH_min,1)
else:
hotplate.color = (1,1-dH/dH_max,1-dH/dH_max)
## r = (val_max - val_min) * (bx-xmin) / (xmax-xmin)
## tr = round(r/step,0) * step
## bn = val_min + (tr)
Tlabel.text="Temperature = " + `round(T)`
# ilabel.text="Ice phase = " + `ice_phase`
if T < 0.0 and T > -272.0:
T = T + dH * dt / (c_ice * ice.mass)
if T > 0.0:
T = 0.0
ice.mass = tot_mass
water.mass = 0.0
elif T > 0.0 and T < 100.0:
T = T + dH * dt / (c_water * water.mass)
if T < 0.0:
T = 0.0
ice.mass = 0.0
water.mass = tot_mass
elif T == 0.0:
#print `ice_phase`
if ice.mass > 0.0 and ice.mass <= tot_mass: #there is some ice in the beaker
T = 0.0
dp = dH * dt / l_fuse #change in phase is positive if heat is added
ice.mass = ice.mass - dp
water.mass = tot_mass - ice.mass
if ice.mass <= 0.0:
ice.mass = 0.0
water.mass = tot_mass
T = T + dH * dt / (c_water * water.mass)
elif ice.mass >= tot_mass:
ice.mass = tot_mass
water.mass = 0.0
T = T + dH * dt / (c_ice * ice.mass)
else: # no ice so either freeze or start warming water
if dH < 0.0: #then freeze
dp = dH * dt / l_fuse #change in phase is positive if heat is added
ice.mass = ice.mass - dp
water.mass = tot_mass - ice.mass
T = 0.0
else: #warm water
ice.mass = 0.0
water.mass = tot_mass
T = T + dH * dt / (c_water * water.mass)
elif T >= 100.0:
if water.mass > 0.0 and water.mass <= tot_mass:
T = 100.0
dp = dH * dt / l_vape
water.mass = water.mass - dp
ice.mass = 0.0
if water.mass <= 0.0:
water.mass = 0.0
ice.mass = 0.0
T = T + dH * dt / (1.0)
elif water.mass >= tot_mass:
ice.mass = 0.0
T = T + dH * dt / (c_water * water.mass)
elif T <= -272.0:
if dH <= 0.0:
T = -272.0
else:
T = T + dH * dt / (c_ice * ice.mass)
fill_beaker(T,ice.mass, water.mass)
Tcurve.plot(pos=(tm, T))
