Transient heat conduction in 1D solid

[nickpestell]

1D solid, 0<x<L, with the following boundary conditions:

The whole solid is at T = T1 at t=0. x = 0 is held constant at initial temperature T1 for all t. There is a constant flow of heat, dQ/dt out of the solid at x = L.
T(0,t) = T1,
T(L,0) = T1,

How do we go about solving the heat equation for such a problem? Does the system reach steady state so that the temperature gradient across the solid is constant in time?

Ultimately, what I really need to be able to work out is the time taken for the system to reach steady state.

Any help would be greatly appreciated.

 

[Chestermiller]

Hi nickpestel. Welcome to Physics Forums!

The way this works is that you have to write down the transient heat conduction and tell us about what kind of attempt you made to solve it. Also, what does your intuition tell you about whether the system will reach steady state? If it does eventually reach steady state, what do you get for that solution? (This should be easy because the time-dependent term drops out)

Chet

 

[nickpestell]

Hi Chet.

I'm not quite sure what you mean by "write down the transient heat conduction". I'm fairly sure the system would reach equilibrium where difference in temperature between the two ends provides a flux through the solid which is the same as that which occurs at x = L. I suppose that further boundary conditions would be that dT/dx = 0 for all x at t = 0 and that dT/dx = 0 for all t at x = 0, as well dT/dx = (T2-T1)/L and dT/dt = 0 for all x and t>ts, where ts is the time taken to reach steady state and T2 is the steady state temperature at x = L. I don't really know how to incorporate the constant flux at x = L.

I'm really unsure how to tackle this problem and I should point out that it is a practical situation which I am trying to model, so I'm unsure if there is a neat solution. I'm sorry if this is unclear but If you could give some direction that would be great!

Thanks
Nick

 

[Chestermiller]

There actually is a neat solution to this problem. But first, for the final steady state solution:

$$T=T_1-\frac{\dot{Q}x}{kAL}$$

where $\dot{Q}$ is the rate of heat removal at x = L, A is the cross sectional area of the sample, and k is the thermal conductivity.

Now. Back to your problem. Suppose that bar were twice as long, running from x = -L to x = +L, and, at the boundary x = -L, rather than removing heat from the bar, you were adding heat at the same rate, such that the heat fluxes at the two boundaries are:

$q = -k\frac{\partial T}{\partial x}=\frac{\dot{Q}}{A}$ at x = +L

$q = -k\frac{\partial T}{\partial x}=\frac{\dot{Q}}{A}$ at x = -L

Because of the symmetry of these boundary conditions, the temperature at x = 0 would not change from its initial value of T₁ over all times. So in the region between x = 0 and x = L, the solution to this problem would be exactly the same as the solution to your problem.

Now you have reduced it to a problem in which the heat flux is specified at the two ends of the bar for all times greater than t = 0, and the heat flux is zero everywhere within the bar at time = 0. We can solve this problem entirely in terms of the heat flux, by taking the partial derivative of the transient heat conduction equation with respect to x, to obtain:

$$\frac{\partial^2 T}{\partial t\partial x}=\kappa\frac{\partial^3 T}{\partial x^3}$$

where κ is the thermal diffusivity.

If we multiply this equation by the thermal conductivity k, this equation is equivalent to
$$\frac{\partial q}{\partial t}=\kappa\frac{\partial^2 q}{\partial x^2}$$

We now have a partial differential equation in terms of the heat flux, in which the heat flux value is specified at the two boundaries, and in which the heat flux is zero throughout the domain initially. This problem is readily solved by separation of variables. You can then integrate the resulting equation for the heat flux with respect to position x to get the temperature.

I realize that this is a very sketchy description, but to describe the solution method in greater detail would require much more space on PF.

 

[nickpestell]

Thanks very much Chet! That was a great help. I have now solved the differential equation, or at least I think I have. I was wondering if you would check my working. I apologize in advance for not using latex but I don't know how.

Firstly I defined a new function:

U(x,t) = q(x,t) - q(x,t>tss), where tss is the time taken for the system to reach steady state. q(x,t>tss) is a constant and is equal to the initial flux through both ends, call it q1.

So U(x,t) = q(x,t) - q1 obeys the same differential equation as q since dq/dt = du/dt and d2q/dx2 = d2u/dx2. The advantage of using U instead of q is that it has the boundary conditions that U(L,t) = q(L,t) - q1 = q1 - q1 = 0, and
U(-L,t) = q(-L,t) - q1 = q1 - q1 = 0.

So the equation to solve is dU/dt = k d2U/dx2. (1)

Assume U is separable:

U(x,t) = X(x)Y(t)

Subs X(x)Y(t) into (1):

XdY/dt = Ykd2X/dx2 = -c^2, where -c^2 is the separation constant.

Solve the X equation:

d2X/dx2 = -Xc^2

X = Asin(cx) + B cos(cx)

Using the boundary conditions to find c gives:

A/B tan(cL) = 0

c = npi/L

Solve the Y equation:

dY/dt = -Ykc^2

Y = Dexp(-c^2kt)

Putting the Y and X part together and absolving constants:

U(x,t) = (Asin(npix/L) + Bcos(npix/L))exp-(npik/L)^2t

By linearity:

U(x,t) = (from n=0 to n=∞)Ʃ(Ansin(npix/L) + Bncos(npix/L))exp(-t(npik/L)^2)

Then use initial condition to get rid of t dependence, and find the Ans and Bns as Fourier coefficients.

U(x,0) = q(x,0) - q1 = q1 - q1 = 0, at x = ± L
= 0 - q1 = -q1 otherwise

The Fourier coefficients are An = 0 and Bn = (-2q1/npi)sin(npi).

So

U(x,t) = (from n=0 to n=∞)Ʃ(-2q1/npi)sin(npi)cos(npix/L)exp(-t(npik/L)^2).

q(x,t) = U(x,t) + q1,

q(x,t) = (from n=0 to n=∞)Ʃ(-2q1/npi)sin(npi)cos(npix/L)exp(-t(npik/L)^2) + q1,

q1 = - (k/L)(T(L,t>tss) - (T(0,t>tss)) (by Fourier's law).

The part of my working which I am most unsure about is working out the Fourier coefficients since U(x,0) is a discontinuous function, although as far as I can tell the contributions at x = ± L have no effect since they are equal to zero so the cos coefficient would be given by

Bn = (from-L to L)(1/L)∫-q1cos(npix/L)dx.

 

[Chestermiller]

As best I can tell, it doesn't look like the boundary conditions on B_n were applied properly. It seems to me that, according to your development B_n should be zero.

I have a suggestion. Before proceeding with the separation of variables, make a simple algebraic substitution: y = x + L. Then the boundary conditions will be at y = 0 and y = 2L. Then the solution will be:

$$U=\sum_1^\infty A_n\sin {\frac{n\pi y}{2L}}\exp (-(\frac{n\pi}{2L})^2\kappa t)$$

The initial condition is U = -q1 at t = 0. To get the Fourier coefficients, you would integrate between y = -2L and y = + 2L, so that you could cover the full period of the sine function. And, of course, in the region y < 0, you would use U = +q1.

 

[nickpestell]

Surly if you change the variable then the initial condition is not just U = -q1 at t = 0 for all y, instead U = 0 at y = 0 and y = 2L, and then U = -q1 elsewhere. This is from the initial condition that the flux is q1 at x= ± L. Does this not make the calculation of the Fourier coefficients very difficult?

 

[Chestermiller]

Surly if you change the variable then the initial condition is not just U = -q1 at t = 0 for all y, instead U = 0 at y = 0 and y = 2L, and then U = -q1 elsewhere. This is from the initial condition that the flux is q1 at x= ± L. Does this not make the calculation of the Fourier coefficients very difficult?

No. I said U is equal to -q1 from 0<y<2L, and U is equal to +q1 from -2L <y<0. But, at the boundary y = 0, U is equal to 0, which is the average of the values from 0<y<2L and -2L <y<0. At y = +2L and y = -2L, U is again zero, because of the periodicity. Check this out in my series solution for U. Basically, initially, U is a square wave, alternating between +q1 and -q1, and we are fitting a Fourier series to this square wave. Of course, from the Fourier series, U is zero at y = 0 and at y = 2L for all times. I hope this makes sense. Just flesh out the entire square wave on a graph (i.e., at higher and lower values of y beyond our range) to see what I'm saying

 

[nickpestell]

Ok, so I can see the square wave and it's periodicity. Am I right in thinking that, since what we are trying to find is U(y,t) from 0 to 2L, function from -2L to 0 is set up in this fashion just for the convenience of working out the Fourier coefficients?

I get An=(2q1/npi)cos(npi).

 

[Chestermiller]

Ok, so I can see the square wave and it's periodicity. Am I right in thinking that, since what we are trying to find is U(y,t) from 0 to 2L, function from -2L to 0 is set up in this fashion just for the convenience of working out the Fourier coefficients?

I get An=(2q1/npi)cos(npi).

Yes. You can also think of it as a "half wave" series. Look up half wave series' in your textbook.

As far as what you got for the coefficient, I haven't had a chance to check it yet, but my initial inclination is that it doesn't look right. I will check it later.

You can also look up tables of half-wave Fourier series for different functions in Abramowitz and Stegan, or Janke and Emde after you've had some experience in doing the integrals.

 

[Chestermiller]

I worked out the Fourier coefficient for the square wave, and it is given by:
$$A_n=-q_1\frac{2}{\pi}\frac{1-(-1)^n}{n}$$
Incidentally, $\cos (n\pi)=(-1)^n$

You can also find the derivation of the Fourier coefficients for this particular square wave in Kreyzig.

 

[nickpestell]

Thank you very much Chet! You're a total legend!