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 / 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 On Mon, Jul 23, 2018 at 11:30 AM Daniel Wheeler <daniel.wheel...@gmail.com> wrote: > On Tue, Jul 17, 2018 at 4:43 PM, Drew Davidson <davidson...@gmail.com> > 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 ] >
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
