Drew - I had found the same sign error in the denominator. You're also correct that the units didn't agree. Multiplying by the face areas is redundant to taking the divergence.
I've updated our notes: https://github.com/usnistgov/fipy/blob/dd3420fb71884d74850051ad2280bff525301824/documentation/USAGE.rst and your latest code: https://github.com/guyer/Convection2DFiPyExample01/commit/ca59e8955319545f9965705c8fadfcbb5abd5951 With these changes, the results look much more plausible. - Jon > On Dec 12, 2018, at 1:22 PM, Drew Davidson <davidson...@gmail.com> wrote: > > Hi Dr. Guyer, > > I looked over your more recent documentation: > https://github.com/usnistgov/fipy/blob/70b72b7abb267c85ada47886b8b3573e1819fffc/documentation/USAGE.rst > > I am embarrassed I did not understand how to choose values for a and b in the > Robin condition before. > > I also cloned your repo to my local computer and ran it: > https://github.com/guyer/Convection2DFiPyExample01 > > When I run the code in your repo, the viewer does show a variation in T with > respect to r, but the size of that variation is extremely tiny. Essentially, > the solution is pinned at the initial temperature. The variation also has the > wrong sign (solid object is warming up), as the ambient air is colder than > the solid object (T_infinity < T_initial), and the solid object should be > cooling down. Maybe this is due to confusion I sowed by having var be T minus > T_infinity; sorry. > > I tried another script in my repo: > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/ConvectionTestProblem2D_01_20181211Version.py > > This script is upgraded in that now, the log file contains information about > the solution (max and min temperatures). I obtain max and min temperatures > using var.value.min() and var.value.max() in order to rule out goof-ups in my > get_solution_along_ray function. > > The solution variable var is now exactly the temperature, in order to get rid > of the previous source of confusion in which var was T minus T_infinity. > > The solution in my script is still pinned at the initial temperature. This > persists even if I crank up the convection coefficient to a huge value (keep > multiplying it by 1000 over and over and still seeing solution pinned at > initial condition). > > I made handwritten notes on top of a pdf I made via Kile Editor of your more > recent documentation: > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage20181211_Annot_CTTC.pdf > > It seems to me that the units of the source terms should match the units of > the terms in the partial differential equation for heat conduction. I thought > that is how it worked from the FiPy examples, but I haven’t gone back and > made sure. I find that for the first source term RobinCoeff*g, the units > don’t seem to match the units of the terms in the PDE for heat conduction (I > didn’t check the 2nd source term since it should be the same situation). I > did this in kind of a hurry so maybe I am goofing up. I have not studied the > finite volume method, and am just going on a hunch. > > It also seemed to me that there might be a minus sign missing in the > RobinCoeff. The solution I get is still pinned at initial condition > regardless of a minus sign or no minus sign there. > > Thanks > > On Fri, Dec 7, 2018 at 6:05 PM Guyer, Jonathan E. Dr. (Fed) via fipy > <fipy@nist.gov> wrote: > Drew - > > > Thanks to your notes, I found a couple of errors in our Robin discussion. As > you noted, the units don't work for the conversion of the divergence of gamma > grad phi to the sum over faces (midway on page 1 of your notes). The left > hand side should be integrated over volume. That line is expressing the > discretization of the divergence theorem. > > I also discovered that the .divergence operator in FiPy doesn't work reliably > on scalars. There's no real reason it should, but I thought I had convinced > myself that it did. As a result, everything in the Robin terms needs to be > multiplied by the surface normal > > Those changes to the FiPy documentation are > [here](https://github.com/usnistgov/fipy/pull/615) for now. > > Beyond addition that, there are some changes necessary to put your boundary > condition into Robin form. The point of the Robin condition is that it ties > the gradient of the field to the value of the field. > So, g isn't h(T - T_inf) / (-k); it's just -h T_inf / (-k). > Likewise, \vec{a} isn't zero. Rather, \hat{n}\cdot\vec{a} = h/(-k). > > Here are the changes I made to your script to reflect these changes: > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/compare/master...guyer:master > > I think I got your normalized temperature accounted for properly, and put the > thermal diffusivity and convection coefficient in the proper places, but it's > worth checking through. > > With these changes, the temperature field is now rotationally symmetric (it > wasn't before, which is why I had to multiply by the normal before taking the > divergence). > > Heat seems to be fluxing in from the outside, so I probably have the sign > wrong. > > I don't have Octave, so I have no idea how this compares with your analytical > solution. > > - Jon > > > On Dec 4, 2018, at 2:34 PM, Drew Davidson <davidson...@gmail.com> wrote: > > > > Hi Dr. Guyer, > > > > First I tried getting rid of the square brackets in > > ConvectionTestProblem2D_01.py (commit ‘changed square brackets to > > parenthesis in convection BC, but get same…’), but results are the same > > (still wrong). > > > > Next: > > As you directed, I took a look at: > > > > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > > > Firefox was showing me raw latex rather than human-readable equations at > > that web address, so I made a version I could read using the kile editor: > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage.pdf > > > > That material is rather challenging for me. I took a stab at it, resulting > > in: > > commit: ‘made a first attempt at a formulation in view of suggestions by > > fipy …’ > > > > ConvectionTestProblem2D_01_2ndTry.py > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/fb47960548650c994eb0c6f990e0db297566e174/ConvectionTestProblem2D_01_2ndTry.py > > > > a few handwritten notes indicating what I was thinking: > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsageConvectionBC.pdf > > > > Here is a bit of the code: > > > > #ref: > > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > > > #warning: avoid confusion between convectionCoeff in that fipy > > documentation, which refers to terms involving "a", and convectionCoeff > > here, which refers to a heat transfer convection coefficient at a boundary > > > > Gamma0=D_thermal > > > > Gamma = FaceVariable(mesh=mesh, value=Gamma0) > > > > mask=surfaceFaces > > > > Gamma.setValue(0., where=mask) > > > > dPf = FaceVariable(mesh=mesh, value=mesh._faceToCellDistanceRatio * > > mesh.cellDistanceVectors) > > > > Af = FaceVariable(mesh=mesh, value=mesh._faceAreas) > > > > #RobinCoeff = (mask * Gamma0 * Af / (dPf.dot(a) + b)).divergence #a is zero > > in our case > > > > b=1. > > > > RobinCoeff = (mask * Gamma0 * Af / b).divergence #a is zero in our case > > > > #eq = (TransientTerm() == DiffusionTerm(coeff=Gamma) + RobinCoeff * g - > > ImplicitSourceTerm(coeff=RobinCoeff * mesh.faceNormals.dot(a))) #a is zero > > in our case > > > > # g in this formulation is -convectionCoeff/k*var, where var=T-T_infinity > > > > eq = (TransientTerm() == DiffusionTerm(coeff=Gamma) + > > ImplicitSourceTerm(RobinCoeff * -convectionCoeff/k)) > > > > > > Now the solution variable remains stuck at the initial condition, as if the > > boundary condition is not being applied. I am sort of out of my depth at > > this point. I was guessing that an ImplicitSourceTerm was the right thing > > to do, since 'g' in the Robin condition depends on the solution variable. > > I did mess around a little in IPython seeing if any terms are coming out > > all zeros. > > > > Again, I put everything at > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01. > > > > Thanks > > > > On Mon, Dec 3, 2018 at 4:21 PM Guyer, Jonathan E. Dr. (Fed) via fipy > > <fipy@nist.gov> wrote: > > Drew - > > > > Apologies for the delayed reply. > > > > There are a couple of issues here: > > > > `-convectionCoeff/k*(var.faceValue-T_infinity)` describes a FiPy > > FaceVariable. > > > > `[-convectionCoeff/k*(var.faceValue-T_infinity)]` is a single element > > Python list that contains a FiPy FaceVariable. > > > > Multiplying that list by other elements has rather unpredictable results. > > > > In short, get rid of the square brackets. We do our best to treat Python > > lists like FiPy vector fields, but there's only so much we can do. A list > > that holds a numpy array is not a vector field. > > > > > > Beyond that, what you've described looks like a Robin boundary condition to > > me. Our best recommendation for Robin conditions is covered at: > > > > > > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > > > Please don't hesitate to ask for more clarification if this doesn't get you > > where you need. > > > > - Jon > > > > > > > > > On Nov 16, 2018, at 11:55 PM, Drew Davidson <davidson...@gmail.com> wrote: > > > > > > Hello, > > > > > > I am stuck in how to correctly apply a simple convection boundary > > > condition in FiPy, in the context of simple transient heat conduction in > > > 2D. > > > > > > Is this correct: > > > var.faceGrad.constrain([-convectionCoeff/k*(var.faceValue-T_infinity)]*mesh.faceNormals,where=surfaceFaces) > > > > > > I have a 2D mesh generated in gmsh. The convection boundary condition is > > > applied to a curved boundary. > > > > > > The code and results are found at: > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/tree/master > > > > > > The project appears to have a heap of files, but It’s a really just a > > > simple 2D problem with a comparison to an analytical solution. The basic > > > script is ConvectionTestProblem2D_01.py. Report20181116.pdf shows current > > > results, which are clearly wrong. > > > > > > Am I correctly applying the convection boundary condition in terms of the > > > FiPy syntax/language for this 2D problem with a gmsh mesh and a curved > > > boundary? > > > > > > Thanks > > > _______________________________________________ > > > 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 ] _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
