Jonathan, Thanks for the reply.
Here is the code I ended up with that recreates the example (what appears to be) accurately: http://pastebin.com/f20733965 At the end of the iterations (~17 hours of "virtual" time) the problem is more or less at steady state. Given the size of the chunk of metal, and comparatively low wattages involved, this seems like an appropriate number. In determining where I needed to fix my units, I think i found an error in the Fipy user manual, with regards to using a Robin condition to simulate a Neumann condition when there is no convection term. The manual states that the user should divide through by the diffusion coefficient, when I think it should actually be multiply through by the diffusion coefficient. Consider the standard "heat convection" boundary condition: k * (dT/dx) = h (T - Tinf) The gradient value is then: dT/dx = g = (h/k) (T - Tinf) with the Robin condition defined as (more or less): R * (dU/dx) + a * U = y and given that a=0 (no convection): R * (dU/dx) = y so: dU/dx = y/R equating dU/dx = dT/dx: y/R = g or y = g*R Q.E.D Concerning the rest of your response, I have some in-line questions/comments below. <-- snip --> > Increasing the time step should be OK, subject to the standard > concerns about accuracy. For the steady-state solution, I think you're > OK, though, as long as the problem converges at all. Since you have > convection at the boundaries, I'm not sure you're free change the > diffusion coefficient without altering the steady-state temperature > profile. With my new, properly formed, problem, if I drop the transient term and do an eqn.solve(T, boundaryConditions=BCs), I do not get the same solution. Is there a different way to solve for the steady solution directly? <-- snip --> > I would have to see what you tried. This works (modulo any adjustments > for density and heat capacity): > > mask = mesh.getFacesBottom() > > outCoeff = (mask * [[0], [-h]]).getDivergence() What does this line ^^^ do? particularly the mask*[[0], [-h]] part? That's really throwing me for loop. > convectionOut = ImplicitSourceTerm(coeff=outCoeff) - outCoeff * > Tinf I will try to reformulate my boundary condition to what you've shown, and see if I get anywhere. Thanks! Cheers, Chris
