Hello fipy listserv users,
I am trying to model a rectangular 3d heat spreader using FiPy. The
problem has a constant heat flux on one side of the domain (neumann
boundary on part of facesFront) and a heat transfer coefficient/biot number
condition on the other side (robin boundary on all of facesBack.) The
equation is simply the steady state heat equation (diffusion equation.) I
tried to use examples I found here and in the examples folder to formulate
this the best I could.
If I use a dirichlet BC on either side I get expected results, but when I
use only neumann and robin conditions the solution seems to blow up.
One thought (that I am having right now) is that I can try to rewrite the
robin BC so that it constrains the variable in terms of its gradient,
instead of constraining the gradient in terms of the variable as it does
now. But I don't know if that is the reason why this isn't working
correctly. There is also an analytic solution to this problem but in a
cylindrical domain, that I could use to see if I am close to the right
answer.
My code is below. You can switch either side to dirichlet conditions by
uncommenting the other two lines where the boundary conditions are set.
Also everything below is in SI units but I can probably non-dimensionalize
the problem or reduce it to a 1D problem if that makes it easier to
troubleshoot as well. The 2D mesh and variable are used only to store the
results for visualization.
Thank you for reading this!
Evan
from fipy import *
Lx=0.060
Ly=0.030
Lz=0.010
nx = 30
ny = 15
nz = 10
dx = Lx/nx
dy = Ly/ny
dz = Lz/nz
L = (Lx+Ly)/2
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny, Lx=Lx, Ly=Ly)
mesh2 = Grid3D(dx=dx, dy=dy, dz=dz, nx=nx, ny=ny, nz=nz, Lx=Lx, Ly=Ly,
Lz=Lz)
#Aluminum Properties
#Conductivity
k=180
#Density
dens=2.7
#Specific Heat
cp=963000
#Problem Parameters
q=5
hTopSurface=400
tinf=25
#Only used with dirichlet conditions
tTopSurface=45
tBottomPlate=30
#Source XY Bounds
Xlbound=Lx/4
Xubound=3*Lx/4
Ylbound=Ly/4
Yubound=3*Ly/4
#Source Area
qa = (Xubound-Xlbound)*(Yubound-Ylbound)
#1D Temperature Rise to try initializing the problem to the right ballpark
dt1d=(q*Lx/(k*qa))
phi = CellVariable(name = "Temperature (DegC)",
mesh = mesh2,
value = tinf+(dt1d/2))
phi2D = CellVariable(name = "Temperature (DegC)",
mesh = mesh,
value = tinf+(dt1d/2))
#diffusivity k/(density*specific heat) - Aluminum 6061
D = (k)/(dens*cp)
#Total Heat Flux
fluxBottomPlate=-q/(k*qa)
X, Y, Z = mesh2.faceCenters
fluxRegion = (mesh2.facesFront & (Y > Ylbound) & (Y < Yubound) & (X >
Xlbound) & (X < Xubound))
topSurface = mesh2.facesBack
phi.faceGrad.constrain(-hTopSurface*(phi.faceValue()-tinf)/k, topSurface)
#phi.constrain(tTopSurface,topSurface)
phi.faceGrad.constrain(fluxBottomPlate, fluxRegion)
#phi.constrain(tBottomPlate, fluxRegion)
eq = DiffusionTerm(coeff=D)
if __name__ == '__main__':
restol = 1e-5
else:
restol = 0.05
res = 1e+10
times = 0
while res > restol:
res = eq.sweep(var=phi)
times = times +1
if __name__ == '__main__':
viewer = Viewer(vars=phi2D)#, datamin=0., datamax=1.)
viewer.plot()
cutPlaneZValue=0.005
x,y = phi2D.getMesh().getCellCenters()
phi2D.setValue(phi((x, y, cutPlaneZValue * numpy.ones(len(x))),
order=1))
viewer.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
