## Abstract

A simple procedure for estimating the uncertainty of estimates of true solutions to problems of deflection, stress concentrations, and force resultants in solid and structural mechanics is introduced. Richardson extrapolation is carried out on a dataset of samples from a sequence of four grids. Simple median-based statistical analysis is used to establish 95% confidence intervals. The procedure leads to simple calculations that deliver reasonably tight estimates of the true solution and confidence intervals.

## Introduction

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:

Of primary concern is the extensive use of under-determined (e.g., two-point) and critically-determined (e.g., three-point) data sets. In the under-determined case, there are infinitely many solutions limiting the usefulness of such estimates and providing ample opportunity for many uncontrolled pathologies (e.g., strong dependence on initial guess, or a strong directional bias); this could be better addressed by increasing the use of regularization methods, although a strategy for specifying the regularization penalties should also be addressed. In the critically-determined case, all of the nonregularized fits will be exact repetitions of each other, and this fact is not taken into consideration in the Robust Verification method.

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.

## 1 Methods and Materials

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.)

### 1.1 Extrapolation From Three Grids.

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

#### Element Size.

The definition of the element size is trivial for uniform grids. For graded meshes, the element size varies from point to point. We assume here that we take one particular graded mesh M0, and produce a series of progressively finer meshes from M0 by scaling the element size as a function of the position $x∈Ω$ by a constant, the refinement factor$α<1$, everywhere in the domain Ω, such that
(1)

For instance, mesh M1 would be produced with element size $h1(x)=αh0(x)$, mesh M2 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 M0 to M1 to M2 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 M1 could be produced with element size $h1(x)=α0h0(x)$, mesh M2 with element size $h2(x)=α1h1(x)=α0α1h0(x)$, and so on.

We may profitably define a relative element size
(2)
Note that this is a constant, independent of the location x. In the following, we shall consider the relative element size a continuous variable, and the extrapolation will be performed with respect to this variable. For any particular mesh j, we may write the relative element size in terms of the refinement factors between successive meshes as
(3)

#### Error Expansion.

Richardson extrapolation is a way of extracting an asymptotic estimate of some quantity of interest from a series of its computed values. If we assume that the true error in the quantity q may be expanded in a Taylor series at element size h =0 (or, alternatively relative element size γ = 0), we may write
(4)
where $qex$ is the unknown true value of the quantity, $qγ$ is the computed approximate value for a particular γ, $Eq(γ)$ is the true error, C is an unknown constant of the leading term $γβ$, with the rate of convergence, $β>0$, again, unknown. Provided C, and β do not depend on γ for small element sizes (this is presumed to hold in the asymptotic range), we might be able to estimate all of the three quantities $qex$, C, and β, if three numerical solutions, $qγi$, are obtained for three different relative element sizes, $γ1<γ2<γ3$. Then, we have a sequence of solutions $qγj$, and so we can write the so-called approximate error as the difference between two successive solutions
(5)
Replacing the approximate solution with the true solution and the true error, $qγj=qex−Eq(γj)$, we obtain
(6)
Substituting from (4), we get
(7)
This equation may be written twice, for two successive approximate errors
(8)
which may be readily solved for C to yield
(9)
This nonlinear equation may be solved for the convergence rate β. Then the constant C follows as:
(10)

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$.

### 1.2 Estimate of the Uncertainty.

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:

The succession of grids should be ideally based on nesting. However, this is typically not the case, and unstructured grids are generated by manipulating some form of mesh control to produce meshes which approximately satisfy (1). This will introduces more or less random error in the target quantity.

The solutions obtained may not (all) be in the asymptotic range. The consequence are errors introduced in the extrapolation procedure, which are not necessarily all random.

Computer-arithmetic errors may introduce errors in the true-solution estimation. Loss of significant digits, or inherent limits of processes iterating to convergence will introduce random errors into the computed numbers that are being piped into the extrapolation. Consequently, the extrapolations will be effected by these random errors.

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 L1 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).

#### Estimate of the True Solution.

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.

Table 1

Extrapolation from four grids

Inputs$qγ1$$qγ2$$qγ3$$qγ4$Extrapolation
$⇒ qex(1)$
$⇒ qex(2)$
$⇒ qex(3)$
$⇒ qex(4)$
Inputs$qγ1$$qγ2$$qγ3$$qγ4$Extrapolation
$⇒ qex(1)$
$⇒ qex(2)$
$⇒ qex(3)$
$⇒ qex(4)$
To estimate the uncertainty of the extrapolated value, we use the robust statistical approach based on the median absolute deviation from the median [15]. First, the median of the four predictions is computed
(11)
Second, median absolute deviation (MAD) from the median is evaluated from
(12)

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].

Finally, the prediction of the true quantity is formed together with its uncertainty interval. Recognizing that $qex(4)$ is computed from the three finest meshes (c.f. Table 1), the prediction from this particular sampling inspires higher confidence than the other three predictions. In addition to $qex(4)$, we have the robust statistic of the median (11). These two quantities, $qex(4)$ and $q̃ex$, can be mixed together in some proportion, for instance, as
(13)

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$.

The estimate $qex$ can be formed with an uncertainty interval from Eq. (12), for instance, as
(14)

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.

#### Estimate of the Convergence Rate.

The four extrapolations $qex(k)$, where $k=1,2,3,4$, required the calculation of four convergence rates, $β(k)$. Clearly, we can use the equivalent of Eq. (11) to estimate the median, $β̃$, compute the median absolute deviation, $β̃mad$ from the equivalent of Eq. (12), and estimate the best guess of the convergence rate $β*$ from the equivalent of Eq. (13). The rate of convergence can then be estimated with a confidence interval of 95% as
(15)

#### Expert Judgment.

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.

## 2 Numerical Examples

The GCI is at this point a defacto standard for reporting convergence, especially in computational fluid dynamics [6,12]. In our notation, it can be defined as
(16)
The expression consists of three terms: Fs is a safety factor, typically taken at the value of 1.25 in well-converging simulations. The second term, the fraction in terms of the computed results, is the normalized approximate error. The last term is a correction factor that converts between the approximate error and the true error. The results are then reported in terms of the GCI as
(17)

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 GCI is a dimensionless indicator for the mesh convergence error relative to the finest-zoned solution. Specifically, when multiplied by the finest-zoned solution, it provides the width of an error band, centered on that solution, within which the exact solution is very likely to be contained.

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.

A related approach applied to stress concentrations had been recently proposed in Ref. [13]. The relative error was defined as
(18)
which may be compared with Eq. (16). The confidence interval may then be formulated as
(19)

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

### 2.1 Stress Concentration in 3D.

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.

The direct stresses σx and σy follow from the radial and angular stress solutions and could be found in Ref. [22] as:
(20)

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$.

Fig. 1
Fig. 1
Close modal

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 2

Stress concentration example. Meshes used, and the computed stress at point A.

Mesh# elements# nodesRelative element size γStress σx (MPa)
020601.000 × 101.96424
11602975.000 × 10–12.35360
2128017852.500 × 10–12.62957
310,240121771.250 × 10–12.79966
481,920895056.250 × 10–22.89550
5655,3606853773.125 × 10–22.94656
Mesh# elements# nodesRelative element size γStress σx (MPa)
020601.000 × 101.96424
11602975.000 × 10–12.35360
2128017852.500 × 10–12.62957
310,240121771.250 × 10–12.79966
481,920895056.250 × 10–22.89550
5655,3606853773.125 × 10–22.94656

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.

Table 3

Stress concentration example. Three successive extrapolations of the stress at point A.

MeshesStress σx (MPa)$2×σxmad$β$2×βmad$
0, 1, 2, 33.098740.167980.661540.12634
1, 2, 3, 43.025470.040270.803540.07961
2, 3, 4, 53.006470.010780.893190.04918
MeshesStress σx (MPa)$2×σxmad$β$2×βmad$
0, 1, 2, 33.098740.167980.661540.12634
1, 2, 3, 43.025470.040270.803540.07961
2, 3, 4, 53.006470.010780.893190.04918

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.

Fig. 2
Fig. 2
Close modal

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.

Fig. 3
Fig. 3
Close modal

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.

Fig. 4
Fig. 4
Close modal

### 2.2 Shell Problem: Scordelis-Lo Roof.

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′$.

Table 4

Scordelis-Lo roof example. Computed deflection at the midedge point for all considered meshes.

Mesh# elements per sideDeflection (ft)
07–0.249830
19–0.266567
212–0.280147
315–0.287211
418–0.291278
523–0.295003
628–0.296972
Mesh# elements per sideDeflection (ft)
07–0.249830
19–0.266567
212–0.280147
315–0.287211
418–0.291278
523–0.295003
628–0.296972
Table 5

Scordelis-Lo roof example. Successive extrapolations of the deflection at the midedge point.

MeshesDeflection w (ft)$2×wmad$β$2×βmad$
0, 1, 2, 3–0.3051330.005511.5020.171
1, 2, 3, 4–0.3026530.001841.6830.090
2, 3, 4, 5–0.3017410.000531.7970.060
3, 4, 5, 6–0.3014270.000231.8630.040
MeshesDeflection w (ft)$2×wmad$β$2×βmad$
0, 1, 2, 3–0.3051330.005511.5020.171
1, 2, 3, 4–0.3026530.001841.6830.090
2, 3, 4, 5–0.3017410.000531.7970.060
3, 4, 5, 6–0.3014270.000231.8630.040

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$.

Fig. 5
Fig. 5
Close modal

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).

Fig. 6
Fig. 6
Close modal

### 2.3 Bending Moments in Equilateral Triangular Plate.

The exact solution of the shear-rigid Kirchhoff model of a simply supported and uniformly loaded equilateral triangular plate is known [26]. In particular, the bending moment resultants at the centroid of the plate (assuming homogeneous, isotropic material) are
(21)

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).

Table 6

Equilateral triangular plate. Bending moment at the centroid. abaqus thin-shell-specialized linear triangular element STRI3.

MeshRelative element size γjMoment $Mx,j$
01.000443750
10.500447393
20.250449892
30.125451133
40.0625451351
MeshRelative element size γjMoment $Mx,j$
01.000443750
10.500447393
20.250449892
30.125451133
40.0625451351

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.

Fig. 7
Fig. 7
Close modal

## 3 Discussion

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.

Table 7

Equilateral triangular plate. Bending moment at the centroid. abaqus thin-shell-specialized quadratic triangular element STRI65.

MeshRelative element size γjMoment $Mx,j$
01.000450601
10.500454784
20.250453153
30.125451558
40.0625451439
MeshRelative element size γjMoment $Mx,j$
01.000450601
10.500454784
20.250453153
30.125451558
40.0625451439
The extrapolation with confidence interval estimation as proposed in this paper is not possible from the first four values. However, we can formulate suitable optimization problems to arrive at some models of the form (4)
(22)
Here f is an objective function expressing the error of the fit of the model. For instance
(23)
for the “infinity norm,” or
(24)
for the “one-norm,” or
(25)

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].

Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal

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.

Fig. 10
Fig. 10
Close modal

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.

## 4 Conclusions

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.

## Acknowledgment

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.

### Appendix: Optimization Solver

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.

Fig. 11
Fig. 11
Close modal

## References

1.
Oberkampf
,
W. L.
, and
Roy
,
C. J.
,
2010
,
Verification and Validation in Scientific Computing
,
Cambridge University Press
, Cambridge, UK.
2.
Oberkampf
,
W. L.
,
DeLand
,
S. M.
,
Rutherford
,
B. M.
,
Diegert
,
K. V.
, and
Alvin
,
K. F.
,
2002
, “
Error and Uncertainty in Modeling and Simulation
,”
Reliab. Eng. Syst. Saf.
,
75
(
3
), pp.
333
357
.10.1016/S0951-8320(01)00120-X
3.
Richardson
,
L. F.
,
1911
, “
IX. The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, With an Application to the Stresses in a Masonry Dam
,”
Philos. Trans. R. Soc. London. Ser. A
,
210
(
459–470
), pp.
307
357
.10.1098/rsta.1911.0009
4.
Editorial,
2008
, “
Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications
,”
ASME J. Fluids Eng.
,
130
(
7
), p.
078001
.10.1115/1.2960953
5.
Roache
,
P. J.
,
1997
, “
Quantification of Uncertainty in Computational Fluid Dynamics
,”
Annu. Rev. Fluid Mech.
,
29
(
1
), pp.
123
160
.10.1146/annurev.fluid.29.1.123
6.
Roache
,
P. J.
,
2016
, “
Verification and Validation in Fluids Engineering: Some Current Issues
,”
ASME J. Fluids Eng.
,
138
(
10
), p.
101205
.10.1115/1.4033979
7.
Stern
,
F.
,
Wilson
,
R. V.
,
Coleman
,
H. W.
, and
Paterson
,
E. G.
,
2001
, “
Comprehensive Approach to Verification and Validation of CFD Simulations–Part 1: Methodology and Procedures
,”
ASME J. Fluids Eng.
,
123
(
4
), pp.
793
802
.10.1115/1.1412235
8.
Xing
,
T.
, and
Stern
,
F.
,
2010
, “
Factors of Safety for Richardson Extrapolation
,”
ASME J. Fluids Eng.
,
132
(
6
), p.
061403
.10.1115/1.4001771
9.
Eca
,
L.
, and
Hoekstra
,
M.
,
2014
, “
A Procedure for the Estimation of the Numerical Uncertainty of CFD Calculations Based on Grid Refinement Studies
,”
J. Comput. Phys.
,
262
, pp.
104
130
.10.1016/j.jcp.2014.01.006
10.
Rider
,
W.
,
Witkowski
,
W.
,
Kamm
,
J. R.
, and
Wildey
,
T.
,
2016
, “
Robust Verification Analysis
,”
J. Comput. Phys.
,
307
, pp.
146
163
.10.1016/j.jcp.2015.11.054
11.
,
G. A.
,
Martin
,
N.
,
Moore
,
C. H.
,
Huang
,
A.
, and
Cartwright
,
K. L.
,
2022
, “
Robust Verification of Stochastic Simulation Codes
,”
J. Comput. Phys.
,
451
, p.
110855
.10.1016/j.jcp.2021.110855
12.
ASME,
2012
, “
An Illustration of the Concepts of Verification and Validation in Computational Solid Mechanics
,” American Society of Mechanical Engineers, New York, Standard No. ASME V&V 10.1.
13.
Beisheim
,
J. R.
,
Sinclair
,
G. B.
, and
Roache
,
P. J.
,
2018
, “
Effective Convergence Checks for Verifying Finite Element Stresses at Three-Dimensional Stress Concentrations
,”
ASME J. Verif., Valid. Uncertainty Quantif.
,
3
(
3
), p.
034501
.10.1115/1.4042515
14.
Eça
,
L.
,
Hoekstra
,
M.
,
Beja Pedro
,
J. F.
, and
de Campos
,
J. A. C. F.
,
2013
, “
On the Characterization of Grid Density in Grid Refinement Studies for Discretization Error Estimation
,”
Int. J. Numer. Methods Fluids
,
72
(
1
), pp.
119
134
.10.1002/fld.3737
15.
Leys
,
C.
,
Ley
,
C.
,
Klein
,
O.
,
Bernard
,
P.
, and
Licata
,
L.
,
2013
, “
Detecting Outliers: Do Not Use Standard Deviation Around the Mean, Use Absolute Deviation Around the Median
,”
J. Exp. Soc. Psychol.
,
49
(
4
), pp.
764
766
.10.1016/j.jesp.2013.03.013
16.
Huber
,
P. J.
,
1981
,
Robust Statistics
,
Wiley
,
New York
.
17.
Rousseeuw
,
P. J.
, and
Croux
,
C.
,
1993
, “
Alternatives to the Median Absolute Deviation
,”
J. Am. Stat. Assoc.
,
88
(
424
), pp.
1273
1283
.10.1080/01621459.1993.10476408
18.
Banks
,
J. W.
,
Aslam
,
T.
, and
Rider
,
W. J.
,
2008
, “
On Sub-Linear Convergence for Linearly Degenerate Waves in Capturing Schemes
,”
J. Comput. Phys.
,
227
(
14
), pp.
6985
7002
.10.1016/j.jcp.2008.04.002
19.
Stein
,
E.
,
de Borst
,
R.
, and
Hughes
,
T. J. R.
,
2004
,
Encyclopedia of Computational Mechanics
, Vol.
3
in Encyclopedia of Computational Mechanics,
Wiley
, New York.
20.
Szabó
,
B.
, and
Babuška
,
I.
,
2021
,
Finite Element Analysis: Method, Verification and Validation
, Wiley Series in Computational Mechanics,
Wiley
, New York.
21.
Sivapuram
,
R.
, and
Krysl
,
P.
,
2018
, “
Improved Recovered Nodal Stress for Mean-Strain Finite Elements
,”
Finite Elem. Anal. Des.
,
146
, pp.
70
83
.10.1016/j.finel.2018.04.005
22.
,
M. H.
,
2014
, “
Chapter 8—Two-Dimensional Problem Solution
,”
Elasticity
, 3rd ed.,
Martin H.
, ed.,
,
Boston
, pp.
159
234
.
23.
Krysl
,
P.
,
2015
, “
Mean-Strain Eight-Node Hexahedron With Optimized Energy-Sampling Stabilization for Large-Strain Deformation
,”
Int. J. Numer. Methods Eng.
,
103
(
9
), pp.
650
670
.10.1002/nme.4907
24.
Krysl
,
P.
, and
Chen
,
J. S.
,
2022
, “
Benchmarking Computational Shell Models
,”
Arch. Comput. Methods Eng.
, pp.
1886
1784
.10.1007/s11831-022-09798-5
25.
Krysl
,
P.
,
2022
, “
Robust Flat-Facet Triangular Shell Finite Element
,”
Int. J. Numer. Methods Eng.
,
123
(
10
), pp.
2399
2423
.10.1002/nme.6944
26.
Timoshenko
,
S.
,
Timoshenko
,
S.
, and
Woinowsky-Krieger
,
S.
,
1959
,
Theory of Plates and Shells
(Engineering Mechanics Series),
McGraw-Hill
, New York.
27.
Anonymous
,
2018
,
ABAQUS/Standard User's Manual
, 6.14 ed.,
Dassault Systemes Simulia
, Providence, RI.