A primary objective of this current work was to develop and test a C++ implementation of Milthorpe's algorithm. The practical reality of this software development is that the development and testing goes in cycles. Some computer code is developed, then it is tested. The result of the testing might show that there was some error in the code development so another cycle of development and testing occurs. Even if the testing shows that the code is working, the next cycle starts with development of the next feature of the code. It is more practical and efficient to develop a complex computer program in stages, adding and testing each feature or section one at a time.
Much of the time spent on this current work was in this development cycle. Features were added one at a time and there were several times where a feature that was initially thought to be working had to be re-developed and re-tested. In a fluid dynamics simulation such as this, it is quite difficult to isolate which calculation might be producing an erroneous result since most calculations are very much interdependent.
There was a reasonably steep learning curve that needed to be tackled to implement Milthorpe's algorithm in the new language. Once CONVECT was written and working to a certain point, however, it was relatively quick to implement further enhancements or bug-fixes. The `building block' approach that an object-oriented programming language offers, meant that rarely did changes need to be made to existing test object classes.
The tuning algorithm described in Section 5.6 was applied to a one dimensional version of the shock tube problem. CONVECT was set to run on the driver gas contamination problem described in Section 3.4, except that only one gas species was used - the test gas. This shock tube problem involved the propagation of just the primary shock. The inlet conditions were set so that they matched the post shock conditions that resulted in no secondary shocks or expansion waves being present. This should theoretically produce a shock that appears as a step function in all flow variables. An error function, based on the theoretical shock location, was minimised as the upwinding parameters were varied. Figures 6.1 and 6.2 show the results before and after a series of tuning runs. Flow quantities have been normalised to show agreement with the ideal step function shape of the shock wave being simulated. The theoretical shock location is also shown in both diagrams. It can be seen that the tuning vastly improved the results, showing greater agreement with the theoretical shock.
Some experimentation was done with the number of cells in the tuning runs. Not surprisingly, it was found that agreement with the theoretical shock improved as the number of cells was increased.
The original upwinding coefficients are shown in Figure 6.3. The final coefficients produced by the tuning process are shown in Figure 6.4.
To verify that the CONVECT code was indeed producing the same results as Milthorpe's original Fortran computer program, the calculation presented by Mallinson and Milthorpe [18] was reproduced (Figure 6.5). Mallinson and Milthorpe were investigating hypersonic laminar boundary layer flow both experimentally and numerically. This included a comparison of heat transfer and pressure data obtained in a free piston shock tunnel to the numerical results produced by Milthorpe's computer program. They showed that there was good agreement between their experiment and the numerical simulation. This was a good test case for demonstrating that CONVECT does indeed implement Milthorpe's algorithm and, indirectly, that the code simulates flow conditions produced in a free piston shock tunnel.
The numerical problem was a calculation over a domain of 500
150
cells. The left hand boundary was set to an inlet condition at the
freestream flow values (Table 6.1). The lower boundary was
set to a symmetry condition for the left most 50 cells and the
remainder was set to a solid boundary condition. The upper and right
hand boundaries were set to the freestream outflow condition. The
results produced by CONVECT are shown in Figure 6.5.
To enable overlays of these results with the original results obtained by Mallinson and Milthorpe [18], the output was scanned (Figure 6.6). Mallinson and Milthorpe stated (pers. comm., October 1999) that their plot was produced with `Gnuplot' [13] and that the velocity contours were generated at automatically determined values of the velocity data. Gnuplot was also used to plot the results obtained by CONVECT (Figure 6.5) with the contours being generated automatically. If the contours of the results produced by CONVECT coincide with the Mallinson and Milthorpe plot, it can be concluded that both simulations were in fact producing the same results, even though we do not know the values at the contours. To show the coincidence of the contour lines the scanned plot was scaled, made transparent using an image manipulation package and overlaid on the CONVECT results. This is shown in Figure 6.7.
The results of running the driver gas contamination problem outlined in Section 3.4 with an early version of CONVECT are shown in this section. These results were first presented at the 21st International Symposium on Shock Waves [20].
At that time it was apparent that the boundary layer produced in the
simulation was much smaller than the theoretical prediction. Some of
the parameters that the CONVECT code accepts in its input file are
the Reynolds numbers for each gas species. This gets translated into
the base viscosity, 
, used in the Sutherland viscosity
calculation. In these early results the Reynolds number was adjusted
so that the boundary layer was a more realistic size.
There was also some concern about the effect of non square cells and
of round-off errors produced by small cell dimensions. In order to
disregard these possible effects, cell dimensions were set so that 
. This meant that the dimensions being simulated were
different from a real shock tube. Table 6.2 shows some of
the parameters used.
| Re driver gas (HeAr) | 427 |
| Re test gas (Air) | 2177 |
| grid dimensions | 500 100 |
| problem dimensions (m) | 75 25 |
Figures 6.8, 6.9 and 6.10 show a velocity contour plot of the problem before and after shock reflection off the end wall. There is quite a sharp shock, which shows the effectiveness of Milthorpe's method to deal with this kind of problem. The boundary layer after the shock is well defined, however there should also be a boundary layer in the driver gas starting at the inlet. It is interesting to note the interaction of the reflected shock with the boundary layer. Shock bifurcation does not occur, but there are changes in the shock angle as it interacts with the boundary layer.
The main reason for implementing a computer program in parallel is to
achieve faster execution. As already stated, the ideal `speedup' for a
parallel computer program is the number of processors it is running
on. In reality, the speedup is often less than that due mainly to
overheads in the parallel implementation. Table 6.3
shows some timings for a typical small problem using 150
50
cells using the parallel code.
| number of processors | time (min) | speedup |
| 2 | 252 | 1.02 |
| 3 | 193 | 1.34 |
| 4 | 140 | 1.85 |
| 7 | 94 | 2.75 |
From this speedup information, it is clear that there are significant overheads for running the parallel code, and these are only overcome by using enough processors. A better indicator of the overheads associated with the parallel implementation is the fact that it took the original single processor version of CONVECT only 133 minutes to run this same problem! It is possible that these overheads would be less dominant on problems with many more cells to calculate over, but there is no timing data for this at the present.
The driver gas contamination problem, whose results appear in Section
6.2, was re-run using the parallel version of CONVECT to
verify that the parallel code worked. This was done first on four Sun
Ultra1 workstations networked together. In that case the computational
domain consisted of 250
50 cells. Figure 6.11 shows
the solution after shock reflection. It can be compared with Figure
6.9, noting that the time at which this output was captured
was a little further into the calculation. The results show a boundary
layer growing from the inlet. This boundary layer appears in Figure
6.11 and not in Figure 6.9 because some bugs in the
diffusion code had been fixed by the time the parallel version of
CONVECT was implemented. Apart from the discrepancy in the boundary
layers, the results are much the same suggesting that the parallel
code is working. The main problem during development of the parallel
code was dealing with the information flow across the parallel
boundaries. The parallel boundaries in this case divide the domain
into four equal slices in the 
direction. No fluctuations or
other numerical artifacts appear in this result, so the parallel
boundaries must be passing information correctly.
To verify that the parallel code could handle different sized problems, the same problem was run on a larger computational domain of 500
100 cells. This was done on a Sun Enterprise 4500
symmetric multi-processor computer using four of its processors. A
contour plot of the calculated velocity is shown in Figure
6.12. It can be seen that changing the number of cells does not
significantly change the results. It does, however, greatly impact on
how long the calculation takes. Increasing the number of cells does
make the solution more accurate, however the smaller timestep required
means that it takes much longer to get to the same stage of the
simulation. The number of calculations and the amount of computer
memory the problem uses also increases with the number of cells,
adding to the calculation time.
A series of runs exploring the inclusion of the artificial viscosity
term (Equation 5.1), and other changes to the way viscosity
was calculated (as outlined in Section 5.5.9), were carried
out. Since the error term is negative, a positive term of the same
form was added. This positive term was multiplied by a factor, 
,
which was varied until the expected boundary layer thickness was
achieved. After some experimenting, a value of 
was found
to produce the best results. Shown here in this section are the
results of a calculation of the driver gas problem outlined in Section
3.4. The computational domain used here was 1500
50 cells. To simplify the problem a little, and to observe the
boundary layer development, just the first 1 meter of the shock tube
was simulated.
Figures 6.13 and 6.14 show the flow variable velocity and pressure after the boundary layer has developed and before the shock has reflected from the end wall. Figure 6.15 shows a combined contour plot of the mass fraction of the driver gas and the pressure (also shown in Figure 6.14). This also illustrates the distance between the primary shock and the contact surface. In these results, the contact surface is quite smeared, leaving only a small amount of uncontaminated test gas behind the primary shock.
|
Figures 6.16 and 6.17 show the flow values of
velocity and pressure at a time after the shock has reflected off the
end wall bringing the test gas at the end of the shock tube to
rest. These results were calculated on a 150
50 cell grid
due to time constraints.
Even though CONVECT was reproducing the results of Milthorpe's original code, there was some concern that there were inaccuracies in the shock tube driver gas contamination problem. A commercial computational fluid dynamics code, GASP [16], was used to provide a comparison. GASP is a very capable computer program, but it takes some time to learn how to use it and apply it to particular problems. Transient problems, such as the driver gas contamination problem, seemed particularly difficult to set up since the documentation only gave examples of solutions to steady state problems.
As a preliminary step towards a full solution of the driver gas contamination problem the following calculation was set up. The flow conditions were taken from Section 3.4. The rest of the options were as follows:
33 grid points.
The `mesh sequencing' feature of GASP was used to improve the accuracy of the calculation. Mesh sequencing involves carrying out the calculation on a coarse grid first then interpolating the results up to a finer grid. In this test case, two levels of grid coarsening were used. The results of this calculation are shown in Figures 6.18 and 6.19. The calculation was completed in only a matter of minutes on a Pentium II 450 running Linux due to the small grid size.
With a successful one dimensional calculation completed, an attempt was made to fully simulate the driver gas contamination problem using GASP. The following options were used:
113 grid points on a grid that becomes finer as
it approaches the wall
This does not seem to be a good simulation of the problem since the level of the velocity (and other flow variables not presented) should have stayed at the inlet condition up to the shock - except in the boundary layer. There seems to be a peak in the value of velocity at the shock near the wall. Figure 6.21 shows that there is a dip in the mass fraction of driver gas behind the contact surface. This result indicates that some air has made its way behind the driver gas. Figure 6.22 shows the same plot as Figure 6.21 with the view point at the downstream end of the shock tube. The mass fraction of the driver gas drops off near the wall. This could be interpreted as a vortex forming where the shock meets the wall, pulling the air into the boundary layer and leaving it behind the contact surface.
)
and pressure, artificial viscosity
