Re: thermal contact between two regions
On Mon, Jul 23, 2018 at 1:06 PM, Drew Davidson wrote: > > 'dx' is actually the thickness of the thermal contact region. It is not the > volume of the cell. If it were a volume, the ratio dx/mesh._cellDistances > would not be dimensionless. I was seeking a thermal contact region of zero > thickness (dx=0), giving a true thermal contact resistance. Sorry, my mistake on that one. > I let this thread drop because the Salazar paper I referenced gave me > expressions for effective thermal conductivity and effective thermal > diffusivity, and this allowed me to arrive at a reasonable-looking FiPy > result for transient thermal behavior. My result still deviates a bit from > the exact solution at steady state, and this deviation persists despite mesh > refinement, but I am busy with other things and have called it good enough > for now. Thanks for following up and sharing your code. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
Hello Dr. Wheeler,
'dx' is actually the thickness of the thermal contact region. It is not
the volume of the cell. If it were a volume, the ratio
dx/mesh._cellDistances would not be dimensionless. I was seeking a thermal
contact region of zero thickness (dx=0), giving a true thermal contact
resistance.
I let this thread drop because the Salazar paper I referenced gave me
expressions for effective thermal conductivity and effective thermal
diffusivity, and this allowed me to arrive at a reasonable-looking FiPy
result for transient thermal behavior. My result still deviates a bit from
the exact solution at steady state, and this deviation persists despite
mesh refinement, but I am busy with other things and have called it good
enough for now.
OK, I won't leave you hanging. Here is the code I was working on. Maybe
you or someone else can spot a defect that would explain why the code does
not exactly replicate the steady state solution. Another thing to do is:
compare the transient result with the transient result obtained another way
(ANSYS etc). I used to have an analytical solution for the transient
laying around, but threw it out in a pile of papers (oops).
#I am trying to look at the script by Wheeler at
#https://www.mail-archive.com/fipy@nist.gov/msg02640.html
#and modify it in view of Guyer as in
#https://www.mail-archive.com/fipy@nist.gov/msg02641.html
#now instead of a steady state problem, do a transient problem
import fipy as fp
import numpy as np
from IPython import embed
import sys
import os
import pygmsh
import meshio
import csv
import matplotlib.pyplot as plt
from shutil import copyfile
#head
def
get_uniform_mesh_with_pygmsh(L_Sample,L_Substrate,cellSize,write_to_msh_file=False):
'''
figure it out in ipython
ref: https://github.com/nschloe/pygmsh/blob/master/test/test_pacman.py
'''
geom = pygmsh.built_in.Geometry()
#geom.add.point? see what the input arguments are in ipython using ?
p1=geom.add_point([0,0,0],cellSize)
p2=geom.add_point([L_Substrate,0,0],cellSize)
p3=geom.add_point([L_Substrate+L_Sample,0,0],cellSize)
l1=geom.add_line(p1,p2)
l2=geom.add_line(p2,p3)
#not much point in doing physical groups since have no way to transmit
this info to 1D fipy mesh
# geom.add_physical_point(p1,label='SubstrateEnd')
# geom.add_physical_point(p2,label='TCR')
# geom.add_physical_point(p3,label='SampleEnd')
# geom.add_physical_line(l1,label='Substrate')
# geom.add_physical_line(l2,label='Sample')
points, cells, point_data, cell_data, field_data =
pygmsh.generate_mesh(geom)
if write_to_msh_file:
mesh=meshio.Mesh(points=points,cells=cells,point_data=point_data,cell_data=cell_data,field_data=field_data)
meshio.write('TCR01.msh',mesh)
mesh=get_fipy_1D_mesh_from_points(points)
return mesh
#head
#head
if __name__=='__main__':
#substrate is copper and sample is pure tin
rho_Sample=6980. #use wolfram alpha to get the density of copper;
kg/m^3
cp_Sample=227. #use wolfram alpha to get the specific heat of copper;
J/kg/K
k_Sample=59.6 #W/m/K
D_Sample=k_Sample/rho_Sample/cp_Sample
L_Sample=2e-6 #m SETTING
rho_Substrate=8960. #use wolfram alpha to get the density of copper;
kg/m^3
cp_Substrate=384.4 #use wolfram alpha to get the specific heat of
copper; J/kg/K
k_Substrate=381.8 #W/m/K
D_Substrate=k_Substrate/rho_Substrate/cp_Substrate
L_Substrate=10*L_Sample #m SETTING
#model an interface thermal resistance between copper and solder, with
the solder being soldered to the copper
R_TIM_unit_area=2e-6 #K/(W/m2) SETTING
TIM_target_dT=1. #K SETTING
HeatFlux=TIM_target_dT/R_TIM_unit_area #W/m^2 heat flux is determined
by TIM_target_dT
#using a uniform mesh:
cellSize=.1*L_Sample #SETTING
mesh=get_uniform_mesh_with_pygmsh(L_Sample,L_Substrate,cellSize,write_to_msh_file=False)
#make sure you like the mesh before doing a run
#sys.exit()
diffCoeff = fp.FaceVariable(mesh=mesh, value=D_Substrate)
diffCoeff.setValue(D_Sample,where=(mesh.faceCenters[0]>=L_Substrate))
xCenters,=mesh.cellCenters
#thermal conductivity of a cell if it's R" were R_TIM_unit_area:
Kcontact=1./R_TIM_unit_area*mesh._cellDistances #type is
numpy.ndarray; [W/m/K]
Kcell=fp.CellVariable(mesh=mesh,value=k_Substrate)
Kcell.setValue(k_Sample,where=(xCenters>=L_Substrate))
Kavg = Kcell.harmonicFaceValue #a FaceVariable
#see 20180717EffectiveProperties.xoj
Keff=Kcontact * Kavg / (Kavg + Kcontact) #W/m/K
rho_times_cp_eff=.5*(rho_Substrate*cp_Substrate+rho_Sample+cp_Sample)
#TODO ask on fipy mailing list as to how to express dAP1/(dAP1+dAP2) in
fipy syntax; the factor of .5 here is only approximate
diffCoeff.setValue(Keff/rho_times_cp_eff, where=(mesh.faceCenters[0] ==
L_Substrate))
v = fp.CellVariable(mesh=mesh)
v.constrain(0, where=mesh.facesLeft)
v.faceGrad.constrain(HeatFlux /
Re: thermal contact between two regions
On Tue, Jul 17, 2018 at 4:43 PM, Drew Davidson wrote: > > I still need FiPy code for: > > > dAP1/(dAP1+dAP2) > > > and > > > dAP2/(dAP1+dAP2) > > > where dAP1 and dAP2 were distances from cell center to cell face for cells > on either side of the interface. If these FiPy expressions are unavailable, > I would think assuming .5 is going to be OK... I think it's outlined in the link that you gave. Kcontact = hc * mesh._cellDistances Kavg = Kcell.harmonicFaceValue K.setValue(Kcontact * Kavg / (Kavg + Kcontact * (1 - dx / mesh._cellDistances)), where=mesh.physicalFaces['thermal contact']) m._cellDistances is the distance between each cell and `dx` is the volume of the cell for a 1D problem. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
Hello, I continued to stare at https://www.mail-archive.com/fipy@nist.gov/msg02641.html. I realized Keff must be effective thermal conductivity as described in: Salazar, Agust n. “On Thermal Diffusivity.” *European Journal of Physics* 24, no. 4 (July 1, 2003): 351–58. https://doi.org/10.1088/0143-0807/24/4/353 . For years I sat right next to grad students whose entire degree was in computing Keff of composite materials via ANSYS/COMSOL, but I still looked at https://www.mail-archive.com/fipy@nist.gov/msg02641.html with incomprehension. I realized all that was being done in https://www.mail-archive.com/fipy@nist.gov/msg02641.html must be to compute Keff of the finite volume which contains the thermal contact resistance, and then to assign thermal conductivity of the finite volume face at the interface to Keff. The Salazar paper also gives equations for effective thermal diffusivity, which is what I was asking for in the current thread. Now I think I mostly have what I need. I still need FiPy code for: dAP1/(dAP1+dAP2) and dAP2/(dAP1+dAP2) where dAP1 and dAP2 were distances from cell center to cell face for cells on either side of the interface. If these FiPy expressions are unavailable, I would think assuming .5 is going to be OK... Thanks On Thu, Jul 12, 2018 at 2:28 PM Daniel Wheeler wrote: > I think you can do this by aligning the mesh with the contact region > and specifying the diffusion coefficient at the contact faces. > > On Thu, Jul 12, 2018 at 12:28 PM, Drew Davidson > wrote: > > Hello, > > > > > > I’m sorry for not being more specific. I was trying to incorporate > > https://www.mail-archive.com/fipy@nist.gov/msg02626.html by reference. > > -- > Daniel Wheeler > > ___ > 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 ]
Re: thermal contact between two regions
I think you can do this by aligning the mesh with the contact region and specifying the diffusion coefficient at the contact faces. On Thu, Jul 12, 2018 at 12:28 PM, Drew Davidson wrote: > Hello, > > > I’m sorry for not being more specific. I was trying to incorporate > https://www.mail-archive.com/fipy@nist.gov/msg02626.html by reference. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: thermal contact between two regions
Hello, I’m sorry for not being more specific. I was trying to incorporate https://www.mail-archive.com/fipy@nist.gov/msg02626.html by reference. Sketch of problem: (1) Region 1 --* thermal contact resistance * -- Region 2 (2) Temperature at (1) is zero. Heat flux at (2) is constant. At time zero, temperature is everywhere zero. ** What is a thermal contact resistance? https://en.wikipedia.org/wiki/Thermal_contact_conductance A dissertation available free online: Thermal contact resistance Author: Mikic, B. B.; Rohsenow, Warren M. Citable URI: http://hdl.handle.net/1721.1/61449 It’s equation (1.1) on page 12 of this dissertation. Thermal contact resistance is directly analogous to an electrical resistor in a simple circuit: (temperature drop across thermal contact resistance) = (heat flux) * (its unit-area thermal resistance) Units: Kelvin = W/m^2 times Kelvin/(W/m^2) Mathematically, the thermal contact resistance has no thickness and no heat capacity (it does not store energy). ** Please look at https://www.mail-archive.com/fipy@nist.gov/msg02641.html. Dr. Guyer and Dr. Wheeler considered the case of steady heat flow (steady state problem with no transient term), and their thermal contact has a finite thickness ‘dx’. They derived an expression for effective thermal conductivity at the thermal contact location, and provided FiPy code. How do things change in a transient problem? What is the effective thermal diffusivity (thermal diffusivity = thermal conductivity/(mass density * specific heat)) at the interface location, given that Region 1 and Region 2 have different thermal conductivities, mass densities, and specific heats? How is it written in FiPy? This is the point at which my knowledge of the finite volume method and FiPy is too weak. A solution where the thermal contact has a finite thickness ‘dx’ is also of interest, since it seems like a thermal contact resistance can be approximated by setting ‘dx’ very small. Trying to get on the same page: thermal diffusivity m^2/s thermal conductivity W/m/K mass density kg/m^3 specific heat J/kg/K Thanks On Thu, Jul 12, 2018 at 9:49 AM Daniel Wheeler wrote: > On Wed, Jul 11, 2018 at 7:00 PM, Drew Davidson > wrote: > > > > 1. Transient heat conduction rather than steady state. Same boundary > > conditions but initial condition can be zero temperature eveywhere. > > Add a TransientTerm to the equation and step through the solution with > time steps and sweeps if necessary. > > > 2. The two regions have differing properties (thermal conductivity, mass > > density, specific heat) > > The coefficients can have spatially varying values in FiPy. Define the > coefficients as face or cell variables. > > > 3. The thermal contact is a true thermal contact resistance (a pure > > resistance); it has zero thickness, and it has zero heat capacity. > > I need to see it defined mathematically to understand how to implement > it in FiPy, but I'm confident it's possible to define it. > > -- > Daniel Wheeler > ___ > 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 ]
Re: thermal contact between two regions
On Wed, Jul 11, 2018 at 7:00 PM, Drew Davidson wrote: > > 1. Transient heat conduction rather than steady state. Same boundary > conditions but initial condition can be zero temperature eveywhere. Add a TransientTerm to the equation and step through the solution with time steps and sweeps if necessary. > 2. The two regions have differing properties (thermal conductivity, mass > density, specific heat) The coefficients can have spatially varying values in FiPy. Define the coefficients as face or cell variables. > 3. The thermal contact is a true thermal contact resistance (a pure > resistance); it has zero thickness, and it has zero heat capacity. I need to see it defined mathematically to understand how to implement it in FiPy, but I'm confident it's possible to define it. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to implement a thermal contact between two regions ?
On May 16, 2013, at 1:23 PM, Daniel Wheeler daniel.wheel...@gmail.com wrote: I believe the correct diffusion coefficient across the face that approximates a thin layer is D = Kc * (1 / (epsilon + (1 - epsilon) * D_ratio)) where D_ratio = Kc / 2 * (1 / K1 + 1 / K2) and epsilon = dx / cellSize Daniel and I talked about this and it's not quite right. I rederived it, considering the possibility of different cell sizes on either side of the boundary (which doesn't change much in the long run). What I get is that: Keff = Kcontact * Kavg / (Kavg + Kcontact * (1 - dx / cellSize)) where Kcontact = hc * cellSize and Kavg = K1 * K2 * cellSize / (K1 * dAP2 + K2 * dAP1) is the harmonic mean of K (dAP1 and dAP2 are the distances from the respective cell centers to the face center. In FiPy, I would write Kcontact = hc * mesh._cellDistances Kavg = Kcell.harmonicFaceValue K.setValue(Kcontact * Kavg / (Kavg + Kcontact * (1 - dx / mesh._cellDistances)), where=mesh.physicalFaces['thermal contact']) For most values of dx, this reduces to hc * mesh._cellDistances, although if hc * mesh._cellDistances is much larger than Kavg, then Keff starts to look like Kavg. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to implement a thermal contact between two regions ?
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.,
Re: how to implement a thermal contact between two regions ?
On May 16, 2013, at 6:39 AM, Frederic Durodie frederic.duro...@googlemail.com wrote: 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. I suspect the problem is that divergence of the flux at the cell just left of the thermal contact is satisfactorily zero, because we are not treating the thermal contact as a flux condition but as a source. The flux condition could be achieved with a ConvectionTerm, but I don't see any way to compose that that retains the jump. A ConvectionTerm applies to the value *at* the face, and you have two values at the face, one on either side. I'll have to think about this. 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). You can resolve this issue by assigning K in two steps: Kcell = fp.CellVariable(name='Conductance', mesh = mesh, value = 0.) Kcell.setValue(K1, where = mesh.physicalCells['Region 1']) Kcell.setValue(K2, where = mesh.physicalCells['Region 2']) K = fp.FaceVariable(name='Conductance', mesh = mesh, value=Kcell.faceValue) K.setValue(0., where=mesh.physicalFaces['thermal contact']) This restores the equivalent of my 1D solution, with a 1 K jump at the interface and zero temperature in Region 1. The solution is also not sensitive to sweeping any longer. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to implement a thermal contact between two regions ?
On Fri, May 10, 2013 at 2:11 AM, Frederic Durodie frederic.duro...@googlemail.com wrote: 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 believe the correct diffusion coefficient across the face that approximates a thin layer is D = Kc * (1 / (epsilon + (1 - epsilon) * D_ratio)) where D_ratio = Kc / 2 * (1 / K1 + 1 / K2) and epsilon = dx / cellSize I implemented this for a 1D problem and it seems to give the correct jump (see http://pastebin.com/qB0XgceL). I set K1 and K2 to just be K in the script to simplify. Cheers -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to implement a thermal contact between two regions ?
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 ]
how to implement a thermal contact between two regions ?
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.0 |
+++++
|0.00500 |0.01000 |0.63490 |1.0 |
+++++
|0.00500 |0.00500 |1.06939 |1.0 |
+++++
|0.01000 |0.00500 |1.37685 |1.0 |
+++++
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
