Download Derivation and Solution of the One-Dimensional Heat Equation and more Summaries Differential Equations in PDF only on Docsity!
The heat equation
Many physical processes are governed by partial differential equations. One such phenomenon is the temperature of a rod. In this chapter, we will examine exactly that.
1 Deriving the heat equation
1.1 What is a partial differential equation?
In physical problems, many variables depend on multiple other variables. For example, the temperature u(x, t) [K] can depend on both position and time. Such variables don’t have normal derivates like du/dt. Instead, they have partial derivates, like ∂u/∂x and ∂u/∂t.
We can set up an equation with multiple partial derivatives. We would then get a partial differential equation (PDE). So a partial differential equation is an equation containing partial derivatives. If a differential equation does not contain partial derivatives, it’s only an ordinary differential equation (ODE).
1.2 Conservation of energy for a one-dimensional rod
Let’s consider a one-dimensional rod of length L. We define the thermal energy density e(x, t) [J/m^3 ] as the energy per unit volume. It depends not only on the position in the rod, but also on time. This is because it can change as time passes by.
There are two reasons why e can vary in time. First, there is the heat flux φ(x, t) [J/m^2 s]. This is the heat flowing to the right through a unit cross-section per unit time. Second, there can be internal heat creation due to heat sources. The amount of heat created is denoted by Q(x, t) [J/m^3 s].
We can now apply the law of conservation of energy to our rod. Let’s examine a thin slice. To be more specific, we examine the rate of change of energy in it. This must be equal to the heat created, plus the heat flowing in, minus the heat flowing out. This gives us
∂e ∂t
∂φ ∂x
+ Q. (1.1)
This equation is called the integral conservation law.
1.3 Deriving the heat equation for a one-dimensional rod
Let’s define u(x, t) [K] as the temperature in the rod. There is a relation between this temperature u, and the variable e. But to find it, we first have to define two other variables.
Let’s define the mass density ρ(x) [kg/m^3 ] as the mass per unit volume. (Usually ρ varies with temperature, and thus also with time. However, these variations are usually small. We thus neglect them.) We also define the specific heat c(x) [j/kg K] as the heat necessary to raise the temperature of a 1kg-mass by 1 Kelvin. (We assume it to be constant in time for the same reasons as for the mass density.)
Using the above definitions, we can derive that
e(x, t) = c(x)ρ(x)u(x, t). (1.2)
This transforms the integral conservation law into
cρ
∂u ∂t
∂φ ∂x
+ Q. (1.3)
Usually Q is known. But u and φ are not. So we still have two unknowns. But you may be wondering, doesn’t the heat flow depend on the temperature as well? In fact, it does. The heat flow depends on temperature differences. The higher these differences, the higher the heat flow. So we can say that
φ(x, t) = −K 0 ∂u ∂x
This equation is called Fourier’s law of heat conduction. The variable K 0 (x) [J/Kms] is called the thermal conductivity. (It is also assumed to be constant for varying temperatures.) The above law now reduces our differential equation into
cρ ∂u ∂t
= K 0
∂^2 u ∂x^2
+ Q. (1.5)
If there are no heat sources (and thus Q = 0), we can rewrite this to
∂u ∂t
= k
∂^2 u ∂x^2
, where k =
K 0
cρ
The important equation above is called the heat equation. By the way, k [m^2 /s] is called the thermal diffusivity.
1.4 Initial and boundary conditions
When solving a partial differential equation, we will need initial and boundary conditions. But what conditions do we exactly need?
If we look at the heat equation, we see that there is only a first time-derivative of u. So we need only one initial condition (IC). (An initial condition is a condition at t = 0.) Usually such a condition takes the form u(x, 0) = f (x).
However, the heat equation contains a second derivative with respect to x. So we will need two bound- ary conditions (BC). (A boundary condition is a condition at a specified position.) These boundary conditions are usually the temperatures at the edges of the rod. So, u(0, t) = T 1 (t) and u(L, t) = T 2 (t).
However, it is also possible to set the heat flow φ (or equivalently ∂u/∂x) at the edges of the rod. We would then have values given for ∂u(0, t)/∂x and ∂u(L, t)/∂x. If the rod is perfectly insulated at its edges, then φ = 0 and thus also ∂u/∂x = 0 at the edges.
It is of course also possible to combine the two possibilities above. In that case, we deal with Newton’s law of cooling. The heat flow then depends on the difference in temperature with respect to a certain reference temperature uB (t). This would give us
−K 0 (0)
∂u ∂x
(0, t) = −H (u(0, t) − uB (t)) and − K 0 (L) ∂u ∂x
(L, t) = H (u(L, t) − uB (t)). (1.7)
Here H is the heat transfer coefficient. Note that if H = 0, we are dealing with an insulated edge. On the other hand, if we have H = ∞, we would have a constant temperature at the edge.
2 Special cases of the heat equation
2.1 Perfect thermal contact
Let’s suppose we have two rods of length L. We can connect them to each other, such that their edges are in contact. So one rod goes from x = 0 to x = L, while the second goes from x = L to x = 2L. We
3 Basic concepts needed to solve the heat equation
It is almost time for us to solve the heat equation. However, before we do that, we will have to look at some other things first.
3.1 Linear operators and linear equations
A linear operator is some operator L for which
L(c 1 u 1 + c 2 u 2 ) = c 1 L(u 1 ) + c 2 L(u 2 ), (3.1)
where c 1 and c 2 are constants. For example, the heat operator
∂ ∂t
− k
∂^2
∂x^2
is a linear operator. A linear equation for u is an equation of the form L(u) = f , with the function f known. If f = 0, we have a linear homogeneous equation. Linear homogeneous equations have a certain advantage. We can apply the principle of superposition to it. Suppose we would have two solutions u 1 and u 2. Then also c 1 u 1 + c 2 u 2 is a solution, for any constants c 1 and c 2.
3.2 Orthogonality
The property of orthogonality comes in very handy when solving heat equations. So let’s examine it. We say that two function f (x) and g(x) are orthogonal on the interval [0, L] if
∫ (^) L
0
f (x)g(x)dx = 0. (3.3)
It can be shown that the functions sin(nπx/L) and sin(mπx/L) (with n and m positive integers) are orthogonal if n 6 = m. And the same goes for cosines. That comes in handy! In fact, the general rules for the interval [0, L] are ∫ (^) L
0
sin nπx L
sin mπx L
dx =
0 if m 6 = n, L/ 2 if m = n.
∫ L
0
cos nπx L
cos mπx L
dx =
0 if m 6 = n, L/ 2 if m = n 6 = 0, L if m = n = 0.
Sometimes, however, we are examining the interval [−L, L]. In this case all values double. So,
∫ (^) L
−L
sin nπx L
sin mπx L
dx =
0 if m 6 = n, L if m = n.
∫ L
−L
cos
nπx L
cos
mπx L
dx =
0 if m 6 = n, L if m = n 6 = 0, 2 L if m = n = 0.
By the way, the functions sin(nπx/L) and cos(mπx/L) are always orthogonal on the interval [−L, L]. So for every n and m we have (^) ∫ L
−L
sin
nπx L
cos
mπx L
dx = 0. (3.8)
4 Solving method for the heat equation
In this part we will present a basic method to solve the heat equation.
4.1 Introducing the method of separation of variables
Let’s try to solve the homogeneous heat equation
∂u ∂t
− k
∂^2 u ∂x^2
Of course, there are also an initial condition u(x, 0) = f (x) and two boundary conditions. To solve this problem, we use the method of separation of variables. According to this method, we assume that we can write u(x, t) as X(x)T (t). Here, the function X(x) only depends on x and T (t) only depends on t. We can now rewrite the above equation to
1 kT
dT dt
X
d^2 X dX^2
We have reduced the PDE to an ODE! It can also be shown that both sides of the above equation equal a certain constant −λ. Here, λ is called the separation constant. We thus get two ordinary differential equations, being d^2 X dx^2
= −λX and
dT dt
= −λkT. (4.3)
Solving the latter one is easy. The solution is
T (t) = ce−λkT^ , (4.4)
where c is a constant. (It depends on the initial conditions.) However, solving the equation for X is a bit more difficult. In fact, we will find that it can only be solved for certain values of λ. These values are called the eigenvalues of the equation. The corresponding solutions for X(x) are the eigenfunctions. Let’s take a look at how we can find them.
4.2 Finding the eigenvalues and the eigenfunctions
We want to solve the ODE d^2 X dx^2 = −λX. (4.5)
We see that X = 0 is a solution. We call this the trivial solution, in which we are not interested. If we ignore this solution, we can distinguish three cases:
- λ < 0. In this case the general solution for X(x) is
X(x) = c 1 e
√ −λx (^) + c 2 e− √ −λx. (4.6)
When applying boundary conditions, we usually only find the trivial solution. Only in certain special cases will there be eigenvalues λ < 0.
- λ = 0. Now we would find the solution X(x) = c 1 x + c 2. Using the boundary conditions, we can solve for c 1 and c 2. Sometimes it turns out that c 1 = c 2 = 0. In this case X(x) is the trivial solution, and λ = 0 is not an eigenvalue. Sometimes, however, X(x) is not the trivial solution. In this case λ = 0 actually is an eigenvalue, with the corresponding eigenfunction X(x).
However, there is one condition we haven’t satisfied yet: the initial condition. And we can use this condition to find the constants Bn. To satisfy the initial condition, we must have
f (x) = u(x, 0) =
∑^ ∞
n=
Bn sin
( (^) nπx L
Now we must apply the property of orthogonality to find the coefficients Bn. To do this, we multiply by sin(mπx/L) and integrate from 0 to L. We then get
∫ (^) L
0
f (x) sin
( (^) mπx L
∑^ ∞
n=
Bn
∫ L
0
sin
( (^) nπx L
sin
( (^) mπx L
Now we can note that every term in the sum on the right drops out (is zero), except for the term with n = m. The right side of the equation thus reduces to BmL/2. It follows that
Bm =
L
∫ L
0
f (x) sin
mπx L
dx. (5.6)
Using this equation, we can find our constants Bn, and thus also our unique solution for u(x, t).
5.2 Second example: Both edges being insulated
In the second example we examine a rod with both its edges insulated. So our boundary conditions are ∂u/∂x(0, t) = 0 and ∂u/∂x(L, t) = 0.
If λ < 0, then we can again show that X(x) = 0 as well. So we only find the trivial solution.
Let’s examine the case where λ = 0. We now do find a non-trivial solution. Once more we have X 0 (x) = c 1 x + c 2. Both boundary conditions imply that c 1 = 0. However, c 2 can be anything. So we have found a non-trivial solution. Thus λ = 0 is an eigenvalue! The corresponding eigenfunction is X 0 (x) = 1. (Remember that we were allowed to ignore constants when examining eigenfunctions.)
Now let’s consider the case λ > 0. This time our solutions will be
λn =
( (^) nπ L
and Xn(x) = cos
( (^) nπx l
The resulting general solution of our PDE (before applying the initial conditions) will then be
u(x, t) = A 0 +
∑^ ∞
n=
An cos
( (^) nπx l
e−λkt. (5.8)
The constants A 0 and An follow from our initial condition. This time we must multiply by cos(nπx/L) and integrate from 0 to L. We then find that
A 0 =
L
∫ L
0
f (x) dx and An =
L
∫ L
0
f (x) cos
nπx L
dx. (5.9)
5.3 Laplace’s equation
We can also extend the method of separation of variables to a two-dimensional plate with width L and height H. However, things get more difficult now. So to make it easier, we only want to find the steady-state solution (with ∂u/∂t = 0). This turns our differential equation into Laplace’s equation, being ∂^2 u ∂x^2
∂^2 u ∂y^2
To demonstrate the solution method, we will show one example. Let’s assume that we have boundary conditions u(0, y) = 0, u(L, y) = f (y), u(x, 0) = 0 and u(x, H) = 0. To apply the method of separation of variables, we assume that u(x, y) = X(x)Y (y). Our differential equation now turns into
1 X
d^2 X dx^2
Y
d^2 Y dy^2
= λ. (5.11)
So we again have two ODEs. For Y (y), we have two rather easy boundary conditions, being Y (0) = 0 and Y (H) = 0. So let’s focus on the y-part. We can solve this using methods we have seen earlier. We find as eigenvalues and eigenfunctions
λn = (nπ/H)^2 and Yn(y) = sin
nπy H
Now let’s find X(x) as well. we now have to solve
d^2 X dx^2
= λX, or equivalently
d^2 X dx^2
( (^) nπ H
X. (5.13)
The solution for this equation is
X(x) = c 1 e
√ λx (^) + c 2 e − √ λx (^) or equivalently X(x) = c 1 cosh^
λx + c 2 sinh
λx. (5.14)
We can use either of the above relations. However, in this case the relation with sinh and cosh is more convenient, since then one term will drop out. So let’s use that one.
One of our boundary conditions is X(0) = 0. If we apply this, we find that c 1 = 0 and the part with cosh will disappear. We thus have as eigenfunctions Xn(x) = sinh
λx. Our general solution for u then becomes
u(x, y) =
∑^ ∞
n=
cn sin nπy H
sinh nπx H
However, we haven’t applied one boundary condition yet, being u(L, y) = f (y). We can use this to find the constants cn. Inserting x = L in the above equation gives
f (y) = u(L, y) =
∑^ ∞
n=
cn sin
nπy H
sinh
nπL H
It can then be derived, using the porperty of orthogonality, that
cn =
H sinh nπLH
∫ H
0
f (y) sin
nπy H
dy. (5.17)
And this completes our solution to the problem.
