Hi,
Thanks for your proposed 1D solution. It bears part of the correct
solution which should look like :
from X=-L/2 the temperature increases to 2.5 deg C at X = 0; then there
is a temperature jump to 3.5 deg C and the temperature increases
linearly to 6 deg C at the X=L/2 boundary.
So I'm a bit puzzled as to why for X from -L/2 to 0 the temperature is
0. : it should increase linearly to 2.5 deg C.
Aside from that the rest of the solution appears just offset by this 2.5
deg C (1 deg C jump at X=0 and a linear increase by 2.5 deg C in region
2)
For your (and other users convenience) I tried to translate your
solution to a 2D formulation. I could not do this in a straight forward
manner as some FiPy methods did not appear to function in the same way
for the 2D problem :
1) I could not make the thermal conductance, K, a Face variable because
if so I could not assign a value for K in a "Physical Surface" of the
mesh for regions 1 and 2
2) but as a consequence of making K a cell variable I could not zero-it
on the "thermal contact" boundary (Physical Line).
Unfortunately the solution it returns is not correct not even the jump
and the linear rise in region 2 come out.
Also I noticed that for this 2D formulation the solution is not
stationary (as in : if I sweep more the solution still varies).
So, in brief : why is the 1D solution wrong in Region 1 and how to
formulate this for a 2D problem ?
Many thanks for any help or insight on these.
"""
test FiPy
"""
import fipy as fp, numpy as np, pylab as pl, os, time
# 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;
// 5 4
// 6 <-------- 5 <-------- 4
// | ^ ^
// | Region 1 | Region 2 |
// 6 | | | 3
// | 7| |
// v | |
// 1 --------> 2 --------> 3
// 1 2
Point(1) = {- L/2, 0, 0, cellSize};
Point(2) = { 0, 0, 0, cellSize}; Line(1) = {1,2};
Point(3) = { L/2, 0, 0, cellSize}; Line(2) = {2,3};
Point(4) = { L/2, W, 0, cellSize}; Line(3) = {3,4};
Point(5) = { 0, W, 0, cellSize}; Line(4) = {4,5};
Point(6) = { -L/2, W, 0, cellSize}; Line(5) = {5,6};
Line(6) = {6,1};
Line(7) = {2,5};
Line Loop(1) = {1,7,5,6};
Line Loop(2) = {2,3,4,-7};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Physical Surface('Region 1') = {1};
Physical Surface('Region 2') = {2};
Physical Line('Boundary 1') = {6};
Physical Line('Boundary 2') = {3};
Physical Line('No transfer bottom') = {1,2};
Physical Line('No transfer top') = {4,5};
Physical Line('thermal contact') = {7};
""" % locals()
# print geo
# produce the mesh
mesh = fp.Gmsh2D(geo)
# physics parameters
# we want a thermal contact with a value of hc = 1E3 :
# 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
# 1 kW/m^2 on boundary 2
Ps = 1E3
print 'delta T = %7.3f (Region 1)' % (Ps/K1*L/2)
print ' + %7.3f (thermal contact)' % (Ps/hc)
print ' + %7.3f (Region 2)' % (Ps/K2*L/2)
print ' = %7.3f (total)' % (Ps/K1*L/2+Ps/hc+Ps/K2*L/2)
# define the thermal conduction over the mesh
# K needs to be a cell variable if one want to set the values for
# Physical Surfaces
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'])
# but then I could not set K to 0. on the 'thermal contact' boundary
# the boundary does not appear in the mesh.physicalFaces.keys()
#
# K.setValue(0., where = mesh.physicalFaces['thermal contact'])
# define the heath transfer coeff. as a variable
hc_coef = fp.FaceVariable(name='thermal contact conductance',
mesh=mesh, value=0., rank=1)
hc_coef.setValue(hc * mesh._faceNormals,
where=mesh.physicalFaces['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
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'])
pl.figure(1)
X, Y = mesh.faceCenters()
faces = mesh.physicalFaces['No transfer bottom']
# equation to solve the heat diffusion problem
eq = fp.DiffusionTerm(coeff=K) +
fp.ImplicitSourceTerm(coeff=hc_coef.divergence)
# solve and plot the intermediate results
for sweep in range(10) :
eq.solve(var=T)
sol = sorted([(x,t) for x, t, face in zip(X,
T.arithmeticFaceValue(), faces)
if face ])
pl.plot([x for x, t in sol], [t for x, t in sol],
'rgbcm'[sweep % 5]+['-','-.',':','--'][(sweep/5)%4],
label = 'sweep %d' % (sweep+1))
# plot the theoretical solution
pl.plot([-L/2,0.,0.,L/2],
[0.,Ps/K1*L/2,Ps/K1*L/2+Ps/hc,Ps/K1*L/2+Ps/hc+Ps/K1*L/2],
'k--', label = 'theory')
pl.xlabel('length along the bottom side [m]')
pl.ylabel('temperature [C]')
pl.legend(loc='best')
pl.show()
On Mon, 2013-05-13 at 12:12 -0400, Jonathan Guyer wrote:
> On May 10, 2013, at 2:11 AM, Frederic Durodie
<frederic.duro...@googlemail.com> wrote:
>
> > 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.
>
>
> I don't know whether the following is the solution you want (it
doesn't look like the solutions from your script at all), but it
satisfies the external boundary conditions and gives exactly 1 K jump at
the internal boundary.
>
> I've solved it on a Grid1D because I didn't feel like ginning up a new
Gmsh, but everything should apply to a Gmsh mesh.
>
> If the flux at the thermal contact is to be -\vec{hc} \Delta T, then
we want to add a term
>
> \nabla\cdot(hc T^+) - \nabla\cdot(hc T^-)
>
> to the diffusion equation and zero out the diffusivity at that
boundary.
>
> Since we don't know the values immediately at the jump condition, we
take the value at the nearest cell centers:
>
> T^+ \nabla\cdot(\vec{hc}) - T^- \nabla\cdot(\vec{hc})
>
> which we can write as an ImplicitSourceTerm and the divergence in the
coefficient takes care of the subtraction.
>
>
>
>
> """
> 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
>
>
> # produce the mesh
>
> mesh = (fp.Grid1D(Lx=L/2., dx=cellSize) + [[-L/2.]]) +
fp.Grid1D(Lx=L/2., dx=cellSize)
>
> X = mesh.faceCenters[0]
>
> # 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.FaceVariable(name='Conductance', mesh = mesh)
> K.setValue(K1, where=X < 0.)
> K.setValue(K2, where=X > 0.)
> K.setValue(0., where=X == 0.)
>
> hc_coeff = fp.FaceVariable(name='thermal contact conductance',
mesh=mesh, value=0., rank=1)
> hc_coeff.setValue(hc * mesh._faceNormals, where=X == 0.)
>
> # 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.facesLeft)
>
> # 1 kW/m^2 on boundary 2
> Ps = 1E3
> T.getFaceGrad().constrain(Ps/K2 * mesh._faceNormals,
> where=mesh.facesRight)
>
> vw = fp.Viewer(vars=T)
>
> # solve the heat diffusion equation
> eq = fp.DiffusionTerm(coeff=K) +
fp.ImplicitSourceTerm(coeff=hc_coeff.divergence)
>
> # ImplicitSourceTerm can only be semi-implicit, so we sweep twice
> for sweep in range(2):
> eq.solve(var=T)
>
> print T([cellSize/2.]) - T([-cellSize/2.])
> print Ps/K2, T.faceGrad[..., mesh.facesRight.value]
>
> # plot the results
>
> fp.Viewer(vars=T).plot()
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
