-
Notifications
You must be signed in to change notification settings - Fork 2
/
planetThermo.py
165 lines (140 loc) · 4.11 KB
/
planetThermo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
"""
Thermodyanmic functions and constants for planetary science
"""
# Load Python dependencies
import numpy as np
######################
# Physical constants #
######################
# Pi
pi = np.pi
# Proton mass
mp = 1.6726e-27
# Molecular masses [kg]
m_h2o = 3.0132e-26
# Molar masses [kg.mole-1]
mw_h2o = 18.015e-3 # Water
# Thermodynamic constants
R = 8.31447 # Universal gas constant [J.mol-1.K-1]
k_B = 1.38064852e-23 # Boltzmann constant [m2.kg.s-2.K-1]
class volatile:
'''
A volatile is a molecular species that undergoes
phase transitions at planetary surface and/or atmospheric
temperatures
'''
#__repr__ object prints out a help string when help is
#invoked on the volatile object or the volatile name is typed
def __repr__(self):
line1 =\
'This volatile object contains information on %s\n'%self.name
line2 = 'Type \"help(volatile)\" for more information\n'
return line1+line2
def __init__(self):
self.name = None #Name of the volatile
########
# ph2o #
########
# -----------------------------------------------
# Equilibrium vapor pressure over water over ice
# Marti and Mauersberger (1993)
# -----------------------------------------------
# Input:
# T = temperature of ice [K]
# Output:
# vapor pressure [Pa]
def ph2o(T):
A = -2663.5
B = 12.537
p = 10**(A/T + B)
return p
########
# pco2 #
########
# --------------------------------------------
# Equilibrium vapor pressure over solid CO2
# Brown and Ziegler (1980)
# --------------------------------------------
# Input:
# T = temperature of solid [K]
# Output:
# vapor pressure [Pa]
def pco2(T):
A0 = 2.13807649*10**1
A1 = -2.57064700*10**3
A2 = -7.78129489*10**4
A3 = 4.32506256*10**6
A4 = -1.20671368*10**8
A5 = 1.34966306*10**9
# Pressure in torr
ptorr = np.exp(A0 + A1/T + A2/T**2 + A3/T**3 + A4/T**4 + A5/T**5)
# Pressure in Pa
p = 133.3223684211*ptorr
return p
########
# th2o #
########
# --------------------------------------------
# Equilibrium H2O frost point at given partial
# pressure
# Marti and Mauersberger (1993)
# --------------------------------------------
# Input:
# p = vapor pressure [Pa]
# Output:
# temperature of solid [K]
def th2o(p):
# Definitions
dT = 0.1 # Temperature precision
Tmin = 20 # Minimum temperature to try
Tmax = 600.0 # Maximum temperature to try
# Temperature and pressure arrays
T_arr = np.arange(Tmin, Tmax, dT)
p_arr = ph2o(T_arr)
# Interpolate to find T
T = np.interp(p, p_arr, T_arr)
return T
########
# tco2 #
########
# --------------------------------------------
# Equilibrium CO2 frost point at given partial
# pressure
# Brown and Ziegler (1980)
# --------------------------------------------
# Input:
# p = vapor pressure [Pa]
# Output:
# temperature of solid [K]
def tco2(p):
# Definitions
dT = 0.1 # Temperature precision
Tmin = 20 # Minimum temperature to try
Tmax = 200.0 # Maximum temperature to try
# Temperature and pressure arrays
T_arr = np.arange(Tmin, Tmax, dT)
p_arr = pco2(T_arr)
# Interpolate to find T
T = np.interp(p, p_arr, T_arr)
return T
######################
# h2oSublimationRate #
######################
# -------------------------------------------------------------------------
# Water ice sublimation rate, in kg.m-2.s-1, based on planar surface
# solution from Estermann (1955): mdot = pv*sqrt(mw_h2o/2piRT), where mh2o
# is the molecular weight of water, R is the universal gas constant, T is
# temperature in Kelvin, and pv is the saturation vapor pressure. This is
# the upper limit on sublimation rate, where the actual sublimation rate
# depends on the difference (p-pv) and the rate of adsorption.
# -------------------------------------------------------------------------
# Input:
# T = temperature of ice [K]
# Output:
# sublimation rate [kg.m-2.s-1]
def h2oSublimationRate(T):
# Equilibrium water vapor pressure over ice
pv = ph2o(T) # [Pa]
# Sublimation rate
mdot = pv * np.sqrt( mw_h2o/(2*pi*R*T) ) # [kg.m-2.s-1]
return mdot