First thank you very much for answering my question! And for working on
Fipy.
If you can, please help me understanding is wrong. My understanding is
that the reference temperature on the Robin boundary condition (tinf)
should serve to "pin down" the solution to a unique value just like a
dirichlet condition would. The 1-D problem here would just simplify to a
"thermal resistance" inside the block of Lz/(k*A) and a convective
resistance outside of 1/(h*A). The maximum temperature would be
q*(Lz/(k*A)+1/(h*A))+tinf. But the minimum temperature in the conducting
domain is still q/(h*A)+tinf. So the minimum temperature (phi) in the
domain is higher than the reference by some amount set by both h and q. I
am thinking of it being "pinned somewhere outside the domain" but I'm not
sure that's a correct understanding.
Actually the number calculated by the other tool I used is close (order of
magnitude) to the 1D number I would get using the above formula.
Trying to get a unique solution is why I tried to specify the robin
condition by constraining phi in my second post. To see if that would fix
the problem you were talking about (I am only constraining faceGrad in my
code) But that didn't seem to work.
To be concise, I will try to write the BCs I am trying to use, though
you'll have to excuse me. I am new to tex:
1. Heat Source q (by "flux" I am just talking about the below equation, a
constant neumann condition, heat per area)
\frac{\delta \phi }{\delta z}\mid_{fluxRegion} = \frac{-q}{kA_{fluxRegion}}
2. Heat Transfer Coefficient h (this is the robin condition, and \phi_{inf}
is a constant)
\frac{\delta \phi }{\delta z}\mid_{topSurface} = \frac{-h}{k}*\left (
\phi_{topSurface} - \phi_{inf}\right )
I also confirmed that the dirichlet cases are giving me correct answers,
which is a reassurance that i'm doing something right. But I can't get the
case with the neumann & robin condition to solve.
I also tried to use another method you wrote about here a few months ago,
instead of setting faceGrad, using a divergence term at the boundary. That
didn't work either. But i'm not sure I did it right. I'm still trying.
Here's what I used:
eq = DiffusionTerm(coeff=D) + (-hTopSurface*(phi.faceValue-tinf) *
mesh2.faceNormals * mesh2.facesBack).divergence
Thanks for the help!
Evan
On Wed, Oct 22, 2014 at 2:51 PM, Guyer, Jonathan E. Dr. <
jonathan.gu...@nist.gov> wrote:
> The diffusion equation, steady-state or otherwise, admits an infinite
> number of solutions in the absence of Dirichlet conditions somewhere on the
> boundary. That's not a numerical issue; it's a property of the equation.
> DiffusionTerm(coeff=D) dictates the curvature, your boundary conditions
> dictate the slopes, but there's absolutely nothing to determine the
> absolute value. The other code may have different default boundary
> conditions. FiPy is no-flux on all boundaries not otherwise specified. I
> would expect other FV codes to do the same, but maybe not. The other code
> may be trying to help you by normalizing the solution in some fashion;
> that's not wrong, but it's not right, either.
>
> Setting aside the absolute value, I see a range in temperature of about
> 0.4, so around half what you get from the other code. This may be because
> you appear to constrain the gradient to be equal to the flux:
>
> phi.faceGrad.constrain(fluxBottomPlate, fluxRegion)
>
> The flux is D \nabla \phi, whereas phi.faceGrad is \nabla \phi.
>
> On Oct 22, 2014, at 1:04 PM, Evan Chenelly <echene...@gmail.com> wrote:
>
> > If anyone trying to read this is wondering what the expected output is,
> I have run the same case in a commercial FV code and the temperature
> profile through that cut plane given the same inputs (q=5, htc=400,
> tinf=25, 30x60x10mm domain, same source size) should go approximately from
> 31.8 to 32.5.
> >
> > I also tried setting the temperature at the boundary in terms of the
> gradient, but that also gave incorrect results This time the temp
> difference was too small instead of too large. Here is the line I used to
> try that:
> >
> > phi.constrain( (-k*phi.faceGrad.numericValue[2]/hTopSurface) + tinf ,
> topSurface)
> >
> > I'm not sure that I have the right index on faceGrad there as well. I
> assume the first index indicates the direction of the gradient? In which
> case I am trying to use the Z or Front-Back direction.
> >
> > I've tried sweeping/iterating the result too but it doesn't seem to
> help. Has anyone tried to use BCs like this before?
> >
> > Thank you!
> >
> > Evan
> >
> > On Tue, Oct 21, 2014 at 2:54 PM, Evan Chenelly <echene...@gmail.com>
> wrote:
> > 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
> > 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 ]
>
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
