Home > Research


Research activities in CERECAM are organised into a number of research programmes which range from those of a fundamental nature to projects having a direct link to industrial and other applications. The research activities are arranged into six broad groupings:

A:  Computational solid and structural mechanics

B:  Numerical analysis and PDEs

C:  Biomechanics

D:  Fluid dynamics

E:  Particluate flow

Details of activities and of student involvement and interactions with outside collaborators follow below:



Strain-gradient plasticity

BD Reddy, F Ebobisse, AT McBride

Collaborators: S Bargmann (University of Wuppertal), ME Gurtin (Carnegie-Mellon University), P Neff (University of Duisburg-Essen), P Steinmann (University of Erlangen-Nuremberg), T Böhlke (Karlsuhe University of Technology)

Students graduated: NJN Richardson (MSc Eng, 2012), T Povall (MSc, 2013), D Gottschalk (PhD, 2014), N Mhlongo (MSc)

Strain-gradient theories of plasticity have been the subject of considerable attention and study for more than three decades, since the early contribution by Aifantis [1]. The interest in such theories lies in their intrinsic ability to introduce a length scale that is absent in classical theories. The length scale together with the non-local nature of the gradient terms allows for a theory that captures the size dependence observed in experiments. Furthermore, the inclusion of gradients of plastic strain makes possible a link at the continuum level between observed size effects and the underlying behaviour of geometrically necessary dislocations. Some representative and important works include those by Fleck and Hutchinson, Gudmundson, Gurtin and Anand, and Nix and Gao [2–6].

There has been sustained activity at Cerecam in the area of strain-gradient plasticity for around a decade. The first works took as a starting point the rate-dependent formulation for strain-gradient plasticity developed by Gurtin, Anand, and co-authors (see for example [5]), and developed a rate-independent model with an associative flow relation and yield criterion in microstress space [ReddyEbobisseMcBride2008, Reddy2011a]. These works also developed variational formulations of the problem and presented conditions for uniqueness and existence of solutions, accounting for both energetic and dissipative gradient effects. A key feature of the variational formulation is that it allows for the flow relation to be formulated in terms of the Cauchy stress, provided that this is done in global form; the local flow relation is expressed in terms of the microstresses, which are indeterminate in the elastic range and which therefore cannot serve as indicators of yield. Closely related work [GurtinReddy2009] showed how the classical isotropic hardening relation may be incorporated equivalently through an irreversible defect energy. The case of single-crystal plasticity is the subject of a companion paper [Reddy2011b] to [Reddy2011a].

Recent work has been concerned with the behaviour of strain-gradient plasticity models under conditions of non-proportional loading. Motivated by the work of Fleck, Willis and Hutchinson [7], the variational structure of the problem has been revisited in [Carstensen-etal2017a,b] with a view to gaining a better understanding of the occurrence of an elastic gap, reported in [7]: that is, an elastic response to a non-proportional change in loading, following loading into the plastic range.

Early computational work on strain-gradient plasticity at Cerecam has focused on the development, analysis and implementation of Discontinuous Galerkin (DG) approximations of the Aifantis model [Djoko-etal2007a,b]. Further elaboration of this approach, including the extension to large-deformation problems, has been reported in [EbobisseMcBrideReddy2008,McBride-Reddy2008,McBride-Reddy2009]. The work [BargmannReddy2011] returned to the viscoplastic Gurtin-Anand model and obtained computational results for oligocrystals.

An important aspect of models of strain-gradient plasticity is on the form of the defect energy included in the free energy. Various forms of this energy have been proposed in the literature. Using the defect energy proposed by Gurtin and Anand, some mathematical analysis was developed in the infinitesimal setting by Ebobisse, Neff and co-authors (see for instance [EbobisseNeff2010] and [EbobisseHackNeff2017]) with the purpose of accommodating the so-called plastic spin. Also, a micromorphic-type formulation of models of strain-gradient plasticity called microcurl was proposed by Forest with the aim of avoiding, for instance, higher order boundary conditions which arise at the interfaces between elastic and plastic parts in the traditional formulation. The mathematical analysis for Forest’s microcurl formulation was considered by [EbobisseNeffForest2017] for both single crystal and polycrystals. Following the work of [Reddy2011b], further mathematical analysis for models of single crystal gradient plasticity was studied in [Ebobisse-Neff-Aifantis2017]. On the other hand, a fourth order model of strain-gradient plasticity based on a defect energy density which is instead a quadratic form in the so-called Kröner’s incompatibility tensor was recently introduced and analysed in [EbobisseNeff2018].

Formulations of self-and latent hardening for single-crystal strain-gradient plasticity have been developed in [GurtinReddy2014], with computational investigations in the context of oligocrystals in [BargmannReddyKlusemann2014]. The stored energy of cold work and related issues in single-crystal thermoplasticity have been investigated theoretically in [AnandGurtinReddy2015] and computationally in [McBrideBargmannReddy2015]. A detailed convergence analysis and implementation was presented in [ReddyWienersWohlmuth2012] for the single-crystal theory with gradients accounted for purely energetically.

There have also been theoretical and computational investigations of a grain-boundary model, due to Gurtin [8] (2008), which accounts for grain misorientation and grain-boundary orientation [Gottschalk-etal2016,McBrideGottschalk-etalTechMech2016]. Related work [Bayerschen-etal2016] explores various slip transmission criteria, with a link to experimental results as appropriate.

Various overviews, some of these based on keynote presentations, provide survey accounts of the developments in strain-gradient plasticity at Cerecam. These include the works [Reddy2013a,b] which include new material on variational formulations of the large-deformation problem for single crystals. The monograph [HanReddy2013] presents a detailed account of the mathematical theory of conventional plasticity, together with material on the gradient theory as developed at Cerecam, as well as associated analyses of convergence of finite element approximations and associated algorithms.

Related work in conventional plasticity. Related work on algorithms for the rate-independent problem of single-crystal plasticity has been reported in [RichardsonMcBrideReddy2011], while a range of hardening relations, including those presented in [GurtinReddy2014] referred to above, are investigated computationally in the context of the large-deformation theory for single crystals. The work [DeLorenzisMcBrideReddy2016] reports on a phase field approach to fracture in conventional single-crystal plasticity. In [GurtinReddy2016] the concept of an intermediate configuration in large-deformation plasticity is revisited; it is shown in this work that what is usually referred to as the intermediate configuration is instead a vector space, and that this intermediate space represents the lattice, when applied to crystals.


  1. Aifantis E., On the microstructural origin of certain inelastic models. Journal of Engineering Materials and Technology 106 (4), 326-330, 1984
  2. Fleck N.A., Hutchinson J.W., Strain gradient plasticity. Advances in Applied Mechanics 33, 295-361, 1997
  3. Fleck N.A., Hutchinson J.W., A reformulation of strain gradient plasticity. Journal of the Mechanics and Physics of Solids 49, 2245-2274, 2001
  4. Gudmundson P., A unified treatment of strain gradient plasticity. Journal of the Mechanics and Physics of Solids 52, 1379-1406, 2004
  5. Gurtin M.E., Anand L., A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: small deformations. Journal of the Mechanics and Physics of Solids 53, 1624-1649, 2005
  6. Nix W.D., Gao H., Indentation size effects in crystalline materials: a law for strain gradient plasticity. Journal of the Mechanics and Physics of Solids 46, 411-425, 1998
  7. Fleck, N.A, Hutchinson J.W. and Willis J.R., Strain-gradient plasticity under non-proportional loading. Proceedings of the Royal Society A 470, 20140267, 2015


High-performance crystal plasticity

M Malahe, BD Reddy

The aim of this project is to implement a framework for computations of polycrystal plasticity formulations that scales efficiently to petascale computers. The framework combines a spectral representation of the grain-level responses implemented on graphics processing units (GPUs) with a homogenised finite element representation implemented on CPUs. The domain is split across computational nodes using the same decomposition for both representations, leading to minimal inter-node communication.

The initial work to achieve this was focused on creating an efficient GPU spectral crystal plasticity solver. The starting point was the solvers by Mihaila et al. [1] and Savage et al. [2], which had achieved a three order of magnitude speedup over CPU-based solvers. The work at Cerecam then introduced additional gains in efficiency primarily through reducing global memory transactions and exposing more instruction-level parallelism in the GPU solver [Malahe2018]. The resulting solver was an additional order of magnitude more efficient in time, and three orders of magnitude more efficient in memory than prior GPU solvers.

The new GPU solver allows for the calculation of tens of millions of grain level responses per second on a desktop, which has driven the question of how to effectively use this new ability to solve engineering problems.  Current work in this direction is on using model predictive control and reinforcement learning to design optimal processing paths for metals to achieve desired textures. These methods require the computation of a large number of possible processing paths to be effective, which is an ideal scenario for using the new solver.

In parallel with the above, work is being conducted on integrating the GPU crystal plasticity solver into a full finite element solver.



Figure 2:  The distribution of grain orientations in oxygen-free high thermal conductivity copper following a simple shear test. Each pole figure shows the distribution for a particular crystal plane. An MRD (Multiples of Random Distribution) of 1.0 indicates the density one would expect from a random distribution.


  1. Mihaila B, Knezevic M, Cardenas A, Three orders of magnitude improved efficiency with high-performance spectral crystal plasticity on GPU platforms. International Journal for numerical Methods in Engineering 97(11), 785–798 (2014). DOI 10.1002/nme.4592
  2. Savage D.J., Knezevic M, Computer implementations of iterative and non-iterative crystal plasticity solvers on high performance graphics hardware. Computational Mechanics 56(4), 677–690 (2015). DOI 10.1007/s00466-015-1194-6


Numerical models for strain-induced crystallization in polymers

BD Reddy

Student: EB Ismail (PhD)

In polymer science, there exists a large family of thermoplastic materials that are present in either amorphous or crystalline solid phases at working temperatures. During unperturbed cooling from the liquid “melt” phase, the amorphous solid phase is usually achieved. A transition from the amorphous to the crystalline phase can occur, however, as a result of strain-induced crystallization. This occurs when the material is deformed at temperatures exceeding the glass transition temperature. In this temperature range the material is still solid, but crystals may grow in the direction of applied stresses.

The amorphous material is often modelled as isotropic and viscoelastic, with flow rules defining the displacement response to applied loading. The crystalline phase is typically anisotropic and elastic. Crystallization occurs progressively, and a mixture of amorphous and crystalline phases is common. 

Several approaches exist for modelling the crystallization process mathematically. The approach followed introduces internal variables that are linked to the deformation of the material and relate to preferred orientations and configurations[1-3]. By tracking strain relative to both the reference and preferred configurations, as well as the crystallization fraction, it is possible to swop phenomenologically the macroscopic behaviour of crystallizing systems. Such models have typically been limited to simple geometric configurations and loading states, where the mathematics remains analytically tractable.

This project aims to develop finite element models within this framework, and associated algorithms that are sufficiently robust for use in applications.


  1. Rajagopal etc.


Numerical simulation of friction welding processes

BD Reddy

Collaborators: AT McBride (University of Glasgow), D Hattingh (Nelson Mandela University)

Student: M Hamed (PhD)

Friction welding is a family of solid state joining processes where friction is used to generate the necessary welding heat.

These processes are used in various important applications where other joining methods are infeasible or produce inferior results. Numerical simulation of friction welding processes can be a useful tool in their deployment in new applications by aiding with selection of welding parameters and with requirement specification for friction welding equipment. In addition, data from friction weld tests can be combined with numerical simulation to analyse material behaviour at the high temperature and deformation rate ranges characteristic of these processes.

In this project, simulation of friction welding processes is carried out through finite element analysis of a large-strain coupled thermo-viscoplasticity model with friction contact. To address the extremely large deformations undergone by the material at the thermo-mechanically affected zone during these processes, an arbitrary-Lagrangian-Eulerian procedure is being implemented into the solver. The ALE method is extended here to account for the consistent projection of the plastic variables as the mesh moves. Additional variables will be incorporated in the model to study the interaction between process parameters and relevant macro and microscopic material properties.

The friction welding model is being developed within the open-source finite element library deal.II.


Micromechanical modelling of advanced heirarchical composites

BD Reddy

Collaborator: S Bargmann, J Wilmers (University of Wuppertal)

Student: E Griffiths (PhD)

The impact of small-scale geometric confinement on deformation mehanisms is the subject of intensive current research in materials science. Nanoporous metals have a micro-structure with an extremely high volume-specific surface content. Due to high local strength and a relatively regular interconnection of the nanoconstituents as well as low mass density, nanoporous metal is filled with polymer, it offers the opportunity to use the material as actuators or sensors.  Actuators translate electrical or chemical energy into mechanical work, resulting in a strong coupling of the mechanical and, in the present case electrical behaviour.  Current work is concerned with the analysis of a representative volume element in the form of a nanoporous gold-polymer composite, with a view to determining its average properties and the nature of the stress strain behaviour within different regions of the ligaments comprising the gold phase under various loading conditions, [GriffithsBargmannReddy2018]



Development, analysis and implementation of new finite elements

BD Reddy, AT McBride

Collaborators: C Carstensen (Humboldt University, Berlin), BP Lamichhane (University of New South Wales), P Wriggers (Leibniz University of Hanover)

Post-doctoral researcher: B Grieshaber

Students graduated: A Chama (PhD, 2014), B Grieshaber (PhD, 2013), D van Huyssteen (MSc Eng)

Current students: F Rasolofoson (PhD), K Penzhorn (PhD), D van Huyssteen (PhD)

Research in this area has been directed towards the development, analysis and implementation of new low-order finite elements that are simple but stable and efficient, particularly in situations involving small parameters such as those that arise in situations of near incompressibility, and for thin domains in which one or more dimensions are much smaller than the others.

The focus has been to undertake detailed analyses and to generate novel approaches, largely within the context of mixed and discontinuous Galerkin formulations. Earlier work, which was based on the use of low-order quadrilaterals, has been extended to broader families of elements. Stokes-stable pairs are an essential feature of the three-field theory, and well-known examples of these pairs have been used as the base on which to build a wide range of stable elements for three-field formulations [LamichhaneReddyMcBride2013, ChamaReddy2013,2014]. The investigation has been further extended to the case of nonlinear elasticity, in which conditions for uniform convergence of the linearized problem are established.

Related work has been carried out on the discontinuous Galerkin (DG) method. The focus has been on the design of quadrilateral elements in two dimensions, and hexahedral in three, which are robust and uniformly convergent in the incompressible limit [GrieshaberMcBrideReddy2015].

An investigation has been carried out of convergence of optimal-order finite element methods for the Bingham flow problem, which is characterized by a variational inequality [Carstensen ReddySchedensack2016].  New results on optimal-order convergence for mixed and nonconforming formulations have been obtained.

A new area of study (with P Wriggers in Hanover) focuses on the virtual element method, a recently developed variant of the finite element method which is able to accommodate arbitrary polygonal elements. The objective will be to extend the existing linear theory to problems of nonlinear solid mechanics. Applications to contact and to nonlinear problems are of particular interest [WriggersRustReddy2016, WriggersReddyRust Hudobivnik2017].

Current work includes the development and analysis of conforming, discontinuous Galerkin and virtual element methods for transversely isotropic elastic materials. As with earlier work, the objective is to develop methods that are stable and uniformly convergent in the incompressible and inextensible limits.


Design and analysis of stable staggered numerical schemes

BD Reddy, AT McBride

Post-doctoral researcher: MF Wakeni (2016 – 2018)

The direct approach for approximating solutions of transient coupled problems is the monolithic method, where all the unknown field variables are solved simultaneously at discrete time points. However, the primary disadvantages of such monolithic schemes include the computational cost of solving the resulting large system of algebraic equations at each time step, and the time step-length restriction which would be imposed for the entire problem in the case of time scales of differing orders of magnitude.

On the other hand, coupled problems can be solved efficiently using staggered approaches (operator-splitting algorithms). They are methods for coupled problems based on the ‘divide-and-conquer’ approach, in the sense that the full coupled problem is split into simpler sub-problems for which appropriate sub-algorithms can be applied to approximate their solutions.

The stability of such staggered algorithms is mainly determined by that of the suboperators in the splitting.

The model coupled problem is that of generalized thermoelasticity with particular attention to the case of finite strains. The generalized model of thermoelasticity is inherently nonlinear and accounts for the second sound phenomenon, a wave mechanism of thermal energy transmission. Its linearization about the reference state (configuration) leads to the wellknown non-classical model of Green and Naghdi. The initial-boundary value problem is a two-way coupling comprising the classical problem of elasticity, which is hyperbolic, and the dominantly hyperbolic generalized heat conduction model.

An operator-split in to mechanical and thermal phases based on isoentropy leads to two stable sub-problems which are dominantly hyperbolic. Time-discontiouous Galerkin finite element (TDG-FE) schemes are shown to be stable, and have been used to develop numerical simulations [WakeniReddyMcBride2016a,b]. Current work includes the design and analysis of PDE-based shock capturing schemes; applying the fully-nonlinear generalized thermoelastic model to the analysis of the thermomechanical response of skin tissue; extending the current TDG-FE method to one which allows discontinuities both in space and time; investigating iterative solvers and preconditioners such as multigrid methods in the framework of staggered TDG-FE schemes; and the formulation of efficient nonlinear solvers such as the Jacobianfree Newton-Krylov (JFNK) method in combination with the staggered TDG-FE schemes for the generalized thermoelastic model at finite strains.


Operator-splitting methods for coupled nonlinear partial differential equations

BD Reddy

Collaborators: JK Djoko, J M-S Lubuma, AR Appadu (University of Pretoria)

Post-doctoral researcher:  HH Gidey (2017 to April 2019)

Nonlinear and coupled nonlinear partial differential equations arise in many fields of science and engineering such as fluid mechanics, gas dynamics, thermodynamics, elasticity, phase separation and material sciences.

This research project is focused on the development and implementation of operator splitting methods, designing multilevel/cost effective time-stepping algorithms for higher order nonlinear scalar pdes and strongly coupled nonlinear pdes. 

An operator splitting approach is applied to find the approximate solution of the 2D convective Cahn-Hilliard equation, where the full model is split into three sub-problems and each sub-equation is solved by existing efficient numerical method. Proofs of the existence of a unique solution, and the stability and convergence of solutions of the splitting methods, are provided [AppaduDjokoGidey2017a,b, AppaduDjokoGideyLubuma2017].

Current work includes the development and analysis of multilevel finite volume methods for the two dimensional incompressible Navier-Stokes equations. The methods are constructed in such a way that some properties of the continuous model are satisfied at the discrete level. 

Further work is concerned with the development and analysis of efficient and robust operator splitting methods and design an effective time-stepping algorithm for multiphase flow models.



Biofluid Mechanics

Hydrogel dispersion in myocardial tissue

MN Ngoepe, T Franz (UCT, Biomedical Engineering)

Collaborators: N Davies (University of Cape Town), D Bezuidenhout (University of Cape Town), S Balabani (University College London), A Passos (University College London)

Current student: A Maluleke (PhD)

Myocardial infarction (MI), a type of cardiovascular disease, affects a significant proportion of people around the world. Traditionally, non-communicable chronic diseases were largely associated with ageing populations in higher income countries. It is now evident that low- to middle-income countries are also affected and in these settings, younger individuals are at high risk.

Currently, interventions for MI prolong the time to heart failure. Regenerative medicine and stem cell therapy have the potential to mitigate the effects of MI and to significantly improve the quality of life for patients. The main drawback with this therapy is that many of the injected cells are lost due to the vigorous motion of the heart. Great effort has been directed towards the development of scaffolds which can be injected alongside stem cells, in an attempt to improve retention and cell engraftment. In some cases, the scaffold alone has been seen to improve heart function.

Our work has focused on polyethylene glycol (PEG) gel, a synthetic hydrogel injected into the heart to improve function [1]. Part of the work has focused on characterising the flow of PEG on the microscale, where the behaviour is likely to deviate from macroscalar flow patterns. Micro particle image velocimetry (μPIV) is used to observe flow behaviour in microchannels, representing the gaps in myocardial tissue. The fluid exhibits Newtonian behaviour at this scale.

An idealised, three-dimensional computational fluid dynamics (CFD) model of PEG flow in microchannels is then developed and validated using the μPIV study. The validated computational model is applied to idealised models and to realistic, microscopy-derived myocardial tissue models.

Figure 3: Flow patterns in different segments of the myocardial tissue model. (A) corresponds to the inlet and outlet profiles of Block 1 in Error! Reference source not found.. The maximum velocity is observed on the outlet. (B) corresponds to Block 2, and high velocities are observed on the inlet and the outlet. (C) corresponds to Block 3 and the maximum velocity is observed on the outlet.


1. Davies N., GoetschK. Ngoepe M., Franz T. and Lecour S., Delivery modes for stem cell therapy.  In Stem Cells and Cardiac Regeneration. Ed. R. Madonna 2016, Springer, Switzerland


Thrombosis in cerebral aneurysms

MN Ngoepe

Collaborators: W-H Ho (UNISA), Y Ventikos (University College London), T Peach (University College London), S Hume (MSc Eng)

Current students: S Hume (PhD), M Taylor (MSc)

Thrombosis, or clotting, is the main underlying condition for a large number of cardiovascular disease cases. Under normal conditions, clotting is a positive physiological feature which ensures that bleeding stops following injury to a blood vessel. In some cases, however, the fine balance sustained during the clotting process is disrupted, leading to the formation of clots that cause blockage of blood vessels, and subsequent morbidity or mortality.

The focus of this work is on the development of computational tools and platforms which can be used to better understand thrombosis in disease cases. At present, the focus is on the development of such tools for the elucidation of clotting in cerebral aneurysms. Cerebral aneurysms are balloon-like bulges which form on blood vessels of the brain as a result of the weakening of the vessel wall layers. Thrombosis is closely linked to cerebral aneurysms and the rupture risk of an aneurysm is influenced by the type of clot that forms in the aneurysm sac [Ngoepe et al, 2018].

The motivation for developing a computational platform for cerebral aneurysm thrombosis is twofold. First, the mechanisms which contribute to the initiation and progression of clotting in cerebral aneurysms are poorly understood. This is mainly because it is difficult to predict when a clot is likely to form in an aneurysm and even more challenging to account for the contribution of key players (e.g. fluid dynamics, blood clotting proteins) during the clotting process. A computational platform is able to give valuable insights, as the different contributors can be included or omitted without catastrophic medical consequences. The second key reason is a computational tool capable of predicting clotting outcome in cerebral aneurysms prior to surgical intervention is highly desirable [NgoepeVentikos2016, PeachNgoepeSprangerZajarias-FainsodVentikos2014]. Given the patient-specific nature of cerebral aneurysms and clotting, being able to predict clotting outcome on a per patient basis and later inferring general population trends would be beneficial for informing treatment approach.


Biothermomechanics of skin

AT McBride (University of Glasgow)

Collaborators:  S Bargmann (University of Wuppertal), G Limbert (University of Southampton)

Student graduated:  D Pond (MSc Eng, 2017)

The skin is the largest organ in the human body. It is the first line of contact with the outside world, being subject to a harsh array of physical loads and environmental factors. The skin also performs numerous physiological tasks such as thermo-regulation, and vitamin D synthesis. The skin undergoes chronological (intrinsic) ageing, whereby there is a general breakdown of tissue function and a decline in mechanical properties. In addition to this, skin undergoes extrinsic forms of ageing through exposure to external factors such as ultraviolet radiation.

The objective of the work is to develop a constitutive model for the ageing of skin that is both capable of capturing the material and structural effects, and is embedded in the rigorous framework of non-linear continuum mechanics. The macroscopic response of the skin is due to microscopic components such as collagen, elastin and the surrounding ground substance and the interaction between them. The ageing process is explored and a firm understanding of the influence of ageing on the substructures is established.

The constitutive model is implemented within a finite element scheme. Published experimental data used to conduct a series of ageing numerical experiments to ascertain the response of the model to changes in material properties associated with ageing. A modified model is then proposed to capture the ageing response of the skin. The key microscopic biophysical processes that underpin ageing are identified, approximated and adapted sufficiently to be of use in the macroscopic continuum model. Aspects of open-system thermodynamics and mixture theory are adapted to the context of ageing in order to capture a continuous ageing response [McBrideBargmannPondLimbert2016, PondMcBrideDavidsReddyLimbert2016].


Computational simulation of bone remodeling during reverse-shoulder arthroplasty

AT McBride

Collaborators:  S Roche (Department of Orthopaedics, UCT), S Sivarasu (Division of Biomecial Engineering, UCT)

Student graduated:  H Liedtke (MSc Eng, 2017)

Bone is a living material. It adapts, in an optimal sense, to loading by changing its density and trabeculae architecture - a process termed remodelling. Implanted orthopaedic devices can significantly alter the loading on the surrounding bone. In addition, these devices rely on bone ingrowth to ensure secure implant fixation. In this work, a computational model that accounts for bone remodelling is developed and used to elucidate the response of bone following a reverse shoulder procedure. The reverse shoulder procedure investigated here is for rotator cuff deficient patients. In this procedure, up to 75% of complications are reported in some clinical series. It is therefore necessary, for the design of successful implants, to understand the loading environment to promote bone growth in the correct areas. The physical process of remodelling is modelled using continuum scale, open system thermodynamics whereby the density of bone evolves isotropically in response to the loading it experiences. The fully-nonlinear continuum theory is solved approximately using the finite element  method.


Computational modelling of the mechanical response of turtle shells

BD Reddy

Collaborator: S Bargmann (University of Wuppertal)

Student: B Alheit (MSc Eng)

A turtle’s shell is a bio-composite material consisting of an underlying layer of interlocking bone plates that are bound together by collagen and covered in a layer of keratin. The constituent components of the shell have poor mechanical properties however, the shell overall acts as a tough and robust material due to its structure. The aim of this project is to develop a computational model of a turtle shell that will allow investigation into how each aspect of the structure of the shell affects the shell’s mechanical response. Thus, gaining insight as to how synthesised composite materials can be structured in order to have more desirable mechanical properties.



Effect of system pressure on leakage from water distribution systems

JE van Zyl, BD Reddy

Collaborator: CRI Clayton (University of Southampton)

Student graduated: E Ssozi (MSc Eng, 2015)

Leakage from water distribution systems is becoming a serious problem worldwide as systems age and water resources are placed under increasing strain. Field and laboratory tests over the past decade have shown that leakage from distribution systems are often substantially more sensitive to changes in pressure that conventional theoretical models predict. Research has shown that this is mainly due to leak areas varying with pressure (in addition to changes in leakage velocity), but that other factors may also play a role. The aim of this project has been to investigate the various mechanisms involved in the pressure-leakage response [SsoziReddyvanZyl2016]. Currently the project is investigating viscoelastic deformation in plastic pipes and the soil leak interaction.


Modelling and analysis of the flow of complex fluids

T Chinyoka

Collaborators: IE Ireka (Fraunhofer Institute for Industrial Mathematics, Germany), OD Makinde (Stellenbosch University)

Graduated students: IE Ireka (PhD, 2015), ZS Nyandeni (MSc, 2017), P Nchupang (BSc Hons, 2018), FNZ Rahantamialisoa (MSc, 2018), A Mavi (MSc, 2019)

Current students: JG Abuga (PhD), Idrees (PhD), A Mavi (PhD), T Chadenga (MSc)

Affiliated Students: N Mudau (MSc), T Hlope (MSc), D Muchiri (MSc) (AIMS, Cape Town)

Ongoing work is concerned with the mathematical modelling, analysis and computational solution of flows of complex (non-Newtonian) fluids under various conditions. Of general interest are investigations of the effects of non-isothermal processes and/or flow non-homogeneities in polymeric (viscoelastic) flows with applications in engineering and biological systems.

The early work on non-isothermal processes in polymeric systems revolved around the investigations on design and efficiency of important industrial applications, such as lubrication [Chinyoka2008]; heat exchanger modelling and analysis [Chinyoka2009a, Chinyoka2009b]; exothermic effects [Chinyoka2008, Chinyoka2010, 2011a,b]; and gravitational effects [Chinyoka et. al. 2013b].

Much of this work has informed the direction of the current and ongoing work on secondary and/or combination processes in the flow of complex fluids; processes such as shear banding [Chinyoka2011a], combined shear banding and lubrication effects [IrekaChinyoka2013], combined shear banding and non-isothermal effects [IrekaChinyoka2016], and dynamics of polymeric foams [Ireka et. al. 2015].

Extensive collaborative studies have been conducted, over the past decade, on both Newtonian and non-Newtonian fluid dynamics. Such studies include particulate flows and pollutant dynamics [MakindeChinyoka2010a,b;ChinyokaMakinde2010a,2013a,b]; convection reaction flows [MakindeChinyoka2010c,2011;2012]

ChinyokaMakinde2010b,2012a,b;2013a,b; Makinde et. al. 2011a,b;Lebelo et. al 2017]; suction-injection and porous media flows [ChinyokaMakinde2010a,2013c,2014b,2015b; Makinde et. al. 2011a, Chinyoka et. al. 2013a]; magnetic effects and electrical conductivity [MakindeChinyoka2010c,2011,2013; ChinyokaMakinde2012b,2013d,2014a,2015a Chinyoka et. al. 2013a].

Fundamental fluid dynamical themes, such as variable viscosity modelling and effects, run through the majority of the studies already listed. A lot of the work with the collaborators led to us being editors of the special issue on New Developments in Fluid Mechanics and Its Engineering Applications [Makinde et. al. 2013].


Optimal design of turbine blades

BD Reddy

Collaborator: G Snedden (CSIR: co-supervisor)

Student graduated: J Bergh (PhD, 2018)

This project concerns the mitigation of the aerodynamic losses associated with the highly three-dimensional (secondary) flows which develop in turbine blade passages. This is important because even relatively small increases in component efficiencies can be shown to have a proportionally larger effect on overall engine thermal efficiency and performance. Recently, non-axisymmetric contouring of the hub and shroud endwalls turbine blade passages have been shown to be an effective method for reducing the strength of the secondary flows[BerghSneddenMeyer2012]. The current investigation employs a surrogate-based optimization (SBO) algorithm coupled to a commercial CFD code to design various sets of blade / contour combinations using a selection of objective functions previously used in industry specifically for this purpose. In addition to the numerical design phase, physical testing of each set of contours will be conducted in a low speed, 1 stage research turbine in order to validate the numerical predictions associated with each design and ultimately determine the relative efficacy of each objective function as a basis for the design of endwall contours within the industrial design process [BerghSnedden2015].



A Mainza, I Govender (University of Kwazulu-Natal)

Post-doctoral researcher: S Bremner

Other researchers: MC Richter (University of Cape Town)

Collaborators:  Paul Cleary (Commonwealth Scientific and Industrial Research Organisation), Marcelo Tavares (COPPETEC), Magnus Evertsson (Chalmers University), Guven Akdogan (Chemical Engineering, Stellenbosch University), Narasimha Mangadoddy (IIT Hyderabad), Dion Weatherly (Julius Kruttschnitt Mineral Research Centre), Q Reynolds (Mintek)

Students graduated (2012 - 2017):  D. Blakemore (MSc), H. Brodner (MSc), A. Morrison (MSc), M. Malahe (MSc), M. Hromnik (MSc), O. Oliwaseun (PhD), M. Elbasher (PhD), S. Bremner (PhD), N. Ngoepe (PhD), P. Thirunavukkarasu (PhD), L. Bbosa (PhD), K. von Daramy (PhD)

Current students:  C. Ndimande (PhD), D. de Klerk (PhD), A. Mabentsela (PhD), T. Latona (PhD)

Student graduated:  T. Povall (PhD, 2018)

Granular Flow Modelling combines the inherent frictional nature of particles with their distinctly fluid-like structure. We focus on two classes of flows: (i) confined flows and (ii) free surface flows.

Tumbling mills, chutes and avalanches are examples of free surface flows while cyclonic separators and stirred mills fall into the class of confined flows. Given the high concentration of solids in industrial flows, the granular flow approximation is unique in its ability to capture the solid-solid interactions within a Navier-Stokes-like framework, typical of fluid flow. The physically valid dissipative mechanisms associated with granular material allows for more realistic quantification of abrasion and attrition in comminution processes while accurately capturing important flow features like free surface shape and velocity field profiles.

Industrial systems such as tumbling mills are typified by rotational and axial flows of rock, steel balls and slurry. The tortuous porous network, created by the non-uniform packing of rock and steel balls, form the channels through which the viscous slurry flows. A multi-pronged approach employing non-invasive nuclear techniques (PEPT, X-ray) and computational modelling (DEM, CFD (principally, Finite Elements and/or Finite Volume Methods), SPH) has formed the key ingredients for mechanistic modelling of such industrial systems.

Biological processes such as breathing are characterised by non-linear deformation of the tissue and musculature in the human upper airway. PEPT and X-ray imaging are employed to elucidate the mechanisms responsible for conditions such as obstructive sleep apnoea. The associated airflow that provides the stimulus for the final deformation is also tracked experimentally using PEPT. These measurements complement coupled FEM-CFD modelling of the upper airway.


Positron emission particle tracking (PEPT)

Positron Emission Particle Tracking (PEPT) is a technique for studying the flow of particulate systems such as tumbling mills in the minerals industry. Initially developed for the medical imaging industry, positron emission tomography has been adapted for engineering applications at the University of Birmingham. The particular value of PEPT is the ability to look deep within the particulate system for extended periods of time, thereby elucidating the in-situ kinematics and dynamics of the flow.

The basic principle of PEPT is based on positron annihilation. A single particle is labelled with a radionuclide that decays via beta-plus decay, resulting in two gamma rays, each of energy 511 keV travelling in exactly opposite directions. Simultaneous detection of the two gamma rays in an array of detectors defines a line along which the annihilation occurred. Detection of a few such events in a very short time interval allows the position of the tracer particle to be triangulated in three dimensions.

Location in space of the tracer particle may be achieved at a frequency up to 1000 Hz with an accuracy which depends on the speed of the tracer particle. PEPT is currently the only non-invasive technique capable of mapping the in-situ flow fields in robust, industrial systems to the level of detail that is demanded for mechanistic modelling. Advances in computing have made numerical modelling of complex flows conceivable. However, simulations are computationally expensive and time intensive. Consequently, numerical modelling work employs simplifying assumptions that ultimately make the computing tasks tractable. The integrity of these assumptions requires validation if they are to gain confidence within industry. PEPT offers detailed validation of the flow field and related parameters.

In 2009, an EXACT3D PET camera that was decommissioned at Hammersmith Hospital, London, was donated to UCT. The scanner had been used for clinical research at Hammersmith Hospital, London, since 1995, and may still be the most sensitive 3D PET scanner in operation today. The scanner is housed at the iThemba LABS cyclotron centre in Cape Town. The situation of the PEPT laboratory at iThemba LABS has all the obvious advantages with respect to radiation handling and licensing. PEPT experiments require positron-emitting radioisotopes which are produced by cyclotron proton beams. iThemba LABS routinely produces radioisotopes for medical PET use and other physics applications and research.


Computational modelling

A wide range of industrial systems depict behaviour at the length scale of interest that is dominated by interactions of discrete material regions. In the minerals industry, for example, there are many processes which may be modelled as discrete particle distributions within a fluid medium. Fluid-particle mixtures constitute a rheologically complex fluid which is transported through a dynamic porous bulk. Mineral processing machines, such as tumbling mills and flotation cells, are presently described mainly by phenomenological models. The preponderance of such models is strongly coupled to the lack of in-situ measurements which are needed to elucidate the physical parameters that constitute the underlying mechanism. A clear understanding of the mechanistic nature is crucial to the efficiency and overall optimisation of mineral plants.

The major advantage of the computational frameworks like the Discrete Element Method (DEM), Finite Element Method (FEM), Finite Volume Method (FVM), and Smooth Particle Hydrodynamics (SPH) is the ability to approximate the mechanical environment with the potential to discern meaningful in-situ measurements and behaviour.

However, the computational demands and lack of sound experimental verification have limited the value of DEM and CFD techniques in many industries. Thus it is necessary to fill the vital gap linking computational results to rigorous experimental data. It is only with validation that any confidence can be given to the predictive capability of such computational tools, especially when the predictive range lies outside the range over which the existing semi-empirical models were developed and tested. In-situ measurement tools like Positron Emission Particle Tracking (PEPT) and X-ray imaging allow a detailed investigation into aspects such as contact models, energy distributions and flow.

Discrete Element Method (DEM)

The Discrete Element Method (DEM) was developed by Cundall and Strack for analysing quasi-static soil mechanics problems. The DEM was revolutionary in that it modelled the individual particles in a discrete system separately and thus represented the first truly discontinuous numerical model. The particles were originally modelled as rigid, circular discs that interacted at contacts only, but DEM has been extended to three dimensional and non-spherical objects by several researchers. Rigid spherical particles as implemented in the commercial software package EDEM developed by DEM Solutions and LIGGGHTS® developed by CFDEM® project are used to perform the simulations of an experimental scale tumbling mill.

Contact between particles are modelled using the soft-contact approach wherein contacting particles are allowed to overlap one another at the point of contact. A contact force law is used to relate the relative overlap between the particles to a restoring contact force. The most commonly used contact force law is the viscous damping model. A stiff numerical spring is applied in the normal direction at the contact point. The undamped restoring spring force in the normal direction is calculated as the product of the relative particle overlap and the spring stiffness according to Hooke's law. A second numerical spring, the shear spring, is applied in a direction orthogonal to the normal spring to simulate the frictional forces at the contact. The maximum frictional force cannot exceed the limiting frictional force governed by the Mohr-Coulomb Law. Energy dissipation to non-frictional mechanisms is modelled using numerical dash-pots that oppose the contact forces.

Newton's second law is used to model the rigid body motion of the particles arising from the contact and gravitational forces acting on the particle. The assumption of rigid particles is valid when the majority of the deformation in the system occurs along interfaces. Thus, the DEM is highly suited for the analysis of granular flow, where the deformation of individual particles has little effect on the mechanical behaviour of the system. The dynamic behaviour of a particle system is analysed by discretising the time domain. The assumption of infinitesimal displacement theory is valid if the time-step is appropriately small. The particle displacements during each time-step are approximated using a conditionally stable, explicit, central difference time-stepping algorithm. The use of an explicit time-stepping algorithm makes the simulation of a system composed of thousands of particles computationally feasible.

Smooth Particle Hydrodynamics (SPH)

Smooth Particle Hydrodynamics (SPH) is a mesh-free method employed typically in function interpolation and the solution of partial differential equations. Despite the name, SPH can be used to solve non-hydrodynamic problems too. Unlike conventional grid-based techniques, in SPH the domain is discretised by introducing SPH particles, which are free to roam about the function domain. These particles need not have any physical interpretation, serving solely as carriers of field variable information and their derivatives. Associated with any SPH technique is a kernel function, in terms of which all functions and their derivatives at the location of an SPH particle can be expressed as weighted sums over those neighbouring SPH particles within a finite support radius.

The power of SPH is realised in problems of high-deformation, where grid-based techniques fail. However, SPH is not without its limitations. Boundary effects and the imposition of boundary conditions are the subject of ongoing research, leading to numerous techniques and correction schemes.



Development of a semi-mechanistic model using two-way coupled DEM-SPH

Stirred media mills are used in fine and ultra-fine grinding applications in mineral processing. This is because they have been reported to be energy efficient in this application. Typical stirred mills found in industry include the Stirred Media Detritor (SMD), Vertimill, the IsaMill and the Higmill.

In instances where stirred milling technologies have been applied, the particle fluid interactions are not yet fully understood. Due to this limitation, models for stirred mill design and process optimisation have not yet matured. The aim of this study is to develop a model that can be used to predict the performance of vertical stirred mills that use pins as the stirring mechanism, the SMD. This study will use the numerical framework of discrete element method (DEM) coupled with smoothed particle hydrodynamics (SPH) to study the particle-fluid interactions and the dominant role the fluid plays in a vertical stirred mill. DEM is used to model the motion of the ceramic media and SPH is used to model the motion of the slurry (water and the ore material). The computational work is performed using the DEM-SPH code developed at CSIRO in Melbourne.  The code can simulate more than 500 000 particles and can solve for both particle-particle and particle-fluid interactions in a two-way couple.

The proposed outcomes of this project include calculating the power draw for the SMD, as well as the energy spectra for the system with fluid and without the fluid. A parametric study involving varying impeller geometry and position will also be completed. Studying 3 different size SMDs allows for the scale up laws from the smallest size to the largest size to be examined.


Multi-directional inhomogeneous granular suspensions

A complete continuum theory of granular flow remains elusive. In the regime of dense (dry) granular flow, Jop et al. (2006) proposed an empirical visco-plastic theory that appears to work for uni-directional flows across a wide range of geometries. Scaling analysis by GDR MiDi (2004) identifies a single dimensionless control parameter that underpins this empirical rheology: the so-called Inertial number. The situation becomes more complicated in the presence of a suspending fluid. Boyer et al. (2011) showed that in the limit that the stokes number is small, i.e. viscous forces dominate, the notion of the inertial number breaks down in favour of the granular viscous number.

The universality of these ideas has yet to be proven across the full phase space of granular flows. In particular, Cortet et al. (2009) already casts significant doubt on the validity of the visco-plastic theory in the context of multi-directional flows within rotating drums using numerical data obtained via the Non-Smooth Contact Dynamics method. Notwithstanding the lack of experimental confirmation of Cortet’s findings, it is reasonable to expect that universality of Boyers theory remains an open question, particularly in the context of multi-directional flows.

This thesis takes genesis from these uncertainties and lack of quantitative measurements. Using numerical simulations - via coupled Discrete Element method (DEM) and Finite Volume Method - and direct validation measurements of granular suspension flows via Positron Emission Particle Tracking (PEPT), we interrogate the theories of Boyer et al. (2011) - and by association Jop et al. (2006) within rotating drums operated in high and low Stokes regimes with a view extending (and improving) the unified friction law of dense suspensions.


Boyer, F., É. Guazzelli, and O. Pouliquen (2011). Physical Review Letters 107.18, p. 188301.

Cortet, P.-P. et al. (2009). European Physical Letters 88, p. 14001.

GDR MiDi (2004). European Physical Journal E: Soft Matter 14, pp. 341–365.

Jop, P., Y. Forterre, and O Pouliquen (2006). Nature 441(8), pp. 727–730.


Granular flow modelling using GPUs

This project aims to improve the computational time required to run a DEM simulation by parallelising it on a GPU. A full investigation on the scientific validity of these simulations is also required.


Positron emission particle tracking of near gravitational material inside a dense media cyclone

An analysis of trajectory and time averaged data of near gravitational material (NGM) representative coal tracer particles within a magnetite medium passing through a 100mm diameter hydrocyclone imaged with PEPT (positron emission particle tracking) at iThemba labs, Cape Town, is performed. The tracking of neutrally buoyant tracers as well as the other NGM allows an empirical study of the approximate dynamical behaviours of both the dense medium as well as the range (size and density) of particles around the separator’s cut size within a partially imaged volume of the cyclone body lying within the field of view of the detector. The study was commissioned as part of a CFD validation study using empirical data from a real-world lab-appropriate system in order to shed light on the inner workings of hydrocyclones and their separation mechanisms. We employ various optimisation and noise reduction techniques within the analysis and follow a detailed investigation of the velocity fields derived from the raw spatial positions with an emphasis on the velocity profile within the forced and free vortex zones.


Investigating the effect of disk geometry on flow structure in the IsaMill

An IsaMill is a horizontally stirred mill used in fine grinding applications in mineral processing. Improving the efficiency of such systems requires and understanding of the dynamics of the flow field and granular rheology within. With this in mind, the effect of the disk geometry (including holes and configuration) on various kinematic properties such as the velocity and volume fraction distributions is investigated using DEM. The shear stress and power dissipation distributions are calculated from the kinematic properties using the constitutive stress model of granular flows developed by Govender and Pathmathas (2016).


Govender, I., and T. Pathmathas. "A positron emission particle tracking investigation of the flow regimes in tumbling mills." Journal of Physics D: Applied Physics 50.3 (2016): 035601.


Numerical modelling of furnace freeeze lining fracture under thermals load in ilmenite smelting operations

This project is concerned with developing a numerical model that can be used by smelting operators to predict furnace freeze lining fracture during transient furnace operations. Such a model will enable the South African smelter operators to reduce occurrence of freeze lining fractures thus reducing the occurrence of molten slag-refractory interaction leading to longer furnace lifespans. This model will have a direct impact on production and hence earnings of South African smelters. This is especially true for ilmenite smelters and aluminium smelters.

The model is based on continuum damage mechanics, making use of the gradient-enhanced non-local continuum damage framework to model damage of the freeze lining. A one way coupling will be used between a heat transport model and mechanical model to provide a complete damage model.


C.T. Jayasundara, R.Y.Yang, B.Y.Guo, A.B.Yu, I.Govender, A.Mainza, A. van der Westhuizen, J.Rubenstein (2011)

G.B Tupper, I. Govender, A.N Mainza, N. Plint (2013)

I. Govender, P.W. Cleary, and A.N Mainza (2013)

M. Narasimha A.N. Mainza, P.N. Holtham, M.S. Powell and M.S Brennan (2014)

I. Govender, P. Thirunavukkarasu, M.C. Richter, D.N. De Klerk (2014)

I. Govender (2016)

G.B Tupper, I. Govender, D. N De Klerk, M. C. Richter, and A.N Mainza (2016)

I. Govender and T. Pathmathas (2017)

I. Govender, M.C. Richter, A.N. Mainza, D.N. De Klerk (2017)

G.B. Tupper, I. Govender, A.N. Mainza (2017)

O. Ogunmodimu, I. Govender, A Mainza, J-P Franzidis (2017)