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 ]
