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 ]
