# Cylinder with heat generation, Separation of variables

1. Jul 30, 2014

### 16universes

I'm having some difficulty setting up a problem. I'm trying to model the temperature of a thermistor connected to a constant current source. The thermistor's resistance varies with temperature, so with a fixed current, I would expect to see the thermistor's temperature to oscillate with time. For simplicity, I'm only interested in the surface temperature $T_s$ as a function of time. I'm trying to derive the result analytically. This is how I am setting it up:

Suppose you have a cylinder with some radius $a$, and length $L$. It has internal heat generation per unit volume $Q(T) = \frac{I^2 R(T)}{\pi a^2 L}$, where $I$ is a fixed current, and $R$ is a resistance that varies with temperature. Assume $R≈-mT+b$ where $T$ is the temperature of the cylinder and $m$ is a constant of proportionality. Assume that heat only flows through the cylinder's side surface (not through it's ends).

At $t=0$, $T = T_0$ everywhere.

Using the heat equation:

$$\rho C \frac{\partial T}{\partial t} = \kappa \nabla^2 T + Q(T)$$

where $\rho$ is the density, $C$ is the heat capacity, $\kappa$ is the thermal conductivity. For simplicity, assume these are all constant.

If we replace $Q(T)$ with $\left(\frac{I^2}{\pi a^2 L}\right) (-mT+b) = -\mu T + \beta$, then:

$$\rho C \frac{\partial T}{\partial t} = \kappa \nabla^2 T - \mu T + \beta$$

My questions are:

1. Have I set this up correctly?

2. What is the appropriate way to solve this for $T(T_s,t)$?

Using separation of variables, you would normally let $T(r,t) = R(r) \Gamma (t)$, but I cannot solve it this way since the $\beta$ term makes it difficult to separate variables.

Any help on this would be much appreciated. Thanks in advance!

Last edited: Jul 30, 2014
2. Jul 30, 2014

### Orodruin

Staff Emeritus
1) Find the eigenfunctions of the spatial operator.
2) The spatial operator is Sturm-Liouville, which leads to an orthogonal basis of the functional space (functions of the spatial coordinates)
3) Expand the solution and $\beta$ (i.e., the constant function) in these eigenfunctions. Note that the expansion coefficients for the solution will depend on time.

This should give you separate ODEs to solve for expansion coefficient.

3. Jul 30, 2014

### 16universes

Don't you have to separate variables before you can find eigenfunctions? Again, the $\beta$ term does not allow me to do this.

Let $T_i = \frac{\partial T}{\partial i}$ and $T_{ii} = \frac{\partial^2 T}{\partial i^2}$ for convenience. Then:

$$\alpha T_t - \kappa \left(T_{rr} + \frac{1}{r} T_r \right) - \mu T + \beta = 0$$

where I've defined $\alpha \equiv \rho C$ to simplify the clutter.

If we seek a solution of the form $T(r,t) = R(r) \Gamma (t)$, we get:

$$\alpha R \Gamma ' - \kappa \Gamma \left(R'' + \frac{1}{r} R' \right) - \mu R \Gamma + \beta = 0$$

Dividing through by $R \Gamma$ gives:

$$\alpha \frac{\Gamma '}{\Gamma} - \frac{\kappa}{R} \left(R'' + \frac{1}{r} R' \right) - \mu + \frac{\beta}{R \Gamma} = 0$$

Unless my understanding is wrong, I cannot find eigenfunctions of $R$ because I cannot separate it from $\Gamma$, due to the $\beta$ term. Where am I going wrong?

Can you show a few steps to get me headed in the right direction? I'm totally lost on this.

4. Jul 30, 2014

### voko

Solve the homogeneous equation first. Then find a trivial solution of the inhomogeneous equation. The general solution of the latter is the sum of these.

5. Jul 30, 2014

### Orodruin

Staff Emeritus
You only need to find the eigenfunctions of the Sturm-Liouville problem. These form an (orthogonal) basis for the function space. You can use this fact to expand the solution in these for an arbitrary time. The resulting solutions will not be on separated form with time and space unless the full PDE allows this.

What voko suggests also works assuming you can find an appropriate particular solution to remove the inhomogeneity, which is not always the case. Expanding in the eigenfunctions of the SL operator always works assuming you are able to solve the resulting ODEs (for whih you have to find particular and homogenous solutions).

6. Jul 30, 2014

### 16universes

You can't do this until the variables are separated. If I am mistaken, please demonstrate how I can do this properly with what I have.

7. Jul 30, 2014

### 16universes

Thanks for the help. I vaguely remember doing this back in grad school. In this case what does the spatial SL equation look like? I thought this comes from the PDE itself. It's not like I can magically form a new PDE for the SL problem. I don't know how to arrive at the SL equation. Can you shed some light on this please?

8. Jul 30, 2014

### 16universes

I should add that I did attempt the SL route earlier. Every resource I found first separated variables, then found the eigenfunctions. I've never seen this process when the variables have not yet been separated.

9. Jul 30, 2014

### Orodruin

Staff Emeritus
So let us take a one-dimensional example to demonstrate the principle. Let us say we want to solve the following PDE
$$\partial_t^2 u - \partial_x^2 u = f(x,t)$$
(I intentionally took the inhomogeneous wave equation instead of the diffusion equation just to spice things up, the principle is the same). Let us further assume the boundary conditions:
$$u(0,t) = u(\pi,t) = 0$$
and the initial conditions
$$u(x,0) = g(x), \quad u_t(x,0) = h(x).$$

The Sturm-Liouville operator in this case is $-\partial_x^2$ with the corresponding (homogeneous) boundary conditions defining the functional space on which it is acting. The eigenfunctions are of the form
$$\phi_n(x) = \sin(nx)$$
where $n$ is an integer due to the boundary conditions. Since the operator is a SL operator, the eigenfunctions are orthogonal (given the correct inner product) and any function of $x$ can be expanded in them (in this case it is just a sinus series). We can do this for the solution $u(x,t)$ as well as the inhomogeneity $f(x,t)$ for each value of $t$:
$$u(x,t) = \sum_{n=1}^\infty u_n(t) \phi_n(x)$$
$$f(x,t) = \sum_{n=1}^\infty f_n(t) \phi_n(x)$$
where $u_n(t)$ and $f_n(t)$ are the expansion coefficients for a given $t$, which may in general depend on $t$. In addition, we can also do this to the initial conditions
$$h(x) = \sum_{n=1}^\infty h_n \phi_n(t), \quad g(x) = \sum_{n=1}^\infty g_n \phi_n(t)$$
which do not depend on $t$ and therefore have expansion coefficients which also do not.

Inserting all of this into the PDE we obtain:
$$\sum_{n=1}^\infty \left( u_n''(t) - n^2 u_n(t) - f_n(t) \right) \phi_n(x) = 0.$$
Since the $\phi_n$ are orthogonal (linearly independent is sufficient), each term of this sum must equal to zero and thus
$$u_n''(t) - n^2 u_n(t) - f_n(t) = 0$$
for all $n$. For each value of $n$, this is a linear ODE which can be solved by standard means of finding particular and homogeneous solutions. The constants in the homogeneous solution will be fixed by the initial conditions, which can be found in an analogous manner.

The eigenfunction property $\mathcal L \phi_n = n^2 \phi_n$ as well as the orthogonality property for eigenfunctions of a SL operator (in this case $\mathcal L = - \partial_x^2$) is what was of importance to arrive at these conclusions.

The SL operator in your case is the (negative of the) Laplace operator: $-\nabla^2$, on which you have to put appropriate (homogeneous) boundary conditions. The no-flow on the cylinder ends is one such condition. To find one for the cylinder surface you will have to relate the flow to the temperature (or just assume a fixed surface temperature - in this case you can homogenize the boundary conditions by shifting the entire solution). Due to the geometry, I suggest looking into eigenfunctions in cylinder coordinates in order to accommodate the boundary conditions.

Edit: typo

10. Jul 30, 2014

### Staff: Mentor

What is the boundary condition at the surface of the cylinder?

Chet

11. Jul 31, 2014

### 16universes

I suppose we can assume a convection boundary and radiative boundary at the surface.

$\partial_x T= h (T_{\infty} - T_s) + \epsilon \sigma (T_{\infty}^4 - T_{s}^4)$

12. Jul 31, 2014

### Orodruin

Staff Emeritus
The T^4 term is slightly problematic since it is non-linear. If you introduce it, you will have to linearize around the stationary solution or try to solve the resulting non-linear problem.

Also the linear term is slightly annoying to tackle, you will have to find the solution to some combination of Bessel functions and their derivatives, which typically will have to be done numerically.

13. Jul 31, 2014

### Staff: Mentor

If that's the case, then there is going to be a long time steady state solution. Do you know how to find that solution? If you find that solution first, solving for the transient portion of the solution will be much simpler.

Chet

14. Jul 31, 2014

### 16universes

I expect the steady state solution to be some complex oscillating function that will probably have to be expressed as an infinite sum. But I'm unsure of the form. The heat generation is proportional to I^2 R. But R varies linearly (approximately) with temperature (it's a thermistor). As the temperature increases, the resistance decreases. Since current is constant, eventually it's temperature will begin to decrease and the cycle repeats. The rate at which temperature increases and decreases will probably not be the same.

I'm suspecting error with how I've set up the problem though. Run constant current through a thermistor and determine the temperature. Is there perhaps a better way to go about this?

It's funny that in grad school they give you complicated problems to solve but I'm finding that setting up the problem appropriately and choosing appropriate boundary conditions that reflect the real world problem is more complicated than actually solving the problem. That's what they should be teaching.

15. Jul 31, 2014

### Orodruin

Staff Emeritus
By definition, the stationary solution does not depend on time.

16. Jul 31, 2014

### 16universes

Yes, I was about to edit that... I don't expect it to reach a steady state due to my expectation that the heat generation will oscillate for all time.

But I came across a section in this fantastic PDE book (PDE's for Scientists and Engineers, Stanley J. Farlow, pg. 58) that touches on a very similar problem, I just need to modify it to include heat generation. I'm convinced that my current setup is incorrect. Perhaps someone has a clever way to modify the problem presented in the book:

The problem has a rod insulated on the ends. It's assumed that heat diffuses in the rod along the axis (taken to be the $x$ direction), and there is convective heat loss on the rod's surface.

They use a function $u(x,t)$ to represent temperature, with the notation

$$u_a = \frac{\partial u}{\partial a}$$

and

$$u_{aa} = \frac{\partial^2 u}{\partial a^2}$$

The problem is as follows

PDE: $u_t = \alpha ^2 u_{xx} - \beta u$

BC's:
$$u(0,t) = 0$$
$$u(1,t) = 0$$

for all $t$. Note that the rod is taken to have a length of 1, with it's ends at $x = 0$ and $x = 1$

IC: $u(x,0) = \phi (x)$

They transform the problem using:

$$u(x,t) = e^{-\beta t} w(x,t)$$

This gives us a new PDE that is easy to solve:

new PDE: $w_t = \alpha ^2 w_{xx}$

We can solve this, then multiply by $e^{-\beta t}$ to get our solution. Essentially, we ignore the convection term and solve, then multiply by the convection factor at the end.

@Orodruin: This sounds similar to the process you mentioned with expanding the spatial solution. I think the explanation in the book cleared up my confusion.

If I can modify this to match my problem, then it shouldn't be too difficult. So the question is now, how do you properly set up the problem? I should have a convection term (I'll ignore the radiation boundary to make it easier), but also have a heat generation term that depends on resistance which varies with temperature. The proper way to write the heat generation is the million dollar question. Any thoughts?

And by the way, thanks to all who've been helping me with this. I don't really have anywhere else to turn. It's all very much appreciated.

17. Jul 31, 2014

### Orodruin

Staff Emeritus
These boundary conditions do not describe insulated ends. Instead they describe ends that are being kept at a constant temperature (0), which will result in a lot of heat flowing out (or in if u is smaller than 0).

The further approximation made is that the object is very thin. So thin that temperatures do not vary with the radius, which is not necessarily the case (if it is you may use something similar). They have also taken the heat loss to be proportional to $u$. These are things you would have to get out of your boundary conditions. You also have a constant source term which will not be removed by the corresponding transformation.

In fact, apart from the missing boundary conditions, I think your first approach to setup the problem was quite reasonable.

18. Jul 31, 2014

### 16universes

You're right, I meant "fixed ends". For my problem, I would require that the derivative at the ends is zero.

The thermistor is tiny indeed. It's probably reasonable to assume that temperatures don't vary with radius. Perhaps this would make the problem easier. I suppose I could also make it easier if I assume spherical geometry. I don't necessarily need the exact solution; just an approximation that models the problem within reason. I expect to see a complicated oscillation for a solution, based on the fact that resistance varies with temperature. If I don't see this, it's likely that I've set it up wrong.

Does an oscillating result sound reasonable?

19. Jul 31, 2014

### 16universes

I wonder if rewriting the problem in terms of resistance $R$ is a better approach. Ultimately, the temperature is determined through the resistance measurement. Then, I won't have that pesky $\beta$ term.

20. Jul 31, 2014

### Staff: Mentor

The problem, as you formulated it, features a heat generation rate that varies only with temperature. If that is what you intended, then your assessment of the long term behavior is incorrect . The long time temperature profile can be obtained simply by setting the time derivative equal to zero and solving the resulting ode. This solution will feature a temperature that is highest at the center of the cylinder, and which decreases monotonously with radial position to the lowest value at the surface.
First solve the problem for constant Q to get your feet wet.

Chet