




























Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Course title Seminar in Engineering Analysis. Analytic and numerical methods applied to the solution of engineering problems at an advanced level. Solution methods are demonstrated on a wide range of engineering topics, including structures, fluids, thermal, thermal energy transport, and mechanical systems.
Typology: Study notes
1 / 36
This page cannot be seen from the preview
Don't miss anything!





























The diffusion equation describes the diffusion of species or energy starting at an initial time, with an initial spatial distribution and progressing over time. The simplest example has one space dimension in addition to time. In this example, time, t, and distance, x, are the independent variables. The solution is obtained for t 0 and 0 x xmax. The dependent variable, species or temperature, is identified as u(x,t) in the equations below. The basic diffusion equation is written as follows.
2
2
x
u t
u
Here, , is a species or 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 initial condition is usually written as follows:
u(x,0) = u 0 (x) [2]
The boundary conditions at x = 0 and x = xmax, can be functions of time in general. These are written as follows (with subscripts L and R denoting left and right).
u(0,t) = uL(t) u(xmax,t) = uR(t) [3]
There is actually a discontinuity in the initial conditions and the boundary conditions. The initial conditions are assumed to hold for all times less than or equal to zero. The boundary conditions hold for t > 0. This discontinuity only occurs at the boundaries. For example, the initial condition may be that u = 0 for all x and the boundary conditions may be that u = 100 at both boundaries for t > 0. This creates a discontinuity at t = 0 for x = 0 and x = xmax. Fortunately, this discontinuity does not affect our ability to solve the problem.
For the simplest cases we will consider both the initial and boundary conditions to be constant. However, in the most general case the boundary conditions may be given by any functional forms for u 0 (x). uL(t), and uR(t)
The solution of the diffusion proceeds by a method known as the separation of variables. In this method we postulate a solution that 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].
2
2 2
2 2
x
X x Tt x
X xTt x
u
t
Tt X x t
X xTt t
u
If we divide the final equation through by the product X(x)T(t), we obtain the following result.
2
x
X x t X x
Tt T t
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 simplify the solution, we choose the constant^1 to be equal to^2. This gives us two ordinary differential equations to solve.
2 2
2 2 ( ) ( )
dx
d X x dt X x
dTt Tt
The first equation becomes ()
Tt dt
dTt
known to be^2 Tt Ae ^ t
2 (^) ( ) . The second differential equation in [7] can be written as
2
2 X x dx
d X x
Ccos(x). Thus, our general solution for u(x,t) = X(x)T(t) becomes
( , ) () ( ) sin( ) cos( ) 1 sin( ) 2 cos( )
2 2
We first consider a simple set of boundary conditions where u(0,t) = 0 and u(xmax,t) = 0. In order to satisfy these boundary conditions on u for all times, we have to have boundary conditions on X(x) that X(0) = X(xmax) = 0. With these homogenous boundary conditions on X(x) and the form of the original differential equation for X(x), the solution
(^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.
integrals of the sine. If we multiply both sides of equation [13] by sin(mx/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 – the upper and lower limits of the Sturm-Liouville problem for X(x) – we get the following result.
max max
max max
0 max
2 1 0 max max
0 max 0 1 max max
0
sin sin sin
( )sin sin sin
x m n
x n
x
n
n
x
dx x
m x dx C x
n x x
m x C
dx x
n x x
m x dx C x
m x u x
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 except for the one where n = m, because the eigenfunctions, sin(nx/xmax), are orthogonal. This leaves only one integral to evaluate. Solving for Cm and evaluating the last integral^4 in equation [14] gives the following result.
max max
max
0 max
0 max 0 max
2
0 max
0 ( ) sin
sin
( )sin x
x
x
m (^) x dx
m x u x x dx x
m x
dx x
m x u x C
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.
As an example, consider the simplest case where u 0 (x) = U 0 , a constant. In this case we find Cm from the usual equation for the integral of sin(ax).
max max max
max 0
max max
0 max 0 max
0 0 max
0 max
cos
sin
( )sin
x x x m x
m x m
x x
dx x
m x x
dx x
m x u x x
m m
x
m m
x x
m x m
x x
Cm 1 cos
cos cos
max
max max
max max max
(^4) 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:
sin 2 4
sin 2 4
sin max max
max max max max 0
max 0 max
2
max max x x
m x m
x x x
m x m
x x dx x
m x
x x
^
Since cos(m) is 1 when m is even and – 1 when m is odd, the factor of 1 – cos(m) is zero when m is even and 2 when m is odd. This gives the final result for Cm shown below.
m even
modd m
Cm 0
Substituting this result into equation [12] gives the following solution to the diffusion equation when u 0 (x) = U 0 , a constant, and the boundary temperatures are zero.
1 , 3 , (^5) max
( ,)^40 sin( )
2
x
n with n
U e x u xt n n
n
(^)
If we substitute the equation for n into the summation terms we get the following result.
1 , 3 , 5
0 max
sin 4 ( ,)
(^2) max
22
n
x
n t
n
x
n x e U uxt
Equation [19] shows that the ratio u(x,t)/U 0 depends on the dimensionless parameters t/(xmax)^2 and x/xmax. Thus we can compute a universal plot of the solution as shown in Figure 1. Because the solution is symmetric about the midpoint x/L = 0.5, only the left half is plotted here.
Although we can obtain reasonable solutions using equation [19] for “large” times, we will not get a correct solution, within the limits of practical computers, if we try to evaluate equation [19] for times near zero. Figure 2 shows the results of using equation [19] to compute u(x,0)/U 0 using various numbers of terms for t = 0. In this case we are trying to use the sine series to evaluate X(x) = 1. We see that taking more terms makes the solution more accurate away from x = 0, but no matter how many terms we take, there are still strong oscillations near x = 0. This result is well known and is called the Gibbs phenomena.
The use of zero boundary conditions to give us a Sturm-Liouville problem appears to be a major limitation on this approach. We can solve a broader class of problems, still using the same approach, by appropriate definitions of variables.
v(x,t) = u(x,t) – uB [20]
This transformation gives v(0) = u(0,t) – uB = uB – uB = 0 and v(xmax) = u(xmax,t) – uB = uB
Since uB is a constant, when we substitute equation [20] into equation [1], we get the same equation in terms of v instead of u. Since the boundary conditions on v are that v = 0 at x = 0 and x = xmax, the differential equation for u is that of equation [1], and boundary conditions for v are the same as those used to derive the solution for u in equation [12]. Accordingly, equation [12] is the solution for v(x,t). When we apply the initial condition that u(x,0) = u 0 (x) to equation [12] for this revised problem, we get the following equation in place of equation [13].
(^1) max
0 (^ ) ( ,^0 ) ( ,^0 ) sin n
B n x
n x u x ux vx u C [21]
Starting with this equation, we can repeat the derivation of equation [15] for Cm, which gives the following result for this case where both x = 0 and x = xmax boundaries have the same common boundary value, uB.
^ ^
max max
0 max
0 0 max max
0 max
cos 1 ( )sin
( ) sin
x B
x m B m
m dx u x
m x u x x
dx x
m x u x u x
Thus, the solution for u(x,t) for the general initial condition, u(x,0) = u 0 (x) and the boundary conditions that u(0,t) = u(xmax,t) = uB is given by the following equation.
^ ^ ^
(^1 )
0 max
( ) sin sin( )
max 2
n
n
t
x u xt uB (^) x u x uB nxdx e x
We can extend the procedure used above to consider the case where the two boundaries are set at two different constant temperatures. That problem is specified by the following equations, where uL and uR are constants.
with ux u x u t uL u x t uR x
u t
u
2 (^ ,^0 )^0 ( ), (^0 ,) , ( max,)
2
The solution to this problem proceeds by defining u(x,t) as the sum of two functions v(x,t) and w(x), both of which satisfy the differential equation in [24].
u(x,t) = v(x,t) + w(x) [25]
Again, the function v(x,t) is defined to have a zero value at x = 0 and x = xmax, to give us the Sturm-Liouville problem, which provides us with a complete set of eigenfunctions. This zero boundary condition on v(x,t) tells us what the boundary conditions on w(x) must be set to satisfy the boundary conditions in equation [24].
w(0) = u(0,t) – v(0,t) = uL – 0 = uL w(xmax) = u(xmax,t) – v(xmax,t) = uR – 0 = uR [25a]
Setting w(0) = uL and w(xmax) = uR means that the solution for u(x,t) proposed in equation [25] satisfies the boundary conditions of equation [24]. Of course the proposed solution must also satisfy the differential equation. Since both v(x,t) and w(x) satisfy the differential equation in [24], their sum, u(x,t) = v(x,t) + w(x), which is the solution that we seek, will also satisfy the differential equation. We then examine the solution to the differential equations for both v(x,t) and w(x).
Since w(x) does not depend on t, we have the following result when we substitute w(x) into the differential equation.
2 2
2
dx
d w x
w t
w
That is, we have a simple, second-order, ordinary differential equation to solve for w.^5 Integrating this equation twice gives the following general solution for w(x).
w(x) = C 1 x + C 2 [27]
Since this equation has to satisfy the boundary conditions that w(0) = uL and w(xmax) = uR, we have the following solution for w(x) that satisfies the differential equation and boundary conditions.
x x
u u w x uL R L max
Since v is defined to be zero at both boundaries, the solution for v(x,t) is the same as the solution originally found for u(x,t) in equation [12]. Thus the solution for u(x,t) = v(x,t) + w(x) is found by adding equation [12] for v(x,t) and equation [28] for w(x) to give.
x x
u u u x
n x u xt v xt wx Ce L R L n
t n n 1 max^ max
( ,) ( , ) ( ) sin
[29a]
Setting t = 0 in this equation gives the initial condition for u(x,t).
x x
u u u x
n x u x v x w C L R L n
n 1 max^ max
( , 0 ) ( , 0 ) ( 0 ) sin
[29b]
Applying the general initial condition that u(x,0) = u 0 (x) gives the following result.
(^5) The w function does not depend on time. If the solution for v(x,t) is such that this function approaches zero for long times, then u(x,t) = v(x,t) + w(x) becomes equal to w(x). Thus, w(x) represents the steady state solution and we can regard the definition of u(x,t) = v(x,t) + x(x,t) as dividing the solution into a transient part that decays to zero and a steady-state part.
0 :^1 ,^3 ,^5 , max
0 : 0 : max^2 ,^4 ,^6 , max
sin
sin
2 max
2 max
n
xn t L
R L
n
xn t L
R L L
L
x
n x e U u n
u u
x
n x e x n
x U u
u u U u
u xt u
In this equation the dimensionless ratio, [u(x,t) – uL] / (U 0 – uL) depends on the dimensionless distance x/xmax, dimensionless time, = αt/(xmax)^2 , and the ratio of initial and boundary conditions, (uR – uL) / (U 0 – uL). A plot of this result for (uR – uL) / (U 0 – uL) = 0.5 is shown in Figure 3. This plot shows that the steady-state solution is a linear profile with the expected boundary values for [u(x,t) – uL] / (U 0 – uL); these are zero at x = 0 and (uR – uL) / (U 0 – uL) = 0.5 at x = xmax. If we changed the value of the parameter (uR – uL) / (U 0 – uL) from 0.5 to some other value, the plot would be similar, but the location of the dimensionless parameter at x = xmax would shift.
Figure 3 – Solution of diffusion equation for nonzero boundaries
t = 0.0001t = 0.
t = 0.
t = 0.005t = 0. t = 0. t = 0. t = 0.
(^) t = 0.
t = 0.
t = 0.
t = 0.15 t = 0.
t = 0. t = 0.
0
0.
0.
0.
0.
0.
0.
0.
0.
0.
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x
(u(x,t) - u
)/(UL
0
- u
)L
(uR - uL) / (U 0 - uL) = 0. t = (time)/(xmax)^2
In a cylindrical coordinate system, 0 ≤ r ≤ R, the diffusion equation has the following form.
r
u r t r r
u
The most general initial and boundary conditions for the radial problem are
u(r,0) = u 0 (r) u/r|r=0,t = 0 u(R,t) = uR(t) [36]
We will first consider the conditions where uR is a constant. We will later see that we will have to use an expansion of Sturm-Liouville eigenfunctions to match the initial condition. Anticipating this, we will obtain a problem with zero boundary conditions by splitting u(r,t) into two parts, as we did in equation [20].
u(r,t) = v(r,t) + uR [37]
The new function, v(r,t) satisfies the differential equation in [35] and has v/r = 0 at r = 0 and v = 0 at r = R. Since uR is a constant, the proposed solution in [37] will satisfy the original differential equation and the boundary conditions from equations [35] and [36]. As before, we seek a solution for v(r,t) using separation of variables. We postulate a solution that is the product of two functions, T(t) a function of time only and P(r) a function of the radial coordinate, r, only. With this assumption, our solution becomes.
v(r,t) = P(r)T(t) [38]
We substitute equation [38] for u in equation [35]. Since P(r) is a function of r only and T(t) is a function of t only, we obtain the following result.
r
Pr r r r
T t r
PrTt r r r r
u r r r
t
T t Pr t
PrTt t
u
If we divide the final equation through by the product P(r)T(t), we obtain the following result.
1 () 2 ( )
r
Pr r t Pr r r
Tt Tt
The left hand side of equation [40] is a function of t only; the right hand side is a function of r only. The only way that this can be correct is if both sides equal a constant. As before, we choose the constant to be equal to^2. This gives us two ordinary differential equations to solve.
The first equation becomes ()
Tt dt
dTt
2 ( ) . The
second differential equation in [40] can be written as ( ) 0
rPr r
Pr r r
We still have to satisfy the initial condition that u(r,0) = u 0 (r). We can do this by using an eigenfunction expansion.
1
m
ur uR u r uR CmJ mr [45]
The values of Cm are found from the general equation for orthogonal eigenfunction expansions which includes a weighting function. Here the weighting function p(r) is equal to r. For any initial condition, u 0 (r), we find the values of Cm from the following equation. (In the second equation below, the integral in the denominator has been evaluated using integral tables and the fact that J 0 (mR) = 0.)
2
0
0 0
0
2 0
0
0 0
rJ r u r u dr
rJ r dr
rJ r u r u dr C m
R m R R m
R m R m
If we consider the simple case where u 0 (r) is a constant equal to U 0 , we have the following result for Cm.
2 1
0
1
2
0 2 1
2
0
0 0
( )
rJ r
U u
J R
rJ r U u dr C m
r R
m r
m R
m
R m R m
0 1 0
0 1
0 2 1
2
1 0
m m
R m m
R m
m
m R m J
U u RJ R
U u R J R
U u C
With this result for Cm our solution for u(r,t) with the constant initial condition is written as shown below after some algebra.
^
1
0 0 0 1 0
0 2
(^2 )
( )
m
m R R
t
m m
R (^) u R
r e J J
U u u rt
m
^ [48]
This equation shows the important variables are r/R and αt/R^2. We can rearrange this equation to get a dimensionless form, plotted in Figure 4.
1
0 0 0 0 1 0
2 (^2 )
( )
m
R m
t
R m m
R R
r e J U u J
u rt u m
^ [49]
Figure 4 – Solutions to radial diffusion with constant initial condition tau = 0. tau = 0.002tau = 0. tau = 0.
tau = 0. tau = 0.
tau = 0.
tau = 0.08tau = 0. tau = 0.
tau = 0.
tau = 0.
tau = 0. tau = 0. tau = 0. tau = 0.
tau = 0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
0.
1.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 r/R
[u(r,t) - u
] / (UR
0
- uR)
In Figure 4, the notation “tau” denotes the dimensionless time αt/R^2. The results in Figure 4 are similar to those in Figure 1 for the rectangular geometry. Here the plot from 0 to R covers half the cylinder and is equivalent to the plot in Figure 1 that covered only half of the one-dimensional region.
Boundary conditions on differential equations can be expressed in terms of the value of u at the boundary, the gradient of u at the boundary, or a combination of the two. The boundary condition which is an equation relating the value of u and its gradient at the boundary is called a mixed boundary condition.
In this example we consider a problem of one dimensional heat conduction. Here we will explicitly use the temperature, T, as the dependent variable (instead of u) and note that α in this case is called the thermal diffusivity. We consider the case where a slab, initially at a constant temperature, T 0 , is placed in a medium at a different temperature, T∞. Convective heat transfer occurs on each side such that the results in the slab are symmetric about the center. In such a case, we can solve for only half the geometry. If we assume that the slab has a thickness of 2L, where –L ≤ x ≤ L, we solve for a region between 0 and L using a symmetry boundary condition that the gradient is zero at the center line. The convection boundary condition at x = L simply states that the conduction heat flux in the solid at the boundary is the same as the convection heat flux leaving the boundary. This is written as follows.
Figure 5 – Eigenvalues for convection problem
-** -
0
5
10
15
20
25
0 2 4 6 8 10 12 14 16 18 20 x
Functions
Line Cotangent Intersections
In this problem, the eigenfunctions are cos(λmx), but the values of λm do not have a simple equation such as they did in previous problems with simpler boundary conditions. Instead the values of λm for this problem, similar to those for the eigenvalues for J 0 , are not evenly spaced and have to be found by numerical solutions. The first few values of λm for hL/k = 1 are^7 λ 1 = 0.860333589, λ 2 = 3.425618459, λ 3 = 6.437298179, λ 4 = 9.529334405, λ 5 = 12.64528722, λ 6 = 15.77128487, and λ 7 = 18.90240996. Changing the value of hL/k would change the slope of the line in Figure 5 and consequently change the intersection points. As m increases, the values of λm increase and the intersections come closer to an integer value of π where the cotangent goes to infinity. For example, λ 7 = 18.90240996 is close to 6 π = 18.85.
Even though the eigenvalues, λm, are not evenly spaced, the part of the separation-of- variables problem that was solved for the ξ dependence is a Sturm-Liouville problem, so we know that we can expand any function, representing the initial condition, in terms of these eigenfunctions.
As usual, the most general solution for Θ is a sum of all the eigenfunctions that satisfy the differential equations and the boundary conditions.
( , ) cos( ) cot( ) 1
2 m m m
m m k
hL C e m^
(^7) You can observe that these values of l are the intersections indicated by the red dots in Figure 5.
If, for the moment, we regard T 0 as some arbitrary measure of a variable initial temperature profile, we can use an eigenfunction expansion to fit any initial condition, Θ(ξ,0) = Θ 0 (ξ) = (T 0 (L) – T)/(T 0 - T∞).
( , 0 ) ( ) cos( ) cot( ) 1
0 m m m
m m k
hL C
We have the usual equation for the eigenfunction expansion, where we can find the integral in the denominator from integral tables.
m m m
m m
m
m m
d
d
d C
cos( )sin( )
2 ( )cos( )
cos( )
( )cos( )
1
0
0 1
0
2
1
0
0 [59]
In particular for a constant initial temperature, T 0 , Θ 0 = 1, and Cm is found as follows.
m m m
m m m m
m
m m
m m m
m m m
d C
cos( )sin( )
2 sin( ) cos( )sin( )
sin( ) 2
cos( )sin( )
2 cos( )
1
0
1
(^0) [60]
Substituting this expression for Cm into equation [57] gives our solution for the dimensionless temperature, Θ, assuming a constant initial temperature.
cot( ) cos( )sin( )
2 sin( ) cos( ) ( , ) 1
2 m m m (^) m m m
m m k
e m hL
[61]
The solution is plotted for hL/k = 1 in Figure 6 and for hL/k = 10 in Figure 7. Both solutions are different from the earlier ones in which the boundary temperature was set to a new value at zero time. Here, the boundary temperature decreases over time because of convection heat transfer with the environment. With hL/k = 10, the convection is ten times as strong as it is with hL/k = 1. This leads to a more rapid change in temperature and a greater temperature gradient (i.e., greater heat transfer) at the x = L boundary.
Consider the diffusion equation with zero gradient boundary conditions shown below.
0 0 , 0 ( ) 0 , ,
2
2
x t x L t
^
1
2
n
t L
n
n
The constants in equation [63] are given by the following equation.
L m
L
0 0
As t → , the exponential terms will vanish and the only term that will be left in the
L f x dx L
uxt C 0
Physically, this result tells us what happens when there is no flux into the region. (Recall that the flux of physical quantities described by the diffusion equation is proportional to the gradient; thus a zero gradient means no flux.) With no flux added to the region, the final potential at all points in the region is simply the average of the initial condition.
Nonzero gradients on two sides of a one-dimensional geometry
This section should be omitted at first reading. Go directly to the heading “Mixed boundary condition in a cylinder” on page 25. If we want to consider two nonzero gradients we have a problem that does not have a steady-state solution unless both gradients are the same. Even if both gradients are the same the usual notion of a solution in terms of two variables, u(x,t) = v(x,t) + w(x) will not have a unique solution. To see this, consider the following problem.
,
0 0 ,
2
2
L x t x Lt
A solution of the form u(x,t) = v(x,t) + w(x) where v(x,t) satisfies the diffusion equation with zero gradient boundary conditions and w(x) satisfies the equation d^2 w/dx^2 = 0 with the boundary conditions that dw/dx = g 0 at x = 0 and dw/dx = gL at x = L will satisfy the differential equation. However, the equation for w cannot satisfy its boundary conditions. To see this, we integrate the equation d^2 w/dx^2 = 0 two times to obtain the general solution w = A + Bx.
Taking the first derivative of this general solution gives du/dx = B. This shows that du/dx is a constant which must be the same at both boundaries so that B = g 0 = gL. However, after we find the value for B, we have no other boundary condition for the constant A. Thus, any solution u = A + g 0 x = A + gLx will satisfy the differential equation and the boundary conditions only when g 0 = gL.
This mathematical result parallels physical reasoning. In the diffusion equation derivation the gradient of u represents a flux. If we have a one dimensional region in which we have a different flux on two sides, there cannot be a steady solution. How can we find a solution to the problem posed in equation [65]?
Since this problem does not admit a steady state solution we can augment our usual division of u(x,t) into a function v(x,t) that satisfies the diffusion equation with zero boundary conditions and a function, w(x), that depends only on the space coordinate. We have to include another function, (t) that provides the long term solution. We write our proposed solution as follows.
u(x,t)=v(x,t)+ w(x) + (t) [66]
Here v satisfies the diffusion equation with zero gradient boundary conditions, and the
function w(x) has the following gradient boundary conditions: 0
g dx
dw x
and
L xL
g dx
dw
. We can show that the solution in equation [66] satisfies the original
gradient boundary conditions.
0 0 0 , 0 , 0 , 0 ,
x t x t x t x t
L L x Lt x Lt x Lt x Lt
, , , ,
The initial condition for v can de derived from the initial condition for u:
v ( x , 0 ) u ( x , 0 )( 0 ) w ( x ) f ( x )( 0 ) w ( x ) [
If we substitute equation [66] into the diffusion equation and note that w(x) is a function of x only and (t) is a function of time only, we obtain the following result.
2 2
2 2
2 2
2
Since v satisfies the diffusion equation, the v terms in the last expression cancel leaving the following relationship between and w.
2
2