Confidence Intervals for Richardson Extrapolation in Solid Mechanics

Models based on partial differential equations are commonplace in science and engineering. Due to the unavailability of analytical solutions in general, such models are converted to discrete forms, and a computer solution is sought. Discretization errors are an unavoidable consequence of the numerical scheme (e.g., finite difference, finite element, or finite volume methods) used to obtain a discrete approximation of the continuous equations. Solution verification is a process by which the discretization error of the quantities of interest is estimated [1]. The goal is to estimate both the error and the uncertainty associated with that estimate [2].

An important characteristic of discretization schemes is the convergence rate, which is the exponent in the power law term relating the numerical truncation error to the value of a parameter associated with the discretization, usually given by the size of the finite element. The factor multiplying this term gives a measure of the overall error of a given scheme. The standard method by which to estimate the error is systematic mesh refinement. Richardson extrapolation (usually traced to Ref. [3]) is currently one of the (if not the) most reliable method for the prediction of numerical uncertainty of the error estimation (c. f. Ref. [4]).

The grid convergence index (GCI) of Roache was perhaps the earliest attempt to quantify the numerical uncertainty associated with inferred convergence parameters [5,6]. Safety factors were used to approach a 95% confidence level. This approach was extended to the correction factor method of Ref. [7]. Xing and Stern later delved deeper, with somewhat pessimistic assessment of the reliability of then current approaches [8].

Heuristics to estimate the numerical uncertainty associated with a set of computed results have been proposed by Eca and Hoekstra [9]. There appears to be an assumption that the underlying numerical scheme has a known theoretical convergence rate, and that this rate is realized in simulations.

The robust verification analysis (RVA) method [10] was developed as an alternative to the GCI method. One of the goals was to allow for verification based on as few result sets as possible. Another goal was to avoid the use of safety factors. A distinguishing characteristic was the calculation of the numerical uncertainty from median-based statistics from a number of constrained-optimization fits the numerical error model. These fits allowed constraints to introduce expert judgment. Various fitting models were used, such as multiple residual norms, various weightings of the data, regularization methods, fits with fixed orders of convergence, and fits for various subsets of the response data. This method also permits (using statistics of multiple fit estimates) the use of the converged result (extrapolation) as the best guess, and estimates the uncertainty about that value, in contrast to the GCI method which estimates the uncertainty about the most resolved simulation value.

The formulation of the RVA approach [10] included regularization methods in order to treat under-determined cases: fits to error models performed on a less-than minimal set of response data (such as a two-point Richardson extrapolation). The under-determined fitting was partially overcome by including regularized optimization methods (Tikhonov, LASSO - least absolute shrinkage and selection operator). Unfortunately, the results of the under-determined fitting would to some degree depend on the regularization employed, and the details of the regularization were not provided in the paper.

The RVA [10] is evidently very powerful. Its sophistication comes with a price tag, though: considerable constrained-optimization software machinery is required, and unless this comes prepackaged, the user also needs to acquire substantial skills to put together the individual software components, and also the statistical sophistication to interpret the results.

There is a newer approach that builds on the RVA [11]. It offers the following critique of the RVA:

The approach called stochastic robust extrapolation error quantification, or StREEQ, was offered as a replacement. However, StREEQ was mainly developed to deal with results from stochastic simulations, even though deterministic simulations were also targeted. The claim was made that it was a robust tool for solution verification of stochastic and deterministic simulation results. The statistical robustness was due to the use of a bootstrap resampling of the response data to increase the extrapolation data pool. The inherent statistical power of the estimation procedure needs to be matched by the statistical know-how on the part of the user.

There GCI approach had been applied to solid and structural mechanics problems in Ref. [12]. A related development reports conservative error estimates based on the notions of Richardson extrapolation in stress-concentration problems [13]. Both of these approaches are compared with the present methodology in our examples.

We offer here a very simple Richardson-extrapolation procedure, which only requires the solution of a single nonlinear equation and trivial median-based statistical processing. The result is an estimate of the true solution and the convergence rate, both with estimates of the 95% confidence interval. In Sec. 1 we first introduce a classical Richardson extrapolation from results obtained using three different grids (meshes). Then we continue by introducing the estimation of the uncertainty of the extrapolation. We identify the sources of the uncertainty, and we propose to obtain a solution on a fourth grid (mesh) in order to build up a dataset of extrapolations to aid in the quantification of uncertainty. The dataset consists of extrapolations based on samples constructed from the four grids. Robust median-based statistics are applied to the dataset to extract the best guess of the true solution, and the median absolute deviation from the median is employed as a measure of the uncertainty of the best guess. In Sec. 2 we illustrate the proposed methodology with solid-mechanics examples: a stress concentration solved with an assumed-strain 3D finite element model, a shell benchmark solved with a shear-gap triangular facet finite element, and a shear-rigid plate problem solved with a commercial finite element code. The computed confidence intervals are discussed in the light of comparisons with established procedures, such as the GCI. Section 3 considers complications due to the failure to obtain data from which extrapolation is possible. Possible extensions are also discussed.

Richardson extrapolation (usually traced in history to Ref. [3]) is currently the most reliable method for the prediction of numerical uncertainty of results from discretized models (c.f. Ref. [4]). As described aptly in [14], the two basic requirements for its applicability are:

1. The grid refinement is defined by a single parameter, the “element size”; and

2. The solutions on the grid are in the asymptotic range, meaning that a single term of a Taylor series expansion dominates the error. (We shall use the terminology “mesh” and “grid” interchangeably in what follows.)

First, we outline the classical Richardson extrapolation from results obtained on three grids (meshes).

For instance, mesh M₁ would be produced with element size $h1(x)=αh0(x)$⁠, mesh M₂ with element size $h2(x)=αh1(x)=α2h0(x)$⁠, and so on.

In the above, the refinement factor was assumed constant when going from mesh M₀ to M₁ to M₂ and so on. Strictly speaking, keeping the refinement factor constant is not necessary. We may choose a different refinement factor for each new mesh. Thus, mesh M₁ could be produced with element size $h1(x)=α0h0(x)$⁠, mesh M₂ with element size $h2(x)=α1h1(x)=α0α1h0(x)$⁠, and so on.

and the estimate of the true error for any element size can then be calculated from Eq. (4). Obviously, with the estimate of the true error, we also have the estimate of the true solution $qex$⁠.

It is common knowledge that the three-grid procedure from Sec. 1.1 cannot be used to quantify the uncertainty of the true solution estimate: four grids are needed [9]. The reason is that there are three unknown parameters, and the solutions from three grids supply just enough equations to solve for the three unknowns. The result of the extrapolation is a number—the estimate of the true solution—without any indication of the uncertainty of this number.

Where does this uncertainty come from? We can identify at least three sources:

Four or more grids have been used in Ref. [9]. The error was estimated by fitting the expansions of the form (4) to the data in the least-squares sense. The selection of the best error estimate was founded on the standard deviation of the fits. The error estimate was formed into an uncertainty using a safety factor. The safety factor was proposed to depend on the observed convergence rate and the magnitude of the standard deviation of the fit.

Robust verification analysis was introduced by Rider et al. [10]. Results from simulations with multiple grids were treated with optimization instead of with linear or nonlinear regression. The measure of central tendency was selected to be the median instead of the mean. This choice was motivated by the median's robustness to outliers [15]. In fact, optimization in the L₁ norm automatically removed outliers from the data. The optimization started from results obtained with fine meshes, and new estimates for course meshes were progressively added and the optimization fit was recomputed for each additional result. From the considerable machinery required for the implementation of the RVA procedure, we derive a readily executable process that does not require anything beyond solving the Eqs. (9), (10), and (4).

One of the main ideas in this work is to form three tuples of inputs to the extrapolation process (9), (10), and (4) by sampling from the results obtained using four meshes. Table 1 shows the sampling scheme. The four meshes with relative element size γ_j result in computed numbers $qγj$⁠. Then the sampling in row one of the table produces the extrapolation $qex(1)$⁠, as indicated by the bullets. The four samplings thus result in four predictions $qex(k)$⁠, where $k=1,2,3,4$⁠. These are subsequently used in the estimation of uncertainty.

Here b is needed to make the estimator consistent for the parameter of interest. In the case of the parameter with a normal distribution, we need to set b = 1.4826 [16,17].

to produce the final estimate of the true solution. Here $0≤W≤1$ is an arbitrary weighting factor. In this paper we always take W = 2/3, meaning we give more weight to the estimate $qex(4)$ then to the median $q̃ex$⁠.

where factor 2 corresponds to the two-standard-deviations confidence interval for normal distributions. It bears emphasis that the estimate (14) does not represent bounds on the true solution, merely a band centered on $qex*$ where the true value can be expected to be contained with a high probability.

The convergence rates computed from Eq. (15) should be subject to expert judgment. In particular, convergence rates which are unreasonably high should lead the analyst to question both the best guess and the confidence interval (15), as well as the best guess of the solution and its confidence interval (14). However, it would be difficult to incorporate expert judgment a priori, as it is sometimes suggested in fluid dynamics [10]: Numerical models in solid mechanics can work correctly, yet converge at unexpected convergence rates. This was reported also in computational fluid dynamics [18], but in solid mechanics models, there are some additional difficulties: various forms of locking (volumetric, shear, Poisson), both for reduced-dimension-models (beams, shells), and for continuum models [19]. Furthermore, various kinds of singularities will affect the convergence rate of the quantity of interest in the primary domain, even when the singularities are present in the secondary domain (outside of the region of interest) [20]. For example, various kinds of boundary conditions, such as hard simple supports for shells, will produce singular force resultants; clamped boundary conditions introduce stress singularities, which even when not inspected directly may pollute the stress results elsewhere; etc.

where $qγ3$ is the solution on the finest grid. The convergence rate in solid and structural mechanics typically cannot be guessed at the theoretical value, which means that it needs to be computed as discussed above. Reference [12] describes the estimation (17) thus

The error-band aspect is consistent with our confidence interval (14). Therefore, in the examples in this section, we compare our approach with the GCI-based estimation.

We will also compare with this alternative estimation, which we will designate by the initials of the authors as BSR.

The following problem was considered in Ref. [21]: An infinite slab 0.1 m thick with a stress-free circular hole of radius $R1=0.1$ m, under far-field unidirectional tensile loading P = 1.0 MPa, was considered. The Young's modulus was E = 2.4 MPa and Poisson's ratio was $ν=0.49995$ (i.e., the material was nearly incompressible). Plane-strain conditions were simulated by fixing the faces parallel to the plane of the slab in the direction normal to them. A region of radius $R2=0.4$ m around the stress-free hole was modeled. Owing to the two-plane symmetry of the model, only a quarter of the geometry was considered for analysis.

where r is the distance from the center of the stress-free hole and θ is the anticlockwise angle with respect to the horizontal X-axis. These stresses were used to form the traction boundary conditions on the outer cylindrical surface.

The target quantity was the direct stress, σ_x at point A shown in Fig. 1. From the analytical solution, one can evaluate $σx|A=3P$⁠.

A high-performance mean-strain hexahedral finite element was used. The element was stabilized (in the sense of correcting the rank deficiency) by setting up a stabilization energy that was sampled either with the full quadrature rule or with the mean-strain quadrature [23]. The element is robust with respect to the near incompressibility of the material. Hence, we can hope to realize rate of convergence close to the theoretical one of 1.0.

The problem was solved using six meshes. Thus we were provided with the opportunity of applying the uncertainty quantification as outlined above three times, as we could form three tuples of four sequences of refined meshes. The coarsest mesh is shown in Fig. 1, and successive refinements have been produced by isotropic bisection, followed by projection onto the cylindrical surfaces. Table 2 shows the meshes used, and the computed stress at point A.

Table 3 presents the three extrapolations performed on the four tuples. For instance, the first row shows extrapolation from the meshes 0, 1, 2, and 3: Extrapolation from $1.96424,2.35360,2.62957$ yields 3.30123 (with convergence rate 0.49659), extrapolation from $1.96424,2.62957,2.79966$ results in the estimate 3.11464 (with convergence rate 0.62294), extrapolation from $1.96424,2.35360,2.79966$ yields 3.18620 (with convergence rate 0.55350) and finally extrapolation from $2.35360,2.62957,2.79966$ gives the estimate 3.07290 (with convergence rate 0.69821). The median of the four predictions is $q̃ex=3.15042$⁠. Thus the best guess according to Eq. (13) is $(2/3)3.07290+(1/3)3.15042=3.09874$⁠. The details concerning the calculation of the uncertainty and an analogous calculation of the convergence rate and its uncertainty are omitted, we just note that the stress is predicted to be $σx=3.09874±0.16798$ MPa.

The other two rows of Table 3 are completed analogously. The convergence rate can be seen to increase with refinement from $0.66154±0.12634$ to 0.89319±0.04918.

Figure 2 illustrates the uncertainty quantification of the predicted stress $σx|A$ that corresponds to the data in Table 3. The predicted stress is shown with a filled dot, and the uncertainty interval is shown with whiskers. As can be seen, the uncertainty interval includes the exact answer of 3.0 MPa.

Now we present the estimation based on the GCI. Note that the six computed results in Table 2 allow for three extrapolations, but four estimates using the GCI. As shown in Fig. 3, the extrapolated value of the stress is above the true value. On the contrary, the calculated values of the stress approach the true value from below. The error band of the GCI contains the true solution, but the error band is approximately five times wider than the confidence interval obtained with our approach.

The estimation based on the formula (19) (BSR) is illustrated in Fig. 4. Similarly to the GCI approach, the true value is approached from below. The error band is somewhat narrower than for the GCI.

The Scordelis-Lo cylindrical roof is one of the original obstacle-course benchmarks and was discussed at length by Krysl and Chen [24]. The target quantity is the deflection of the midpoint of the free, straight, edge. The robust flat-facet triangular shear-flexible finite element [25] produced the limit value of the deflection $−0.30148′$ (in the direction of the load), extrapolated from extremely refined meshes (512 element edges for the finest mesh).

The meshes with a given number of elements per side were used (c.f. Table 4), i.e., a total of seven meshes were used, which made it possible to use the extrapolation procedure with uncertainty quantification for a grand total of four estimates: four tuples of four meshes were available. The results are summarized in Table 5. The estimate in the last row is just 0.02% off of the reference value of $−0.30148′$⁠.

This shell is known to respond to the applied loading with a mixed membrane-bending stress state. Given that both the membrane and bending+shear mechanics of the element will be involved, it isn't clear what the rate of convergence should be. It is interesting to note that the evaluated convergence rate varies noticeably from the coarsest meshes to the finest (refer to Table 5). Clearly, until the estimation procedure begins to work with reasonably fine meshes (let us say mesh 4 and higher), there can be no doubt that the asymptotic range had not been reached. Yet, the estimates of the limit solution and its associated uncertainty are quite reasonable, and the 95% confidence interval includes the reference solution. This is also shown in Fig. 5, which graphically summarizes the estimates of the limit deflection (filled solid dots), and the uncertainty intervals of $±2×wmad$⁠.

In order to satisfy the assumptions of the GCI calculation, the solution was recalculated using a sequence of grids related by nesting, i.e., the meshes were generated by bisection in both axial and circumferential direction. Meshes 0 had 4 element edges along each direction. Estimation using the GCI is shown in Fig. 6. Clearly, the error band reflects the distance between the solution on the finest mesh for each triplet and the true solution. The width of the error band is much larger than our confidence interval (Fig. 5).

where ν is the Poisson ratio, q is the magnitude of the uniform distributed loading, L is the length of the side of the triangle. Here we consider a very thin plate (with thickness $0.05 mm$⁠), and the above parameters of $L=500 mm, ν=0.3, q=100 N·mm−2$⁠. The target value is the bending moment resultant at the centroid of $Mx,exact=451389 N·mm/mm$⁠.

The problem was solved using the abaqus finite element program [27]. The thin-shell-specialized three node triangular finite element with a discrete Kirchhoff-Love constraint, STRI3, was used with unstructured meshes. Mesh 0 was constructed of elements of size $60 mm$⁠, and subsequent (unstructured) meshes were produced by halving the element size. The moment resultants were evaluated at the node placed at the centroid (refer to Table 6).

The extrapolation estimate of the true solution from the first four grids resulted as $Mx=452624±1832 N·mm/mm$⁠, and from the grids 1, 2, 3, and 4, we obtained $Mx=451446±376 N·mm/mm$⁠. Both estimates contained the analytical solution (refer to Fig. 7). For the convergence rate we obtained $β=0.923±0.288$⁠, and $β=2.198±0.868$⁠, respectively. The convergence rate should be theoretically close to 2.0 in the asymptotic range.

The procedure outlined above is predicated on the data being suitable for extrapolation. Essentially that means that the calculated results approach a limit value monotonically, either from below or from above, and the successive approximate errors (i.e., differences of successive solutions) decrease in magnitude with each refinement. If these conditions are not met, the extrapolation formulas cannot be applied, in particular, Eq. (9) will not have a solution.

Recall that we have postulated the availability of four solutions obtained from four meshes. It is possible that the data in this dataset is not suitable for extrapolation for one or two of them for extrapolations. This would result in a reduced-size dataset, which could still be used to calculate the confidence interval (14). However, we should probably be wary of putting too much faith into an estimate based on so little data. An alternative would be to use the methods of the RVA method [10].

Since Ref. 17 does not provide sufficient information to reproduce their calculations, we furnish here examples of computations that arrive at error models in the form of Eq. (4) by optimization approaches. The equilateral plate of Sec. 2.3 will be used as an illustration. The quadratic triangle STRI65 has been shown to be too flexible in thin shell applications [24]. In fact, one might argue that the finite element cannot be relied upon. We illustrate the problems here with the values of the centroid bending moment resultant computed with the STRI65 triangular element starting from mesh 0 of element size $50 mm$⁠, and then halving the element size four times to produce the data in Table 7. Clearly, the element does not converge monotonically, even though the last four results do.

for the “two-norm” of the error of the model fit. The nonsmooth optimization problems could be solved with the evolutionary optimization solver widely available in Excel (refer to the Appendix).

Figure 8 shows the data points from the first four rows of Table 7 with markers and the three models with curves: the dashed line is the model for the objective function (23), the dotted line is the model for the objective function (24), and the dash-dotted line is the model for the objective function (25). Further optimization models could be derived by weighting the results on the finer meshes more strongly than the results from the coarse meshes [10].

Figure 9 shows the models for the same objective functions calculated for the entire range of five simulation results. Evidently, the predictions of the limit value are heavily dependent on the formulation of the objective function. Neither model could predict the analytical solution very well.

The finest four meshes produced results that appear to be compatible with extrapolation (refer to Fig. 10). The scatter of the predicted true values of the bending moment is substantial, and again the predictions depend quite strongly on the type of the objective function. However, in this case, our procedure was applicable, and the results were: $Mx=451229±1754 N·mm/mm, β=2.977±2.210$⁠. While the prediction of the bending moment resultant was spot on, the uncertainty was quite large. That uncertainty was also reflected in the convergence rate, which was quite high, with a very broad confidence interval.

A tentative conclusion may be that it is better to have data suitable for extrapolation and to use the proposed procedure, rather than to have a general optimization procedure for using data unsuitable for extrapolation.

The Richardson extrapolation procedure was extended to include quantification of the uncertainty of the estimate of the true solution. The estimate of the true solution is computed from a database of extrapolated values extracted from four three-tuples of calculated solutions based on a sequence of four progressively refined meshes. The database can be processed using median statistics to extract the best guess of the solution and the median absolute deviation as a measure of the uncertainty. The estimator produces the best guess of the solution and a confidence interval (or, alternatively, an error band). It was shown here that the confidence interval contained the true solution, at least in our examples, but no bounds can be guaranteed in general. The estimator does not require substantial software infrastructure, as it relies purely on the solution of a single nonlinear equation, easily treated with a handheld calculator. The precondition of the procedure is a monotonically converging sequence of computed solutions. We showed in the discussion section that it may be possible to accommodate data unsuitable for extrapolation with a procedure based on optimization. However, preference should be given to the acquisition of data suitable for extrapolation.

William J. Rider, Sandia National Laboratories, Albuquerque, is thanked for clarifications concerning the paper cited herein. The anonymous reviewers provided valuable comments, and their efforts are sincerely appreciated.

Figure 11 shows the setup of the spreadsheet for the optimization model (23). The objective function is defined as MAX(ABS($B$11-C2:C6-$C$11*POWER(B2:B6, $D$11))). Note that the evolutionary solver will not be able to converge without reasonably tight constraints on the “changing variables” in cells B11, C11, D11. While it is reasonably easy to guess the ranges of the true value and the convergence rate, guessing the range of the constant C11 can be tricky.