Steady, nonlinear heat equation with pure neumann BCs

mathematics
2

I believe there are two mistakes in the formulation of the weak problem. First, there is a sign error in the boundary term. If we derive a weak problem via integration by parts, we get the following:

\int_\Omega -\nabla\cdot(k\nabla u)v\,d\mathbf{x} = \int_\Omega fv\,d\mathbf{x}
\Rightarrow \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} ~-~\int_{\partial\Omega}(k\nabla u\cdot\mathbf{n})v\,d\mathbf{s}~=~\int_\Omega f v\,d\mathbf{x}
\Rightarrow \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} ~+~\int_{\partial\Omega}h(u - u_\text{air})v\,d\mathbf{s}~=~\int_\Omega f v\,d\mathbf{x}\text{ ,}

under the assumption that the desired BC is

-k\nabla u\cdot\mathbf{n} = h(u - u_\text{air})\text{ ,}

which differs from the equation in the post by a factor of k on the left. (I have followed the notation/conventions used in the code.) Notice that the boundary term enters with a + sign, unlike the code, where it is subtracted.

Second, with a Robin BC like this, there is no need to introduce an auxiliary mechanism to pin down an arbitrary constant. Unlike the pure Neumann demo, if u satisfies the weak problem, u + C (for some nonzero constant C) does not satisfy it, because the boundary term would change.

Another minor point, comparing the code and equations in the post: If k varies in space, the weak form term \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} corresponds to the strong form term -\nabla\cdot(k\nabla u) (as in the derivation above), not -k\Delta u; the first one is what you want, based on deriving the PDE from a physical conservation law.

2 Likes