How to Solve the Heat Equation with Neumann Boundary Conditions?

[johan_ekh]

Hi all,
I have tried to solve the heat equation with a Fourier-Bessel approach but I fail to implement the boundary condition, which is a Neumann condition. Every textbook that I have available treats the corresponding Dirichlet problem but not the Neumann one. Below I have tried to summarize the problem and my findings so far in the hope that someone could show me how to solve the problem with Neumann boundary conditions. Please be aware that my skills in solving PDEs with analytical methods are limited.

The problem is to calculate the temperature in a unit 2D disc as a function of time, starting with a known temperature distribution at time 0 (t=0) and with heat exiting through the boundary of the circle. Thus,

$$\frac{\partial u}{\partial t} = a \left[ \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial \varphi^2} \right]$$

with boundary condition (BC)

$$\frac{\partial u}{\partial n} +h(u-u_c)=0$$

where $$u_c$$ is the temperature of the surrounding material, and initial condition (IC)

$$u(r,\varphi,t=0)=T_0$$

where we can assume that $$T_0$$ is constant.

Separation of variables leads to 3 differential equations. We assume now also that the problem is axisymmetric leaving only 2 diff. eqs. according to

$$T'(t)=-a \lambda^2 T(t)$$
$$r^2 R''(r) + r R'(r) + \lambda^2 r^2 R(r) - \mu R(r) =0$$

Solution to the first equation will be on the form

$$T(t)=A e^{-a \lambda^2 t}$$

and by making the variable substitution $$x=\lambda r$$ the second equation turns in the well known Bessel equation with solutions according to

$$X=c_1 J_0(x)$$

where $$J_0(x)$$ is a Bessel function of first kind and zero order. Thus, the solution is on the form

$$u(r,t)=cJ_0(\lambda r) e^{-a \lambda^2 t}$$

I believe that this is a standard result in meny textbooks. Standard procedure in these books is then to proceed by explaining how to determine what values on $$\lambda$$ are valid and how the coefficients $$c_k$$ in the Bessel function are identified. This is pretty straight forward in the case of a fixed temperature on the boundary. But how is it done with the BC stated above? If someone could show me the steps, or point me to relevant literature, I would be very thankful.


Johan Ekh

 

[dhris]

You are most of the way there, you just need to plug your proposed solution into the boundary condition. First of all, since we're dealing with the unit disk, the normal direction is r, so n=r. Secondly, you'll want to rid yourself of the pesky $$u_c$$ in the boundary condition. You can either make the substitution

$$v=u-u_c$$

or just realize that you can add $$u_c$$ to your solution and it still solves the equation, so actually you have

$$u=c J_0(\lambda r) e^{-a \lambda ^2 t} + u_c$$ .

(I'm assuming you have solved the equations correctly. I take it that $$\mu$$ is the eigenvalue associated with the axial dependence and you set it to 0?) Anyway, now take your solution and plug it into your B.C.

$$\frac{\partial u}{\partial r}+h(u-u_c)=c \lambda J_0'(\lambda r) e^{-a \lambda ^2 t}+h c J_0(\lambda r) e^{-a \lambda ^2 t}$$

Evaluating at r=1 and canceling the time dependence gives
$$\lambda J_0'(\lambda)+h J_0(\lambda)=0$$

Unfortunately, that's your condition for $$\lambda$$. You can use the identity,

$$J_0'(\lambda)=-J_1(\lambda)$$

which is found here:http://www.math.sfu.ca/~cbm/aands/page_361.htm
to get

$$\lambda J_1(\lambda)=h J_0(\lambda)$$.

This has many solutions $$\lambda_n$$ but they are not expressible in terms of elementary functions as far as I know, especially not with arbitrary h.

So now you have

$$u=\sum_n c_n J_0(\lambda_n r) e^{-a \lambda_n ^2 t} + u_c$$

and you want to know the $$c_n$$ from the initial condition. That is you need

$$\sum_n c_n J_0(\lambda_n r) + u_c =T_0$$ (*)

For this, you need some kind of orthogonality for the Bessel functions, which is found here:
http://www.math.sfu.ca/~cbm/aands/page_485.htm
Direct your attention to 11.4.5. Yours is the case a=h and b=1 (and $$\nu$$=0 of course). Now you multiply both sides of (*) by $$r J_0(\lambda_m r)$$, integrate from 0 to 1 and use the orthogonality condition to find the coefficients. I will leave that to you, though.

 

[johan_ekh]

Thank you dhris,
it works exactly as I want. I will continue with non-constant initial conditions and also a heat source term. Thus I might still run into some problems! In that case I will continue this thread. If I succeed in solving the problem completely I will give the solution in this thread in case anyone is interested.Johan