Hello,
I have looked through a lot of documentation and emails but I am not sure I
have correctly constructed the equation for heat flow in a conductor where
the source term and the coefficient of the diffusion term depend on the
dependent variable.
\frac{\partial u}{\partial t} = D(u)\frac{\partial^2 u}{\partial x^2} + S(u)
[image: Inline image 2]
eq = fp.TransientTerm(coeff=1.0) == \
fp.DiffusionTerm(coeff=(k[areas[0]] / (C[areas[0]] *
hx.DENSITY_CU))) + \
fp.ImplicitSourceTerm(coeff=((r[areas[0]] * I ** 2) / (areas[0] **
2 * hx.DENSITY_CU * C[areas[0]])))
I have several sample areas that I use, hence the indexing of the
cellVariables, r, k, and C. I declare/initialize them like this:
C = fp.CellVariable(value=hx.specific_heat2(u.value), name='Specific Heat',
mesh=mesh)
Should this be a faceVariable? These temperature dependent variables are
updated in my solution loop. The loop algorithm is:
Outer Loop: while elapsed time < final time:
1 updateOld
2 solve (from Richards Equation email)
3 update r, k, and C face variables
4 update/show plot
5 Inner Loop: while residue > tolerance:
1 sweep
2 update r, k, and C face variables
End Inner Loop
6 update elapsed time
End Outer Loop
I performed an analysis to choose the area of my conductor and I shouldn't
see a rise above 493 K for 1500 seconds with a cross-sectional area of
0.196 cm^2. The simulations I've done in FiPy produce very large
temperatures using the parameters from my earlier analysis. I have attached
the code used to make the following graph.
[image: Inline image 3]\
Thank you in advance for any help.
Cheers,
Chris
__author__ = 'Chris Jones'
__maintainer__ = 'Chris Jones'
__email__ = '[email protected]'
__status__ = 'Development'
"""
heat_xfer.py
Uses numpy 1.9.1 and written for python 3.4
Functions for thermal conductivity of copper and aluminum. Expansion
coefficient,
resistivity and specific heat for copper.
"""
import numpy as np
DENSITY_CU = 8.96 # g / cm^3
INTR_RESIST_CU = 1.71e-6 # Ohm * cm
ALPHA_CU = 3.86e-3 # 1/K ... Should this be 0.00386?6.8e-3
def therm_cond_cu(t, rrr):
"""
Returns thermal conductivity of OFHC copper for
a given t in Kelvin and purity ratio rrr.
Returned values have units of W/(cm K).
Valid between 4 K and 300 K only.
t must be an integer or numpy array
:rtype : float
Valid for the following rrr values: 50, 100, 150
Source
http://cryogenics.nist.gov/MPropsMAY/OFHC%20Copper/OFHC_Copper_rev.htm
"""
a = {50: 1.8743, 100: 2.2154, 150: 2.3797}
b = {50: -0.41538, 100: -0.47461, 150: -0.4918}
c = {50: -0.6018, 100: -0.88068, 150: -0.98615}
d = {50: 0.13294, 100: 0.13871, 150: 0.13942}
e = {50: 0.26426, 100: 0.29505, 150: 0.30475}
f = {50: -0.0219, 100: -0.02043, 150: -0.019713}
g = {50: -0.051276, 100: -0.04831, 150: -0.046897}
h = {50: 0.0014871, 100: 0.001281, 150: 0.0011969}
i = {50: 0.003723, 100: 0.003207, 150: 0.0029988}
# if rrr not in (50, 100, 150):
# return 0
#
# if type(t) == np.ndarray and (t.min() < 4 or t.max() > 300):
# return 0
#
# if (type(t) == int or type(t) == float) and (t < 4 or t > 300):
# return 0
return (1 / 100) * 10 ** ((a[rrr] + c[rrr] * t ** 0.5 + e[rrr] * t + g[rrr]
* t ** 1.5 + i[rrr] * t ** 2) /
(1 + b[rrr] * t ** 0.5 + d[rrr] * t + f[rrr] * t
** 1.5 + h[rrr] * t ** 2))
def therm_cond_al(t):
"""
Returns thermal conductivity of 6061-T6 Aluminum for
a given t in Kelvin.
Returned values have units of W/(cm K).
Valid between 4 K and 300 K only.
t must be an integer or numpy array
:rtype : float
Source http://cryogenics.nist.gov/Papers/Cryo_Materials.pdf
"""
a = 0.07918
b = 1.09570
c = -0.07277
d = 0.08084
e = 0.02803
f = -0.09464
g = 0.04179
h = -0.00571
# if type(t) == np.ndarray and (t.min() < 4 or t.max() > 300):
# return 0
#
# if (type(t) == int or type(t) == float) and (t < 4 or t > 300):
# return 0
return (1 / 100) * 10 ** (a + b * np.log10(t) + c * np.log10(t) ** 2 + d *
np.log10(t) ** 3
+ e * np.log10(t) ** 4 + f * np.log10(t) ** 5 + g *
np.log10(t) ** 6 + h * np.log10(t) ** 7)
def specific_heat1(t):
"""
Returns specific heat of OFHC copper for a given t in Kelvin.
Returned values have units of J/(g*K).
Valid between 4 K and 300 K only.
Source
http://cryogenics.nist.gov/MPropsMAY/OFHC%20Copper/OFHC_Copper_rev.htm
"""
# Polynomial Coefficients
a = -1.91844
b = -0.15973
c = 8.61013
d = -18.996
e = 21.9661
f = -12.7328
g = 3.54322
h = -0.3797
# Polynomial
return (1 / 1000) * 10 ** (a + b * np.log10(t) + c * np.log10(t) ** 2 + d *
np.log10(t) ** 3
+ e * np.log10(t) ** 4 + f * np.log10(t) ** 5 + g *
np.log10(t) ** 6 + h * np.log10(t) ** 7)
def expansion_coeff(t):
"""
Returns expansion coefficient of OFHC copper for a given t in Kelvin.
Returned values have units of 1/K.
Valid between 4 K and 300 K only.
Source
http://cryogenics.nist.gov/MPropsMAY/OFHC%20Copper/OFHC_Copper_rev.htm
"""
# Polynomial Coefficients
a = -17.9081289
b = 67.131914
c = -118.809316
d = 109.9845997
e = -53.8696089
f = 13.30247491
g = -1.30843441
# if not 4 <= t.all() <= 300:
# return 0
# Polynomial
return (1 / 1000) * 10 ** (a + b * np.log10(t) + c * np.log10(t) ** 2 + d *
np.log10(t) ** 3
+ e * np.log10(t) ** 4 + f * np.log10(t) ** 5 + g *
np.log10(t) ** 6)
def resistivity(t, r0=INTR_RESIST_CU, a0=ALPHA_CU):
"""
Returns resistivity at a temperature t for a given temperature coefficient
a0
and a given intrinsic resistivity r0.
Return units are Ohm * cm
:param t: K
:param r0: 1.71e-6 Ohm * cm
:param a0: 3.86e-3 1/K
:return: Ohm * cm
"""
return r0 * (a0 * t + 1)
def d_ratio(t, rrr):
"""
Returns the ratio of k to c. Used in nonlinear heat diffusion.
:param t:
:param rrr:
:return:
"""
return therm_cond_cu(t, rrr) / specific_heat1(t)
def specific_heat2(b):
a = []
c = []
u_min = 0.0
temp = b / 100.0
for u in temp:
if 0.29 <= u < 0.4473:
a.extend([1.528, 15.929, 47.676, -89.7484])
u_min = 0.29
elif u < 0.692:
a.extend([4.864, 24.266, 5.3237, -29.051])
u_min = 0.4473
elif u < 1.3946:
a.extend([10.695, 21.653, -16.003, 5.2425])
u_min = 0.692
elif u < 2.0:
a.extend([19.827, 6.93, -4.9524, 1.8736])
u_min = 1.3946
elif u < 3.3:
a.extend([22.623, 2.9936, -1.5496, 0.3724])
u_min = 2.0
elif u < 12.37:
a.extend([24.714, 0.8526, -0.09737, 0.00873])
u_min = 3.3
else:
raise ValueError('Temperature must be between 29 and 1237 Kelvin.')
u -= u_min
c_J_per_mol_K = a[0] + a[1] * u + a[2] * u ** 2 + a[3] * u ** 3
c_J_per_gram_K = c_J_per_mol_K / 63.54 # Atomic Weight Cu = 65.34 g/mol
c.extend([c_J_per_gram_K])
a = []
return c__author__ = 'Chris Jones'
import sys
import fipy as fp
import matplotlib.pyplot as plt
sys.path.insert(0, 'C:/Users/Chris
Jones/Documents/PythonProjects/OFHC_Therm_Cond')
import heat_xfer as hx
areas = [10., 5.] # cm^2. First value in areas used for transient analysis
steadyState = 1
transient = 0
finalTime = 150.0
dt = 1e-1
tol = 1e-8
elapsedTime = 0.0
sweeps = 0
maxSweeps = 100
# Our working constants
I = 200 # Electric current in Amps
# a = 0.19635 # Conductor Area in cm^2
# Setup and create a grid
nx = 100 # cm - number of cells
lx = 50 # cm - length of rod
mesh = fp.Grid1D(nx=nx, Lx=lx)
m = mesh.cellCenters()[0]
# Our initial boundary conditions
valueLeft = 40.0 # K
valueRight = 293.0 # K
u = {}
C = {}
k = {}
r = {}
p = {}
l = {}
for a in areas:
u[a] = fp.CellVariable(mesh=mesh, name='Area:{:.3f}'.format(a), hasOld=1)
u[a].value = fp.numerix.logspace(fp.numerix.log10(valueLeft),
fp.numerix.log10(valueRight), num=u[a].shape[0])
C[a] = fp.CellVariable(value=hx.specific_heat2(u[a].value), name='Specific
Heat', mesh=mesh)
k[a] = fp.CellVariable(value=hx.therm_cond_cu(u[a].value, 150),
name='Thermal Conductivity', mesh=mesh)
r[a] = fp.CellVariable(value=hx.resistivity(u[a].value), mesh=mesh)
u[a].constrain(value=valueLeft, where=mesh.facesLeft)
if steadyState:
for a in areas:
(fp.DiffusionTerm(coeff=(k[a] / (C[a] * hx.DENSITY_CU)))
+ fp.ImplicitSourceTerm(coeff=((r[a] * I ** 2)
/ (a ** 2 * hx.DENSITY_CU *
C[a])))).solve(var=u[a])
p[a], = plt.plot(m, u[a].value)
l[a] = r'Area:{:.3f} $cm^2$'.format(a)
plt.legend(([x for l, x in p.items()]), ([x for l, x in l.items()]),
loc='upper left')
plt.xlabel('Position along Conductor (cm)')
plt.ylabel('Degrees (K)')
plt.title('Steady State Temperature Distributions by Conductor Area')
plt.grid(True)
plt.show()
# End Steady State
if transient:
eq = fp.TransientTerm(coeff=1.0) == \
fp.DiffusionTerm(coeff=(k[areas[0]] / (C[areas[0]] * hx.DENSITY_CU)))
+ \
fp.ImplicitSourceTerm(coeff=((r[areas[0]] * I ** 2) / (areas[0] ** 2 *
hx.DENSITY_CU * C[areas[0]])))
fig = plt.figure()
ax = fig.gca()
fig.show()
elapsedTime = 0.0
# Outer Loop - Time Step
while elapsedTime <= finalTime:
res = 1
# print('Before UpdateOld: {}'.format(u[areas[0]].value[10]))
u[areas[0]].updateOld()
# print(' After UpdateOld: {}'.format(u[areas[0]].value[10]))
eq.solve(var=u[areas[0]], dt=dt)
print('Solve: {:.3f}'.format(u[areas[0]].value[nx - 1]))
C[areas[0]].value = hx.specific_heat2(u[areas[0]].value)
k[areas[0]].value = hx.therm_cond_cu(u[areas[0]].value, 150)
r[areas[0]].value = hx.resistivity(u[areas[0]].value)
ax.cla()
ax.plot(m, u[areas[0]].value)
ax.grid(True)
fig.canvas.draw()
# Inner loop - Refine Time Step
while res > tol:
try:
res = eq.sweep(var=u[areas[0]], dt=dt)
C[areas[0]].value = hx.specific_heat2(u[areas[0]].value)
k[areas[0]].value = hx.therm_cond_cu(u[areas[0]].value, 150)
r[areas[0]].value = hx.resistivity(u[areas[0]].value)
except RuntimeError:
print('Exception occurred. Ending simulation.')
print('Total time: {:.1f}s'.format(elapsedTime))
print('nx={:.3f}_dt={:.5f}s_res<={}'
.format(nx, dt, res))
quit()
sweeps += 1
if sweeps >= maxSweeps:
raise Exception('Maximum sweeps reached.')
print('Elapsed Time: {:.3f}'.format(elapsedTime))
print('Sweeps {}: {:.3f}'.format(sweeps, u[areas[0]].value[nx - 1]))
sweeps = 0
elapsedTime += dt_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
