Goal oriented error estimation in multi-scale shell element finite element problems

A major challenge with modern aircraft design is the occurrence of structural features of varied length scales. Structural stiffness can be accurately represented using homogenisation, however aspects such as the onset of failure may require information on more refined length scale for both metallic and composite components. This work considers the errors encountered in the coarse global models due to the mesh size and how these are propagated into detailed local sub-models. The error is calculated by a goal oriented error estimator, formulated by solving dual problems and Zienkiewicz-Zhu smooth field recovery. Specifically, the novel concept of this work is applying the goal oriented error estimator to shell elements and propagating this error field into the continuum sub-model. This methodology is tested on a simplified aluminium beam section with four different local feature designs, thereby illustrating the sensitivity to various local features with a common global setting. The simulations show that when the feature models only contained holes on the flange section, there was little sensitivity of the von Mises stress to the design modifications. However, when holes were added to the webbing section, there were large stress concentrations that predicted yielding. Despite this increase in nominal stress, the maximum error does not significantly change. However, the error field does change near the holes. A Monte Carlo simulation utilising marginal distributions is performed to show the robustness of the multi-scale analysis to uncertainty in the global error estimation as would be expected in experimental measurements. This shows a trade-off between Saint-Venant’s principle of the applied loading and stress concentrations on the feature model when investigating the response variance.


Introduction
An understanding and evaluation of error magnitudes and bounds is vitally important in almost any engineering problem. It is not possible to assign a quantifiable level of confidence against model predictions without this information, resulting in conservatism and waste. The present work is concerned with practical error approximation in an industry relevant analysis problems. Before discussing the particulars of the present work, it is important to recognise the different types of error and analysis types [1].
Approximation error is concerned with evaluating the inaccuracies which are inherent to the discretization methods that are required in order to approximate the solutions to mathematical models. Modelling error, on the other hand, is concerned with how well an abstract model approximates physical phenomena in the real world. Estimates of approximation error may be termed a priori estimates or a posteriori estimates [2]. The former utilises the problem definition and the discretization to estimate error, whereas the latter uses the solutions approximation itself to estimate error. A distinction should also be made between the approximation of error bounds and the estimation of error itself [3]. The former can be guaranteed but may well be inaccurate (error bounds may be large, for example), whereas the latter can, in general, not be guaranteed but will be constrained. As highlighted by Grätch and Blathe error estimates should possess several properties [3]. Clearly, error estimates should be accurate (i.e. predicted error is comparable to the actual, true error). A good error estimate should also: asymptotically tend to zero as the mesh density increases (i.e. as the approximate solutions approaches reality), produce tight bounds for the error, be computationally simple and inexpensive, be robust enough to be applicable in a wide range of foreseeable (potentially non-linear) applications, and inform mesh refinements such that approximate solution processes can be optimised. Grätch and Blathe acknowledge that, at present, no one estimator can satisfy all of these requirements [3].
Considering the preliminary definitions, the focus of this work can be outlined. This work is concerned with the application of "goal orientated" a posteriori error estimates of approximation error. More specifically, goal orientated error estimates (GOEEs) are determined for coarse mesh approximations of a system using shell elements formulations. Errors in translational and rotational degrees of freedom (DoFs) at the shell mid plane are then propagated to continuum element realisations of the detailed sub regions. A challenge that the aerospace industry is currently facing is performing accurate simulations of modern designs (particularly for quantities such as failure detection) and quantifying the sensitivity of simulation predictions in order to avoid overly conservative assumptions afterwards, such as large factors of safety. One of the main reasons for this challenge is the structural components within the design being based on a variety of length scales. Design features such as fillets, laminate layers, small holes for instrumentation, and others would typically require the use of hyper-refined models that are infeasible for a realistic system. System level computations, like stiffness, can be accurately calculated via homogenisation. However, several desired computations, like the onset of failure, require explicit information of the design features on the varied length scale. Multi-scale analysis accounts for the multiple structural length scales, wherein DoFs values are passed from the coarse global model to a localised detailed representation. The present work considers a common response to this challenge, wherein large global structures/assemblies are modelled using shell elements (approximating stiffness and allowing for the estimation of translation and rotation DoFs) and finer detailed regions are modelled using continuum elements (these models are driven by DoF results for the coarse shell model and are implemented to estimate local stresses and strains). This type of scale/realisation bridging is commonplace in industry and practical stress analysis but has not been extensively explored in the literature, at least in relation to error estimation. One of the major issues with using a multi-scale modelling technique is the confidence in the use of multiple meshes and mesh sizes due to the difference in scale. However, in most cases, the large global model mesh refinement is limited due to the computational effort to perform the required simulations on these large models. This inflexibility of the global mesh leads to the confidence issue that this work addresses with both the classification of global meshing error and how this error propagates into the multi-scale simulation's failure analysis. The present work utilises GOEE in shell element global models and propagates uncertainty in driving DoFs to continuum element sub-models, such that uncertainty in stress/strain quantities (and any other dependent terms) can be evaluated using standard Monte-Carlo sampling approaches.
Approximations of errors in elliptical partial differential equations (PDEs) has received a good deal of attention in the literature [4][5][6]. This is encouraging for the stress analyst, as the equations that govern linear elasticity (both the displacement based Navier-Lamé formulation and the stress based Beltrami-Mitchell formulation) take elliptical forms. In many cases, error measures are utilised to drive mesh refinement through, for example, polytree decomposition algorithms [5,[7][8][9][10]. It is worth noting here the wide range of problems and discretization approaches that have been considered in the error estimation literature. In addition to "conventional" finite element analysis (FEA), error estimates have been derived for boundary element methods [8,9,11], immersed surface methods [12], multi-grid and composite FEA methods [7,13,14], extended finite elements (XFEM) [2], and stress singularity problems [15]. The work of Larson and Runesson deserves particular note here as it is concerned with error estimation in multiscale problems [7]. Importantly, "seamless" scale bridging was achieved through the development on a single error estimator on the macro scale that drives mesh refinement at all scales. Numerous reviews of a priori and a posteriori error estimators are available in the literature [2,3]. In most cases, error estimators are categorised as energy norm based (including element residual and subdomain residual methods) or recovery based estimators. The latter includes the well know Babuska and Rheinboldt estimator, the Kelly, Gage, Zienkiewicz and Babuska estimator, and the Zienkiewicz-Zhu patch recovery technique [2]. In the present work GOEEs are utilised, rather than providing a general indication of approximation error, as they allow for the quantification of error in a specific quantity of interest (QoI) [1,3]. GOEE can be traced back to the 1990's [16] through the work of Prudhomme, Oden, and Ainsworth, [1,2,17], Ladevéze [18], and Bathe [3]. In the interest of clarity, readers should note that QoIs in the present work are DoFs (translational and rotational) in shell element models. If meaningful error estimates of these can be derived, it is straightforward to sample from the resulting distributions and propagate uncertainties to continuum element sub-models.
Multi-scale modelling methods are a group of powerful techniques that can answer vitally important questions in many engineering sectors, particularly aerospace. Fundamentally, the multi-scale methods allow for many length and time scales to be incorporated in solving material/structural analysis problems, reducing the computational costs associated with a traditional refined simulation method [19,20]. An example of this is the consideration of microscopic behaviour in crystal plasticity to the macroscopic response of full sized engineering structures. Multi-scale methods can, in a sense, be viewed as a set of approaches which homogenise the heterogeneity observed in all materials at some length scale. These can equally be applied to metallic material as to composite materials, both of which are ubiquitous in aerospace structures. It is well known that the importance of retaining the local features depends on a variety of factors, such as the relationship between micro and macro structures, domain of interest, and the behaviour of interest. Multi-scale methods allow for the efficient introduction of additional information required to describe these effects and solve related problems. Multi-scale modelling techniques have been applied to a wide variety of fields, including aerospace [21][22][23], marine [24], and civil [25].
Many different multi-scale methods can be found in the literature and several characterisation schemes have been proposed. Weinan considers two different categorization types for multi-scale problems: Type A in which there are local defects, singularities etc. that require a locally micro-scale model in an otherwise coarser global model, and Type B problems which have micro-scale features throughout, and require fine-scale modelling everywhere, for instance within some form of computational homogenisation [20]. Other authors, such as Geers in [26], have classified these same categories but under different names. In [26], Type A methods are called Hierarchical methods while Type B methods are called Concurrent methods. There is inevitably some overlap between the techniques applicable to Type A and Type B problems, for example Kim and Swan used adaptive refinement of voxel meshes of representative volume elements within their numerical homogenisation approaches [27].
Type A problems are typically solved using mature, classical multi-scale modelling approaches such as sub-modelling [28], domain decomposition [29,30] and local mesh refinement including adaptive mesh refinement. Of particular note here are the procedures used within the CleanSky2 project MARQUESS, wherein a pre-computed database of solutions enables a single sub-model to be applied to multiple instances of a recurring feature on a global finite element model of a composite component [31]. Type A methods are the primary focus of the present work. A posteriori error estimates have been applied to type A problems in the work of Tirvaudey et al. [32], wherein a weighted residual based goal orientated estimate is applied in non-intrusive sub-modelling problems. Error contributions due to model, discretization, and convergence sources were evaluated.
An example of a Type B problem is the multi-scale modelling of complex fibre architectures, such as 3D woven textile composites. Many techniques can be identified for tackling this type of problem. Computational homogenisation, for example, attempts to determine equivalent properties via representative volume elements within a periodic displacement field. The Finite Element Squared (FE2) approach, on the other hand, utilises a fine mesh discretization linked to the Gauss points of a coarser mesh. The multi-scale finite element method (MsFEM) uses a fine mesh substructure to replace each element in the global model [20,33]. Examples from the recent literature include Liang et al.'s use of voxel models (generated using the well-known TexGen software) to analyse woven textile composites with each element being assigned to a particular material component [34]. Shi et al. use a three-scale model involving representative volume elements at the micro and meso-scales to model the fracture of braided composites [35]. Liu et al. use a variant of a voxel Finite Element (FE) mesh termed the inhomogeneous FE mesh to model a woven textile composite, with the material varying from integration point to integration point rather than with material boundaries being assumed to follow the mesh [36]. A posteriori error estimates for MsFEM problems have been developed in the work of Chung and Chamoin [37,38].
While there are many methods to perform this multi-scale modelling, a method widely used in industry utilizes pre-computed sub-models [39] and is particularly well suited to simulating local features, such as the ones of interest in this work. Specifically, this method utilizes a bottom-up and top-down approach to identify local failure locations [39][40][41]. The methodology is particularly useful, especially for this study, due to the superposition principle used in the sub-modelling procedure. This typically takes the displacement from the global model and converts it into the stress field of the sub-model that can be converted into a specified failure criterion.
Error in this multi-scale formulation is non-obvious with many researchers studying how to quantify and propagate the error. For this work, the study of error focuses on a Goal Oriented Error Estimator (GOEE) approach in order to represent the uncertainty in a physically relevant quantity that can be propagated into an error on a failure criterion [42][43][44][45]. This failure criterion is mainly used to identify areas of high interest in the refined sub-models called hot-spot identification. These areas are used to give a higher resolution of high importance sub-areas in the model (such as fillets, small holes, or multi-layered materials).
The work presented in this paper introduces a novel use of GOEE into this multiscale methodology using shell element formulations. While the use of shell elements is exceedingly common in the aerospace industry for large designs, this application of multi-scale propagation of a GOEE field is novel. The chosen GOEE approach focuses on the error associated with the mesh on the coarse global model with a focus on the spatial distribution of errors. Specifically, this approach propagates the error field into the refined sub-model and determines the effect of this error on the failure criterion (such as von Mises stress). The novelty of this work is in the use of an error field as opposed to single point measurements and the multi-scale propagation of the error. Previous work in this methodology presented in [41] focuses on the calculation of the GOEE spatial distribution. This work expands this calculation into a multi-scale analysis for isotropic, metallic materials. The methodology is expected to be valid for more complex materials, such as composite materials, commonly used in the aerospace industry.

GOEE methodology
The work presented in this paper uses a GOEE approach that utilizes a dual formulation, the Zienkiewicz-Zhu (ZZ) smoothing recovery for the strain, and a GOEE definition that incorporates these methods, expanding the work in [45]. This GOEE estimates the error due to the mesh in the global system. The other aspect of this paper focuses on multiscale propagation of this GOEE field into a refined local feature, discussed in further detail in "Application to multi-scale GOEE propagation" Section. To better illustrate this GOEE approach and multi-scale propagation, a general workflow is presented in Fig. 1. This illustrates how some methods are used multiple times, such as the ZZ recovery. As a note of nomenclature, the use of the term "primal solution" refers to a traditional FE evaluation utilizing the coarse global mesh with applied Boundary Conditions (BCs) and forces resulting in the nodal displacements. This section is split into three sections that describes the main formulations used in this work, which are further referenced in In "Global GOEE approach" Section, the definition of the dual formulation, how it is calculated, the ZZ recovery, and the GOEE formulations are expressed. "Custom element formulation" Section explains the custom shell elements used in this work to mimic the ABAQUS S4 shell element using python, but with the ability to output the required quantities used in the dual formulation. Finally, "Specifics for using the GOEE" Section describes how this specific implementation makes adjustments on the original definition of this methodology due the multi-scale nature and the material used in the test system.

Global GOEE approach
The first method used in the GOEE approach utilized in this work is the introduction of a dual formulation to characterize a Quantity of Interest (QoI) such as stress at a location. Utilizing this formulation requires an additional FE evaluation to calculate the QoI that is used to estimate the error [46,47]. The first evaluation is the standard primal formulation. This is a traditional FE evaluation with applied forces resulting in nodal displacements. The system considered in this work is static, so only the stiffness matrix is required for these calculations, but this methodology is not limited to static simulations. Additional information (in addition to the stiffness matrix and displacement) from the primal solution is required to formulate the GOEE estimate and is discussed later in this section. For any given system, only one primal solution is required. In the current implementation, the stiffness matrix from the primal solution is also used in the dual equations of motions to eliminate the need to regenerate the system multiple times thus reducing the computational time.
The dual formulation, used in quantifying the error, is heavily dependent on the selection of the QoI. This formulation requires the generalized displacement due to a generalized forcing vector. For a static system, this can be expressed as: where [K ] is the stiffness matrix, {Z i } is the generalized displacement corresponding to the DoF's contribution to the ith QoI, and {Q i } is the generalized forcing vector that depends on the QoI and location. The stiffness matrix is the same as the primal stiffness in the cases presented in this paper, while the generalized forcing vector is dependent on the type of QoI (average stress, max displacement, etc.) and the location of interest.
To calculate {Q i }, the QoI must be written in the FE setting as: where {u} is the displacement field from the primal solution. In this work, the average displacement in a region is considered due to the utilization in the multi-scale propagation workflow. As an example, the average displacement in the y-direction, Eq. 2 is expressed as:ū where the value of the Q vector is calculated as: with N y being a matrix of the y-component of the C 0 elemental shape functions, î is a pointer matrix to identify the DoFs associated with the location within the integration (this is 24 by the number of DoFs where each row is comprised of zeros with a single unity entry), and 0i is the total domain being considered. This integration is typically decomposed by the element boundaries since both N y and î are dependent on the which element is being evaluated. This formulation is useful if a specific region is desired. For this work, average QoI measurements are of interest due to the multi-scale propagation. Taking the displacement at a specific point would require information at an exact location, which can introduce errors associated with interpolation. To eliminate this issue, a sharp Gaussian distribution is applied centered at the location of interest. This results in, with the addition of Gauss quadrature integration, the calculation of the Q vector in Equation 4 is approximated via: where k is the Gauss quadrature location with a total of N int in the domain, |J k | is the determinate of the Jacobian matrix, W k is the Gauss quadrature weight for the integration point, andŴ k is a spatial distribution weight (via the narrow Gaussian). This formulation is used since the information from the FE (Gaussian weights, Jacobian, etc.) is known at the Gaussian integration locations. With this formulation, there is one normalization constraint of: due to the applied spatial weighting function.
The spatial weights (Ŵ k ) in the analysis is based on the Euclidean distance from the target location. In this paper, these weights are based on a Gaussian distribution via: where l is a user-defined length, x i is the location of the centre of the distribution, x k is the location of the Gaussian integration point, and a is a normalization factor to ensure Eq. 6 is valid. Once the Q vector is determined for a specific QoI, then the FE analysis can be recomputed in order to determine the dual solution {Z i }.
The main calculation of this GOEE approach is based on the difference between the discontinuous and smoothed strain fields. To create a smoothed field, ZZ recovery is used. ZZ recovery creates a piece-wise continuous field [48,49] of the directional strains. For the shell element considered, the shape functions are C 0 continuous, such that the derivative (strain) is not continuous across element boundaries. In order to make the strain continuous, the ZZ recovery method is used to create a smooth strain field [50,51].
To perform the ZZ recovery, the strain must be known at specific locations. For the use in FE, these locations are specified as the Gauss integration locations (similar to the calculation of the stiffness matrix). Once these values are known, then the ZZ recovery can be performed for the location of each node in each component of stress individually. For each node, the ZZ recovery averages the known strain values at the nearby integration points. This work utilizes the Nearest Neighbour's (NN) approach, which requires the information about the nodal connectivity (which nodes are in each element). The NN approach then takes the nearest integration point location for each element that node is connected to and averages them based on distance. Note that this information is easily available in the mesh data that describes each element as a collection of nodes.
Once the nodal values of the recovered strain field are determined, a surface is applied to the system. For ease-of-use, the C 0 element shape functions are used as the basis function for this surface. The equation for this surface is defined as: where { * } e is the recovered strains at the nodal values for element e, the superscript * represents a recovered value, and N e (x) are the element shape functions.
Once the ZZ recovery is performed on both the primal and dual strain fields, the GOEE is calculated. This formulation of the GOEE closely follows the work in [52] with slight differences on how the dual problem results are incorporated. The calculation is based on a modification of the energy norm. For the GOEE approach used in this work, the dual problems are implemented into this energy norm to be calculated as: with * u and * z i being the ZZ recovered strain for the primal and dual problem respectively, u h is the displacement determined by the primal FE problem, z hi is the generalized displacement from the dual FE problem (in units of the force normalized QoI), is the domain of interest, A : C is the double dot product of the two tensor quantities, C is the material constitutive tensor, and ∇ is a vector operator of derivatives. In simplified terms, the difference ( * − ∇u h ) is a measure of the difference between the discontinuous strain field and the smoothed recovered strain field.

Custom element formulation
Commercial finite element codes, such as ABAQUS [53], are widely used for the structural analysis of components and structures within industrial applications. However, for the purposes of this work, additional information is required to calculate the error estimator that is not returned by ABAQUS. This information is used within the element formulation (to create the stiffness matrix) but is not stored or returnable to the user. To overcome this, a custom, general purpose, iso-parametric, flat shell element has been developed for the simulation of 3D components within Python.
The custom shell element, named Q4, is used to store all the information used in the calculation. In the GOEE evaluations, values such as the strain-displacement matrices, integration point locations, and other quantities are required. For a typical ABAQUS analysis, this information is either not stored in memory or is not easily available. Ideally, this analysis could be performed using ABAQUS's S4 element if all the information was stored, but due to the black-box nature of the ABAQUS element formulation, this is not possible.
The proposed, degenerate continuum, element formulation, is implemented using the Mixed Interpolation of Tensorial Components (MITC) approach, first attributed, for a quadrilateral, four node elements (MITC4), to [54,55]. The element formulation uses a bilinear interpolation with 4 nodes, with 2 × 2 Gauss points for the bending and membrane integration with a single point contribution for shear, derived from the Reissner-Mindlin theory assuming linear material/geometric properties and small strains. The MITC4 formulation, developed to reduced shearing locking, has five DoF per mid-surface node, as the rotational stiffness about the z-axis, often termed the drilling DoF, is neglected due to the thin nature of the elements. To allow for three dimensional system and subsequent non-planar global coordinates, the element is implemented with a sixth DoF, but is constrained as a boundary condition. For the purpose of this paper, as the mid-surface is assumed as flat (i.e. there is no curvature to the elements and all nodes are co-planar), the shell element can essentially be classified by the superposition of plane stress, plate bending and shear stress, where the effects are assessed independently [56]. The internal energy for this formulation is defined for each element by a linear geometric interpolation scheme throughout the element, expressed as: where σ α and α are defined for the corresponding bending, membrane and shear components {α} for each element domain e , and σ α · α is the tensor dot product of the stress and strain. The linear-elastic stress-strain relations are defined for a homogeneous isotropic material as: where α is the applied strain and the material matrix C α is defined by the constitutive equation for plane stress/plate bending as: where E and ν are the material Young's modulus and Poisson's ration, t is the shell thickness, which is constant over the shell, and G is the shear modulus given as: where κ is an additional classical shear correction factor. The coefficient is applied to take into account the thickness variation at the surface rather than the theoretically defined constant distribution for the transverse shear stress across the thickness. In accordance with [57], k is usually taken equal to 5/6 for a homogeneous, isotropic, rectangular section with no curvature. The generalised strain-displacements for bending, membrane and shear are independently interpolated in the local coordinates by: where N n (ξ , η) are the shape functions of a standard bi-linear four node element and u n is the nodal deflection for node n. With strains computed from displacements from the localised strain matrices for the element via: To account for elements within three-dimensional coordinate space, a transformation of nodal displacements and forces from local to global Cartesian coordinate system is performed separately. For each strain-displacement matrix, this rotation is performed by: where L e α is the transformation matrix, which for flat elements, is constant for all element nodes and is defined using the following expression: where: λ x x is the dot product of the axes x and x etc. [56]. In Eq. 18, the matrix is a despite the The element stiffness matrix is therefore obtained by numerical integration for each element by: The vector of nodal forces, which are equivalent to distributed forces P are then calculated as: wheret is the surface traction and e are the Dirichlet BCs. This selective integration for both the stiffness matrix and force vector, using a classical shell theory, is considered a simple procedure for avoiding shear locking of the element.

Element formulation Verification
To evaluate the performance of the custom element implementation, two benchmark studies have been investigated, the Cook's trapezoidal skew beam and the hemispherical shell with an 18 o hole. The purpose of the verification is to assess the accuracy of the shell element (Q4) by comparing it to a referenced benchmark solution and the widely used S4 element within ABAQUS. The Cook-trapezoidal beam, proposed in [58], is used to assess the in-plane membrane performance when loaded in shear, under moderate distortion. The standard skew beam test, represented in Fig. 2a, is defined as a tapered beam, clamped on the left edge and subjected to a uni-axial traction load. The structure has a thickness of h = 1.0, and plane stress, material properties: Young's modulus E = 3×10 7 Pa and Poisson's ratio ν = 1/3. Where the loading, given as P = 1.0, is specified as a uniformly distributed shear load across the right end edge of the beam. taken from [59] and [60] respectively, using a refined numerical model. With the stresses at point A and B calculated as: where σ and τ are the principle stresses. The displacement and stress results show that although the Q4 element is not as accurate for coarse meshes, both elements are sensitive to distortion in membrane deformation problems. With Q4 requiring a finer mesh to be able to converge to the reference solution. Verification results by [61], using the same mesh densities, show that the measured vertical tip displacement for the Q4 element is equivalent to the MITC4 element results they present. The verification of both displacement and stress shows that the Q4 element is valid for natural values (displacement) and its derivative (stress/strain), which are both used in the GOEE approach in this work. The hemispherical shell with 18 o hole problem, proposed by [62], and represented in Fig. 2b, is investigated to evaluate the elements performance under in-extensional bending deformations and rigid body rotations, normal to the shell surface. The, doublecurved, hemisphere, with a radius R = 10.0, is defined as a thin shell, thickness h = 0.04, with an 18 o open hole at top and in plane stress material properties: Young's modulus E = 6.825×10 7 Pa and Poisson's ratio ν = 0.3. Where the loading is defined as two pairs of opposite radial concentrated loads P = 1.0. Utilizing axial symmetry, one quarter of the structure, corresponding to the region ABCD, is modelled with symmetrical boundary conditions along edges AC and BC. As with the previous benchmark, the same element mesh densities (N ×N ) are evaluated, with uniform, structured mesh patterns. A radial displacement coincident value at point A, U REF (A) = 0.094, is used as the reference solution. [63]. The normalized convergence displacement (U A /U REF (A) ) is illustrated in Fig. 3c, comparing Q4 and S4 against the referenced solution. One interesting aspect to note is the increased accuracy, even with a course mesh, compared to the Cook's skew beam. These two validation cases demonstrates that the Q4 element used in this work is comparable to other shell elements commonly used. Currently the proposed formulation is suitable for the analysis of thin shells of arbitrary shapes and has been verified against well known membrane and bending benchmark problems. With the latter agreeing well with the referenced solution for lower mesh densities. For membrane dominated problems, although its overall accuracy, under in plane shear deformation, is reduced, it can converge to an accurate solution when finer meshes are used. Minimising mesh distortion. To increase overall accuracy, a shell element formulation, with a second order geometric interpolation scheme, would need to be adopted. Where the formulation for bending and membrane contributions interact and the cannot be treated independently using the current super position approach [64]. This is a subject of ongoing research. However, it should be noted that the example problems implemented in the current work avoids this complexity by implementing only thin, flat elements, with a uniform mesh.

Specifics for using the GOEE
The GOEE approach described in "Global GOEE approach" Section uses a single dual problem with a single material constitutive matrix. In this work, these two details are modified to match the desired analysis. The first detail discussed is the single material constitutive matrix. This approach works perfectly when the element is constructed with a single material constitutive matrix. However, this work (as described in "Custom element formulation" Section) separates the construction of the stiffness matrix into three components of bending, membrane, and shear, thus creating three material constitutive matrices.
Due to the calculation of the stiffness matrix being comprised of bending, membrane, and shear components, the GOEE is also comprised of three components. The GOEE is calculated for each component then summed as: where GOEE α is the GOEE calculated with the bending, membrane, and shear material constitutive matrix. One aspect that can be noticed in "Custom element formulation" Section is that the material constitutive matrices do not span all six DoF of the system. To account for this, Eq. 9 does not use all six DoF for each component, only the DoF used for each component is utilized for the GOEE calculation. For example, the shear component only uses the in-plane rotational DoF. This calculation is done on a local coordinate system then transformed back into the global coordinate system. The second main modification is the use of multiple dual problems. In the initial approach, there is only one dual problem being performed. However, in this work, multiple dual problems are used since it takes the error at multiple points (boundary DoFs for the sub-model). To account for this, multiple GOEE are also computed. For each dual problem, there is a corresponding GOEE value. These values are stored and used in a method that is described in "Application to multi-scale GOEE propagation" Section.

Application to multi-scale GOEE propagation
The GOEE approach described in this paper is used to describe the meshing error in the global model in a single QoI, such as the displacement in a specific direction at a specified location (then iterated for each location and direction). To incorporate the multi-scale aspect of this work, a pre-computed sub-modelling approach (also called feature modelling) is used to take the information from the global model and estimate the quantities (specifically the element centroid stresses) within the feature model. This methodology is based on the generation of pre-computed unit normalized solutions for the sub-model. The pre-computed solution is used in a linear superposition calculation to approximate the loading on the sub-model.
The multi-scale methodology used in this work uses a two stage multi-scale approach. This involves using two models of the local feature to generate a matrix of results for linear superposition calculation. For clarity, the three models used in this multi-scale analysis are: 1. The global model of the entire system. This is typically constructed using only shell elements. These elements are assumed to be quadrilateral, but current work is being performed to expand this methodology to triangular and a mixture of triangular and quadrilateral elements. In the current implementation, this model is solved within python. 2. The shell surrogate model of the local feature. This also utilizes shell elements but includes physical features such as holes and a more refined mesh. For this work, these are S4R ABAQUS shell elements.
3. The feature model of the local feature. This uses the same region as the shell surrogate model but is expanded into full 3D elements. The mesh for the continuum elements does not necessarily correspond to the mesh of the shell surrogate due to the use of ABAQUS's shell-to-solid sub-modelling function. This model uses C3D8R ABAQUS continuum elements.
In addition to these models, two sets of nodes are defined within these models. Driving nodes are the nodes of the global model that define the boundary of the feature. The other set is the driven nodes comprising of the node set for the shell surrogate model on the boundary. These two sets of nodes overlap in physical space, thus there is a 1-to-1 mapping of driving and driven nodal location. In practice, this is not a requirement if an interpolation function is used. Using an interpolation function can cause discontinuities and discrepancies with the energy transfer between global and feature models, so this work enforces the boundary of the shell surrogate model to have the same nodal location as the global model while the interior contains a more refined mesh. This will be shown for the demonstration system in "Demonstration system" Section.
In order to calculate the pre-computed solution, unit deflection is applied to each driven DoF then propagated using the shell-to-solid procedure in ABAQUS. Simply, within the shell surrogate model, a single driven DoF are individually displaced by unity and the other driven DoF are kept to zero displacement; this is looped for each driven DoF. This can be quite expensive, but it is assumed that the shell surrogate models are relatively small compared to the global model. One reason this method was chosen is due to the industrial use cases that have repeated features (such as fillet/bonds) that would use the same pre-computed results but applied at different locations. This will have different macro-level results applied to the same pre-computed solution. After this calculation, the feature model directional stresses are calculated using the shell surrogate model as a basemodel and the feature model as the sub-model within ABAQUS utilizing the shell-to-solid sub-modelling option. After the shell-to-solid analysis within ABAQUS, a report of the directional stresses at the centroid of each element is written and converted into a large In order to calculate the feature model stress from the pre-computed solution, the [M] f matrix is used for each element via: where σ xy is the directional stress of the element centroid, M xj is the directional stress field of the feature model due to a unit displacement of driven DoF j in the shell surrogate model with a total of N DR , and U j is the displacement of the driving DoF in the global model that corresponds to the driven DoF j. One aspect to note, is that for this work, each driven/driving node has 6 DoF with 3 displacement and 3 rotational DoF. Each DoF for each driven node is utilized in the matrix [M] with the knowledge that one of the rotational DoF is zero due to the element formulation. Currently, the driving DoF are taken and applied directly to the feature model. Future work will incorporate an interpolation methodology to lessen the assumption that the global and shell surrogate models have the same boundary mesh. Equation 24 is useful for the calculation of the directional stresses. However, this work is focused on the error propagation. Through a simple variation analysis, the error on the feature model stress can be determined via the error on the displacement calculated via the GOEE. This analysis generates the equation with δσ xy representing the expected error on the stress, [M] f is the element pre-computed solution, and GOEE j being the calculated GOEE for the driving DoF j. To calculate the error on a failure criterion, the error of the stress components are propagated into the feature model. Then the failure criterion is performed for both the nominal values and the error adjusted values, and the difference is the reported error on the failure criterion. The main work presented in this paper assumes that the GOEE is deterministic (representing an offset) since a dual problem is performed for each driving DoF. A second analysis is also performed to quantify some of the robustness of this multi-scale methodology. This assigns a distribution to the driving DoF and performs a Monte-Carlo simulation to give an estimate on the robustness of the various feature models. The robustness measures are the 95% confidence interval of the von Mises stress at the location of maximum stress and the standard deviation field. In this paper, each driving DoF distribution is treated as independent. With the variances seen in these results, this assumption does not lead to large differences between adjacent nodes. Additionally, the robustness measure presented in this work is focused mainly to show how variability propagates into the feature model. If this approach is used for decision making, additional work might be required to determine a correlated multi-dimensional distribution for the Monte-Carlo sampling to ensure a smooth displacement field and prevent the possibility of local failures such a crack initiation.
It should be noted here that, with the a posteriori error estimates of driving DoFs in hand, a user has many options of how to take samples such that distributions in resultant quantities of interest (von Mises stresses in this work) can be developed. It is relatively straightforward, although at a larger computational cost, to develop a covariance matrix between error estimate observation points. This would allow conditional distributions in DoFs to be estimated, rather than the simple marginal approximations used here. The development of multivariate distributions will be the focus of a future publication, however the emphasis of the present work is the implementation of GOEE methods in shell element problems with a view to propagating uncertainties. The ease of implementation and the simplicity of the marginal sampling method are the main motivations for its application in the present work, however readers are encouraged to note that more robust alternatives can be readily implemented.
To perform this second simulation, the nominal displacement at the driving nodes are displaced by a normal distribution with a mean of the GOEE value calculated via the dual problems. The calculation of this distribution utilizes the limits of the error described in [65] where the limits of the error is based on a scaling factor of the true error. For this analysis, it is assumed that the true error is on the order of the average GOEE of the boundary nodes. In this simulation, the standard deviation of the driving nodes is based on a random number between zero and one and the mean GOEE for that direction among all the driving nodes. This random number is assigned per node from a standard uniform distribution. For node i in direction k, the displacement DoF is distributed as: where U j0 and GOEE i are is the nominal displacement and GOEE respectively generated from the previous analysis, Uni n is the random number from a standard uniform distribution for node n, and Mean j is the average GOEE for all the driving nodes in the direction of DoF j. It is noted that GOEE i corresponds to driving DoF j. In the Monte Carlo simulation, this displacement distribution is sampled and applied through Equation 25 to calculate the error of the directional stresses, then converted into the error of the von Mises stress (by subtracting the adjusted von Mises stress to the nominal von Mises). One aspect to note is that since the nodal value is based on a random number, the exact values shown in this analysis are illustrative to show the trends of how the randomness of the driving nodes in the global model affect the variability of the feature model.

Demonstration system
This work is demonstrated on a beam section to explore this methodology with a variety of local feature models, which can be seen in Fig. 4. The global system is comprised of aircraft relevant Aluminium-2024 and has a constant thickness of 1 mm with the red dots identifying the driving nodes of the system and where the feature is located. Although the multi-scale methodology is designed for composite material, this system uses Aluminium. This should be considered as a preliminary step to ensure that this methodology is reasonable and provides satisfactory results. Composite materials are currently being studied to be presented in future work. For this system, there are 61 driving nodes with a total of 366 DoF. This results in 366 dual simulations being performed. Additionally in Fig. 4, the BC of fixing the displacements of one end and the applied point-force of 400 N, located at the other end, is shown in orange. One of the main focuses of this work is to compare possible local features as a study of the sensitivity of this methodology to the feature model design. In total, there are four local features used, for which the shell surrogate models of each design are pictured in Fig. 5 with fully compatible boundary meshes to the global model. One aspect to note in Fig. 5 is that the features contain holes that are not modelled in the global model. This is due to the size of these holes being smaller than one element in the global model. However, since stress concentration are expected near holes, this multi-scale modelling approach is utilized in order to predict the stress near these holes with the propagated GOEE used to quantify the confidence of the stress given the course global mesh. The various feature models are based on the initial design (Fig. 5a) that contains one hole with a diameter of 10 mm and one with a diameter of 5 mm, both in the section called the flange. This work also makes modification of these holes as well as the addition of holes to the other section called the webbing.
Within the various feature designs, one design modifies the holes on the initial design, and two designs add holes to the webbing. Figure 5b expand the small hole into a 30 mm long slot with the same diameter (5 mm). In addition to modifying the already existing holes, two modifications are added to the webbing while maintaining the nominal holes in the flange. The first design contains a 10 mm hole in Fig. 5c, and the second replaces the single large hole with three individual 5 mm holes in Fig. 5d.
Another aspect to note is the mesh on the models in Fig. 5. This mesh enforces the boundary to be compliant with the global model, such that the driven and driving nodes correspond to the same physical locations. However, due to the utilization of the shellto-solid use in the sub-modelling technique, there is no requirement for the meshes of the shell surrogate and feature model to be identical (with the only requirement being related to the number of elements in the through thickness). So, for plotting the results for the local feature models, simplified, yet more refined, feature models are used. This refined mesh creates a large [M] matrix but is still feasible on a fairly standard desktop. As a test, the several feature models meshes were created while keeping the shell surrogate mesh constant. Results from these tests (not presented in this work) showed no noticeable differences since the main calculation occurs in the shell surrogate model.

Results
The focus of this work is based on the multi-scale propagation of meshing error in the global model. To give a reference to the errors experienced in this system, the global model GOEE field is presented in Fig. 6. These results are the error in the displacement in the bending direction (in meters). One thing to note is that similar fields are also calculated for each of the six directions. Additionally, the propagation analysis only requires the GOEE at the driving nodes, but the GOEE is computed for the entire global model to show the overall trend of the error. Overall, the errors are small compared to the nominal deflection, so the global model mesh is well defined for this loading condition. This gives an estimate of 0.062 mm for a 8 mm deflection.
Since the GOEE used in this work is based on the size of the global mesh, it is expected that a more refined mesh produces less error. While this has not been validated for the sub-modelling propagation, it has been validated on the global model. The results in [41] demonstrate the GOEE methodology in a single-scale analysis with various meshes. It was seen that the error decreases as the mesh become more refined.
For comparing the propagated GOEE, there are three simulations for each local feature. The first simulation does not use the GOEE but is used as a baseline. This is the von Mises stress of the feature. As a baseline, this gives a good indicator on possible failure and a simple comparison between designs. The second simulation is the error on the von Mises stress for each feature. This is calculated via a perturbation of the propagated GOEE for the components of stress, then the computed von Mises stress is subtracted from the nominal stress profile calculated from the first simulation. The final simulation is the robustness testing by adding a stochastic displacement of the driving DoF, propagated via a Monte Carlo sampling. This Monte Carlo sampling shows two results with the first result being the standard deviation field on the webbing, and the second being the 95% confidence interval at the location of maximum von Mises stress found in the first simulation. The standard deviation field gives a qualitative understanding on how the error spatially affects the feature model, while the confidence interval gives a quantitative estimate on the probability of failure.

Initial design
One of the key aspects of the initial design is the two different size holes in the flange. Adding these holes that are not part of the global model introduce stress concentrations around the holes. To demonstrate whether there are any effects due to the stress concentrations, the results for the first two simulation for the initial design are shown in Fig. 7 for both the nominal von Mises stress (in MPa) in Fig. 7a and the error of von Mises stress (also in MPa) in Fig. 7b.
The first aspect to investigate is the nominal stress in the feature model, seen in Fig. 7a. This shows a maximum stress of 217.9 MPa at the top of the webbing. It is interesting to denote the shape of the stress field. There is a gradient along the webbing with almost zero stress in the flange, with a small increase near the holes. This low stress on the flange makes sense due to the loading condition and the use of shell elements. In a shell with a small thickness, such as the ABAQUS S4R, there is very little through thickness stiffness, so the major component of stress for these elements on the flange is due to stretching as opposed to bending. In the custom user element, it is assumed to be a thin shell so there is zero through thickness stiffness. For this loading condition, the global bending is much larger than the in-plane stretching so the stress on the webbing, which experiences this bending, is much larger than in the flange. As a verification, although not shown in this work, a refined global model that contains the initial design feature was analysed. Since this is a simple and reasonably small system, this is still feasible as a validation study; This would not be feasible for a realistic system. The results from this simulation show nearly identical von Mises stress found in the feature region as well as the same field trend. There are some slight differences (less than 1% of the maximum stress) due to the stiffness of a continuum element compared to a shell element, but the differences in the stress are very small. This validation gives confidence in the multi-scale methodology used in this work. While the nominal stress gives a baseline understanding of the system, the main purpose of this work is for the GOEE propagated into the feature. The error of the von Mises stress can be seen in Fig. 7b. This is calculated by applying the GOEE to the directional components of the stress, taking the von Mises stress of the adjusted stress then taking the difference from the nominal von Mises stress. A major result from this analysis is that the maximum error on the stress is 3.5 MPa and is located at the bottom of the webbing at one edge with other zones of larger error also along the edge of the webbing. The location of the errors on the edge is plausible since the error is based on the global model and propagated through the edges. When looking solely at the maximum stress, the error is less than 2%, but the largest error is not located at the region of high stress.
The final simulation is the study of robustness in the multi-scale approach via the Monte Carlo sampling resulting in the standard deviation field and the 95% confidence interval of the von Mises stress at the location of maximum stress. This location is at the top of the webbing section. The confidence interval for this location is 220.0 ± 8.8 MPa. This produces a greater than 95% confidence that the maximum von Mises stress is below the yield criterion. In addition to this specific interval, the field for the standard deviation of the Monte Carlo samples for the webbing is shown in Fig. 8. One of the most interesting aspect of this standard deviation is that it is localized to the edges with a maximum

Slotted hole design
The second design expands the small hole in the initial design into a 30 mm long slot. This extension of the circular hole into a slot is thought of as a design that can accommodate more sensors on the design for testing and prototyping, for example adding a thermocouple near a pitot tube. The nominal stress and propagated error can be seen in Fig. 9. Both fields are very similar to the initial design feature with even the same maximum values.
In addition to the same maximum values, the slotted feature also closely matches the general field. The main reason for this is the fact that the flange experiences low stress due to the loading and the use of thin elements. This same effect can be seen for the propagated GOEE as the error in the von Mises stress. Both the nominal stress and the error fields look identical to the fields for the initial design.
Due to the low stresses on the flange, there is little difference between this design and the initial design. This is also experienced in the Monte Carlo simulation. With the Monte Carlo simulation, the results were nearly identical to the initial design results (within the error of the Monte Carlo simulation). This is primarily since no modifications were made to the areas where there is high stress. Because there are no identifiable differences, the slotted design results for the Monte Carlo simulation are not presented.

Single webbing hole design
While the first two designs contain holes only in the flange section, the last two contain holes in both the flange and webbing section to get a better understanding regarding these design changes in areas of higher stress. This specific design puts a single 10 mm hole in the webbing section. The reason for the location of this hole is based on the knowledge gained from the previous two designs. They showed that the holes on the flange do not make a large difference in any of the metrics used in this work for this specific loading case. To better test the methodology, the next two designs utilize holes near the areas  Fig. 10.
The first thing to observe in Fig. 10 is the addition of the webbing profile to better identify any stress concentration due to the additional hole. This is primarily used to highlight the location of the concentrations and as a comparison between the single and triple hole features. For the von Mises stress in Fig. 10a, the maximum stress is 289.7 MPa. The maximum stress for this single hole slightly exceeds the yield criterion for this type of aluminium of 265 MPa. With this stress value, it is believed that local plastic deformation would occur, but not cause a full failure. It is important to note that this simulation only accounts for linear-elastic material properties. This is not a requirement for the multi-scale methodology, just for the current implementation. In the current analysis, only linear-elastic behaviour is modelled, thus these results when applied to a physical system is only valid before yielding. Further nonlinear behaviour (such as plasticity) can be estimated in this analysis by a series of linear moduli, however this would complicate the generation and utilization of the [M] matrix. This would require multiple [M] matrices and a custom interpolation algorithm between these multiple matrices. Doing this is theoretically possible but is beyond the scope of this work.
In addition to the increase in the maximum value, the stress field is also significantly different compared to the initial design. The difference is apparent in two main areas. First is the location of the maximum stress. In the initial design, the maximum occurred on the upper edge, while the single hole has it located on the hole via a stress concentration, occurring at the top of the hole. The second is an almost swirling effect since the side of the hole has low stress that diffuses into the stress field. This is believed to be due partially to the size of the feature. It is believed that if the edge was expanded, there would be a more localized field adjustment as opposed to the mixing/swirling affect seen in Fig. 10a.
The other main aspect of Fig. 10 is the error in the von Mises stress in Fig. 10b. One interesting aspect is that the maximum error is 3.4 MPa is nearly identical to the initial design. The location of this stress, however, is not solely located on the boundary. There is also an area of large error at the hole. One interesting aspect is the comparison between the stress concentration and the location of large error on the hole. There appears to be a rotation between these locations. This rotation seems to be about 37.5 • , while the exact rotation is difficult to measure due to the FE approximation.
The final analysis is the Monte Carlo simulation to show the robustness of this design to uncertainty. For the location of maximum stress, the 95% confidence interval found is 295.3 ± 8.6 MPa. While the nominal stress is larger than for the initial design, the variance is nearly the same. To better understand how these two designs, have similar variances, the standard deviation field from this simulation is shown in Fig. 11. In general, this has the same features as the initial design, where the largest standard deviation is near the boundary with a maximum of 20.3 MPa compared to 20.7 MPa. However, one difference is that there is a local area of high error near the hole. This localization increases the variation in the confidence interval but is still small compared to the maximum. The fact that the initial design and the single hole design have similar variances is believed to be a trade-off between the Saint-Venant's principle of being far away from the loading and the nature of the stress concentration.

Triple webbing holes design
The final design presented in this work is a modification of the single webbing hole. This replaces the single large hole with three 5 mm holes. The inclusion of a smaller hole is predicted to introduce a larger stress concentration, while using multiple holes is used to investigate interactions between the fields/stress concentrations created by these multiple holes. The nominal stress and error of the von Mises stress results from this analysis can be seen in Fig. 12.
The first aspect of Fig. 12a that is noticeable is the maximum von Mises stress. This feature has a maximum von Mises stress of 315.1 MPa that greatly exceeds the yield criterion. By looking at the webbing, the maximum stress occurs at the top of one of the holes. While not modelled in this simulation, this large stress is expected to be near the ultimate strength and cause a crack to form, resulting in a global failure compared to the local failure seen with the single hole design. If this design was presented to an analyst, it is expected that this design would be rejected. Despite this rejection, there is a large amount of information that can be gathered from this design.
One of the other aspects of interest in Fig. 12b is the interaction between the holes. The holes are 5 mm in diameter with the centres distanced 10 mm apart. In general, these holes did not affect each other. In the lowest hole, the same swirling field as the single hole (shown in "Single webbing hole design" Section) is present with the same predicted reason based on the feature size.
The error of the von Mises stress in Fig. 12b shows some interesting aspects. Firstly, the maximum error is smaller than the other designs, although by a small amount. The cause of this decrease is unknown but shows an interesting trend where the larger the stress concentration, the smaller the error concentration. Some future work is planned to identify if this trend is accurate for the method or only this specific model. The stress concentration is solely geometrically dependent while the error concentration is both geometric and numerically dependent since the error is based on the mesh and is described in terms of a GOEE. Additionally, the same rotation of the stress and error concentration as the single hole is seen.
The final results for this design are from the Monte Carlo simulation. For the location of maximum stress, the 95% confidence interval is 319.0 ± 17.9 MPa. The interval has about double the spread compared to the other designs. This increase is especially interesting considering the standard deviation field shown in Fig. 13. For the entire design, the maximum standard deviation is 19.7 MPa compared to 20.7 MPa in the initial design. Despite this overall smaller variance, the maximum stress has a larger variance; this is due to the increased stress concentration in addition to the location being closer to the top edge.

Discussion
The results from all of these feature models show that the global model is well-defined and introduces very little error into the feature model. These results highlight several aspects; the first being that the total error is small compared to the nominal stress. This instructs the analyst that this modelling techniques is well defined for this loading, even at the stress concentration due to the holes. Another aspect is that the maximum error and the maximum stress do not necessarily coincide in physical space. For example in the initial design, the maximum stress occurs at the top of the webbing while the maximum error occurs near the intersection of the webbing and flange. While this is not a major issue, it is important to know since typical mesh refinement is performed at the stress concentration but might not be the location that has the largest stress value considering the effect of the error. An exact reason for this is currently unknown and is being studied, but the error does follow the expected pattern of being at the boundary of the feature model. The error measurement is based on the global mesh, so it follows logic that the error will be largest at locations that coincide with the global mesh.
One of the most interesting findings in this work is the differences between the feature designs that contain webbing holes. The initial design showed that the area of high stress was on the webbing. For designs that contain holes/stress concentrations in that region, the nominal stress increases to and past the von Mises yield criterion. Despite this increase in stress, the GOEE propagated into the feature model did not show many differences. The major difference is based on the GOEE field around the webbing holes. While the maximum error is nearly identical for the various designs, the single hole and triple holes have a GOEE concentration at the holes similar to the stress concentration. The one aspect that is different between the error and stress concentrations is that there is a rotation around the hole. It is difficult to specifically identify the rotation due to the FE resolution, but it appears to be approximately 37.5 • .
The other main finding is the comparison of the robustness of the various feature models. This is done with a Monte Carlo sampling of stochastic driving DoF. There are two major results from this simulation. The first is that the largest standard deviation of the feature model is located near the driven nodes. This showcases Saint-Venant's principle that the farther away from the loading conditions, the smaller the effects of variation are upon the stress. The second result is the effect of the feature design on the 95% confidence interval at the location of maximum von Mises stress. While the initial design and the single hole have approximately the same confidence interval, this is contributed to by a few factors. The first factor is the location itself. For the initial design, the maximum location is on the top part of the webbing near the driven nodes, while the single hole design is located near the hole. With the location near the hole, there is a trade-off between the decrease due to the distance from the edge and the increase due to the stress concentration. This trade-off is dominated by the stress concentration for the triple hole feature design due to the small size of the holes, thus producing the largest confidence interval and the largest von Mises stress.

Conclusions
The work presented in this paper introduces a novel approach to propagate the uncertainty from the global model using a custom element formulation into the feature model by utilizing dual problems and multi-scale modelling. This provides information to the design engineer about the failure of complex areas without the need to perform large mesh refinements that are computationally non-feasible. To demonstrate this methodology, this work takes a simplified beam section and evaluates four local feature designs. The global model experiences near yielding stress via the von Mises criterion, while the initial design, which contains two holes in the flange section, does not experience any yielding. In addition to the stress, the error in the von Mises stress showed little error near the maximum stress with the location of maximum error near the intersection of the webbing and flange sections of the feature. Future planned work is to introduce a Bayesian response surface through a Gaussian process to reduce the total number of performed dual problems, thus reducing the computational evaluation time with little to no decrease in accuracy.
In addition to the initial local feature design, three additional designs are tested in order to determine the sensitivity of the multi-scale modelling procedure and the error propagation. One design makes a modification to the initial holes, while the other two add additional holes to the webbing section. The nominal hole modification did not show any difference in the stress distribution and error field compared to the initial design. This is mainly due to the applied loading since the elements do not have through thickness stiffness due to the thin shell assumption. To fully test this methodology, the other two features contain holes near the location of maximum stress. These holes increase the maximum stress to near and past the yielding criterion resulting in one feature having gross yielding suggesting a failed design if this system would be considered for an industrial case. In addition, confidence in these results is quantified via the error and show high confidence in these results.
To test the robustness of the multi-scale methodology, some variance is added to the driving nodes to represent the use of sensors in the analysis that introduce errors not accounted for using the GOEE (such as instrumentation error, manufacturing tolerances, etc.). Using this variability, the 95% confidence interval for the location of maximum von Mises stress via a Monte Carlo simulation and the standard deviation field are compared between the different feature designs. These results show that the variability decreases as the distance from the driven nodes increases, showcasing Saint-Venant's principle. However, one aspect observed is that the greater the stress concentration, the larger the confidence interval despite not being located near the driven nodes where the largest variance is identified.