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/***@nist.gov/msg02640.html
#and modify it in view of Guyer as in
#https://www.mail-archive.com/***@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 / k_Sample, where=mesh.facesRight)
#follow examples at
https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html#module-examples.diffusion.mesh1D
eq1=fp.TransientTerm()==fp.DiffusionTerm(coeff=diffCoeff) #implicit
dtExplicitLimit_Substrate = cellSize**2 / (2 * D_Substrate)
dtExplicitLimit_Sample = cellSize**2 / (2 * D_Sample)
dtExplicitLimit_InterfaceCell=cellSize**2/(2*Keff.value.mean()/rho_times_cp_eff)
#play with Keff in ipython; .mean seems close enough to what you want
dt=10*min([dtExplicitLimit_Substrate,dtExplicitLimit_Sample,dtExplicitLimit_InterfaceCell])
#SETTING
#viewer=fp.Viewer(vars=(v,),datamin=0,datamax=10.)
viewer=fp.Viewer(vars=(v,))
TJumpAtInterfaceExpected=HeatFlux*R_TIM_unit_area
TJumpInBarNoInterfaceExpected=HeatFlux*(L_Substrate/k_Substrate+L_Sample/k_Sample)
TJumpInBarExpected=TJumpAtInterfaceExpected+TJumpInBarNoInterfaceExpected
blurb='''
at steady state:
expected T jump at interface is %.2E
expected T jump in bar with no interface is %.2E
expected total T jump in bar is %.2E
hit return to continue
''' %
(TJumpAtInterfaceExpected,TJumpInBarNoInterfaceExpected,TJumpInBarExpected)
#run with --inline if you have weave working
steps = 15000 #SETTING
plotEvery=100 #SETTING
saveDataEvery=200 #SETTING
#TODO need code to stop when at thermal steady state
times=[]
rod_dTs=[]
outFilenameBase='20180717UniformMesh04' #SETTING
outFilename=outFilenameBase+'.csv'
assert (not os.path.exists(outFilename)), '%s already exists on disk;
this script requires that you manually update outFilenameBase'
with open(outFilename,'w') as csvfile:
csvWriter = csv.writer(csvfile, delimiter=' ',quotechar='|',
quoting=csv.QUOTE_MINIMAL)
for step in range(steps):
eq1.solve(var=v,dt=dt)
if step%plotEvery==0:
print 'step is %d of %d\n' % (step,steps)
viewer.plot()
if step%saveDataEvery==0:
#saving the temperature drop across the entire rod vs time:
times.append(step*dt)
rod_dTs.append(v.value[-1]-v.value[0])
csvWriter.writerow([step*dt,v.value[-1]-v.value[0]]) #TODO
am looking at values at cell centers; really want face values
raw_input(blurb)
fig, ax = plt.subplots()
ax.plot(times,rod_dTs)
ax.set(xlabel='time (s)', ylabel='temperature drop across entire rod,
K')
ax.grid()
fig.savefig(outFilenameBase+'.png')
plt.show()
raw_input('hit enter to continue')
#make a copy of this script with a new name so you have a good record
of what produced the result
#copyfile(src,dst)
copyfile('ModifyWheelerInViewOfGuyer20130517TransientVersionImplicit03.py',outFilenameBase+'.py')
#don't forget to run with --inline if you have weave installed and
working
Post by Drew DavidsonPost by Drew DavidsondAP1/(dAP1+dAP2)
and
dAP2/(dAP1+dAP2)
where dAP1 and dAP2 were distances from cell center to cell face for
cells
Post by Drew Davidsonon either side of the interface. If these FiPy expressions are
unavailable,
Post by Drew DavidsonI 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
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]