Cylinder with heat generation, Separation of variables

[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!

 

[Orodruin]

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.

 

[16universes]

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.

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.

 

[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.

 

[Orodruin]

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).

 

[16universes]

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.

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.

 

[16universes]

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).

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?

 

[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.

 

[Orodruin]

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

 

[Chestermiller]

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

Chet

 

[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)$

 

[Orodruin]

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.

 

[Chestermiller]

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)$

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

 

[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.

 

[Orodruin]

I expect the steady state solution to be some complex oscillating function that will probably have to be expressed as an infinite sum.

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

 

[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.

 

[Orodruin]

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$

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.

 

[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?

 

[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.

 

[Chestermiller]

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.

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

 

[Orodruin]

The oscillating solution is not a reasonable solution to the problem you have set up. Your solution is based on heat conduction and will therefore be of exponential character.

Furthermore, if you assume that heat conductivity is large (which is equivalent to saying that the object is thin and takes the same temperature across), then you will have the same temperature in the entire thing (the solution cannot depend on z as neither your initial values or source does (and the constant function is an eigenfunction to the Laplace operator). You could then just as well formulate the problem as an ODE directly:
$$
\partial_t T = - \kappa T + \beta
$$
for some constants $\kappa, \beta$. If radiation is a significant contribution to cooling, you will also have a term for heat loss through radiation on the RHS, subtracting something proportional to $T^4$ (and changing $\beta$) and the area of the object.

The stationary solution to this equation is given by letting the LHS = 0 yielding
$$
\kappa T = \beta
$$
or solving a 4th order equation if you add the radiation term.
If you linearize the problem around this solution you will find a solution that exponentially approaches the stationary solution.

 

[16universes]

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 .
Chet

Take a thermistor and run constant current through it. What long term behavior would you expect to see for the temperature? That's the problem I'm faced with.

Now, I've set it up according to how I view the problem. Everyone has commented on how to solve the problem I've set up, but no one has commented on how I've set up the problem. Think about how YOU would expect temperature to behave. Is an oscillation of some kind not reasonable? If so, what should we expect to see? I think this is an important question to answer, since the model should reflect our expectations within reason.

 

[16universes]

The oscillating solution is not a reasonable solution to the problem you have set up. Your solution is based on heat conduction and will therefore be of exponential character.

I do agree with this. Which leads me to believe that I have set it up incorrectly. Before setting up any PDE, in your opinion, how should the system behave? Disregard any assumptions I have made. You have a thermistor and run constant current through it. What is its behavior?

 

[Orodruin]

Why do you expect the solution to be oscillating? I would expect it to heat up until equilibrium is reached between radiated and added energy.

I did not comment on how you set up the problem because it was relatively reasonable apart from the missing boundary conditions. As I just said, if the object is small in comparison to the heat conduction you do not need a spatial part of the equation.

 

[16universes]

After thinking about it more, my oscillation assumption was incorrect. You are absolutely right. In the steady state, heat can only be generated at the rate it can be dissipated.

 

[Chestermiller]

What's the L/D ratio of the cylinder? If it's large, you can neglect axial heat conduction. That's what I was assuming when I suggested that the long term differential equation reduces to an ode.

You were saying that this is nothing like the heat transfer problems you encountered in graduate school. In my judgement, this is only a minor modification to problems we had in graduate school, involving a heat generation rate that depends on temperature, rather than a constant heat generation rate.

Chet

 

[16universes]

I would say L/D >> 1

I never dealt with heat transfer. Scary, I know. Jackson E&M was the majority of my experience with PDEs. I didn't get much out of it because most of the time the material was way over my head. If I can't visualize the problem and mathematics, I struggle to understand it.

 

[Chestermiller]

With L/D large, the heat flow out the ends of the cylinder is going to be negligible compared to the radial heat flow. So you can neglect the axial derivative terms in the Laplacian. So the Laplacian reduces to the radial derivative terms.

At long times, the temperature will be a function only of radial position. Show us how you would solve for that limit, please.

Chet

 

[AlephZero]

After thinking about it more, my oscillation assumption was incorrect. You are absolutely right. In the steady state, heat can only be generated at the rate it can be dissipated.

There will always be a steady state solution, but whether it is stable or not might depend on whether you have a NTC or PTC thermistor. Possibly, a small perturbation from the stable solution would give you "thermal runaway", i.e. the solution would diverge to infinity.

But even if the steady state solution is unstable, that doesn't necessarily mean it oscillates. As somebody else said, the differential equation for linear thermal conduction is first order. To get oscillations, you need some more terms in the equation from the nonlinear behavior of the thermistor, and you need the higher order system to have zero damping, otherwise any oscillations will die out, or never even start.

The thermistor is tiny indeed. It's probably reasonable to assume that temperatures don't vary with radius.

If you are looking for non-steady-state solutions, I wouldn't make that assumption. In real life the temperature does vary with position inside the thermistor. Heat is generated everywhere, so there must be a temperature gradient to conduct the heat to the surface.

Of course you can get an approximate steady state solution by assuming the boundary temperature does not vary with position and ignoring the internal temperature gradient, but that is probably throwing out the effects you are looking for (if the exist at all).

 

[16universes]

With L/D large, the heat flow out the ends of the cylinder is going to be negligible compared to the radial heat flow. So you can neglect the axial derivative terms in the Laplacian. So the Laplacian reduces to the radial derivative terms.

At long times, the temperature will be a function only of radial position. Show us how you would solve for that limit, please.

Chet

I had to put this problem on hold. I will post my approach/solution in the not so distant future when I get some time to work through it.

 

[16universes]

First, I would like to thank everyone's help with this problem. I spent some time thinking about it, and realized that I probably made it more complex than it really was. After discussing the problem with my advisor, he suggested the following solution:

Let $T$ be the temperature of the thermistor, and $R = R_0 + \alpha T$ be the resistance of the thermistor as a function of its temperature. Note that $\alpha < 0$. Heat generated by the thermistor at constant current is $I^2 (R_0 + \alpha T)$. Assume heat loss due to the airflow is $h_c (T - T_0)$, where $h_c$ is the heat transfer coefficient.

Combining everything together:

$$C_p \frac{dT}{dt} = I^2 (R_0 + \alpha T) - h_c (T-T_0)$$

where C_p is the heat capacity of the thermistor. If we define $m \equiv I^2 \alpha - h_c$ and $b \equiv I^2 R_0 + h_c T_0$, then this can be rewritten as

$$C_p \frac{dT}{dt} = mT + b$$

Note that $m < 0$ and $b > 0$.

This has the following solution:

$$T(t) = k_1 e^{\frac{m}{C_p}t} - \frac{b}{m}$$

Again, noting that $m<0$ and $b>0$, the temperature asymptotically approaches a positive value. Furthermore, according to www.engineeringtoolbox.com, $h_c$ can be approximated as $h_c \approx 10.45 - v + 10v^{1/2}$, so we can also determine $T$ as a function of airflow velocity.

I'm happy with this solution. I was way overthinking it, but this was quite elegant. What do you guys think?

 

[Orodruin]

$$C_p \frac{dT}{dt} = mT + b$$
...
I'm happy with this solution. I was way overthinking it, but this was quite elegant. What do you guys think?

...
You could then just as well formulate the problem as an ODE directly:
$$
\partial_t T = - \kappa T + \beta
$$
for some constants $\kappa, \beta$.
...
The stationary solution to this equation is given by letting the LHS = 0 yielding
$$
\kappa T = \beta
$$
or solving a 4th order equation if you add the radiation term.
If you linearize the problem around this solution you will find a solution that exponentially approaches the stationary solution.

This is what I think ;)

 

[16universes]

Ha I apologize. I was still caught up "seeing" the problem differently. It took some time to change how I was looking at it. But yes, you were correct. Thanks