Frequency Response Analysis of Perforated Shells with Uncertain Materials and Damage

In this paper, we give an overview of the issues one must consider when designing methods for vibration based health monitoring systems for perforated thin shells especially in relation to frequency response analysis. In particular, we allow either the material parameters or the structure or both to be random. The numerical experiments are computed using the standard high order finite element method with stochastic collocation for the cases with random material and Monte Carlo for those with damaged or random structures. The results display a wide range of responses over the experimental configurations. In perforated shell structures, the internal boundary layers can play an important role especially when damage is allowed within the penetration patterns. The computational methodology advocated here can be used to build statistical databases that are necessary for development of probabilistic damage identification methods.


Introduction
Vibration based structural health monitoring remains a vibrant research area with an extensive body of literature. A review article from 2004 by Carden and others [1] alone cites over 150 publications. A more recent (2011) review article is by Fan and Qiao [2] concentrating on damage identification methods. To use the language of inverse problems, modelling of damage identification methods requires the forward part, i.e., simulation of the mechanical systems we are interested in, and the inverse part or identification of the structure or changes within. In this paper, we give an overview of the issues one must consider when designing methods for vibration based health monitoring systems for perforated thin shells, especially in relation to frequency response analysis.
Engineering constructions made of perforated materials are challenging from the computational mechanics point of view. Some examples of relevant industrial constructions are drying drums and Trommel screens. However, the literature on perforated materials is somewhat sparse. The introduction in Martikka et al. [3] serves as an excellent overview of the complexities of perforated shell designs-our interest-and their applications.
It is natural to seek a simplified problem via homogenization, a process where the material properties of a non-perforated structure are adjusted to match those of the perforated one. Recently, Jhung and his collaborators have published a number of papers within the nuclear reactor design context, see [4][5][6]. For a detailed and exhaustive case study for a single cylindrical shell submerged in fluid, see [7]. In this paper, however, we do not advocate such an approach. The homogenization of shell structures under loading and especially in dynamic settings is not well understood. The main culprit is that this class of engineering problems is a parameter-dependent problem, where the parameter is the (dimensionless) thickness. It is characteristic of such problems that there are local features in the solution that vary in size and amplitude depending on the parameter and one can always define a critical value of the parameter where the local features become dominant, and homogenization fails. In frequency response analysis, there is an added complication. Let us consider shells of revolution. It follows immediately that in free vibration the eigenmodes must have integer valued wave numbers in the angular part of the mode. This wave number in turn depends on the parameters of the problem, including the shell geometry and the material parameters; see [8]. Thus, successful homogenization would require matching of modes which would be a more difficult problem than simply associating the first ten modes, say, in both cases. In addition, the dependence on the shell geometry means that it is not feasible to homogenize a perforated plate (a sheet of metal) and then ignore the effects of the curvature changes.
Here, we are interested in the frequency response of perforated shells when either the material parameters such as Young's modulus or the structure or both are random in a well-defined way. The structural randomness is introduced via a simple damage model where the perforation patterns are perturbed with uniform intensity over the whole domain. We study two structures, a cylindrical shell which is of practical importance, and an example of different shell geometry, a perforated hyperbolic shell (cooling tower), which is more of an academic interest but serves to underline the central issues. It follows naturally that inclusion of uncertainty in the simulations leads to distribution of frequency responses that require statistical analysis [9].
The main contribution of this work can be divided into two interconnected parts: First, through understanding and analysis of the model problems, we can design a set of experiments that illustrate the expected phenomena, and, second, the experiments can be computed and the results analyzed within the standard stochastic framework. We see that in multi-parametric problems the responses can be widely varying as well. In perforated shell structures, the internal boundary layers can play an important role, especially when damage is allowed within the penetration patterns. The computational methodology advocated here can be used to build statistical databases that are necessary for development of probabilistic damage identification methods. The important observation is that at least in some cases the distributions of quantities of interest in such settings can be recovered and hopefully utilized in the design of, e.g., damage identification methods.
All simulations are computed with an hp-version of the finite element method [10,11], using standard elements. This choice is based on our desire to avoid any issues related to numerical locking.
The rest of the paper is structured as follows: In Section 2, we discuss the deterministic shell eigenproblems and issues specific to perforated structures; Section 3 sets the reference results for the free vibration of our model structures. The important concept of mixing of modes is demonstrated by comparing results over perforated and non-perforated shell geometries; frequency response analysis under uncertainty is briefly outlined in Section 4; our extensive set of numerical experiments is the topic of Section 5; Finally, conclusions are drawn in Section 6.

Preliminaries
In this section, the basic characteristics of perforated shell problems are introduced. Much of the material has been covered in the literature before, but here the notations have been normalized to conform to discussion in this paper. We assume that the reader is familiar with the finite element method and do not cover it here at all. It is also assumed that the existence of internal boundary layers in structures is at least accepted.
The shell problem is discussed as an eigenproblem, since it is the underlying problem despite the eventual formulation of the actual frequency response problem below. The dimensionally reduced shell model and the effect of the shell geometry are outlined in this context.
Special issues related to perforated domains are the penetration patterns, here with the damage model, and in the shell context the boundary layers. All three are briefly covered in this section as well.

Shell Eigenproblems
Assuming a time harmonic displacement field, the free vibration problem for a general shell leads to the following abstract eigenvalue problem: Find u ∈ R 3 and ω 2 ∈ R such that Above, u = {u, v, w} represents the shell displacement field, while ω 2 represents the square of the eigenfrequency. In the abstract setting, S and M are differential operators representing deformation energy and inertia, respectively. In the discrete setting, they refer to corresponding stiffness and mass matrices.

Shell Geometry
In the following, we denote x = (x 1 , x 2 , x 3 ) ∈ R 3 and let a profile function x 2 = f (x 1 ) revolve about the x 1 -axis. Naturally, the profile function defines the local radius of the shell. Shells of revolution can now formally be characterized as domains in R 3 of type where D is a (mid)surface of revolution, n(x) is the unit normal to D, and d(x) = d 1 (x) + d 0 (x) is the thickness of the shell measured in the direction of the normal. An illustration of such configurations is given in Figure 1. Throughout this paper, we consider parabolic and hyperbolic shells with profile functions f P (x 1 ) = 1 and f H (x 1 ) = 1 + (x 1 /π) 2 , respectively. This classification of surfaces follows from Gaussian theory of surfaces, where it holds that for general profile functions f (x 1 ) = 0 for parabolic, and f (x 1 ) > 0 for hyperbolic surfaces. Let us have d(x) = d (constant) and define principal curvature coordinates, where only four parameters, the radii of principal curvature R 1 , R 2 , and the so-called Lamé parameters, A 1 , A 2 , which relate coordinates changes to arc lengths, are needed to specify the curvature and the metric on D. The displacement vector field of the midsurface u = {u, v, w} can be interpreted as projections to directions where Ψ(x 1 , x 2 ) is a suitable parametrization of the surface of revolution, e 1 , e 2 are the unit tangent vectors along the principal curvature lines, and e 3 is the unit normal. In other words, Assuming that the shell profile is given by f (x 1 ), the geometry parameters are and We can now define the dimensionless thickness t as t = d/R, where R ∼ R 2 is the unit length. In the following we assume that R ∼ 1 and thus use the two thickness concepts denoted by d and t interchangeably.

Reduction in Thickness
The free vibration problem for a general shell (1) in the case of shell of revolution with constant thickness t leads to the following eigenvalue problem: Find u(t) and ω 2 (t) ∈ R such that Again, u(t) represents the shell displacement field, while ω 2 (t) represents the square of the eigenfrequency. The differential operators A m , A s , and A b account for membrane, shear, and bending potential energies, respectively, and are independent of t. Finally, M(t) is the inertia operator, which in this case can be split into the sum M(t) = t M l + t 3 M r , with M l (displacements) and M r (rotations) independent of t. Let us next consider the variational formulation of problem (7). Accordingly, we introduce the space V of admissible displacements, and consider the problem: where a m (·, ·), a s (·, ·), a b (·, ·) and m(t; ·, ·) are the bilinear forms associated with the operators A m , A s , A b and M(t), respectively. Obviously, the space V and the three bilinear forms depend on the chosen shell model (see for instance [12]).

Two-Dimensional Model
Our two-dimensional shell model is the so-called Reissner-Naghdi model [12], where the transverse deflections are approximated with low-order polynomials. The resulting vector field has five components u = (u, v, w, θ, ψ), where the first three are the standard displacements and the latter two are the rotations in the axial and angular directions, respectively. Here, we adopt the convention that the computational domain D ⊂ R 2 is given by the surface parametrisation and the axial/angular coordinates are denoted by x and y: Deformation energy A(u, u) is divided into bending, membrane, and shear energies, denoted by subscripts B, M, and S, respectively: Notice that we can safely cancel out one power of t.
Bending, membrane, and shear energies are given as [13] where ν is the Poisson ratio (constant), and E(x, y) is the Young's modulus with scaling 1/(12(1 − ν 2 )). The symmetric 2D strains, bending (κ), membrane (β), and shear (ρ), are defined as In our experiments, the shell profiles are functions of x only, but not necessarily constant, hence the strains are presented in a general form. The stiffness matrix S is obtained after integration and assembly. Here, the mass matrix M is the standard 2D one with density σ except for the scaling of the rotation components and surface differential A 1 (x, y)A 2 (x, y).

Perforated Domains
Perforated domains are characterized by the penetration patterns which in turn depend on the underlying manufacturing processes and the related hole coverage, typically given as a percentage. Here, we consider standard structured patterns and extend the concepts to quasi-random patterns. We also consider a simple randomized damage model.

Penetration Patterns
The quantity used to characterize perforated sheets of metal is the ligament efficiency η. Let us assume that the holes are ellipses with a, b as the horizontal and perpendicular semi-axis, and the separation of the centres is P x and P y , respectively. Following [5,14,15], we define horizontal and perpendicular ligament efficiency, denoting them η x , η y , respectively. For regular arrays of holes and for triangular arrays, allowing for alternating layers, For circular holes, the radius r = a = b, of course, and further if the pattern is regular η = η x = η y . Both pattern types are illustrated in Figure 2. Notice that the triangular pattern in the figure has a tighter packing than that implied by Equation (23). An intriguing option is to replace a regular pattern with one induced by a quasi-random process. Notice that quasi-random means that the process is deterministic which has implications to manufacturing. One of the simplest sequences of points generated by such a process is the Halton sequence [16]. Definition 1 (Halton). Let i ∈ N 0 and basis b ∈ N be such that b ≥ 2. Then, the inverse representation function φ b (i) is given by Let q 1 , q 2 , . . . , q s be the s first primes. The s-dimensional Halton sequence is The choice of the primes is crucial, and, in the context of this paper, only two have to be chosen. The effect of this selection is shown in Figure 3, where the local structures of two sequences are compared. As can be seen from the definition of the Halton sequences, in 2D, the points populate the unit square [0, 1] 2 and have to be mapped to the actual computational domain. We do not attempt to define averaging ligament efficiencies for such penetration patterns. However, a useful concept used to relate different kinds of patterns is the cell size , which is the size of the cell containing exactly one hole on average. In the case of a regular n × n-grid on a shell of revolution over (x, y) ∈ [−π, π] × [0, 2π], we have = 2π/n. In the limit as → 0, we expect regular and quasi-regular penetration patterns to converge.

Damage Model
Since we are committed to simulating the systems with faithful geometry description, we can similarly adopt a destructive damage model, where cracks between two holes are modelled simply by removing the area between the respective holes. This is illustrated in Figure 4. The model is straightforward: first, a sample X is chosen with intensity µ from the set of holes C. Every hole x i ∈ X is paired with one of its four nearest neighbors, and together they form a damaged region d k . As can be seen in Figure 4b,d, it can happen that, for both holes in a pair d k = (x i , x j ), it holds that x i ∈ X and x j ∈ X, i = j, and thus the damage can spread over a larger region. Thus, one must ensure that the domain remains simply connected.

Layer Chart over Representative Cells
Even though the shell problems are elliptic partial differential equations, the solutions may exhibit local features, formally boundary layers, with characteristic length scales much larger than the chosen cell size . Therefore, it is necessary to review the possible boundary layers that may occur in different scenarios.
For our discussion, it is useful to define the concept of boundary layer generators first (see [17]). The layer generators are independent of the length scale of the problem under consideration.
Definition 2 (Layer Generator). The subset of the domain from which the boundary layer decays exponentially, is called the layer generator. Formally, the layer generator is of measure zero.
In perforated structures, the hole boundaries are natural layer generators. The boundary layers are not concentrated only around the holes but can emanate along the characteristics of the shell surface. Elliptic, parabolic, and hyperbolic structures each possess a distinctive set of layer deformations. The layer structure is classically assumed to be an exponential solution to the homogeneous Euler equations of the shell problem. In Pitkäranta, Matache, and Schwab [17], it is shown using the Ansatz u(s, r) = Ue λs e ikr (24) that solutions with Reλ < 0 such that the characteristic lengths L = 1/Reλ → 0 are of the form L ∼ d 1/n , where n ∈ {1, 2, 3, 4}. Here, s is the coordinate orthogonal to the layer generator and r the one along it. From these, the layer with n = 2 is present in all geometries, whereas layers with n = 3 and n = 4 are present only in hyperbolic and parabolic geometries, respectively, and thus play a role here as well. The case n = 1, i.e., the shortest one arises from a shear deformation and can be captured by the shallow shell formulation [17], if necessary. In the 3D visualizations for the very thin case (Figures 5b and 6b), distinct features of the solution, that is, the boundary layers, are clearly visible. In Figure 7, the layer chart is given both for the parabolic and hyperbolic representative cells, with the exception of the shortest one, which is present on every boundary.    The characteristic length scales are indicated in terms of the dimensionless thickness. In the hyperbolic case for the non-axial layers, the directions depend on the actual geometry, in this case, the decay is in the orthogonal direction to the one indicated by the arrow.

Free Vibration of Reference Configurations
In this section, we establish the reference configurations on which the numerical experiments below are based and explain the rationale for these experiments. As already given above, the profile functions are f P (x 1 ) = 1 and f H (x 1 ) = 1 + (x 1 /π) 2 , for the parabolic and hyperbolic cases, respectively.
In Table 1, the first 10 eigenpairs are tabulated for both non-perforated shells with t = 1/100. In the two other tables (Tables 2 and 3), the first 10 eigenpairs are listed for perforated configurations, including quasi-random ones. The eigenmodes are characterized using axial (K x ) and angular (K y ) wave numbers, where the fractional axial wave numbers refer to symmetric modes that is, the pair (K x , K y ) = (1/2, 4) refers to mode ϕ(x, y) ∼ cos(x/2) cos(4y) or ϕ(x, y) ∼ cos(x/2) sin(4y) due to symmetry. The numbers are universal, that is, no material parameters are included except for the Poisson ratio ν = 1/3. For non-perforated shells, every eigenvalue (disregarding torsional modes) is a double one due to symmetry. Theoretically, the lowest mode can have an order four, but this occurs only at critical values of t which is not the case here. Even for perforated shells, the eigenvalues are observed to appear in clusters of two.
The discrepancies between the eigenvalues in a cluster are a simple measure of homogeneity. In the Halton cases, the discrepancies are larger than in the regular case but still remarkably small.

Remark 1.
Even though the regular perforation pattern results in a symmetric configuration, there exist clusters with relatively large discrepancy. If the pattern is a n × n-grid, then the modes with angular wave number n, for instance, belong to such a cluster.
Two further observations: first, the reference eigenvalues cannot be simply scaled to match the perforated ones, and, second, the modes are observed out of order. For homogenization, it would be necessary to perform minimization over a set of perforations to get the eigenvalues reasonably approximated, but similar process for the eigenmodes cannot be performed unless the subspaces are the same. At the moment, we are not aware of any efficient way short of computing the whole eigenspaces to accomplish this. Based on these results, we do not attempt to select modes in the frequency response a priori, but compute the results using full systems regardless of the cost.

Frequency Response Analysis under Uncertainty
In this section, we consider the stochastic version of the shell eigenproblem, which is obtained by introducing uncertainties in the material properties or geometry of the shell. We discuss key properties of the problem and present a framework for frequency response analysis under uncertainty. A stochastic collocation algorithm combined with statistical sampling is proposed for computing the expected value as well as higher moments for the quantity of interest.

The Stochastic Eigenproblem
As before, let D ⊂ R 2 denote the computational domain. Assume that the Young's modulus is a random field expressed in a series of the form where M ∈ N and ξ = (ξ 1 , ξ 2 , . . . , ξ M ) is a vector of mutually independent random variables. For simplicity, we assume here that each ξ m is uniformly distributed in some interval of R and, therefore, after suitable rescaling, the random vector ξ takes values in the domain Γ := [−1, 1] M . A parametrization of the form (25) is standard in the literature and may be seen to result from a Karhunen-Loéve expansion of the underlying random field E(x, y, ξ), see e.g., [18]. We assume that E(x, y, ξ) is strictly uniformly positive and uniformly bounded, i.e., there exist E min , E max > 0 such that Given an integrable function v = v(ξ), we denote by its expectation taken over the parameter space Γ. Assuming the stochastic model (25) for the Young's modulus, the eigenpairs of the free vibration problem now depend on ξ ∈ Γ. The stochastic version of the eigenvalue problem is obtained simply by replacing the differential operators in Equation (7) with their stochastic counterparts and assuming that the resulting equation holds for every realization of ξ ∈ Γ. In variational form, the problem reads: find functions u(t, ·) : Γ → V and ω 2 (t, ·) : Γ → R such that, for all ξ ∈ Γ, we have ta m (ξ; u(t, ξ),v) + ta s (ξ; u(t, ξ), v) + t 3 a b (ξ; u(t, ξ), v) Here, a m (ξ; ·, ·), a s (ξ; ·, ·) and a b (ξ; ·, ·) are stochastic equivalents of the deterministic bilinear forms in Equation (8).
After integration and assembly, we obtain the stochastic stiffness matrix S(ξ), which can now be expressed as Here, each matrix S (m) corresponds to a term in the series (25). Observe that the mass matrix M does not depend on the parameter vector ξ ∈ Γ.

Remark 2.
From the free vibration examples of the previous section, it should be clear that special care must be taken in matching of eigenmodes. When the problem is cast into the stochastic setting, the eigenvalues are bound to cross within the parameter space, and therefore it makes sense to match eigenmodes according to their axial and angular waveunumbers K x and K y rather than by simple increasing enumeration. This has previously been illustrated in [8,9].

Frequency Response Analysis
In frequency response analysis, the idea is to study the excitation or applied force in the frequency domain. Therefore, the uncertainty in the eigenproblem is directly translated to frequency response. Given ξ ∈ Γ the procedure is as follows [19]: We start from the equation of motion for the system where M is the mass matrix, C the viscous damping matrix, S(ξ) the stochastic stiffness matrix, f the (deterministic) force vector and v(ξ) the displacement vector. The matrices M and S(ξ) are defined as above, whereas the viscous damping matrix is taken to be a constant diagonal matrix C = 2δI, where δ = 1/2000. The choice of the viscous damping matrix is simply for convenience. There are many problem-dependent choices such as Rayleigh damping where C would be some linear combination of M and S(ξ).
In the case of harmonic excitation, a steady-state solution is sought and the force and the corresponding response can be expressed as harmonic functions as Taking the first and second derivatives of Equation (28) and substituting using (29) leads to thus reducing to a linear system of equations: Any quantity of interest can then be derived from the solutionv(ξ, ω). In the following, we consider the mechanical energy defined by the natural energy norm as the quantity of interest in matrix notation This is a mathematically reasonable choice, yet in practice it may be difficult to measure. We have ruled out the use of maximal deflections since the effects of the boundary layers should then be taken into account.

Stochastic Collocation
We consider an anisotropic Smolyak-type collocation operator defined with respect to a finite multi-index set (see e.g., [20][21][22][23]). In principle, a simpler formulation would also suffice for the purpose of our numerical experiments. However, we wish to keep the general form of the operator here, since it allows for efficient computation even in high-dimensional parameter spaces.
Let L p be the univariate Legendre polynomial of degree p. Denote by {χ Now, let A ⊂ N M 0 be a finite set of multi-indices. For α, β ∈ A, we write α ≤ β if α m ≤ β m for all m = 1, . . . , M. We assume that A is monotone In the following, sense: The sparse collocation operator is defined as where I (m) −1 := 0. This may be rewritten in a computationally more convenient form The collocation points are of the form for some γ ∈ N M 0 such that γ ≤ α ∈ A. Similarly, we define the tensorized quadrature weights w Statistics, such as the expected value and variance, for the collocated solution may now be computed by applying the quadrature rule on the terms in (34).
The accuracy of the collocated approximation is ultimately determined by the smoothness of the solution as well as the choice of the multi-index set A ⊂ N M 0 . We refer to [21,23] for a detailed analysis.

Computing Statistics under Hybrid Uncertainty
Stochastic collocation is an efficient method for resolving the effects of the random material model. However, due to the nature of our damage model, the effects of the random deformations to the actual computational domain are easier to resolve using traditional statistical sampling, in other words the Monte Carlo method. In Algorithm 1, we have outlined a hybrid strategy for computing the expected value for some quantity of interest Φ(v) derived from the solutionv =v(ξ), when both the material parameters and the computational domain are random. Higher moments may be computed in a similar way.

Algorithm 1
Let S ∈ N and A ⊂ N M 0 . For j = 1, 2, . . . , S do (1) Construct a computational domain D ⊂ R 2 using the damage model from Section 2.2.2.
Return the sample average as the approximate expected value.

Numerical Experiments: Frequency Response
In order to reduce combinatorial complexity and thus focus the message, let us fix some parameters over the whole set of experiments. As throughout this paper, the representative parabolic profile function is f P (x 1 ) = 1 and the corresponding hyperbolic one f H (x 1 ) = 1 + (x 1 /π) 2 . In all cases, we let x 1 ∈ [−π, π] so that the 2D computational domain is Material constants adopted for all simulations are: E = 2.069 × 10 11 MPa, ν = 1/3, and ρ = 7868 kg m −3 , unless otherwise specified. All numerical experiments are solved with the hp-FEM using standard elements at p = 3 or 4. We select the 20 × 20 regular penetration pattern as the representative configuration for both geometries, that is, let = π/10 with ligament efficiency η = 0.61, and hole coverage of 12%. This corresponds to a hole diameter 2r where the radius r = √ 3π/50 ≈ 0.061, with (20 × 20 × πr 2 )/(4π 2 ) = 0.12 exactly. The corresponding Halton grid H(7, 37) is used in one of the experiments. The dimensionless thickness is constant t = 1/100, which is a realistic value in engineering context. The loads or excitation functions act only on the transverse direction w and are simply Fourier modes in the angular direction and constant in the axial direction, i.e.,f (x, y) = C cos(Ky) N, where C = 1000 is a scaling parameter, and K is the angular wave number (We used K y above). The random Young's modulus depends only on the axial direction and has the form where we choose E 0 (x) = E (constant), and E m (x) = √ λ m E sin(mx) with λ m = 1/(m + 1) 4 . Finally, in the damage model, the intensity µ is constant over the whole domain.
The experiments progress in a natural order: first, we only consider the effect of the excitation (or loading) in a deterministic setting, second, the material parameter, Young's modulus is assumed to be random, and the responses are computed at the collocation points, third, the damage model is introduced with deterministic materials, and finally the damage model is combined with material uncertainty to complete the sequence so that the set of experiments spans the different possibilities discussed above.
On average, the meshes have 6000 triangles corresponding to 150,000 d.o.f. at p = 3 and 250,000 at p = 4. One full solve at p = 3 took roughly 2 min on Aalto cluster with Intel Xeon processors (varied models) (Santa Clara, CA, USA). In other words, one full experiment with uncertain material parameters required one hour, and the last experiment with 100 such instances 100 h.

Effect of Loading
First, we study the effect of the angular wave number in the loading in a purely deterministic setting. We choose three cases, K = 1/2, 4, and 20, and let p = 4. The last one is chosen to illustrate the special case where the related cluster has the largest discrepancy. In the example, this is reflected in differing responses to trigonometric basis functions with the same wave number.
The angular frequency range is for the first two cases ω = 2π f Hz, f = 5, 10, . . . , 200, and for K = 20, f = 100, 200, . . . , 2000, where f denotes the ordinary frequency. The responses are shown in Figures 8 and 9, respectively. Indeed, for K = 20, we observe the curious effect that, despite rotational symmetry, the frequency responses for cosine and sine-type loadings are different. This does not occur at K = 1/2 or 4. The reason for this is illustrated in Figure 10; with regular penetration patterns, it is possible that either one, cosine, say, has maximal amplitude exactly at the holes, whereas the other one has very small ones. Thus, the observed energies can have a difference of an order of magnitude. In the following, this effect does not play a role since we are content to use K = 4 only. However, this effect is real and its importance depends on the overall configuration, of course.  The frequency response graphs show typical resonance patterns. The responses clearly depend not only on the loading, but also on the shell geometry. The observed maximal energies differ considerably, and, for a fixed frequency interval with K = 1/2, the highest one is hyperbolic, but with K = 4 parabolic. Corresponding energies have an order of magnitude difference (see Figure 9).

Effect of the Random Material Parameter
In the second experiment, we let the Young's modulus be random. We let M = 3, and use full three point Gaussian quadrature in all stochastic dimensions. Thus, the associated collocation operator is given by (33)  The analysis of the results can now be performed on two levels, either by collocation or by sample statistics. Indeed, these results are very close, but of course not the same. See Figures 11 and 12. We can take one step further and compute higher moments and even plot them. In Figures 13 and 14, we show two solutions (expected solutions) with highest energies and their respective variances, one on a regular grid and another one on the Halton grid. Notice that for both geometries the variances in the regular grid cases reflect the fact that the random material parameter depends only on the axial direction.
In Figure 14c, the response on the Halton grid is shown in the hyperbolic case. Interestingly, the irregular grid leads to a response that has an approximate angular wave number K H = 7, even though the excitation has K = 4. This is an example where it is very difficult to predict the number of modes for the subspace if one wanted to apply some reduced basis method. In Table 3, we see that in the hyperbolic case the lowest modes are in a tight cluster and therefore it is likelier than in the parabolic case that perturbations result in higher order responses.

Effect of Local Low Intensity Damage
Next, we only perturb the penetration pattern using the damage model described above. Here, two interesting events can take place. It can happen that the damaged model leads to a response that is a close approximation of one the eigenmodes. An example of this is shown in Figure 15b,c. On the other hand, the damage may be close to the boundary and the associated boundary layer can result in local high energy feature in the response.
In Figures 15a and 16a, the standard response profiles are shown. Again, a rare, but significant response is indicated with a large standard deviation (i.e., large variance). In Figure 15a, we see that the important frequencies are f = 40 and f = 110. We interpret the differences in the responses at the two frequencies so that it is likelier that there is a response at f = 40, but at f = 110 a very energetic response is possible but rare.
Of particular interest is the hyperbolic case of Figure 16c. Even though the loading is constant in the axial direction and periodic in the angular, the transverse deflection clearly shows internal layers (refer to Figure 7 above). The damaged holes perturb the regular pattern and generate internal layers which travel along the characteristics of the surface. Since the shells are surfaces of revolution, the characteristic lines align themselves close to the boundary. In fact, a similar phenomenon occurs for the parabolic case in Figure 15c but is obscured by the effect of the axially aligned loading.

Hybrid Uncertainty
Finally, we consider the full model where both the material parameter and the damage are included in the uncertainty quantification. Since the damage model does not preserve the symmetries, it follows that the variances of the highlighted solutions do vary over the whole domain and the structure of for instance Figure 13a is not preserved. In the spirit of Algorithm 1, we simulate S = 100 damage models each with a random material parameter as above leading to 2900 solutions per geometry type. The frequency responses and highlighted cases are given in Figures 17 and 18. Again, the internal boundary layers are clearly visible in the transverse deflections, and even in the variance plots as well! In Figure 19, the energy distributions at the frequencies where the most energetic cases occur are shown. The distributions are exponential, which makes design of reliable damage detection algorithms challenging. However, this underlines the fact that comprehensive statistical understanding of the problem class is required, and systematic use of the stochastic finite element method and related techniques is well motivated.

Conclusions
For perforated shell structures, frequency response analysis is complicated by the wide range of parameters affecting the results. Indeed, shell is not referred to as the primadonna of structures for nothing. For instance, one has to have a good a priori understanding of the boundary layer structure, since, for sufficiently thin shells, the local features dominate and the damage detection algorithm has to be selected or tuned accordingly.
In general cases, any successful homogenization procedure should be able to match the subspaces between the perforated and non-perforated configurations. At the moment, solutions to this problem are at least as complex as direct simulations, which is the approach chosen here. Adding material uncertainties and simple damage models makes this problem even harder to resolve. In specific cases where the relevant parameters and their ranges can be reliably identified, homogenization can be accomplished with less effort.
Nevertheless, the simulations indicate that statistical and/or probabilistic modeling techniques are likely to be the best candidates for widely applicable damage detection methods in this context. Through reliable computations, it is possible to establish good approximations to distributions of the quantities of interest. The computational methodology advocated here can be used to build statistical databases that are necessary for development of probabilistic damage identification algorithms.