EGN 5456 Homework 1 Fall 1998
Numerical Solution of Heat Conduction in a Bar

Problem Statement

Find the temperature distribution T in a bar of length $\ell$ with a heat conduction coefficient $\kappa$.The ends of the bar are kept at zero degrees for all time. However, initially the center of the bar is at a larger temperature $T_{\mbox{max}}$. You can assume that this initial temperature decreases linearly from the middle to zero at each end. Plot the temperature distribution T for times $\kappa t/\ell^2 =$ 0, 0.125, 0.25, 0.375, and 0.5 in a single graph against x.

Analytical Solution

Solve the problem analytically.

• Find the separation of variables solution as a Fourier sine series. Rewrite this solution as a complex Fourier series, just to see that it can be done.
• Use the Fourier sine series solution to create a function subprogram T_ana that takes x, t, $\kappa$, and $\ell$ as arguments and returns the exact temperature at that point and time. Some points to consider:
  • The first time this subprogram is used, it should read the value of $T_{\mbox{max}}$ from an input file. But do not keep reading it in.
  • Pay particular attention to the best place to terminate the Fourier series. Sum as many terms as needed to get the most accurate solution, but no more than that.
  • Avoid underflow, which may impact performance.
  • Since the series converges slowly for time t=0, you may want to do this time separately.
• Use your function in a computer program to plot the exact solution T(x) in a single graph for times $\kappa t/\ell^2 =$ 0, 0.125, 0.25, 0.375, and 0.5. Use a plotting software package of your choice.

Numerical Solution

Solve the problem numerically using a finite difference method.

• Discretize the bar into J+1 x-locations spaced $\Delta x = \ell/ I$ apart by defining

  $$x_j = {j\over J} \ell \mbox{\ for } j=0, 1, 2, \ldots, J.$$

• Similarly, restrict the times to N+1 discrete times spaced $\Delta t = t_{\mbox{max}}/N$ apart by defining

  $$t^n = {n\over N} t_{\hbox{max}} \mbox{\ for } n=0, 1, 2, \ldots, N.$$

  Here $t_{\mbox{max}}$ is the last time you want to compute.
• You need to determine all (J+1)(N+1) temperature values $T^n_j \equiv T(j\Delta x, n\Delta t)$and then plot the ones corresponding to the requested times. Do that by first initializing the initial values T⁰_j as given. (It is best to do this using your exact solution function; then you can compute a different initial condition using exactly the same program).
• Next compute the values for later times using the formulae

  $${T^{n+1}_j -T^n_j\over \Delta t} = \kappa {T^n_{j+1} - 2 T^n_j + T^n_{j-1}\over (\Delta x)^2}$$

  (the left hand side approximates $\partial T/ \partial t$, the right hand side $\partial^2 T/\partial x^2$)This formulae can be used first for n=0 to compute the values of T¹_j for all $j = 1, 2, \ldots, J-1$.The end values T¹₀ and T¹_J can be found from the given Dirichlet boundary condition. (Again it is neatest to do this using your exact solution function). Next, now that all Tⁿ_j values are known for n=1, the values of T²_j can be computed for all values of $j=0,\ldots,J$, and so on.
• You will need to check at each time level whether it needs to be output.
• Try a value for J such as 16 or 32.
• Try 3 different values for N. The first value should give (as accurately as possible) $\Delta t = 0.25 \Delta x^2/\kappa$,the second $\Delta t = 0.5 \Delta x^2/\kappa$,and the third $\Delta t = \Delta x^2/\kappa$.
• For each of the three cases, make a separate plot of T versus x for times $\kappa t/\ell^2 =$ 0, 0.125, 0.25, 0.375, and 0.5. In these graphs, present the exact solution by a continuous line and the computed values T_jⁿ by symbols. (You may want to use more than J points to plot the exact solution accurately. About 50 points is a good choice.)
• Comment on the differences in the results.
• Follow the recommended programming techniques. This includes reducing the storage needed for the temperature values significantly from a straightforward two-dimensional array temp(j,n).
  Click here