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 <[email protected]> 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 <[email protected]> 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 > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
