Summary
MatLab code was produced to model the heat flow through a thermal tile on a space shuttle during re-entry. The model predicted the one-dimensional temperature distribution with time. The outside boundary is determined by a Dirichlet boundary and the inside is a Neumann boundary condition. The results produced by the code were in line with expected results and allowed appropriate tile thicknesses to be determined ranging from . The assumptions, limitations and effects on the model are highlighted. Extensions were made to generalise the problem, while simultaneously improving user experience. Four finite difference methods were compared and Crank-Nicolson elected due to its convergence and second order approximation in both space and time.
Introduction
Experimental trials are often unfeasible early in the design stage. In the case of the space shuttle the tile thickness must be determined to protect the internal structure, as well as the crew, without adding unnecessary weight. Modelling the heat flow through the tile with a finite differencing method is necessary to have a thorough understanding of temperatures throughout the tile.
Heat flow through a solid is governed by Fourier’s equation , which in the one-dimensional case simplifies to . This parabolic PDE has a wide range of applications throughout engineering and finances, such as particle diffusion. It is therefore useful to have efficient and accurate finite difference approximations to this form of PDE.
The aim of the investigation was to produce MatLab code to accurately model the temperature throughout a tile with time employing different numerical methods and to consequently elect the most effective method. The code is to be used to determine the optimal thickness of the tile.
A method is stable when truncation and round-off errors decay as the method moves to the next step, assessed by Von Neumann stability analysis.
It is consistent if the errors become small by refining the mesh size, and consequently allowing the discretised operator to approach the original differential operator . In the limit this becomes:
Convergence occurs when the discretised solution approaches the exact solution of the PDE as the mesh is refined to in the limiting case:. Commonly demonstrated using Lax’s Equivalence Theorem.
[1] [2] [3]
Results
The stability, consistency and convergence of the four different methods was assessed for different time and spatial steps. This was achieved by plotting the solution against the step size. (Appendix 7: stabletx.m)
For forward differencing, the Neumann Stability parameter was plotted against space and time. For the solution to be stable this value must be below [5]. This is plotted in Figure 1 instead of final temperature, as forward differencing is a consistent method (Appendix 1: Consistency Proofs). There is a region of stability for large spatial steps and small time steps. Large time steps have less of an influence on as per the small coefficient .
Similarly, for Dufort-Frankel the final internal temperature was plotted against mesh size as the method is unconditionally stable, but inconsistent (Appendix 1: Consistency Proofs). The solution appears to be accurate however, for larger spatial steps, as the surface plateaus.
The backward differencing method can be shown to be both consistent and unconditionally stable (Appendix 1: Consistency Proofs). In Figure 3 the end temperature is plotted for various spatial and time steps. There appears to be a linear error in time and a less dominant quadratic error in space. The method is significantly more accurate for all mesh sizes compared to previous methods.
The Crank-Nicolson method is stable and consistent and therefore convergent (Appendix 1: Consistency Proofs). In Figure 4 the final internal temperature is plotted against mesh size. The method is insignificantly affected by large time steps and has a second order spatial error. It produced the most stable solution as there is comparatively little change in with mesh size.
A 2D representation for all four methods was produced to assess the accuracy of the solution as well as determine a range of acceptable time and spatial steps, shown in Figure 5 and Figure 6. It was assumed that a 1% error margin is sufficiently accurate. Step sizes were selected for Crank-Nicolson, by observing temperature stays constant for , which was reduced to . A 1% error is introduced for , which was reduced to . Reductions were made as a safety factor.
The code (Appendix 2: shuttle.m) was produced for four separate methods: Forward differencing, Dufort-Frankel, Backward differencing and Crank – Nicolson. Forward differencing and Dufort-Frankel are an explicit finite differencing method, while backward differencing and Crank-Nicolson are explicit methods.
Figure 7 shows the one-dimensional heat flow through tile 597, with 50mm thickness over an interval of 4000 seconds. The Crank-Nicolson method was used to produce the plot. The right boundary is the temperature profile of the external temperature. The internal temperature is modelled with a Neumann-boundary and reaches a maximum internal temperature of approximately 219oC. The code can also be run for all other tile numbers with any of the four methods.
The tile thickness can be determined by setting a maximum internal temperature. A temperature of 100oC is suggested to prevent loss of structural strength of aluminium [5]. Maximum internal temperature decreases with thickness and occurs later in time, as shown in the temperature profiles in Figure 8. (Appendix 8: thicknessplot.m)
The optimal tile thickness was produced for each tile, using the shooting method and setting a maximum temperature of 100oC. For this the Crank-Nicolson method was used for a period of 4000 seconds and the results are recorded in Table 1. (Appendix 4: thicknessopt.m)
Tile
468
480
502
590
597
711
730
850
Thickness
0.0885
0.0904
0.0842
0.0848
0.0870
0.0904
0.0890
0.0897
Table 1: Optimal Tile Thicknesses
Discussion
Observations
The results produced are as expected and very accurate given an appropriate and . The time derivative is usually the most significant source of error, arising from the first order approximations of the time derivative. The second order time approximation of Dufort-Frankel and Crank-Nicolson produces accurate results for an appropriate . This allows for an accurate approximation to be produced for a large time step, reducing computation.
Forward differencing is limited as it requires larger spatial steps to make the method stable, which sacrifices resolution of the mesh. Dufort-Frankel has an additional error arising from the introduction of the wave equation, which is the critical error . This means must grow slower than to minimise the effect. Backward differencing has a linear divergence in time and a slight parabolic divergence in space . To produce accurate solutions a small should be selected while can be moderately large to improve computational efficiency.
Crank-Nicolson is the best method as it is unconditionally stable and consistent implying via Lax’s equivalence theorem that the method converges. Additionally, it is a second order method in both space and time , implying smaller errors than with other methods. It also has the most consistent temperature output. The method for this application also appears to be unaffected by ‘spikiness’ due to spurious oscillations. It is largely unaffected by time step. The spatial step is the critical error, which can be reduced to improve accuracy, although this is not necessary as only a 0.5oC change was observed over a range of 1 to 5.5 mm. It is however necessary that both and are relatively small in the frame of reference to remain accurate [6].
Assumptions
Several assumptions were made throughout the modelling process to simplify the problem. Furthermore, logical assumptions were made to make the problem both more comprehensible and easier to use the graphical user interface (GUI). These simplifications; however, result in deviations in the model from the actual solution.
The real problem is three dimensional. The assumption was made that the problem can be approximated with a one-dimensional method. This was done as very little information is known about dimensions of tiles and temperature distributions across the entire surface of the shuttle. A one-dimensional model is sufficient and provides accurate results given the inherent limitations of the model. A higher dimensional model would not have added to the quality of the solution, as the edges of the tile would be modelled with a Neumann boundary, as a constant surface temperature would have to be assumed across the tile.
The outside boundary was assumed to be a known temperature after the shuttle has landed, whereas it would most likely be a Robin boundary, as the hot tile will transfer heat to the surrounding air by means of radiation and convection. This however, was not included in the model, as little information is known about the boundary. Additionally, this would have added complexity in the noncritical region, as heat is no longer being added.
The internal boundary was modelled with a Neumann boundary, which implies perfect isolation and no heat flow across the boundary. This is unrealistic, as the aluminium structure would increase in temperature as it is in contact with a heat source. This means the internal temperature and consequently optimal tile thickness were overestimated.
The material is assumed to be eutectic. Inclusions or changes cannot me assessed in the model. This is a reasonable assumption, as the tiles are manufactured to a high level of precision. Therefore, any significantly anisotropic tiles will most likely be discarded.
Extensions
Several extensions were made to the task, to improve the user experience and allow a more general problem to be solved. The first extension performed was an automatic scan function was added to determine the temperatures on the outside boundary (Appendix 3: plottemp_exp.m). This provides more data points for the model and reduces the workload of the user. Furthermore, the code was generalised to allow any tile to be scanned.
To optimise tile thickness, the shooting method was employed to find a tile thickness where the maximum internal temperature does not exceed a pre-specified value with a default value of 100oC. This produces a precise result without having to resort to a graph. (Appendix 4: thicknessopt.m)
As the external boundary may be jagged a smoothing function was introduced to make the temperature profile smoother. If unsmoothed data is used unnatural noise is introduced to the solution, as small errors in the initial data results in oscillatory behaviour in mesh analysis. Smoothing allows the profile to appear more realistic and removes significant oscillations from the stability plot, such that only oscillations arising from the method itself are observed.
Additionally, for a more in-depth analysis of the effects on the solution by time and spatial steps surface plots were produced of the solution for mesh size. This allows regions of convergence to be determined visually, meaning appropriate time and spatial steps can be selected for each method. This also allows the user to explore the errors associated with each method, for example a first order approximation will show a linear divergence from the exact solution whereas a second order approximation will show a less dramatic parabolic divergence.
To further improve the user experience, a GUI was produced (Appendix 9: shuttlerunner.mlapp). The user can select any tile and method on the input tab. Furthermore, the user can elect to either determine the temperature profile for given parameters in the ‘Standard Run’ tab, optimise tile thickness in ‘Thickness Optimisation’ and explore stability of the methods in the ‘Stability’ tab. All graphs are displayed in the ‘Output’ tab. When the forward method is used, a stability panel appears which determines whether the input parameters provide a stable solution and suggests new time and space steps to make the model stable.
Figure 9: GUI demonstration
Additional extensions were made to increase the robustness of the code to account for user error. One such addition is the error message which appears if no tile is selected and an attempt is made to run the code. Furthermore, a warning appears, when a computationally expensive function is being run to alarm the user to wait.
Conclusion
Crank-Nicolson was determined to be the most accurate method due to the second order approximations in space and time. It is also convergent suggesting that the smaller the mesh the more accurate the approximation. It was observed that the method is relatively unaffected by time step which can be increased to reduce the computational load.
The tile thickness can be optimised by setting a maximum internal temperature and using the shooting method. The assumptions made in the model were largely accurate besides the inner boundary condition, which overestimates the temperature, as it assumes no heat flow. Tile thickness was overestimated which may be used as a safety factor in the design process to account for inclusions, anisotropic properties or minor damage.