Developing skills for neutronic modelling of nuclear power reactors in South Africa

In recent years, due to economic and social infrastructure development and growth, South Africa has been facing growth in energy demand. Addressing this demand includes building more coal power stations, however with attention paid in designing them to reduce greenhouse gas emissions. A second response is to deploy more power sources using renewable and nuclear energy. The South African government has plans to add about 9.6 GW of nuclear energy to the electricity grid. Accepting that South Africa will seek an international vendor or vendors to supply nuclear plants, a certain degree of localisation of manufacture and operation should be planned. One localisation task that can be actively pursued is reactor analysis, including criticality, burnup, shielding and accident analysis of the reactor. Such development of expertise will support both economic and safety aspects of building and running a nuclear reactor. With this in mind, neutronic analysis of the VVER 1000 reactor was initiated. The government’s intention to build a new fleet of reactors means it is important that the VVER1000 reactor be included in studies done by the reactor analysis group at the School of Mechanical and Nuclear Engineering at the North-West University. The analysis was performed using MCNP6 for the cold zero power state at the beginning of cycle with the specifications obtained from the open literature. The input file was generated using the in-house code NWURCS. To ensure accuracy and precision of the results produced by the MCNP6 code, convergence studies of the MCNP6 models were carried out. Once a satisfactory model was obtained, the critical reactor state was calculated by adjusting the boron concentration in the water. Furthermore, the control rod worth, reactivity coefficients and beff were also calculated and are reported in this paper.


Introduction
In recent years, due to economic and social infrastructure development and growth, South Africa has been facing growth in energy demand (Winkler, 2006).Given that the major source of power is coal power stations, however this demand will result in the depletion of coal.Coupled to this problem are the greenhouse gas emissions caused by the use of coal, and this has led to two outcomes.Firstly, new coal-powered plants are designed and built with features to reduce such emissions (SA-Advocacy, 2015).Secondly, South Africa is interested in expanding its total energy mix to include more environmentally friendly power sources.Considering the finite coal resources, these alternative power sources must also be sustainable.One such power source is nuclear energy, which is believed to have a capability to supply base-load electricity.The South African government has plans to add about 9.6 GW of nuclear energy to the electricity grid, as stated in the Integrated Resource Plan for electricity, IRP2010 (Department of Energy (DOE), 2013).
Producing power from nuclear energy is not new to South Africa, which has successfully been operating the Koeberg nuclear power station near Cape Town for over 30 years.Koeberg consists of two French designed pressurised water reactors (PWRs), each of 930 MWe, and contributes about 6% of the country's electricity (Stott, 2013).A choice for South Africa is therefore to expand its current nuclear capacity with PWRs.
Currently, South Africa is unable to produce its own nuclear power plants, which requires years of experience, since both safety and economics must be considered in depth.Given this, it is accepted that South Africa will seek an international vendor or vendors, with a certain degree of localisation of manufacture and operation.The degree of localisation will in all likelihood depend upon agreements made between the various stakeholders, but certain aspects of localisation should be mandatory, with particular reference to safety considerations.Given the nature of nuclear power, it would be unwise to build a fleet of nuclear power plants without developing adequate skills, knowledge and know-how in terms of safety of nuclear power plants.In 2014 the government signed memoranda of understanding with several countries experienced in nuclear energy, with Russia being among the prospective countries that would assist in implementation of the new nuclear build.The VVER-1000 type reactor (pressurised water reactor of Russian design) is therefore chosen as a relevant reactor to investigate.
One localisation task that can be actively pursued is reactor analysis, including criticality, burnup, shielding and accident analysis of the reactor.These studies can be classified into two broad groups: experimental and plant operational data analysis, and calculation analysis.Whilst the first option would be the natural choice, it becomes prohibitively expensive to study all possible scenarios for any given reactor, and it has become standard to perform high fidelity calculations with computers to characterise the various scenarios.Experimental and plant data are then used to validate the calculation models.There are many ways to characterise computational studies.One way is to separate the calculation space into the reactor core and the balance of the plant.When the core is studied, one would first characterise the core in terms of the neutronics, fuel performance, thermal-hydraulics and burn-up of so called 'steady states'.Thereafter, transient states can be characterised, to model both operational transients and accident scenarios.The reactor analysis group at the School of Mechanical and Nuclear Engineering at the North-West University has initiated such studies, and neutronic calculations on a steady state fresh core of a chosen PWR, the VVER 1000, are reported in this paper.

Aspects of reactor physics studied
In a nuclear reactor, the nuclear fission reaction is the principal source of nuclear energy.It is the reaction between a neutron and a heavy nucleus (e.g. 235U, 239 Pu) and it results in the production of a highly unstable compound nucleus, which splits into fragments accompanied by the emission of a few neutrons and a large amount of energy.The most important purpose of the reactor operation is to achieve a controllable sustainable nuclear chain reaction.A nuclear chain reaction is a series of nuclear fissions, each initiated by a neutron released in a preceding fission.A sustainable chain reaction means that the ratio between the numbers of neutrons that are produced and the number of neutrons that are absorbed or leave the system through the outer boundary must remain constant as time progresses.
The criticality of the nuclear reactor is characterised by k eff (multiplication factor), which is an eigenvalue to the transport equation (Stacey, 2007).When k eff = 1, the system is considered critical.With k eff > 1, the system is supercritical and the number of fissions in the chain increases with time.With k eff < 1, the system is subcritical and the chain reaction will not be sustained, since the generation of neutrons will be terminated.The reactivity of the reactor core can be defined by Equation 1 (Stacey, 2007): and when the reactor is critical, r = 0.In a typical PWR the reactor is kept critical by first adding boron to the moderator.As the reactor cycle progresses from the beginning to the end of cycle, the concentration of boron in the moderator is gradually depleted to compensate for the depletion of the 235 U.In addition, the build-up of other actinides and fission products in the core also affects the rate of boron dilution.Control rods are also present in the reactor core, arranged into banks (clusters), and entire banks are moved when required.
One main purpose of the control rods is to shut down the reactor in an emergency, so it is important that they can do so from any operational or accident state of the reactor when the rods can still be inserted into the reactor core.This ability is measured by evaluating the shut-down margin, which, in terms of International Atomic Energy Agency regulations, should be larger than 5100 pcm for PWRs (IAEA, 2003).The shut-down margin is the measure of the control rod worth, which is the change in reactivity due to the control rod movement.The shut-down margin will differ depending on the state of the reactor.For example, the hot full power state will have a different value from the hot zero power state, due to the difference in temperatures of the two states.The control rod worth analysis is divided into two types: (a) the differential worth (pcm/cm), which is the reactivity change that is brought about by the movement of a rod, and (b) the integral rod worth (pcm) at a given position, which is the summation of the entire differential rod worth up to the point of withdrawal (Burn, 1988).The reactivity r of the system as defined above can change as a function of certain parameters, some important ones being the fuel temperature, moderator temperature, and void fraction of the moderator (which for a PWR is borated water).In a reactor transient involving an increase in power, the core temperatures will increase, together with the void fraction.In a properly designed reactor, however, the overall effect on the reactivity should be a decrease in reactivity with the increase in temperature, since this will assist in terminating an undesirable power excursion.The effect is normally analysed by the definition of the reactivity coefficient of a parameter as shown in Equation 2. The reactivity coefficients that are considered most for PWRs are the Doppler coefficient (DC) and the moderator temperature coefficient (MTC).The Doppler coefficient is associated with the fuel temperature.As an example of this effect, consider a low enriched 235 U core.The resonance absorption of 238 U will increase with temperature, and this decreases the reactivity.By defining DC as: (2) With Dr negative and DT ƒ positive, a negative reactivity coefficient results.
The MTC is the measure of the change in reactivity due to a change in the temperature of the coolant.When the temperature increases, the den-sity decreases, with an accompanying decrease in the boron concentration in the moderator for a typical PWR.If this is the only effect present, then the reactivity of the core will increase.A related effect is that of increasing the void fraction.This will also result in a decrease in the boron concentration for PWRs, with the potential for increasing the reactivity.There is also a decrease in the number of hydrogen isotopes with an increase in the void factor, however, so the opposing effect of the decrease in moderation will also be present.It is noted that the PWRs operate such that the heat transfer at the fuel pin wall above a certain height of the active core occurs in the nucleate boiling regime.Therefore, the vapour phase will exist, although localised at and near the wall surface, leading to the definition of a void fraction.The limits on reactivity feedback coefficients in all reactor states for PWRs for the DC and the MTC are as follows: -4.90 < DC < -2.90, and -70.0 < MTC < 0.00.
Another important parameter of the reactor is the delayed neutron fraction b eff .During the fission process, about 99.35% of the total neutrons will be produced immediately at the time that fission reaction occurs, and these are termed prompt neutrons.However, there will be some neutrons that are produced later than the prompt neutrons due to decay and reactions of the fission products.About twothirds of the delayed neutrons originates from lighter fission product nuclides, the remainder from heavier isotopes (Neeb, 1997).The fraction of the number of delayed neutrons to the total number of neutrons produced is defined as b eff .This quantity is important when analysing the time (transient) behaviour of the reactor.b eff is limited to a fraction 0.0067.For a VVER-1000, the limit is at 0.0074 according to the FSAR of the VVER-1000 (BNPP-1, cited in Jahanbin & Malmir, 2011).There are also other important parameters that must be considered from both safety and economic considerations, such as fuel burn-up.Since these are not presented in this paper and are being pursued in current studies, they will be discussed in appropriate later reports and papers.

Computational tools
In this paper, the computer codes MCNP6 (Monte Carlo N-Particle code, version 6) and NWURCS (North-West University Reactor Code Suite) were used for the neutronics simulations, and are briefly described in the following sections.

The MCNP6
The MCNP is a Monte-Carlo general-purpose continuous, generalised-geometry computer code designed to track many particle types over broad ranges of energies.The latest version of MCNP is MCNP6, which combines the capabilities of MCNP5 and MCNPX, together with new building capabilities (Pelowitz, 2013) The output from the code includes flux tallies and the multiplication eigenvalue k, (k eff or k inf ) .To perform the neutronics analysis, MCNP6 requires an input file that contains information that describes the specific geometry, the materials of the medium, a selection of the cross-section evaluation, the type of particles to be transported, the geometry of the source, and the type of tallies (such as fission power).In order for MCNP6 results to be applicable, source convergence, k convergence and tally convergence tests must be done initially.
MCNP6 is considered to be a high fidelity code, since accurate modelling of both geometry and materials are possible, together with a continuous energy treatment of the nuclear parameters.Models based on MCNP6 are therefore often used as benchmark models for benchmark activities (Perin, et al., 2015).Since no full core MCNP6 models of the VVER-1000 were found in the open literature, this study therefore helps to fill that gap in modelling the VVER-1000 reactor.

NWURCS
Computer codes have been developed that are used as 'code wrappers', meaning that they are used to generate input files for the main code being used and also to extract data from the output file generated by the main code (Nuttin, et al., 2013;Aghaie, et al., 2012).The NWURCS suite of codes is such a code, used to generate the MCNP6 input file used for neutronics simulations (Naicker, et al., 2015).The suite of codes can also be used to create steady state RELAP5 input files for the reactor core and for coupling between flux and temperature calculations, viz.coupling of MCNP6 and RELAP5.The suite of codes can also be used to extract data from an MCNP6 or RELAP5 output file, making the data analysis easier for the user.Because the code was recently developed, verification of the input files generated by the code for MCNP6 and RELAP5 are necessary.This paper reports on some of the verification tasks that have been completed on NWURCS.

The VVER reactor
The PWRS are light water reactors, used in the large majority of the world's nuclear power plants.There are about 437 PWR nuclear power currently in operation globally, with about 52 VVER-1000s among the operational PWRs (Katona, 2011).The VVER-1000 is a Russian design, generation III reactor producing 1000 MW electric power output and it incorporates international safety standards (Gidropress, 2008).A generation III reactor is a generation II reactor with evolutionary design improvements in the areas of fuel technology, modularised construction, safety systems, and standardised designs (Goldberg & Rosner, 2011).The design chosen for the study was of the VVER-1000 type.
A schematic of the primary loop of the VVER-1000 reactor, which contains four loops, is shown in Figure 1.The reactor has horizontal steam generators while other PWRs typically have vertical steam generators.One of the advantages of the horizontal steam generator is the moderate velocity of the medium in the secondary loop which prevents any danger of vibrations of the heat-exchanger tubes and damage by foreign objects.
The reactor core of the VVER-1000 consists of 163 hexagonal fuel assemblies (FAs) placed in a lattice of hexagonal symmetry shown in Figure 2.This makes it different from other PWR cores, which have square fuel assemblies, with the comparison being shown in Figure 3.
According to the benchmark definition of (Lotsch, et al., 2010) includes a reactor pressure vessel, which is the pressure boundary of the reactor core and its internals, as seen in Figure 2.
The fuel assembly (FA) as shown in Figure 4 contains 312 fuel rods, 18 guide tubes and one instrumentation tube.It is different in geometry from other PWR fuel assemblies.The fuel pin of the VVER-1000 reactor also has a different structure from those of other PWRs.It has a central hole in its fuel pellet as seen in Figure 5, which is filled with helium, its purpose being to provide lower centre temperatures and a free volume to allow any released fission gas to expand into and thus reduce internal pressure (IAEA, 2006).The fuel pin material is uranium dioxide (UO 2 ).The 235 U of the UO 2 has an enrichment of up to 4%.Other assemblies contain fuel pins that have 5% of gadolinium-oxide (Gd 2 O 3 ) integral burnable poison in their compositions (Lotsch, et al., 2010).The use of gadolinium (Gd) burnable poison allows for a reduction in the quantity of the initial boric acid concentration in the water.Low boric acid concentration helps to ensure a negative moderator temperature coefficient of reactivity (Allen, 2003).
Table 1 gives a summary of the VVER-1000 specification.These specifications were not available from a single source; they were formulated from Lotsch, et al. (2010), Pazirandeh, et al. (2011) and Kosourov, et al. (2003).

NWURCS verification
Development of the NWURCS in-house code was started in 2012.The code is used to generate the input files for MCNP as previously mentioned.This input file contains specification of the reactor geometry which includes the surfaces that defines the various parts of the reactor, the material in the reactor, together with references to the nuclear parameters contained in the nuclear data libraries, and general parameters, such as the number of source points per cycle.This NWURCS code is therefore useful because it reduces the amount of time it would take for an input file to be written and also helps to reduce human error.Verification of the input model is therefore crucial to ensure that what the user wants to model is exactly what is being generated by NWURCS.It is also important in terms of reporting any computer bugs present in the code.Verification was done in the ways discussed below.It is important to note that verification is an ongoing process, and that further studies are planned to continue with the verification effort.
The first method of verification was a rigorous inspection of the lines of the input file.try is constructed in terms of nested arrays.The fuel pins and guide tubes comprise an array which defines the fuel assembly which, in turn, comprises an array which defines the core.One needs further to add to the modelling of the heterogeneity in terms of fuel pins (different material compositions), temperature and burnup.This geometric and material definitions were, therefore, closely inspected to ensure consistent nesting and material assignment.
It is noted that a full core input contains about 300 000 lines of input data, so each line was not inspected, but the symmetry used in NWURCS (for example the nesting outlined above) allowed for samples of three from each type of array to be examined.Further details regarding the input data is given in Nyalunga (2016).
The second method of verification was by visual inspection of plots of the reactor core.This can be seen in Figure 6, where the positions of the fuel pins within the fuel assemblies and then the subsequent positioning of the fuel assemblies within the reactor core can be seen.Given that this nesting was modelled using arrays embedded within arrays, it was important to verify that the arrays were assigned properly.The axial definition of the reactor core was also visually inspected, as shown in Figure 7.One can see that the various components are positioned properly within the layout.The figures as mentioned verifies that what was defined by the user as the reactor system was produced.
Varying various input parameters to ensure that the physics behaviour was as expected allowed for a third verification method.Two such cases were the effect of boron concentration in the water, and the effect of control rod motion on the multiplication constant k.A fourth verification check was to ensure that a full core calculation produced the same multiplication constant as a 1/6 core calculation.The full core and 1/6 core calculations yielded values of k eff = 1.00353 ± 0.00003 and k eff = 1.00366 ± 0.00003.The slight difference (of 15pcm) as seen in the two results are expected due to the statistical nature of the MCNP calculation.
All of the above checks showed that NWURCS generated an input model according to the user's intent.

Convergence tests
Due to the statistical nature of the MCNP calculation, various convergence tests must be satisfied in order to use the results: source convergence, k convergence and tally convergence.
The statistical accuracy of the MCNP6 calculation is controlled by three parameters, viz., the number of source points, the total number of cycles to be completed, and the number of cycles to be skipped before beginning tally accumulation and statistically averaging the multiplication constant.At the start of the calculation, the user needs to define the fission source distribution.At the end of the first cycle, the neutrons that are produced as a result of the fission process are used as the source definition for the second cycle.This process continues as the cycles progress, and the source will be progressively spread over all the cells that contain fissile material.Statistical counting should only commence once the source is satisfactorily distributed over all fissile cells.This convergence of the source distribution is monitored by defining the Shannon entropy, as in Equation 3.
(3) Where P j is the number of source points for each cell or mesh j (Brown, 2006).When the source is well distributed, the Shannon entropy should be statistically constant as a function of the generation cycle.
One therefore finds that one has to perform a preliminary calculation monitoring the Shannon entropy to determine the number of cycles that must be skipped.The Shannon entropy is automatically produced during the MCNP calculation and the output indicates whether enough cycles have been skipped.It is also convenient to assess H src by plotting H src as a function of the cycle number.To test this convergence, three cases were run, with sources points of 400 000, 800 000 and 1000 000.
From Figure 8 it can be concluded that for the sources points of 400000, 800000 and 1000000, about 116, 123 and 122 cycles respectively, are required to converge the fission source.This means that at least that many cycles should be skipped before tallying any results.
Once the source has been converged, the k values calculated thereafter for each cycle is used to determine the average value of k together with its statistical error.It has been shown by Naicker, et al. (2016) (where the Monte Carlo code KENOVI instead of MCNP6 was used) that the product of the total number of active cycles counted and the number of fission source points determines the accuracy to which k can be calculated.One strategy to obtain the required accuracy in k would be to first determine the total number of source points and the number of skipped cycles required to obtain an adequately converged source.Thereafter, the required accuracy in k can be achieved by controlling the number of active cycles.
The final set of convergence tests that can be performed are on the tallies, and needs to be carried out only should any tally be required.This set of convergence testing places a heavier burden on the computational resources and calculation time since more cycles are required in order for the tests to be passed.There are 10 tally tests (or checks as referred to in the MCNP manual (Pelowitz, 2013)) that must be passed.A sample listing of these tally tests as obtained from a MCNP output file is shown in Figure 9, and a complete description of these tallies are given by Pelowitz.The statistics required for each tally check to pass are not the same, with tally check 10 being the most expensive.However, one can see the effect of tally convergence by studying  the power profile, as shown in Figure 10.Since the reactor was assumed to be at a constant temperature of 300 K, and with the fuel assemblies displaying axial symmetry, an axial symmetrical profile was expected.Therefore, if one intends to use the power profiles further (for example in thermal-hydraulic coupling calculations), a converged power distribution is required.

Reactor characterisation
The criticality of the system was reached without inserting any control rods in the reactor core, but by dilution of boric acid in the reactor core.Figure 11 shows the graph used to predict the critical boron by linear interpolation.The critical boric acid concentration was found to be equal to 1328.59 ppm at k eff = 1.00003 ± 0.00003 with the reactor temperature being 300 K.
The Doppler and moderator temperature reactivity coefficient (DC and MTC respectively) were examined with a single fuel assembly model considered.To calculate the DC, a criticality calculation was done by changing only the fuel temperature in all the fuel pins by 10 K.The reactivity change was   then calculated using Equation 2, where r was calculated using Equation 1.The DC was found to be -3.52 pcm/K.A similar calculation was done to calculate the MTC.The moderator temperature throughout the fuel assembly was changed by 5 K using the approximation in Equation 4.
(4) A value of 3.5 pcm/K was obtained.This positive value for the MTC can be of concern, since a negative MTC is desirable because of its self-regulating effect.A possible reason for the slightly positive value could be that the evaluation was carried out at 300 K.It is envisioned that further studies will involve a definition of the reactor temperatures in terms of coupled neutronic-thermal hydraulic calculations, and this would most probably yield a negative value for the MTC.
The next set of results presented is that of the control rod worth.A control rod in this reactor is composed of two axial sections, as shown in Figure 12.The material at the tip of the rod (i.e. the end that is inserted into the reactor) is composed of Dy 2 O 3 •TiO 2 of thickness 30 cm.The rest of the control rod is made up of B 4 C.The Dy 2 O 3 •TiO 2 is included because of its high dimensional and structural stability and it occupies the part of the control rod with the highest radiation dose (Risovany et al., 2000).To calculate the control rod worth and to also test the effect of Dy 2 O 3 •TiO 2 , two sets of calculations were done, the first with Dy 2 O 3 •TiO 2 and the second without Dy 2 O 3 •TiO 2 .The calculations were performed for a FA model.Equation 5 was used to determine the rod worth.

Rod worth =
(5) where x is the distance a control rod is moved, Dx is the rod displacement between two successive calculations in one set, and Dr is the corresponding change in the reactivity.The integral rod worth calculations are shown in Figure 13, where the integral rod worth is the cumulative value of the rod worth.For most of the insertion depth, the rod worths are the same for both B 4 C and B 4 C + Dy 2 O 3 •TiO 2 .However, between insertion fractional depth 0.7499 and 0.87499, the rod worth for the control rod with B 4 C is seen to be larger than the rod worth of the rod with B 4 C +Dy 2 O 3 •TiO 2 , with the maximum difference calculated to be 550.68 pcm at 0.87499 fractional insertion depth.Control rod worth curves for an EPR design are shown in Figure 14.The results were seen to be in good agreement with these results.
The final result to be shown is that of the effective delayed neutron fraction b eff .The delayed neutron fraction is calculated by switching the TOTNU card off in the MCNP6 input file.In this way only prompt neutrons are modelled.b eff is then calculated using Equation 6 (Michalek, et al., 2003): Where k p is the calculation using prompt neutrons only and k t is the criticality constant calculated with both prompt and delayed neutrons.With the values for k t and k p calculated as 1.00353 ± 0.00003 and 0.99657 ± 0.00003 respectively, the effective delayed neutron fraction was then obtained as 0.00694 ± 0.00004.This is 9.0% smaller than the value of 0.00761 listed in Jahanbin and Malmir (2011) and 6.6% smaller than the value of 0.0074 recommended in the FSAR of the VVER-1000 (BNPP-1, cited in Jahanbin and Malmir).However, it is larger than the value of 0.0067 listed in Stacey (2007).

Future work
Having successfully built an input model for the VVER 1000-type reactor, the next step will be to couple this model with a RELAP5 model.The RELAP5 calculation will characterise the temperature of the various materials in the reactor core, thus resulting in the hot full power state at the beginning of core being established.Thereafter, the timedependent progression of the core towards the end of core state will be calculated using the burnup feature in MCNP6.Due to the limit on computational resources placed by these types of calculations, burn-up will be carried out on individual fuel assemblies.The full definition of the core states will then be calculated using MCNP6 critical calculations, with the input files for these states written using features in NWURCS.These features are presently being developed.

Conclusions
The prospect of a nuclear new build in South Africa brings with it a number of key challenges that should be addressed, including localisation.It has been identified that nuclear reactor analysis is one   k p k t field that can have a strong national localisation.This will support both economic and safety aspects of building and running a nuclear reactor.With this in mind, a light water reactor was chosen for reactor analysis.The reactor chosen was a VVER-1000 type reactor.The neutronic analysis was performed using MCNP6 for the cold zero power state at the beginning of cycle with input specifications obtained from open literature.The input file was generated using the in-house code NWURCS and some aspects of the verification of NWURCS were presented.The model was constructed so that the convergence tests were passed.Thereafter, the critical reactor state was obtained by calculating the critical boron concentration.Control rod worth, reactivity coefficients and were calculated.The results show that these parameters are within acceptable limits.However, this task also need to be carried out at hot full power states at the beginning of cycle.This will be carried out by coupling the MCNP6 calculation with a RELAP5 calculation.
The geome-68 Journal of Energy in Southern Africa • Vol 27 No 4 • November 2016

70
Journal of Energy in Southern Africa • Vol 27 No 4 • November 2016

Figure 6 :
Figure 6: VVER-1000 core layout showing the positions of the fuel pins within the fuel assembly, and the position of the fuel assembly within the reactor core.

Figure 7 :
Figure 7: Reactor pressure vessel and internals showing the layout of the various components within the reactor pressure vessel.

Figure 8 :
Figure 8: Full core model fission source convergence.

Figure 9 :
Figure 9: Ten statistical tests as they appear in MCNP6 output.

Figure 11 :
Figure 11: Calculation of the critical boron concentration.

74
Journal of Energy in Southern Africa • Vol 27 No 4 • November 2016