Homework 9: Implicit Solution of a Transient Conduction Problem

Part A

Time to move closer to real world applications.  Convert your Homework 8 program to use  a banded direct method such as dgbfa/dgbsl  to provide solution of the system of linear equations rather than the dgefa/dgesl full matrix solver.  Read the comments in dgbfa  carefully to understand how to construct your sparse data storage for the matrix.   You can also consult my solution to HW3 for an example of implementation of dgbfa and dgbsl.  (In case you haven't figured it out yet, with the right selection of variable values in DataM.f90, my HW4 sample is fully capable of transient calculations.)

Part B

Now upgrade your program from homework 8, so that it will run the same problem with Fully Explicit, Fully Implicit, or Crank-Nicholson methods.  As a first step in the homework, derive the stability limit for the 2-D fully explicit method, using a von Neumann (Fourier) approach.  Provide me with your derivation and install lines in your program that will not let the time step exceed this limit if a fully explicit method has been selected.  To demonstrate that your program is running correctly do the following:

  1. Run all three methods on the 3x3 grid for a 5 second transient and  a time step of 0.0001 seconds.  Plot the time history of the central volume's temperature for all three methods on Figure 1 of your report.  The curves should overlay.
  2. Make four more runs for the fully implicit method with time step sizes of 0.05, 0.1, 0.2, and 0.4 seconds..  Plot the time histories together with the original 0.0001 case on Figure 2 of your report.  Limit your temperature axis between 300K and 650K.
  3. Make four more runs for the Crank-Nicholson method with time step sizes of 0.05, 0.1, 0.2, and 0.4 seconds.  Plot these time histories together with the original 0.0001 case on Figure 3. Again, limit your temperature axis between 300K and 650K.
  4. For the explicit method make runs to 0.8 second with time steps of 0.001, 0.002, and 0.004 seconds on a 9x9 mesh (the stability limit is about 4.5e-3).  At 0.8 seconds apply Richard Extrapolation analysis of time discretization error to all points on your mesh, but do not include any points for which you would generate the logrithm of a negative number when calculating the order of accuracy.  Provide me with the estimates of mean absolute error and average order of accuracy associated with the lowest time step size.  Also provide me with minimum and maximum absolute error and minimum and maximum order of accuracy from all valid values calculated on the mesh.  If you extract data to a spreadsheet or other analysis package be certain to carry enough significant figures to do a meaningful error analysis.
  5. Repeat the analysis in the previous item using a fully implicit method and timesteps of 0.01, 0.02, and 0.04
  6. Repeat the item 4 analysis using a fully implicit  method and time steps of 0.1, 0.2, and 0.4
  7. Repeat the item 4 analysis using a Crank Nicholson method and time steps of 0.01, 0.02, and 0.04
  8. Repeat the item 4 analysis using a Crank Nicholson method and time steps of 0.1, 0.2, and 0.4


Sample Solution:

Plots:  crankNic.pdf and implicit.pdf