How to Solve the Heat Equation with Neumann Boundary Conditions?

johan_ekh
Messages
3
Reaction score
0
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,

<br /> \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] <br />

with boundary condition (BC)

<br /> \frac{\partial u}{\partial n} +h(u-u_c)=0<br />

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

<br /> u(r,\varphi,t=0)=T_0<br />

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&#039;(t)=-a \lambda^2 T(t)
r^2 R&#039;&#039;(r) + r R&#039;(r) + \lambda^2 r^2 R(r) - \mu R(r) =0

Solution to the first equation will be on the form

<br /> T(t)=A e^{-a \lambda^2 t}<br />

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

<br /> X=c_1 J_0(x)<br />

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

<br /> u(r,t)=cJ_0(\lambda r) e^{-a \lambda^2 t}<br />

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
 
Physics news on Phys.org
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&#039;(\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&#039;(\lambda)+h J_0(\lambda)=0

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

J_0&#039;(\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.
 
Last edited:
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
 
Are there any good visualization tutorials, written or video, that show graphically how separation of variables works? I particularly have the time-independent Schrodinger Equation in mind. There are hundreds of demonstrations out there which essentially distill to copies of one another. However I am trying to visualize in my mind how this process looks graphically - for example plotting t on one axis and x on the other for f(x,t). I have seen other good visual representations of...
Back
Top