Problem Definition

This assignment explores solution of the heat conductance equation (Eqn. 1), to model the temperature distribution along a thin, solid, homogeneous bar, which is insulated along its length.

=                                                                                    (1)

The bar is of length L = 100cm, and is made from aluminium giving a value of κ = 0.835cm2s − 1  – as depicted in Fig. 1(a).  During a manufacturing process, the entire bar is heated to 500◦ C, with both ends then plunged into ice (which can be assumed to remain at 0◦ C throughout the process). The manufacturing team need to know the temperature distribution along the length of the bar after a period of 10minutes (600s).  The target of this assignment is to develop a numerical solution to the given problem, and assess its performance.

The assignment explores solution of this partial differential equation using a numerical approach based on finite difference approximations of derivatives.  It is suggested to think of the problem on a two-dimensional grid (as shown in Fig. 1(b)): with location along the length of the bar x, on the horizontal axis; and time t on the vertical axis.  A spatial grid with four internal nodes will be used, with nodes spaced evenly along the length of the bar with inter-node distance ∆x = 20cm. The time grid will be formed with a step size of ∆t = 100s. The filled nodes

(a)

(b)

Figure 1: Problem description: (a) schematic of problem domain including auxiliary conditions (b) finite difference grid over which a solution will be evaluated.

represent locations where auxiliary conditions are applied, while the unfilled nodes represent locations at which a numerical solution will be evaluated. The assignment explores two solution strategies: in Part 1, an explicit solution method will be investigated; while in Part 2, the implicit Crank-Nicholson method is implemented. The analytical solution to this problem can be evaluated via Eqn. 2 (noting that n takes odd values only in the summation). This will help provide a reference solution against which to compare the numerical solutions developed within the assignment.

T(x,t) =    sin exp − t

for n = 1, 3, 5,..., ∞                                  (2)

Report Submission (5 marks)

You should submit a PDF report demonstrating your solutions to the tasks below.  Effort should be spent pro- portional to the marks awarded for each task. An additional 5 marks will be awarded based on the quality of the report, giving an overall total of 35 marks for the assignment. The report should contain concise technical answers expressing your findings and reasoning on the given tasks, and include high-quality figures plotting numerical so- lutions as instructed.  The 5 additional marks will be awarded according to the criteria in Tab. 1.  Note that no code is required to be submitted with this assignment, therefore the report should contain all information pertaining to solution of the tasks. The report should be no more than 8 pages in length, however references in a bibliography are not included in this limit.

Part 1: Explicit Solution of Heat Conductance Equation (10 marks)

Solution of the parabolic heat equation will first be explored using an explicit finite difference method.  Finite difference approximations will be substituted into Eqn. 1, and used to develop an update scheme for each interior node of the bar.  The explicit nature of this update scheme means it is possible to compute the solution at each node at a time in the future based on the present value of the solution at the node and its neighbours.

1. Construct the finite difference equations: use a central difference approximation for the second order spatial derivative, and forward difference approximation for the first order time derivative. The report should contain a description of the update equations using the same naming conventions as in the computational molecule from Point 2.

2. Draw the computational molecule for a generic interior node based on the explicit finite difference scheme. A copy of the diagram should be included in the final report, along with a brief description.

3. Implement the finite difference solution scheme in software (language of your choice) and solve for temperature T(x,t), along the length of the bar, over the range 0 < t ≤ 600s, using a grid with spacing ∆x = 20cm, ∆t = 100s (as seen in Fig. 1(b)). Discuss the main steps of the implementation in the report, including how boundary conditions are applied for the explicit solution.

4. Create two plots displaying model results:  (1) in the same graph, plot the temperature distribution along the bar at 100s intervals; and (2) on a separate graph plot the evolution of temperature with time at x = 20cm, including the analytical solution to enable calculation of the error. Additionally tabulate the resulting temperature at x = 20cm, when t = 600s. Analyse and discuss these results in the report text.

5. Adjust the solver to use ∆t = 50s, keeping ∆x = 20cm, and re-run the model and capture results.  Now adjust the model such that ∆t = 100s, but ∆x = 10cm, recording results as before. Add results to the table and second plot from Point 4, analysing how the model responds to these changes, and discuss findings in the report text.

Part 2: Implicit Solution of Heat Conductance Equation (10 marks)

Explicit solution methods can encounter problems related to stability and accuracy.  Therefore the second part of this assignment will implement and analyse an implicit numerical solution to the heat conductance equation (Eqn. 1), in order to understand how accuracy and stability can be improved, at the expense of computational effort. Here the Crank Nicholson method will be employed, which evaluates the finite difference approximations at the mid-point of the time step, producing a solution which is second order accurate in both space and time.  The drawback of this method is that unlike the explicit method from Part 1, the solution at a particular location can no longer be evaluated at a future point in time using known information. Instead a system of equations is assembled, and subsequently solved simultaneously via the application of boundary conditions.

1. Construct the finite difference equations using the Crank-Nicholson method. This constructs approximations of the derivatives at the step midpoint (see Week 8 notes for further details).  The report should contain a concise description of the system of equations, and how they can be solved to compute the temperature distribution at each node for a particular time.

2. Draw the computational molecule for a generic interior node based on the implicit Crank-Nicholson method. A copy of the diagram should be included in the final report, along with a brief description.

3. Implement the Crank-Nicholson solution method in software (using a language or your choice), and solve for temperature T(x,t), along the length of the bar, over the range 0 < t ≤ 600s, using a grid with spacing ∆x = 20cm, ∆t = 100s. Discuss the main steps of the implementation in the report, including how boundary conditions are applied for the implicit system of equations, and the chosen solution algorithm.

4. Create two plots displaying model results: (1) on the same graph, plot the temperature distribution along the length of the bar at 100s intervals; and (2) on a separate graph plot the evolution of temperature with time at x = 20cm, also plot the analytical solution to enable calculation of the error.  Additionally, tabulate the resulting temperature at x = 20cm, when t = 600s. Analyse and discuss these results in the report text

5. Adjust the implicit solver to use ∆t = 50s, keeping ∆x = 20cm, and re-run the model and capture results. Now adjust the model such that ∆t = 100s, but ∆x = 10cm, recording results as before. Add results to the table and second plot from Point 4, analysing how the model responds to these changes, and discuss findings in the report text

Part 3: Modelling the Initial Temperature Change (10 marks)

The manufacturing team are now interested in modelling the temperature of the bar, shortly after it is plunged into ice. It is therefore requested to provide accurate results for how the temperature of the bar varies over the first minute (60s) of cooling. Them team have also noticed that the initial temperature distribution is non-uniform, and instead T(0,x) is described by Eqn. 3.

T(0,x) = −xk(x − L) + c       k = 0.1, c = 400                                                 (3)

Develop a solution method (e.g. based on Part I or Part II), which provides an accurate and efficient solution for T(t,x) over the time period 0 < t ≤ 60s, with these updated initial conditions.  You should assume there is no analytical solution for this updated problem description, and therefore look for alternative methods to justify the accuracy of your proposed solution, and the trade-off against computational effort. Findings should be written up in the report, with appropriate plots/tables to support your proposed approach.