Solution of the Diffusion Equation



[pic]

|College of Engineering and Computer Science

Mechanical Engineering Department

Mechanical Engineering 501B

Seminar in Engineering Analysis | |

| |Spring 2009 Class: 14443 Instructor: Larry Caretto |

February 9 Homework Solutions

Kreyszig, page 561 problems 5, 13, 14 and 27 – Finding the temperature in a one dimensional solid as a function of time.

These problems involve the solution of the diffusion equation with various boundary conditions and initial conditions. All problems start with the same separation of variables process which is described below.

[pic] [1]

Here, α = k/ρc, is the thermal diffusion coefficient with dimensions of length squared over time. The initial condition specifies the value of u at all values of x at t = 0. This condition is written as follows:

u(x,0) = f(x) [2]

In problem 5 we have f(x) = sin(0.4πx) and in problem 13, f(x) = x. The boundary conditions at x = 0 and x = xmax, can be functions of time in general. These are written as follows.

u(0,t) = uL(t) u(xmax,t) = uR(t) [3]

In problem 5 we have uL = uR = 0 and in problems 13 and 14 we have ∂u/∂x = 0 at x = 0 and x = xmax.

We use the basic separation of variables approach in which we assume that the solution, u(x,t), is the product of two functions, T(t) a function of time only and X(x) a function of the distance x only. With this assumption, our solution becomes.

u(x,t) = X(x)T(t) [4]

We do not know, in advance, if this solution will work. However, we assume that it will and we substitute it for u in equation [1]. Since X(x) is a function of x only and T(t) is a function of t only, we obtain the following result when we substitute equation [4] into equation [1].

[pic] [5]

If we divide the final equation through by the product αX(x)T(t), we obtain the following result.

[pic] [6]

The left hand side of equation [6] is a function of t only; the right hand side is a function of x only. The only way that this can be correct is if both sides equal a constant. This also shows that the separation of variables solution works. In order to simply the solution, we choose the constant[1] to be equal to −λ2. This gives us two ordinary differential equations to solve.

[pic] [7]

The first equation becomes[pic]. The general solution to this equation is known to be[2] [pic]. The second differential equation in [7] can be written as [pic]. The general solution to this equation is X(x) = Bsin(λx) + Ccos(λx). Thus, our general solution for u(x,t) = X(x)T(t) becomes (after defining AB as C1 and AC as C2)

[pic] [8]

At this point we have to consider the initial conditions, boundary conditions and property values for each problem separately. Equation [8] will be the starting point for the analysis of all four problems.

Problem 5 – Find the temperature, u(x,t) in a bar of silver (length = 10 cm, constant cross section of area 1 cm2, density 10.6 g/cm3, thermal conductivity 1.04 cal/(s cm oC), specific heat 0.056 cal/(g oC) that is perfectly insulated laterally, whose ends are kept at temperature 0oC and whose initial temperature (in oC) is sin(0.4πx).

In this case the boundary conditions are u(0,t) = 0 and u(xmax,t) = 0. If we substitute the boundary condition at x = 0 into equation [8], get the following result because sin(0) = 0 and cos(0) = 1.

[pic] [9]

Equation [9] will be satisfied for all t if C2 = 0. With C2 = 0, we can apply the solution in equation [8] to the boundary condition at x = xmax.

[pic] [10]

Equation [10] can only be satisfied if the sine term is zero. This will be true only if λxmax is an integral times π. If n denotes an integer, we must have

[pic] [11]

Since any integral value of n gives a solution to the original differential equations, with the boundary conditions that u = 0 at both boundaries, the most general solution is one that is a sum of all possible solutions, each multiplied by a different constant, Cn. This general solution is written as follows

[pic] [12]

We still have to satisfy the initial condition that u(x,0) = f(x). Substituting this condition for t = 0 into equation [12] gives

[pic] [13]

Since the ordinary differential equation for X(x) forms a Sturm-Liouville problem, the solutions Cnsin(nπx/xmax) are eigenfunction solutions to the problem in the interval 0 ≤ x ≤ xmax. Thus we can expand any initial condition function, f(x) in terms of these eigenfunctions. We can obtain the usual equation for eigenfunction expansion coefficients of Sturm-Liouville solutions by using the orthogonality relationships for integrals of the sine. If we multiply both sides of equation [13] by sin(mπx/xmax), where m is another integer, and integrate both sides of the resulting equation from a lower limit of zero to an upper limit of xmax, we get the following result.

[pic] [14]

In the second row of equation [14] we reverse the order of summation and integration because these operations commute. We then recognize that the integrals in the summation all vanish unless m = n, leaving only this integral to evaluate. Solving for Cm and evaluating the last integral[3] in equation [14] gives the following result.

[pic] [15]

For any initial condition, then, we can perform the integral on the right hand side of equation [15] to compute the values of Cm and substitute the result into equation [12] to compute u(x,t). We recognize that equation [15] is the usual equation for eigenfunction expansions in Sturm-Liouville problems.

In this problem we have f(x) = sin(0.4πx). No dimensions are given for the constant 0.4 in this equation, but we know that it must have units of inverse length to make the argument of the sine a dimensionless quantity. Since other length data are given in units of cm, we assume that the units value multiplying πx is 0.4 cm-1; since xmax = 10 cm, the value of 0.4 cm-1 = 4/xmax and we can write this initial condition as f(x) = sin[(4)πx/xmax]. If we substitute f(x) = sin[(4)πx/xmax] into equation [15] we find Cm as follows.

[pic] [16]

Because sin(mπx/xmax) is an orthogonal function in the region 0 ≤ x ≤ xmax, the only time that this integral is non zero is when m = 4. We have just shown3 that the integral of sin2(mπx/xmax) between x = 0 and x = xmax = xmax/2 so that equation [16] tells us that Cm = 1 when m = 4 and Cm = 0 for all other values of m.

[pic] [17]

Substituting this result into equation [12] gives the following solution to the diffusion equation when f(x) sin(πx/xmax) and the boundary temperatures are zero.

[pic] [18]

For the data in this problem we can compute the thermal diffusivity as follows.

[pic] [19]

The solution to the problem uses cm as the unit for length.

Problem 13 – Adiabatic means no heat exchange with the neighborhood, because the bas is completely insulated, also at the ends. Physical information: The heat flux at the ends of a bar is found to be proportional to ∂u/∂x there. Show that for the completely insulated bar, ux(0,t) = ux(L,t) = 0, u(x,0) = f(x) and separation of variables gives the following solution with An given by (2) in section 11.3.

[pic]

Again we use equation [8] as the starting point.

[pic] [8]

Taking the first derivative of this equation with respect to x gives.

[pic] [20]

The boundary condition that the ∂u/∂x = 0 at x = 0 gives the following result.

[pic] [21]

We can satisfy this condition by having C1 = 0. Setting C1 to zero in equation [20] and applying the boundary condition that ∂u/∂x = 0 at x = L gives.

[pic] [22]

We can satisfy this condition is the sine term is zero, which it will be if λL is an integer multiple of π.

[pic] [23]

As usual, we take the general solution to be the sum of all solutions that satisfy the differential equation and the boundary conditions. Note that when the eigenfunction is cos(nπx/L) we do have a nonzero eigenfunction for n = 0, so we include this in our sum. However, this eigenfunction is simply cos(0) = 1. We will write this term separately outside the sum.

[pic] [24]

As usual we find the coefficients, An, by using an eigenfunction expansion. A specific example of this is given in the next problem. The development in the next problem solution shows that the equations for these coefficients match those of equation (2) in section 11.3.

The solution gives u → A0 as t → ∞. That is, the steady-state solution is a constant temperature everywhere in the region. This seems correct. For the adiabatic boundary conditions, there is no way for heat to leave the region. The initial temperature profile will head to heat transfer in the x direction so long as there is a temperature difference. At steady state, there is no mechanism to maintain any temperature difference so that the solution approaches a constant.

Problem 14 – Find the temperature in the bar of problem 13 with L = π, α = c2 = 1 and f(x) = x.

Here we use the general solution in equation [24] that satisfies the differential equation and the boundary conditions and apply it to the given initial condition. Setting L = π in equation [24] gives λn = n and the following result for our general solution.

[pic] [25]

As usual we set t = 0 so that u(x,t = 0) becomes the initial condition, f(x) = x.

[pic] [26]

If we multiply both sides of this equation by cos(mx) and integrate the result from 0 to π, we obtain.

[pic] [27]

We have to consider the case of m = 0, in which case the eigenfunction is 1, separately. This gives.

[pic] [28]

The integral [pic] Thus, all the integrals on the right, except for the integral of A0dx are zero. (This is consistent with the notion that the set of eigenfunctions ym = cos(mπx/L) for the cosine solution, which starts at n = 0, gives y0 = 1; since the y0 eigenfunction must be orthogonal to all other eigenfunctions, we expect that the integral of any eigenfunction, other than y0, will be zero.) This gives the following result for A0.

[pic] [29]

Returning to equation [27], we find Am for m ≥ 1, using the result for the orthogonality of the cosine function to give.

[pic] [30]

Using integral tables gives the following result for the integral in the numerator.

[pic] [31]

The integral in the denominator can also be found using integral tables as

[pic] [32]

Using these results and noting that (cos(mπ) – 1) equals 2 for odd m and zero for even m gives Am as follows.

[pic] [30]

Substituting this result for Am (m ≥ 1) and the previous result for A0 into equation [25] gives.

[pic] [31]

Setting α = 1 (given) gives the following result for the first few terms in the series.

[pic] [32]

Problem 27 – (Non-homogenous heat equation) Show that the problem consisting of

ut – αuxx = Ne-ax

and (2), (3) can be reduced to a problem for the homogenous heat equation by setting u(x,t) = v(x,t) + w(x) and determining w so that v satisfies the homogenous equation and the conditions v(0,t) =v(L,t) = 0, v(x,0) = f(x) – w(x). (The term Ne-ax may represent heat loss due to radioactive decay in the bar.)

The differential equation and the boundary conditions given in the problem statement are shown below.

[pic] [33]

Substituting the proposed solution, u(x,t) = v(x,t) + w(x) into the diffusion equation gives, after some manipulation, the result shown in equation [34]. In the first step we drop the ∂w/∂t derivative since w is a function of x only. In the second step we drop the v terms since they are exactly the diffusion equation and v satisfies that partial differential equation.

[pic] [34]

Since w is a function of x only, the first and last terms in [34] give us an ordinary differential equation that we can solve for w. Writing this equation and integrating it two times gives.

[pic] [35]

[pic] [36]

[pic] [37]

We can combine the solution definition, u(x,t) = v(x,t) + w(x), with the initial condition, u(x,t) = f(x), to get an initial condition for v.

v(x,0) = u(x,0) – w(x) = f(x) – w(x) [38]

We would have to use this initial condition for v(x,0) in developing any eigenfunction expansion to match the initial conditions. The constants of integration in equation [37] can be found by matching boundary temperatures. If u(0,t) = w(0) = uA and u(L,t) = w(L) = uB we have the following results form equation [37].

[pic] [39]

[pic] [40]

Solving equation [40] gives C1 as follows.

[pic] [41]

Substituting the results of equations [40] and [41] for C1 and C2 into equation [37] gives.

[pic] [42]

A rearrangement gives

[pic] [42]

The solution for v(x,t) which satisfies the diffusion equation with zero boundary conditions has already been found; it is given by equation [12].

[pic] [12]

The coefficients Cn are found by the following formula for eigenfunction expansion coefficients. This is the same formula as equation [15] with f(x) replaced by f(x) – w(x).

[pic] [43]

Substituting equation [42] for w(x) gives the equation necessary to compute the expansion coefficients, Cm.

[pic] [44]

In this problem we are given uA = uB = 0 which makes only a slight simplification in the final result.

-----------------------

[1] The choice of –λ2 for the constant as opposed to just plain λ comes from experience. Choosing the constant to have this form now gives a more convenient result later. If we chose the constant to be simply λ, we would obtain the same result, but the expression of the constant would be awkward.

[2] As usual, you can confirm that this solution satisfies the differential equation by substituting the solution into the differential equation.

[3] Using a standard integral table, and the fact that the sine of zero and the sine of an integral multiple of π are zero, we find the following result:

[pic]

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download