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 ]
