Dear Evan, dear all,
This is my first post here, so first of all I'd like to thank the developers for sharing their great work. I've only been using FiPy for a few weeks, I'm not an expert, but I'd like to share my experience which I think is relevant to this post. As a preliminary step before tackling more complex problems, I've been testing a simple model of thermal diffusion with internal heat generation, both in steady-state and transient, to be benchmarked with analytic solutions or with a commercial FEM software. Like the original poster, I also had difficulties implementing a convection boundary condition, or in general a BC that is a function of the variable solved for. These are my main remarks: 1. It seems to me that when such a BC is implemented using for example phi.faceGrad.constrain(-hh * (phi.faceValue - T0)/k, mesh.facesRight), the value of the variable is not updated properly throughout the calculation. At least in the cases I was testing. The workaround I came up with, after much trial and error, is to place this call within a sweeping loop, which seems to work and yield the expected results. This could be possibly a FiPy issue. 2. Even in presence of a Dirichlet condition on part of the boundary, when a convection boundary is implemented the solution could diverge. I've been using the standard solver, which with pysparse is a PCG solver. This can be dealt with with a proper choice of a coefficient of underrelaxation. 3. In Evan's case there are no source or transient terms and in this case I am not so sure of the proper choice for the coefficient of the diffusion term; as the coefficient is constant, it could be canceled out of the equation which reduces to a Laplace equation. I would personally put either 1.0 or the thermal conductivity k, rather then the thermal diffusivity k/rho.cp. Maybe someone can comment on this. My remark is that with Evan's figures, D ~7e-5 and a tolerance 1e-5, the calculation would stop at the first iteration without reaching convergence. It seems to work with a tolerance 1e-9. 4. With the above changes to the code I get a global temperature range from 30.8C to 31.4C, which is still lower than the reported expected result. I hope this helps. Any comments are welcome! Roberto ________________________________ From: [email protected] [[email protected]] on behalf of Evan Chenelly [[email protected]] Sent: 23 October 2014 00:53 To: [email protected] Subject: Re: Fipy gives incorrect solution when using neumann and robin boundary conditions, but no dirichlet 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. <[email protected]<mailto:[email protected]>> 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 <[email protected]<mailto:[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]<mailto:[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]<mailto:[email protected]> > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected]<mailto:[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 ]
