Dear FiPy users and developers,
could you help me with how to implement a thermal contact, hc
[W/(m^2.K)], between two regions : so there is like a discontinuity in
the temperature between the two regions.
In some cases I could implement it as a thin layer, dx [m], with a given
thermal diffusion coefficient, lambda [W/(m.K)], such that hc =
lambda/dx.
However in some cases this distords the geometry somewhat and moreover
the results seem to depend rather strongly on the details of the mesh.
So, I was wondering if that contact could be included somehow in a sort
of BC.
In the example below I have a thermal contact of hc = 1 kW/(m^2.K)
between two slabs (of length L/2) of good thermal conductors K = 200
W/(m.K). On side is held to 0 C while on the other side a power flux Ps
= 1 kW/m^2 is applied.
The temperature jump theoretically should be Ps/hc = 1 K (deg C) for
stationary conditions.
I made a little table of the temperature jump estimated with FiPy vs.
the meshing size (cellSize) and the width, dx, of the layer that
simulates the thermal contact, hc.
+------------+------------+-------------------------+
| | | delta T |
| cellSize/L | dx/L | |
| | | FiPy | Theory |
+------------+------------+------------+------------+
| 0.01000 | 0.01000 | 1.05961 | 1.00000 |
+------------+------------+------------+------------+
| 0.00500 | 0.01000 | 0.63490 | 1.00000 |
+------------+------------+------------+------------+
| 0.00500 | 0.00500 | 1.06939 | 1.00000 |
+------------+------------+------------+------------+
| 0.01000 | 0.00500 | 1.37685 | 1.00000 |
+------------+------------+------------+------------+
"""
test FiPy
"""
import fipy as fp, numpy as np, pylab as pl
# geometry parameters
L = 1. # lenght of the slab
W = 0.4 # width of the slab
dx = 0.005 * L # thickness of the thermal layer
cellSize = 0.01 * L # cellSize function of problem
# define the geometry
geo = """
cellSize = %(cellSize)g;
L = %(L)g;
W = %(W)g;
dx = %(dx)g;
// 7 6 5
// 8 <-------- 7 <------ 6 <-------- 5
// | ^ | ^
// | Region 1 | thermal | Region 2 |
// 8 | | contact | | 4
// | 10 | | 11 |
// v | v |
// 1 --------> 2 ------> 3 --------> 4
// 1 2 3
Point(1) = {- L/2, 0, 0, cellSize};
Point(2) = {-dx/2, 0, 0, cellSize}; Line(1) = {1,2};
Point(3) = { dx/2, 0, 0, cellSize}; Line(2) = {2,3};
Point(4) = { L/2, 0, 0, cellSize}; Line(3) = {3,4};
Point(5) = { L/2, W, 0, cellSize}; Line(4) = {4,5};
Point(6) = { dx/2, W, 0, cellSize}; Line(5) = {5,6};
Point(7) = {-dx/2, W, 0, cellSize}; Line(6) = {6,7};
Point(8) = { -L/2, W, 0, cellSize}; Line(7) = {7,8};
Line(8) = {8,1};
Line(10) = {2,7};
Line(11) = {6,3};
Line Loop(1) = {1,10,7,8};
Line Loop(2) = {3,4,5,11};
Line Loop(3) = {2,-11,6,-10};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Physical Surface('Region 1') = {1};
Physical Surface('Region 2') = {2};
Physical Surface('thermal contact') = {3};
Physical Line('Boundary 1') = {8};
Physical Line('Boundary 2') = {4};
Physical Line('No transfer bottom') = {1,2,3};
Physical Line('No transfer top') = {5,6,7};
Physical Line('common 1') = {10};
Physical Line('common 2') = {11};
""" % locals()
# produce the mesh
mesh = fp.Gmsh2D(geo)
# physics parameters
# we want a thermal contact with a value of hc = 1E4 :
# so the equivalent thermal conduction in the region of the
# thermal contact which is dx wide should be Kc = hc * dx
hc = 1E3 # [W/(m^2.K)] thermal contact conductance
Kc = hc * dx # [W/(m.K) thermal conduction of the material representing
the contact
K1 = 200. # [W/(m.K) thermal conduction in region 1
K2 = 200. # [W/(m.K) thermal conduction in region 2
# define the thermal conduction over the mesh
K = fp.CellVariable(name='Conductance', mesh = mesh, value = 0.)
K.setValue(K1, where = mesh.physicalCells['Region 1'])
K.setValue(K2, where = mesh.physicalCells['Region 2'])
K.setValue(Kc, where = mesh.physicalCells['thermal contact'])
# define the variable T and boundary conditions
T = fp.CellVariable(name='Temperature', mesh = mesh, value = 0.)
# fixed temperature for boundary 1
T.constrain(0., where = mesh.physicalFaces['Boundary 1'])
# 1 kW/m^2 on boundary 2
Ps = 1E3
T.getFaceGrad().constrain(Ps/K2 * mesh._faceNormals,
where=mesh.physicalFaces['Boundary 2'])
# no heat exchange with the sides
T.getFaceGrad().constrain(0. * mesh._faceNormals,
where=mesh.physicalFaces['No transfer top'])
T.getFaceGrad().constrain(0. * mesh._faceNormals,
where=mesh.physicalFaces['No transfer bottom'])
# solve the heat diffusion equation
fp.DiffusionTerm(coeff=K).solve(var=T)
# plot the results
X, Y = mesh.faceCenters()
Ts = T.arithmeticFaceValue()
pl.figure(1)
faces = mesh.physicalFaces['No transfer bottom']
for x, t in [(x,t) for x, t, face in zip(X, Ts, faces) if face ] :
pl.plot(x, t, 'r.')
pl.xlabel('length along the bottom side [m]')
pl.ylabel('temperature [C]')
pl.figure(2)
faces = mesh.physicalFaces['common 1']
n1, avg1 = 0, 0.
for y, t in [(y,t) for y, t, face in zip(Y, Ts, faces) if face ] :
pl.plot(y, t, 'r.')
n1 += 1
avg1 += t
avg1 /= n1
pl.axhline(avg1, color='r', linestyle=':')
faces = mesh.physicalFaces['common 2']
n2, avg2 = 0, 0.
for y, t in [(y,t) for y, t, face in zip(Y, Ts, faces) if face ] :
pl.plot(y, t, 'b.')
n2 += 1
avg2 += t
avg2 /= n2
pl.axhline(avg2, color='b', linestyle=':')
pl.xlabel('length along the contact side 1 (red) and 2 (blue) [m]')
pl.ylabel('temperature [C]')
print '+------------+------------+-------------------------+'
print '| | | delta T |'
print '| cellSize/L | dx/L | |'
print '| | | FiPy | Theory |'
print '+------------+------------+------------+------------+'
print ('| %10.5f | %10.5f | %10.5f | %10.5f |'%
(cellSize/L, dx/L, avg2 - avg1, Ps/hc))
print '+------------+------------+------------+------------+'
print
pl.show()
Many thanks.
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
