Advanced Search
RSS LinkedIn Twitter

Journal Archive

Johnson Matthey Technol. Rev., 2015, 59, (2), 74


How Good is Your Model?

“Just because the results are in colour, it doesn’t mean they are right” (1)

  • By Hugh Stitt*, Michele Marigo and Sam Wilkinson
  • Johnson Matthey Technology Centre,
  • PO Box 1, Belasis Avenue, Billingham, TS23 1LB, UK
  • Tony Dixon
  • Worcester Polytechnic Institute,
  • 100 Institute Road, Worcester, Massachusetts 01609, USA
  • *Email:

Article Synopsis

Models, which underpin all chemical engineering design work, vary widely in their complexity, ranging from traditional dimensionless number correlations through to modern computer based techniques such as computational fluid dynamics (CFD) and discrete element method (DEM). Industrial users require confidence in a model under the conditions it is to be applied in order to use it for design purposes and this can be a reason for slow acceptance of new techniques. This paper explores the validity of models and their validation using a variety of examples from heat transfer, reaction kinetics as well as particle and fluid flow, considering both traditional and modern computational-based approaches.

The examples highlight that when comparing models to experimental data the mathematical form of the equations can contribute to an apparently good ‘fit’ while the actual adjustable parameter values can be poorly predicted; residuals or least squares alone are not a reliable indicator of quality of model fit or of model discrimination. When fitting models to experimental data, confidence in the adjustable parameter values is essential. A finite set of experimental data can fit many different models and often with many sets of parameter values. Not all of these models are of course useful for design. For that purpose it needs to be founded upon the real physics of the system and the adjustable parameters represent real quantities which can be measured, computed or estimated independently. The examples show also the importance of validating a model against more than one output parameter; instances are shown where a too simplistic validation exercise can be misleading.

This paper shows therefore across a range of modelling approaches and applications that extreme care is required when validating a model. Models require validation under the conditions they are to be applied and against more than one output parameter, using appropriate data across appropriate scales and the paper encourages the practice of validating models in order to better persuade industry to adopt more advanced modelling approaches in the future.

1. Introduction

While the way engineers work has transformed massively over the last 50 years, the tools used have not progressed as much as sometimes seems apparent. Despite the extensive use of computational tools to solve the basic design equations, many of the inputs that are used in these calculations are still derived from traditional dimensionless correlations; this is especially true of heat and mass transfer coefficients. In using these correlations we should however be aware of how they represent the data.

Consider heat transfer data for fluidised beds correlated using classical dimensionless number relations. Dimensionless numbers are a critical part of chemical engineering analysis and each one holds its importance because it represents the ratio of two critical energy terms, forces or velocities. Heat transfer data analysis and correlation has classically included in Equations (i)(iv):

Reynolds No.: Re = ρdV/μ (i)

Nusselt No.: Nu = hd/K (ii)

Prandtl No.: Pr = Cp μ / K (iii)

Grashof No.: Gr = g β(TST)D3/v2 (iv)

Rowe (2) takes a heat transfer data set and plots the data on a series of log-log plots, all of which appear to represent a reasonable correlation with a few data points apparently as outliers; a selection of the correlations is shown in Figure 1. This all seems perfectly satisfactory until the raw data set is revealed as being random numbers! To explore how this can be so requires an examination of the structure of dimensionless correlations.

Fig. 1.

Dimensionless number correlations of heat transfer data. St = Stanton number; Re = Reynolds number; Gr = Grashof number; Nu = Nusselt number (Reproduced from Rowe (2))

Dimensionless number correlations of heat transfer data. St = Stanton number; Re = Reynolds number; Gr = Grashof number; Nu = Nusselt number (Reproduced from Rowe (2))

A typical heat transfer correlation is of the form (Equation (v)):

Nu = aRebPrc (v)

or expanding the dimensionless numbers in Equation (vi):

It is noticeable that a key parameter, d, the characteristic dimension, appears on both the right side and left side of the equation (as indeed does the thermal conductivity, K). This means that the correlation is in fact, generalising the variables, yz vs. xz and also plotted on log-log axes. Unless care is taken therefore this approach can be misleading.

The above should not be taken as an indictment of dimensionless correlations. They are an invaluable method of representing data and their relationship with the key independent variables. Nor indeed should it be taken to indicate that they are the only offenders. From the world of reaction kinetics, the Langmuir-Hinshelwood expression format is frequently used to represent adsorption influenced reaction rate data in Equation (vii):

Generalising this equation leads to the observation (3) that this is effectively (Equations (viii) and (ix)):

or rearranging y(1 + bx)q = axp (ix)

As with the heat transfer example, the expression features the same variable on both sides of the equation: it is thus equivalent to Equation (x):

y · xxc (x)

It is hardly surprising therefore that this equation fits many data sets satisfactorily.

To re-emphasise, as with dimensionless correlations, the use of Langmuir-Hinshelwood type kinetic expressions is not wrong; they are based on a specific model of adsorption and reaction phenomena. Rather, users should be aware of the mathematical robustness that arises from the forms of these equations.

Given the typical accepted accuracy of dimensionless transport coefficient expressions, variously cited at ±20% to ±30%, engineers over the years have taken this into account by the judicious use of so-called ‘design margins’. This is what lies behind the common practice of adding 20% to the answer to provide a counterweight to the possibility of under-design due to the inaccuracy of some of the design equations or the input parameters. If this is done by three successive designers then the degree of potential overdesign is now 1.23 = 1.75; a 75% overdesign. If a common cost scaling exponent of 0.7 is applied to this then the cost inflation of the designed plant is 1.750.7 ≈ 1.5; the plant cost is inflated by 50% (4). Use of inaccurate input data and subsequent use of design margins to account for that uncertainty has a significant impact on the final capital cost of the plant.

The approach continues to be used and trusted due to the large amount of knowledge and expertise built up over 50 years in industry. New modelling approaches need to provide us with a similar amount of confidence. That requires validation and cross referencing.

Advanced modelling techniques such as CFD are widely used across many industries (5) such as combustion (6), aerospace design (7) and environmental hydraulics (8, 9). Within the chemical and process industries, for example, it is used in furnace design (10, 11). To extend use into other and new applications requires confidence in the model results.

The industrial designer will want to use the model for scale up and design purposes. It must therefore give the correct answer where it is to be used; not just under conditions where it was derived and tested. This paper will explore the question: how valid does a model need to be and how and at what scale should it be validated?

2. Model Validation Examples in Reaction Kinetics

The accuracy of modelled or fitted kinetic parameters has been much debated due to the robust structure of kinetic expressions, as noted in Equations (viii)–(x). To demonstrate, two key, seminal papers will be reviewed below and their observations highlighted.

2.1 Effect of Adsorption Model

Corma et al. addressed the issue of the quality of the adsorption model and its impact on a kinetic model (12). They selected nineteen adsorption influenced reaction rate data sets from the literature, all of which had been previously modelled by reaction kinetics expressions incorporating a Langmuir isotherm, even where a non-ideality was known to exist. The Langmuir isotherm assumes that the adsorption energy is independent of surface coverage, viz. it is constant. The alternative adsorption isotherm models allow a reduction in adsorption energy with increasing coverage:

  • Langmuir – adsorption energy independent of surface coverage,

  • Temkin – linear decrease with increased coverage,

  • Freundlich – logarithmic decrease with increased coverage.

It is a moot point as to whether these dependencies are real or whether they are simply a reflection of multiple adsorption sites of differing adsorption energies. That notwithstanding, the nineteen data sets were all refitted to reaction rate equations based on the three alternate adsorption models. The outcome is presented in Figure 2.

Fig. 2.

Effect of adsorption isotherm on kinetic model error (Based on data presented by Corma et al. (12))

Effect of adsorption isotherm on kinetic model error (Based on data presented by Corma et al. (12))

The quality of fit is surprisingly good for most cases, remembering that in all cases experimental or data error will probably exceed ±5%. Not only is the same level of fitting obtained in most cases but also few cases actually give a poor fit. This approach therefore simply does not allow model discrimination or specifically in this case discrimination of the appropriate adsorption model.

The authors observe that the effect of the distribution of site energies on the global kinetics is weak (12). Reviewing the same paper Keil observes that “the simple fitting of an [kinetic] equation hides significant features of reactions” (13). Berger et al. note that the mathematical forms of the fitting equations are too robust to enable discrimination (3).

2.2 Reactor Simulation using a Kinetic Model

Another example from reaction kinetics explores the effect of kinetic model inaccuracies or errors on the results of a consequent reactor model (14). The authors used an idealised theoretical reaction model to generate ‘experimental’ data from virtual, in silico experiments. Different kinetic models are fitted to the ‘data’ and these models are then used to simulate a reactor and the reactor model output results compared.

The selected test system was methanol synthesis (Equation (xi)), simplified by neglecting the presence of CO2 in the reaction model.

CO + 2H2 ⇔ CH3OH (xi)

A mechanistic kinetic model, the ‘true’ kinetics were derived from an elementary step and thermodynamic analysis. In silico experiments using this ‘true’ kinetic model were carried out on a statistical experimental design (27 ‘experiments’) with flow, temp, total and partial pressures as prime variables and included outlier experiments. A 5% random error was applied to the results. Nineteen different kinetic models were fitted to these data and all gave a high quality fit based on residual least squares (R2 values all in excess of 95%, many exceeding 99%). These kinetic models were then used in the same reactor design simulation: a steam raising converter with a specified shell side pressure (hence shell side temperature: 210°C). The results were compared using the predicted methanol production rate (Figure 3) and temperature profiles (Figure 4).

Fig. 3.

Predicted methanol production rate based on different kinetic models (Based on data presented by Berty et al. (14))

Predicted methanol production rate based on different kinetic models (Based on data presented by Berty et al. (14))

Fig. 4.

Reactor temperature profiles calculated based on the different kinetic models: (a) acceptable temperature profiles; (b) runaway predicted (Based on data presented by Berty et al. (14))

Reactor temperature profiles calculated based on the different kinetic models: (a) acceptable temperature profiles; (b) runaway predicted (Based on data presented by Berty et al. (14))

The predicted methanol production rates show a wide divergence; the reactor simulation based on the ‘true’ kinetics gave a production rate of 35.5 te h–1. The explanation comes from Figure 4(a), which shows also a wide disparity in the predicted temperature profiles. Figure 4(b) shows that some of the kinetic models in fact predicted a runaway. It should be noted that the shell side temperature was set uncomfortably close to the temperature at which the true kinetics predicted runaway (≈212°C) which makes the results especially sensitive to minor differences in the kinetic model and rate predictions. Looking at this from a different viewing point, were the shell side temperature slightly higher a number of kinetic models with a high quality statistical fit (R2 ≈ 99%) would fail to predict a thermal runaway.

2.3 Discussion

Quality of fit (residual squares: R2) is a poor criterion for model discrimination. The best predictive model in this study was in fact far from best in terms of fit. That is, R2 alone is inadequate and inaccuracies not highlighted by this simple quality of fit criterion may lead to gross misrepresentation of reactor behaviour. Residual squares values indicating accuracies that exceed the experimental data confidence are meaningless statistically; more on this later.

It has long been accepted and is well documented in the reaction engineering literature that a kinetic model must be tested thoroughly before reactor design is fixed and traditionally that has required independent data, measured on a different type of laboratory reactor or pilot unit (15). This is true in other fields as well.

Considering specifically the statistical model fitting exercise, a key issue is that simple univariate model validation is not adequate. In many cases the inherent parameter estimation is an ill-posed mathematical problem: there are more parameters than there are independent plus measured variables. The evaluation of quality of fit should extend to consideration of the statistical quality of the parameter estimates and a check for parameter cross correlation.

More generally, arising also from the ill-posedness of the mathematical problem, fitting a model against one objective function can lead to a poor model. It does not validate the entire model and it does not support other model outputs. Different aspects of a model therefore require independent validation parameterisation.

3. Assessing Parameter Estimate Quality

Fundamentally, fitted model parameters should act as estimates of physically meaningful values. To achieve this, it is of critical importance that the fitted parameters are statistically significant. A common approach is to assess uncertainty in parameter estimates using a 95% confidence interval. The value of this is influenced by correlation between observations, noise in the data and degrees of freedom in the estimation process (16). When a confidence interval is greater than the estimated value, the fitting parameter can be seen as indiscriminate from zero and discounted.

The field of reaction kinetics has utilised more systematic and robust assessment processes. Parameter sensitivity on a local and global basis can be analysed via the Jacobian matrix (i.e. the impact parameter value perturbations have on model responses). Low sensitivity parameters can be systematically discounted from the fitting process (1719). The significance of this with respect to model responses is then checked by the statistical F-test. Cross correlation between parameters is also assessed with these methods, the presence of which impedes reliable parameter estimates.

Principal component analysis (PCA) can also be applied to the sensitivity values of the fitted parameters. PCA provides a ranked series of eigenvalues which contain contributions from the fitting parameters. The most important parameters will have a strong contribution to the largest eigenvalues. Early examples include application to formaldehyde oxidation kinetics (20) and more recently there has been application in micro-kinetic (21) and molecular models (22). Elsewhere, methods to understand parameter range validity (23) and the impact of noise (24) have been incorporated into these approaches.

3.1 Case Study – Paste Extrusion Modelling using the Benbow-Bridgwater Equation

This example, from non-reaction kinetics, will exemplify the importance of considering parameter estimate quality in equation fitting. The Benbow-Bridgwater equation to predict extrusion pressure drop (P, MPa) is described in Figure 5 (25) and contains five input variables (Vext, D0, D, N and L) and six fitting parameters (σ0, α, m, τ0, β and n). As fitting parameters outnumber input variables the problem may be ill-conditioned. While not necessarily ill-posed (where, strictly, the fitted parameters outnumber the quantity of data) this is clearly a situation that requires good structuring of the data set, relevance of the parameters and attention to the independence of the fitted parameters. As an illustration, with a well-designed data set a good fit of a univariate polynomial or other multi-parameter f(x) can be achieved. By contrast with data inadequacy even a simple linear or power law fit can yield significant cross correlation and thus poor confidence intervals of the two fitting parameters. The original version of the Benbow-Bridgwater equation features four parameters, each of which has a distinct physical relevance. The additional parameters, exponents m and n, have however often been used for so-called ‘enhancement of fit’.

Fig. 5.

Schematic of paste extrusion and the associated Benbow-Bridgwater equation

Schematic of paste extrusion and the associated Benbow-Bridgwater equation

To assess the six-parameter model, a 20 point ceramic extrusion dataset featuring five extrusion velocities (Vext) at four different die lengths (D) is utilised. All observations are repeated twice. Fitting the six-parameter model produced a high R2 (0.997), however, as Table I shows, the parameter estimates are not satisfactory. A large confidence interval is seen in β; α could not be estimated and β and n were strongly cross-correlated with one another.

Table I

Parameter Estimates and Confidence Intervals for Four- and Six-parameter Modelsa

Six fitting parametersFour fitting parameters
ParameterEstimated value (2 sig. fig.)95% Confidence interval (2 sig. fig.)Estimated value (2 sig. fig.)95% Confidence interval (2 sig. fig.)
σ0 0.41 ± 0.043 0.44 ± 0.047
α 87 (?) Indeterminate 12 ± 5.7
m 0.72 ± 0.056 1.0 FIXED
τ0 0.0060 ± 0.0020 0.0090 ± 0.0010
β 8.8 ± 6.9 1.8 ± 0.10
n 0.78 ± 0.11 1.0 FIXED

[i] aParameter estimation carried out using Athena Visual Studio v.14.2, Stewart and Associates Engineering Software, Inc

Additional sum-of-squares (model vs. experimental error residual) analysis revealed that the presence of a fifth and sixth parameter in the model simply described experimental errors between the repeats. A t-test of m and n with reference to a value of 1 as the null hypothesis showed both to be insignificant as fitting parameters. Table I reveals that fixing m and n greatly improves the estimation quality of the other parameters, with minimal compromise to fit (R2 = 0.996). A final note on this topic comes from the original authors, “This six parameter fit to data is very frequently inappropriate or not warranted in view of the reliability of the primary data or the accuracy of the prediction required” (26). This important caveat is frequently overlooked.

4. Simulation of Catalyst Loading

How catalyst pellets load into a reactor is critical to the subsequent performance of the reactor. This is especially true of tubular reactors with low tube to particle diameter ratios, typical of steam reforming, oxidation and hydrogenation reactors. To study this experimentally is difficult as it requires forming of the pellets and different pellet shapes if that is a parameter of the study, as well as the characterisation of the bed structure. The latter tends to need advanced and expensive measurement techniques, such as magnetic resonance imaging (MRI) (27) or X-ray micro-tomography (XMT) (28). In view of these difficulties the potential to do this computationally was explored. At the time DEM codes were not sufficiently well developed for them to cope with complex particle shape so an alternative approach was used.

DigiPac® is a particle packing algorithm developed at the University of Leeds, UK, (29, 30) and appeared promising for the simulation of tube loading. In the code, space is discretised (Cartesian basis) and particles are represented by a coherent collection of voxels. The key components are a particle contact detection algorithm with particle behaviour on contact based on a Monte Carlo algorithm which itself has three parameters (rebounding, falling, rotating) which are assigned median values (0–1) on an empirical basis to fit measured validation data.

Preliminary attempts to validate this code against experimental MRI data indicated that although visually the results appeared comparable, it badly under-predicted packing density (31). Constraining the range of the Monte Carlo random inputs led to satisfactory prediction of packing density. Detailed consideration of the actual bed structure, using pellet orientation data, indicated however that the detailed structure predicted was far from correct.

Further improvements to the model included the incorporation of a Hertz-Mindlin contact model – allowing the inclusion of friction into the overall particle contact considerations. Re-validation of the model versions against XMT data showed that this final version not only predicted the packing density correctly but also gave a good representation of the bed structure (characterised again using pellet orientation distributions), Figure 6 (32). Both poured and tapped bed densities and structures were adequately represented.

Fig. 6.

Packing algorithm validation – packing density and orientation (Reprinted with permission from (32). Copyright (2009) American Chemical Society)

Packing algorithm validation – packing density and orientation (Reprinted with permission from (32). Copyright (2009) American Chemical Society)

Summarising, had validation been based only on packing density then that validation would have been false. Only by carrying out a model validation using adequate and sufficient objectives was model verification achieved.

5. Parameterisation and Validation of Discrete Element Method Particle Flow Models

Particle flow simulation using DEM is becoming more popular and more frequently used due to advances in both commercial and open-source codes coupled to the rapid advances in computing power. One of the issues impeding its implementation for industrial design is the selection of the key particle properties to be used as input parameters. Traditionally powder properties have been measured on bulk powders, whereas the DEM code is specifically based on single particles.

5.1 Discrete Element Method Input Parameters

There are many parameters for a standard Hertz-Mindlin-type contact model: shape, size, density, coefficients of static and rolling friction (particle-particle, particle-wall), coefficient of restitution, Poisson ratio, Young's and shear moduli. There are two primary schools of thought on how they should be selected:

  1. Use results of bulk experiments and ‘fit’ parameters to achieve the correct results using, for example the “sandpile” test (33), discharge from a flat bottomed hopper (33, 34), a shear tester (34) or a powder rheometer (35).

  2. Measure the parameters directly using single particles, for example the direct measurement of sliding friction (36), rolling friction (37) and coefficient of restitution (38, 39).

How best these values should be measured, estimated or ‘calibrated’ is a matter of much debate within the DEM community and on which the present authors have recently presented a much broader discussion elsewhere (40).

In the example to be presented here the input parameters were based as far as possible on directly measured values, some of which were taken from the literature.

5.2 Discrete Element Method Simulation and Validation

The Turbula® Mixer is a commercial bench scale powder blender commonly used in the pharmaceutical industry in formulation development. It has a complex motion including rotational, gyrational and axial motion. This complex motion was successfully imported into the commercial EDEM DEM code available from DEM Solutions Ltd, UK, and the representation assured by comparison with positron emission particle tracking (PEPT) data (41). Qualitative comparison with MRI based measurements for powder blending also looked promising (42). A more detailed, quantitative validation is however required.

A detailed comparison of DEM predicted velocity distributions and those measured using PEPT was carried out. A global summary of the results is presented in Figure 7 (43). The overall comparison is reasonable, although a systematic under-prediction of the velocities by about 10% is observed. It is noted that a remarkably similar systematic error in velocity prediction against PEPT data is obtained in an independent study of a paddle mixer (44). These studies appear therefore to validate the DEM model to within a tolerance that is generally acceptable. Does this infer though that predictions of powder blending will also be correct?

Fig. 7.

Comparison of PEPT and DEM velocity data (Based on data in (46))

Comparison of PEPT and DEM velocity data (Based on data in (46))

Axial and radial dispersion coefficients were calculated from the raw velocity data for both DEM and PEPT (43). The comparison of the model and experimental dispersion coefficient values as a function of the Turbula speed is shown in Figure 8 (note that the y-axis scales are different on the axial and radial coefficient plots). DEM was found to over-predict dispersion coefficients by a factor of approximately two.

Fig. 8.

Comparison of PEPT and DEM dispersion coefficients (Redrawn based on data in (43))

Comparison of PEPT and DEM dispersion coefficients (Redrawn based on data in (43))

While comfort may be drawn from the fact that DEM correctly predicts the regime change at approximately 46 rpm, it should be noted that this can be predicted also using simple dimensionless number (Froude number) based analysis (45). Despite significant sensitivity work on the DEM input parameters, it is still not clear why the error in prediction of mixing related parameters such as dispersion coefficient is so poor (46).

It is tempting to blame the offset on the experimental method. However a detailed comparison of PEPT data with results from the well-established and cross-validated particle imaging velocimetry (PIV) method has been reported for studies in turbulent liquid mixing (47). These studies have indicated an excellent agreement overall, but noted a discrepancy at the impeller tip where the highest velocities and rapid direction changes prevail. The origin of these errors is the data averaging techniques used in the reconstruction algorithm that have subsequently been refined (48). The conclusion is that while at the higher Turbula speeds the PEPT may misrepresent the full trajectory of the particles, that under-prediction is smaller than the error observed in the dispersion coefficient prediction. It also would not explain why the error is almost independent of drum speed.

In the context of this communication it is simply important to note that while validation based on a simple velocity analysis indicates that the model could be used predictively, albeit with the caveat that average velocities appear to be over-estimated by around 10%, it does not accurately predict mixer performance. The validation of DEM models needs to be against data that are closely related to the final application.

6. Validation of Computational Fluid Dynamics Models

CFD is now widely used in many fields, including reaction engineering. There has been a steady growth in its application to the ‘high fidelity’ modelling of flow, heat transfer and reaction in packed beds where the packing shapes are fully represented and meshed in detail (4954), for example Figure 9, as opposed to a simplified ‘porous media’ type representation.

Fig. 9.

CFD of flow through a packed tube showing: (a) flow pathlines; (b) velocity vectors; and (c) pressure contours (Reprinted with permission from (56). Copyright (2009) American Chemical Society)

CFD of flow through a packed tube showing: (a) flow pathlines; (b) velocity vectors; and (c) pressure contours (Reprinted with permission from (56). Copyright (2009) American Chemical Society)

First validations of detailed CFD were made against packed tube heat transfer data (55). Simulations of heat transfer for an idealised DT/dP = 2 structured bed of 44 packed spheres compared favourably with experimental radial temperature profiles. Single-parameter validation work in the literature includes: measured pressure drop data (54, 56, 57), heat transfer correlations (58) and classic reaction engineering (pseudo-homogeneous 1D) model results (59). All of these examples attempt to validate the model against a single model output parameter. They thus allow us to state only that the simulation matches the given interrogated output.

The particle processing simulation examples above suggest that one simulation output which matches experimental data cannot be taken to imply that all other results are correct or indeed that the model components are adequate. Full validation of a CFD model requires detailed and multi-objective interrogation of model components. These are detailed below.

6.1 Fluid Flow Fields and Turbulence Model

The published packed bed simulations predominantly use the k-ϵ turbulence closure model. Given the large specific surface area, extent of flow curvature and high particle Reynolds numbers (ReP) it is questionable whether the inherent assumption of isotropic turbulence is valid. A number of studies have compared MRI flow fields with both CFD (60, 61) and with Lattice Bolzmann direct numerical simulation (DNS) simulations (6264). This is however for liquid (aqueous) laminar (or at best transitional) flow and does not help validate the selection of the turbulence closure. The impact of different turbulence closures has been evaluated by comparing the predicted heat transfer with established heat transfer correlations (65, 66), concluding that under their conditions of low ReP (transitional flow) there is no benefit in a more complex model over a single equation form such as that of Spalart Almaras.

Gas-surface contact and (momentum and heat) transfer is a critical aspect of packed tube reactor simulation. To evaluate how best to model this requires a reversion to basics. Comparative Reynolds-averaged Navier-Stokes (RANS) CFD and fundamental DNS based on employing large eddy simulations (LES) of gas flow over a single ‘sphere in a box’ at high Re have recently been carried out to ascertain the correct turbulence model (67) and the results checked against published experimental data. The results show that flow and transport can be accurately calculated using a RANS method with shear-stress transport (SST) k–ω closure provided that the mesh at the particle surface is fine enough and covers most of the boundary layer.

6.2 Particle Contact Points

There are a large number of particle-particle and particle-wall contact points in a packed bed. This is especially a problem where a spherical packing is used. These glancing contact points lead to intractable meshing problems and it is thus necessary to ‘modify’ them to avoid a near-zero contact angle. Contact point representation is vital to different aspects of packed column simulation. It affects packing voidage (hence pressure drop), particle contact area (so conductive heat transfer) and particle spacing (so near surface flows).

A number of strategies are used in the literature:

  • Shrink particles slightly (gap) (55)

  • Expand particles slightly (overlap) (65)

  • Expanded contact (bridge) (68, 69)

  • Flatten the touching curved surfaces (caps) (70)

Dixon, Nijemeisland and Stitt classify and compare the four alternative approaches as shown in Figure 10 (71). They report and compare RANS simulation results for all four approaches and note that all approaches have a significant effect on the simulation results. The effects differ however from one strategy to another. Global approaches significantly affect voidage and thus pressure drop. The removal of the contact point (gaps and caps) distorts particle-particle conductive heat transfer. Overall, bridges appear the best overall solution – but not always. This has been confirmed, but a sensitivity to the diameter of the bridge also noted (72).

Fig. 10.

Classification of contact point modelling approaches (Reproduced from (71) by permission of Elsevier)

Classification of contact point modelling approaches (Reproduced from (71) by permission of Elsevier)

6.3 Heat Transfer Validation

Packed bed heat transfer includes convective heat transfer from wall-to-particles via the gas and conductive heat transfer via particle contact points.

Initially, the conduction model in the CFD was validated against a fundamental analysis based on well-established literature methods (73, 74). The ‘zero flow’ CFD results, incorporating only solid phase conduction, showed good correspondence with the theoretical analysis for a range of input particle thermal conductivities (75) inferring that the conduction model is correct.

In the knowledge that conduction is correctly represented, validation of the heat transfer was made by comparison of experimental heat transfer results with analogous CFD simulations. Heat transfer experiments were carried out in a 98 mm diameter by 0.6 m packed tube at 2200<ReP<27,000 with a heating jacket. The experimental data and comparative CFD results are shown in Figure 11 for the tube packed with ceramic spheres (N = DT/dP = 7.44). The simulation results are in reasonable agreement with the experimental data (76). Given that all other elements of the model have been independently validated, this demonstrates the validity of the convection heat transfer model. Alternatively the convective aspect of gas-particle transfer has been validated by comparison with theoretical or empirical models for heat (77) and mass transfer respectively (78).

Fig. 11.

Comparison of CFD and experimental heat transfer data (Reproduced from Dixon et al. (76))

Comparison of CFD and experimental heat transfer data (Reproduced from Dixon et al. (76))

6.4 Including Intra-pellet Diffusion and the Reaction

Guardo et al. have validated an intra-particle diffusion and reaction model against experimental data for different sized catalyst particles (79). For the steam reforming case, the above sections report the successive validation of key elements of the overall CFD model for flow and heat transfer. Attention can now therefore turn to assessing the inclusion of intra-particle diffusion and reaction terms in an overall reaction model (80).

Intraparticle diffusion was modelled assuming a uniform porosity and Fickian effective diffusivity evaluated based on the dusty gas model assuming pressure variation in the pellet is small compared to total external pressure (81) and user defined scalars for the species balances. The reaction rate terms were input by user-defined code based on the steam reforming kinetics of Hou and Hughes (82).

For validation of this last model component a more structured pellet arrangement was used (83). The set up used for the experiments and simulation is shown in Figure 12 and used a ‘string’ of commercial pellets (Johnson Matthey Mini-Q 57-4). The reaction zone was six catalyst pellets, with a further six inert pellets in the feed and exit zone. The holes in the pellets were used effectively as thermo-wells and thus flow was only on the quasi-annular exterior.

Fig. 12.

Validation of CFD model including reaction and diffusion: (a) experimental set up; (b) graphical results from CFD (Reproduced from (83) by permission of Elsevier)

Validation of CFD model including reaction and diffusion: (a) experimental set up; (b) graphical results from CFD (Reproduced from (83) by permission of Elsevier)

Temperature profile measurements of a pellet string under methane steam reforming reaction conditions are shown in Figure 13. The simulation results are in reasonable agreement with the experiments, with an error in the order of only 2°C–5°C for the temperature predictions. Given that there are no adjustable parameters in this model other than those in the kinetic expression this seems to be a remarkably good result. This is a result, of course, of the fact that individual model components have been independently modified and validated and that in this final step only the diffusion and reaction terms represent any uncertainty.

Fig. 13.

Comparison of CFD to experiment for: (a) methane conversion at different inlet gas temperatures, methane gas feed; (b) inside active particle temperature (K) for the different inlet gas temperatures at S:C ratio = 4.5, methane gas feed (Redrawn based on data in (83))

Comparison of CFD to experiment for: (a) methane conversion at different inlet gas temperatures, methane gas feed; (b) inside active particle temperature (K) for the different inlet gas temperatures at S:C ratio = 4.5, methane gas feed (Redrawn based on data in (83))

7. Conclusions

Advanced modelling techniques such as CFD and DEM are transforming the way chemical engineers carry out process and equipment research and evaluation. They provide massively enhanced levels of information and understanding on many complex process operations. In order for them to achieve widespread use as primary design tools it is essential to achieve the same levels of confidence that exist in traditional design methods and approaches.

This paper has addressed the level and scope of validation required for models, considering some traditional modelling approaches as well as modern computationally intensive models. In considering kinetic and rheology models, the examples have highlighted that when fitting models to experimental data, confidence in the adjustable parameter values is essential if that constituent model is to be reliable when applied to a design. Reliance on residuals (least squares) alone will not reliably parameterise or validate a model.

In developing or fitting a model one should also bear in mind the sagacious words of the physicist and statistician George Box (84). Often recalled only for the aphorism “All models are wrong; some are useful”, the extended version reminds us that we cannot always improve a model by the addition of parameters and terms: “Since all models are wrong the scientist cannot obtain a “correct” one by excessive elaboration. On the contrary following William of Occam he should seek an economical description of natural phenomena. Just as the ability to devise simple but evocative models is the signature of the great scientist so over-elaboration and over-parameterization is often the mark of mediocrity.”

In the specific context of this paper it is noted that a finite set of experimental data can fit many different models and often with many sets of adjustable parameter values. Not all of these models will be useful for design. A model derived from the real physics of the system, where the parameters represent real quantities which can potentially be measured, computed or estimated independently from the data set, is much more likely to be useful for design. Furthermore, in accordance with Box's observations above, the addition of parameters with limited physical meaning to achieve an enhancement of fit is not beneficial to the predictive quality of the model or indeed to the estimation of the true physical fitting parameters.

The above observations are even more important for models of increased complexity, typified by discrete particle and CFD models. Two examples of particle flow modelling were presented, demonstrating that although the respective models matched well the initial data set more detailed consideration showed that other facets or model outputs were very poorly predicted.

Single-objective validation is not enough. It is necessary to evaluate the veracity of the underpinning model components and closures. Detailed validation of detailed models can be carried out using appropriate data and at appropriate scales and this is to be encouraged if the industry is to adopt successful modelling approaches in future.

This paper is based on a Keynote Lecture given at the 9th European Congress in Chemical Engineering, Den Haag, Netherlands, 21st–24th April 2013.



  1.  M. P. Dudukovic, F. Larachi and P. L. Mills, Chem. Eng. Sci., 1999, 54, (13–14), 1975 LINK
  2.  P. N. Rowe, ‘The Correlation of Engineering Data’, The Chemical Engineer, March 1963, p. CE 69. Republished: The Chemical Engineer, May 2001, p. 54
  3.  R. J. Berger, E. H. Stitt, G. B. Marin, F. Kapteijn and J. A. Moulijn, CatTech, 2001, 5, (1), 36 LINK
  4.  E. H. Stitt, ‘Past and Future in Multiphase Reactor Design: A View from Industry’, 4th International Symposium on Catalysis in Multiphase Reactors, CAMURE-4, Lausanne, Switzerland, 22nd–25th September, 2002
  5.  A. E. Holdo, Int. J. Appl. Mech. Eng., 2002, 7, (1), 85
  6.  P. J. Stopford, Appl. Math. Modelling, 2002, 26, (2), 351 LINK
  7.  K. Fujii, Prog. Aerosp. Sci., 2005, 41, (6), 455 LINK
  8.  S. N. Lane, K. F. Bradbrook, K. S. Richards, P. A. Biron and A. G. Roy, Geomorphology, 1999, 29, (1–2), 1 LINK
  9.  “Computational Fluid Dynamics: Applications in Environmental Hydraulics”, eds. P. D. Bates, S. N. Lane and R. I. Ferguson, John Wiley & Sons Ltd, West Sussex, UK, 2005 LINK
  10.  “Computational Fluid Dynamics in Industrial Combustion”, eds. C. E. Baukal Jr., V. Gershtein and X. J. Li, Industrial Combustion, CRC Press LLC, Florida, USA, 2001
  11.  A. Habibi, B. Merci and G. J. Heynderickx, Computers Chem. Eng., 2007, 31, (11), 1389 LINK
  12.  A. Corma, F. Llopis, J. B. Monton and S. W. Weller, Chem. Eng. Sci., 1988, 43, (4), 785 LINK
  13.  F. J. Keil, Stud. Surf. Sci. Catal., 2001, 133, 41 LINK
  14.  J. M. Berty, S. Lee, F. Szeifert and J. B. Cropley, Chem. Eng. Commun., 1989, 76, (1), 9 LINK
  15.  F. M. Dautzenberg, ‘Ten Guidelines for Catalyst Testing’, in “Characterization and Catalyst Development: An Interactive Approach”, eds. S. A. Bradley, M. J. Gattuso and R. J. Bertolacini, ACS Symposium Series, Vol. 411, American Chemical Society, Washington DC, USA, 1989 LINK
  16.  J. Sjöblom, ‘Parameter Estimation in Heterogeneous Catalysis’, PhD Thesis, Department of Chemical and Biological Engineering, Chalmers University of Technology, Sweden, 2009 LINK
  17.  J. Jansson, ‘Studies of Catalytic Low-temperature CO Oxidation over Cobalt Oxide and Related Transition Metal Oxides’, PhD Thesis, Chalmers University of Technology, Sweden, 2002
  18.  A. S. Quiney and Y. Schuurman, Chem. Eng. Sci., 2007, 62, (18–20), 5026 LINK
  19.  S. K. Wilkinson, ‘Reaction Kinetics in Formulated Industrial Catalysts’, Eng. D. Thesis, School of Chemical Engineering, University of Birmingham, UK, 2014 LINK
  20.  S. Vajda, P. Valko and T. Turányi, Int. J. Chem. Kinet., 1985, 17, (1), 55 LINK
  21.  A. B. Mhadeshwar and D. G. Vlachos, Catal. Today, 2005, 105, (1), 162 LINK
  22.  S. Raimondeau, P. Aghalayam, A. B. Mhadeshwar and D. G. Vlachos, Ind. Eng. Chem. Res., 2003, 42, (6), 1174 LINK
  23.  J. Song, G. Stephanopoulos and W. H. Green, Chem. Eng. Sci., 2002, 57, (21), 4475 LINK
  24.  G. Meinrath, C. Ekberg, A. Landgren and J. O. Liljenzin, Talanta, 2000, 51, (2), 231 LINK
  25.  J. Benbow and J. Bridgwater, “Paste Flow and Extrusion”, Oxford Series on Advanced Manufacturing, Clarendon Press, Oxford, UK, 1993
  26.  J. J. Benbow, S. H. Jazayeri and J. Bridgwater, Powder. Tech., 1991, 65, (1–3), 393 LINK
  27.  A. J. Sederman and L. F. Gladden, Chem. Eng. Sci., 2001, 56, (8), 2615 LINK
  28.  D. Toye, P. Marchot, M. Crine and G. L’Homme, Meas. Sci. Technol., 1996, 7, (3), 436 LINK
  29.  X. Jia and R. A. Williams, Powder Tech., 2001, 120, (3), 175 LINK
  30.  R. Caulkin, A. Ahmad, M. Fairweather, X. Jia and R. A. Williams, Comput. Chem. Eng., 2009, 33, (1), 10 LINK
  31.  C. Xu, X. Jia, R. A. Williams, E. H. Stitt, M. Nijemeisland, S. El-Bachir, A. J. Sederman and L. F. Gladden, Comput. Model. Eng. Sci., 2008, 23, (2), 117 LINK
  32.  R. Caulkin, X. Jia, C. Xu, M. Fairweather, R. A. Williams, H. Stitt, M. Nijemeisland, S. Aferka, M. Crine, A. Léonard, D. Toye and P. Marchot, Ind. Eng. Chem. Res., 2009, 48, (1), 202 LINK
  33.  Y. Li, Y. Xu and C. Thornton, Powder Technol., 2005, 160, (3), 219 LINK
  34.  ‘The Beginning of a New Era in Design: Calibrated Discrete Element Modelling’, Australian Bulk Handling Review, 2011, 16, (6), p. 14
  35.  R. Bharadwaj, W. R. Ketterhagen and B. C. Hancock, Chem. Eng. Sci., 2010, 65, (21), 5747 LINK
  36.  B. C. Hancock, N. Mojica, K. St. John-Green, J. A. Elliott and R. Bharadwaj, Int. J. Pharmaceutics, 2010, 384, (1–2), 39 LINK
  37.  W. R. Ketterhagen, R. Bharadwaj and B. C. Hancock, Int. J. Pharmaceutics, 2010, 392, (1–2), 107 LINK
  38.  R. Bharadwaj, C. Smith and B. C. Hancock, Int. J. Pharmaceutics, 2010, 402, (1–2), 50 LINK
  39.  H. Dong and M. H. Moys, Min. Eng., 2003, 16, (6), 543 LINK
  40.  M. Marigo and E. H. Stitt, KONA Powder Part. J., 2015, 32, 236 LINK
  41.  M. Marigo, D. L. Cairns, M. Davies, M. Cook, A. Ingram and E. H. Stitt, Comput. Model. Eng. Sci., 2010, 59, (3), 217 LINK
  42.  M. Marigo, D. L. Cairns, M. Davies, A. Ingram and E. H. Stitt, Powder Technol., 2011, 212, (1), 17 LINK
  43.  M. Marigo, M. Davies, T. Leadbeater, D. L. Cairns, A. Ingram and E. H. Stitt, Int. J. Pharmaceutics, 2013, 446, (1–2), 46 LINK
  44.  A. Hassanpour, H. Tan, A. Bayly, P. Gopalkrishnan, B. Ng and M. Ghadiri, Powder Technol., 2011, 206, (1–2), 189 LINK
  45.  C. Mayer-Laigle, C. Gatumel and H. Berthiaux, Chem. Eng. Res. Des., 2015, 95, 248 LINK
  46.  M. Marigo, ‘Discrete Element Method Modelling of Complex Granular Motion in Mixing Vessels: Evaluation and Validation’, Eng.D. Thesis, School of Engineering, University of Birmingham, UK, 2012 LINK
  47.  P. Pianko-Oprych, A. W. Nienow and M. Barigou, Chem. Eng. Sci., 2009, 64, (23), 4955 LINK
  48.  F. Chiti, S. Bakalis, W. Bujalski, M. Barigou, A. Eaglesham and A. W. Nienow, Chem. Eng. Res. Des., 2011, 89, (10), 1947 LINK
  49.  A. G. Dixon, M. Nijemeisland and E. H. Stitt, Adv. Chem. Eng., 2006, 31, 307 LINK
  50.  R. K. Reddy and J. B. Joshi, Chem. Eng. Res. Des., 2008, 86, (5), 444 LINK
  51.  T. Atmakidis and E. Y. Kenig, Chem. Eng. J., 2009, 155, (1), 404 LINK
  52.  P. Magnico, AIChE J., 2009, 55, (4), 849 LINK
  53.  F. Augier, F. Idoux and J. Y. Delenne, Chem. Eng. Sci., 2010, 65, (3), 1055 LINK
  54.  M. J. Baker and G. R. Tabor, Comput. Chem. Eng., 2010, 34, (6), 878 LINK
  55.  A. G. Dixon and M. Nijemeisland, Ind. Eng. Chem. Res., 2001, 40, (23), 5246 LINK
  56.  H. Bai, J. Theuerkauf, P. A. Gillis and P. M. Witt, Ind. Eng. Chem. Res., 2009, 48, (8), 4060 LINK
  57.  A. Guardo, M. Coussirat, M. A. Larrayoz, F. Recasens and E. Egusquiza, Ind. Eng. Chem Res., 2004, 43, (22), 7049 LINK
  58.  A. Guardo, M. Coussirat, F. Recasens, M. A. Larrayoz and X. Escaler, Chem. Eng. Sci., 2006, 61, (13), 4341 LINK
  59.  M. Kuroki, S. Ookawara, D. Street and K. Ogawa, ‘High-fidelity CFD Modeling of Particle-to-fluid Heat Transfer in Packed Bed Reactors’, Proceedings of European Congress of Chemical Engineering (ECCE6), Paper 1102, Copenhagen, Denmark, 16th–20th September, 2007 LINK
  60.  P. R. Gunjal, V. V. Ranade and R. V. Chaudhari, AIChE J., 2005, 51, (2), 365 LINK
  61.  D. J. Robbins, M. S. El-Bachir, L. F. Gladden, R. S. Cant and E. von Harbou, AIChE J., 2012, 58, (12), 3904 LINK
  62.  B. Manz, L. F. Gladden and P. B. Warren, AIChE J., 1999, 45, (9), 1845 LINK
  63.  M. D. Mantle, A. J. Sederman and L. F. Gladden, Chem. Eng. Sci., 2001, 56, (2), 523 LINK
  64.  H. Freund, J. Bauer, T. Zeiser and G. Emig, Ind. Eng. Chem. Res., 2005, 44, (16), 6423 LINK
  65.  A. Guardo, M. Coussirat, M. A. Larrayoz, F. Recasens and E. Egusquiza, Chem. Eng. Sci., 2005, 60, (6), 1733 LINK
  66.  M. Coussirat, A. Guardo, B. Mateos and E. Egusquiza, Chem. Eng. Sci., 2007, 62, (23), 6897 LINK
  67.  A. G. Dixon, M. E. Taskin, M. Nijemeisland and E. H. Stitt, Comput. Chem. Eng., 2011, 35, (7), 1171 LINK
  68.  S. Ookawara, M. Kuroki, D. Street and K. Ogawa, ‘High-fidelity DEM-CFD Modeling of Packed Bed Reactors for Process Intensification’, Proceedings of European Congress of Chemical Engineering (ECCE6), Paper 1105, Copenhagen, Denmark, 16th–20th September, 2007
  69.  M. Kuroki, S. Ookawara and K. Ogawa, J. Chem. Eng. Jpn, 2009, 42, Suppl. 1, s73 LINK
  70.  T. Eppinger, K. Seidler and M. Kraume, Chem. Eng. J., 2011, 166, (1), 324 LINK
  71.  A. G. Dixon, M. Nijemeisland and E. H. Stitt, Comput. Chem. Eng., 2013, 48, 135 LINK
  72.  S. S. Bu, J. Yang, M. Zhou, S. Y. Li, Q. W. Wang and Z. X. Guo, Nucl. Eng. Des., 2014, 270, 21 LINK
  73.  P. Zehner and E. U. Schlünder, Chem. Ing. Technik, 1970, 42, (14), 933 LINK
  74.  P. Zehner and E. U. Schlünder, Chem. Ing. Technik, 1972, 44, (23), 1303 LINK
  75.  A. G. Dixon, A. K. Gurnon, M. Nijemeisland and E. H. Stitt, Int. Commun. Heat Mass Trans., 2013, 42, 1 LINK
  76.  A. G. Dixon, G. Walls, H. Stanness, M. Nijemeisland and E. H. Stitt, Chem. Eng. J., 2012, 200–202, 344 LINK
  77.  A. Guardo, M. Coussirat, F. Recasens, M. A. Larrayoz and X. Escaler, Chem. Eng. Sci., 2006, 61, (13), 4341 LINK
  78.  F. S. Mirhashemi, S. H. Hashemabadi and S. Noroozi, Int. Commun. Heat Mass Trans., 2011, 38, (8), 1148 LINK
  79.  A. Guardo, M. Casanovas, E. Ramírez, F. Recasens, I. Magaña, D. Martínez and M. A. Larrayoz, Chem. Eng. Sci., 2007, 62, (18–20), 5054 LINK
  80.  A. G. Dixon, M. E. Taskin, M. Nijemeisland and E. H. Stitt, Ind. Eng. Chem. Res., 2010, 49, (19), 9012 LINK
  81.  R. H. Hite and R. Jackson, Chem. Eng. Sci., 1977, 32, (7), 703 LINK
  82.  K. Hou and R. Hughes, Chem. Eng. J., 2001, 82, (1–3), 311 LINK
  83.  M. Behnam, A. G. Dixon, P. M. Wright, M. Nijemeisland and E. H. Stitt, Chem. Eng. J., 2012, 207–208, 690 LINK
  84.  G. E. P. Box, J. Am. Stat. Assoc., 1976, 71, (356), 791 LINK

The Authors

E. Hugh Stitt is a Scientific Consultant at Johnson Matthey Technology Centre, Billingham, UK. He is a Visiting Professor at the University of Birmingham, UK, and at Queen's University Belfast, Fellow of the Institution of Chemical Engineers and a Fellow of the Royal Academy of Engineering. He has 25 years of industrial research experience across a variety of themes related to catalytic reaction engineering and catalysts with over 100 refereed publications.

Michele Marigo is a Principal Scientist at Johnson Matthey Technology Centre. He obtained an undergraduate degree with Master's in Mechanical Engineering from University of Padua, Italy, and doctorate in Chemical Engineering (EngD) from University of Birmingham. Michele's expertise includes materials science and particle engineering, discrete element modelling (DEM) and finite element modelling (FEM).

Sam K. Wilkinson is a Senior Scientist at Johnson Matthey Technology Centre. He obtained his MEng and his EngD, both in Chemical Engineering from the University of Birmingham. Sam's specialist areas include developing steady state and transient (micro) kinetics for business applications and in new research in addition to multi-scale modelling, statistics, parameter estimation and design of experiments.

Anthony G. Dixon is a Chemical Engineering Professor at Worcester Polytechnic Institute, Massachusetts, USA. His research is in fixed bed and membrane reactors, mathematical modelling and CFD. He has co-authored over 100 refereed journal papers. He is a Fellow of American Institute of Chemical Engineers (AIChE). He has edited a recent volume of Advances in Chemical Engineering and is Executive Editor for Reaction Engineering and Catalysis for Chemical Engineering Science.

Read more from this issue »