Low-Rank Approximation of Frequency Response Analysis of Perforated Cylinders under Uncertainty

Abstract

Frequency response analysis under uncertainty is computationally expensive. Low-rank approximation techniques can significantly reduce the solution times. Thin perforated cylinders, as with all shells, have specific features affecting the approximation error. There exists a rich thickness-dependent boundary layer structure, leading to local features becoming dominant as the thickness tends to zero. Related to boundary layers, there is also a connection between eigenmodes and the perforation patterns. The Krylov subspace approach for proportionally damped systems with uncertain Young’s modulus is compared with the full system, and via numerical experiments, it is shown that the relative accuracy of the low-rank approximation of perforated shells measured in energy depends on the dimensionless thickness. In the context of frequency response analysis, it then becomes possible that, at some critical thicknesses, the most energetic response within the observed frequency range is not identified correctly. The reference structure used in the experiments is a trommel screen with a non-regular perforation pattern with two different perforation zones. The low-rank approximation scheme is shown to be feasible in computational asymptotic analysis of trommel designs when the proportional damping model is used.

1. Introduction

Cylindrical shell structures are common in engineering. Within computational mechanics, they belong to the class of thin solids, the class which arguably remains computationally one of the most challenging ones. Here, the focus is on perforated cylinders—more specifically, trommel screens and their frequency response. Perforated cylinders are used, for instance, in mass dryers and as screening devices, and are naturally subjected to large loads, which are often of impact type or otherwise concentrated [1,2]. Ideally, one would like to be able to analyse perforated structures as homogenised, where the effect of the perforation pattern is transferred to material parameters, the so-called effective material parameters [3,4,5].

The structural models used in shell analysis are dimensionally reduced, where the thickness is treated as a (dimensionless) parameter. Somewhat surprisingly, there is an intricate interplay between the perforation patterns and the thickness, which sets limits to homogenisation both under static loading and in dynamical eigenmode setting; see [6,7]. In this study, no attempts at homogenisation are made, and all simulations are carried out in exact geometry.

Frequency response analysis of thin solids is a multiparametric problem. In the structures considered here, uncertainty is included in the material properties; hence, the parameters can be divided into the stochastic parameters defining the material properties (Young’s modulus), and the deterministic parameters, i.e., the dimensionless thickness and the frequency considered. From the point of view of the analysis, no changes are necessary even if the perforation patterns were randomly perturbed as well. Many of these issues have been addressed in a previous study [8]. Modern stochastic analysis of multiparametric problems leads to a large number of simulations. Low-rank approximation can be a highly effective way of reducing the computational complexity. The Krylov subspace approach pioneered by Freund is the one considered here [9]. Naturally, it becomes necessary to establish limits within the parameter space when such approximations are applicable in the problem class of interest.

Frequency response analysis is of interest in many engineering fields, and hence the body of literature is not only vast but also scattered over a large number of journals. On structural health monitoring alone, two review articles cover nearly 300 publications [10,11]. In this paper, the underlying mathematical models are taken to be linear, and the materials isotropic in order to simplify the discussion. As examples of more complicated configurations in related problems, Wang et al. consider vibration composite beams under random loads [12] and Vu et al. corrugated cylinders and nonlinear buckling [13]. A very recent paper by Fu et al. [14] provides references to papers on complex vibration problems and related uncertainty quantification. A well-known paper by Worden et al. [15] highlights specific problematic issues in the uncertainty quantification of nonlinear systems. In the context of collocation methods such as those used here, the fundamental question is whether the interpolation property of the system is preserved. It is reasonable to expect that, within linear elasticity, this is never violated.

The two foundational papers in introducing the Krylov subspace method in this context are by Eid et al. [16] and Han [17]. Three more widely cited papers building on the connection with the proportionally damped systems are Johnson and Wojtkiewicz [18], Wu et al. [19], and Sepahvand and Marburg [20], where also material uncertainties are included in a similar way as here. Wu et al. have continued with a series of papers [21,22]. Moreover, Delissen et al. [23] have studied this approach in the context of topology optimisation, where also the computational complexity can increase significantly.

The main theoretical contribution of this work is that, through carefully designed experiments, it is shown that the relative accuracy of the low-rank approximation of perforated shells depends on the dimensionless thickness. This is a consequence of the boundary layers—that is, local features becoming dominant as the thickness tends to zero. This is precisely the same phenomenon which affects homogenisation. Determining the limits of low-rank approximation requires extensive and possibly expensive simulations. Here, it is shown how computational asymptotic analysis becomes feasible if the damping models and the related low-rank approximation scheme are chosen in the computationally most advantageous way.

The reference problem is the trommel screen. Trommel screens are devices used in the grading of raw materials and solid waste classification. Consider a perforated cylinder as in Figure 1, but tilted so that the smaller holes are higher, and the whole cylinder rotating. As the raw materials are fed into the trommel from the top, the perforation structure separates objects of different sizes. What makes this device an interesting example is that the perforation patterns are not regular either in size or distribution. For shells of revolution, the eigenmodes have integer-valued wave numbers in the angular direction if the material properties are isotropic, and since this condition is not satisfied, the frequency responses of trommel screens can be more varied. The quantity of interest is the total energy of the response, since there are no guarantees for the convergence of local features such as pointwise displacement. In Figure 1, some of the effects of the deterministic parameters are illustrated.

The rest of the paper is as follows: First, the shell eigenproblem and the related frequency response problem are reviewed also in the stochastic setting; next, some recent results on the free vibration of perforated results are discussed, including the lack of predictive power in the general perforation setting. In Section 4, the low-rank approximation approach is outlined, including the powerful preconditioner necessary for the iterative construction of the reduced basis; finally, an extensive set of numerical experiments is discussed, before the conclusions in Section 6.

2. Preliminaries

This section introduces both the shell eigenproblems and the collocation scheme, which is used to compute the frequency responses in the cases where the material parameters, e.g., Young’s modulus, are random. The discussion below follows closely that in [8].

2.1. Shell Eigenproblems

The free vibration problem for a general shell under the assumption of a time harmonic displacement field leads to the following abstract eigenvalue problem: Find $\mathbf{u}\in {\mathbb{R}}^3$ and ${\omega }^2\in \mathbb{R}$ such that

$$\left\{\begin{array}{c}\mathbf{S}\mathbf{u}={\omega }^2\mathbf{M}\mathbf{u}\hfill \\ +\mathrm{boundary}\mathrm{conditions}.\hfill \end{array}\right.$$

(1)

Above, $\mathbf{u}=\{u,v,w\}$ represents the shell displacement field, while ${\omega }^2$ represents the square of the eigenfrequency. In the abstract setting, $\mathbf{S}$ and $\mathbf{M}$ are differential operators representing deformation energy and inertia, respectively. In the discrete setting, they refer to corresponding stiffness and mass matrices.

Reduction in Thickness

In the case of a shell of revolution with constant thickness t, the free vibration problem for a general shell (1) can be stated as the following eigenvalue problem: Find $\mathbf{u}\left(t\right)$ and ${\omega }^2\left(t\right)\in \mathbb{R}$ such that

$$\left\{\begin{array}{c}t\mathbf{A}{}_M\left(t\right)+t{\mathbf{A}}_S\mathbf{u}\left(t\right)+t^3{\mathbf{A}}_B\left(t\right)={\omega }^2\left(t\right)\mathbf{M}\left(t\right)\mathbf{u}\left(t\right)\hfill \\ +\mathrm{boundary}\mathrm{conditions}.\hfill \end{array}\right.$$

(2)

Again, $\mathbf{u}\left(t\right)$ represents the shell displacement field, while ${\omega }^2\left(t\right)$ represents the square of the eigenfrequency. The differential operators ${\mathbf{A}}_M$, ${\mathbf{A}}_S$ and ${\mathbf{A}}_B$ account for membrane, shear, and bending potential energies, respectively, and are independent of t. Finally, $\mathbf{M}\left(t\right)$ is the inertia operator, which in this case can be split into the sum $\mathbf{M}\left(t\right)=t{M}^l+t^3{M}^r$, with $M^l$ (displacements) and $M^r$ (rotations) independent of t.

The finite element method is based on the variational formulation of problem (2). Accordingly, when the space V of admissible displacements is introduced, the problem becomes: Find $(\mathbf{u}\left(t\right),{\omega }^2\left(t\right))\in V\times \mathbb{R}$ such that

$$\begin{array}{c}\hfill t{a}_m(\mathbf{u}\left(t\right),\mathbf{v})+t{a}_s(\mathbf{u}\left(t\right),\mathbf{v})+t^3{a}_b(\mathbf{u}\left(t\right),\mathbf{v})=\\ \hfill {\omega }^2\left(t\right)m(t;\mathbf{u}\left(t\right),\mathbf{v})\forall \mathbf{v}\in V,\end{array}$$

(3)

where $a_m(·,·)$, $a_s(·,·)$, $a_b(·,·)$ and $m(t;·,·)$ are the bilinear forms associated with the operators ${\mathbf{A}}_M$, ${\mathbf{A}}_S$, ${\mathbf{A}}_B$ and $\mathbf{M}\left(t\right)$, respectively. Obviously, the space V and the three bilinear forms depend on the chosen shell model (see, for instance, [24]).

The energy norm $\left|\right||·|\left|\right|$ is defined in a natural way in terms of the deformation energy $\mathcal{E}\left(\mathbf{u}\right):=\left|\right|\left|\mathbf{u}\right|\left|\right|:=\sqrt{\mathbf{A}(\mathbf{u},\mathbf{u})}.$ Similarly, for (squared) bending, membrane, and shear energies, e.g., $\mathcal{B}\left(\mathbf{u}\right):=t^2{\mathbf{A}}_B(\mathbf{u},\mathbf{u})$. The shell model used throughout the paper, the Reissner–Naghdi model, is outlined in Appendix A.

2.2. Frequency Response Analysis under Uncertainty

The stochastic version of the shell eigenproblem is obtained by introducing uncertainties in the material properties or geometry of the shell. The key properties of the problem and a framework for frequency response analysis under uncertainty are outlined in this section. For computing the expected value as well as higher moments of the quantity of interest, a stochastic collocation algorithm is proposed.

2.2.1. The Stochastic Eigenproblem

Let $D\subset {\mathbb{R}}^2$ denote the computational domain (see Appendix A). Assume further that the Young’s modulus is a random field expressed in a series of the form

$$E(x,y,\xi )=E_0(x,y)+\sum _{m=1}^M{E}_m(x,y){\xi }_m,(x,y)\in D,$$

(4)

where $M\in \mathbb{N}$ and $\xi =({\xi }_1,{\xi }_2,\dots ,{\xi }_M)$ is a vector of mutually independent random variables. For simplicity, it is assumed here that each ${\xi }_m$ is uniformly distributed in some interval of $\mathbb{R}$ and, therefore, after suitable rescaling, the random vector $\xi$ takes values in the domain $\mathrm{\Gamma }:={[-1,1]}^M$. A parametrisation of the form (4) can be seen as derived from a Karhunen–Loéve expansion of the underlying random field $E(x,y,\xi )$; see, e.g., [25], in which case the functions $E_m(x,y)$ are the eigenfunctions of some covariance function $C(x,y)$. These eigenfunctions are known explicitly only in some simple geometric settings, but can also be resolved numerically as eigenvectors of the sampled covariance matrix; see, e.g., Chapter 8 in Ref. [26] and [25,27]. The series has to be truncated for numerical computations, of course. In the sequel, M terms are used in applications.

The standard assumption is that $E(x,y,\xi )$ is strictly uniformly positive and uniformly bounded, i.e., there exist $E_{min},E_{max}>0$ such that

$$E_{min}\le \underset{(x,y)\in D}{\mathrm{ess}\inf }E(x,y,\xi )\le \underset{(x,y)\in D}{\mathrm{ess}\sup }E(x,y,\xi )\le E_{max}\forall \xi \in \mathrm{\Gamma }.$$

(5)

Given an integrable function $v=v\left(\xi \right)$, denote by ${\mathbb{E}}_{\Gamma }\left[v\right]=2^{-M}{\int }_{\mathrm{\Gamma }}v\left(\xi \right)d\xi$ its expectation taken over the parameter space $\mathrm{\Gamma }$.

The uncertainty is introduced to the problem via material properties—more specifically, by representing the Young’s modulus using the stochastic model (4) above. As a consequence, the eigenpairs of the free vibration problem depend on $\xi \in \mathrm{\Gamma }$. Obtaining the stochastic version of the eigenvalue problem becomes simple: one only has to replace the differential operators in (2) with their stochastic counterparts and assume that the resulting equation holds for every realisation of $\xi \in \mathrm{\Gamma }$. In variational form, the problem is stated as: Find functions $\mathbf{u}$ and ${\omega }^2(t,·):\mathrm{\Gamma }\to \mathbb{R}$ such that for all $\xi \in \mathrm{\Gamma }$

$$\begin{array}{cc}\hfill t{a}_m(\xi ;\mathbf{u}(t,\xi ),\mathbf{v})& +t{a}_s(\xi ;\mathbf{u}(t,\xi ),\mathbf{v})+t^3{a}_b(\xi ;\mathbf{u}(t,\xi ),\mathbf{v})\hfill \\ & ={\omega }^2(t,\xi )m(t;\mathbf{u}(t,\xi ),\mathbf{v})\forall \mathbf{v}\in V.\hfill \end{array}$$

(6)

Here, $a_m(\xi ;·,·)$, $a_s(\xi ;·,·)$ and $a_b(\xi ;·,·)$ are stochastic equivalents of the deterministic bilinear forms in (3).

The stochastic stiffness matrix $\mathbf{S}\left(\xi \right)$ is obtained after integration and assembly. Using (4), one has

$$\mathbf{S}\left(\xi \right)={\mathbf{S}}^{\left(0\right)}+\sum _{m=1}^M{\mathbf{S}}^{\left(m\right)}{\xi }_m,\xi \in \mathrm{\Gamma },$$

(7)

where each matrix ${\mathbf{S}}^{\left(m\right)}$ corresponds to a term in the series. Notice that there is no dependence on the parameter vector $\xi \in \mathrm{\Gamma }$ in the mass matrix M.

2.2.2. Frequency Response Analysis

In frequency response analysis, the idea is to study the effects of excitation or applied force in the frequency domain. It is clear that the uncertainty in the eigenproblem is strongly connected to the frequency response.

The starting point is the equation of motion for the system [28]:

$$\mathbf{M}\ddot{\mathbf{v}}\left(\xi \right)+\mathbf{C}\dot{\mathbf{v}}\left(\xi \right)+\mathbf{S}\left(\xi \right)\mathbf{v}\left(\xi \right)=\mathbf{f},$$

(8)

where, given $\xi \in \mathrm{\Gamma }$, $\mathbf{M}$ is the mass matrix, $\mathbf{C}$ the viscous damping matrix, $\mathbf{S}\left(\xi \right)$ the stochastic stiffness matrix, $\mathbf{f}$ the (deterministic) force vector, and $\mathbf{v}\left(\xi \right)$ the displacement vector. The matrices $\mathbf{M}$ and $\mathbf{S}\left(\xi \right)$ are defined as above. The viscous damping matrix $\mathbf{C}$ is a linear combination of $\mathbf{M}$ and $\mathbf{S}\left(\xi \right)$. This choice is often referred to as proportional Rayleigh damping, and is central to the low-rank approximation scheme considered below.

In the case of harmonic excitation, a steady-state solution is sought. The angular frequency is $\omega =2\pi f$, where f is the ordinary frequency. The force and the corresponding response have harmonic function representations as

$$\begin{array}{c}\mathbf{f}=\widehat{\mathbf{f}}\left(\omega \right)e^{\mathrm{i}\omega t},\\ \mathbf{v}\left(\xi \right)=\widehat{\mathbf{v}}(\xi ,\omega )e^{\mathrm{i}\omega t}.\end{array}$$

(9)

Taking the first and second derivatives of Equation (8) and substituting using (9) leads to

$$-{\omega }^2\mathbf{M}\widehat{\mathbf{v}}(\xi ,\omega )e^{\mathrm{i}\omega t}+\mathrm{i}\omega \mathbf{C}\hat{\mathbf{v}}(\xi ,\omega )e^{\mathrm{i}\omega t}+\mathbf{S}\left(\xi \right)\widehat{\mathbf{v}}(\xi ,\omega )e^{\mathrm{i}\omega t}=\widehat{\mathbf{f}}\left(\omega \right)e^{\mathrm{i}\omega t},$$

(10)

finally reducing to a linear system of equations:

$$(-{\omega }^2\mathbf{M}+\mathrm{i}\omega \mathbf{C}+\mathbf{S}\left(\xi \right))\widehat{\mathbf{v}}(\xi ,\omega )=\widehat{\mathbf{f}}\left(\omega \right).$$

(11)

In principle, any quantity of interest can be derived from the solution $\widehat{\mathbf{v}}(\xi ,\omega )$. Often, in engineering practice, maximal deflections at specific locations are used. For shells, the presence of boundary layers makes it difficult to guarantee satisfactory convergence. A mathematically reasonable choice is to take the mechanical energy defined by the natural energy norm as the quantity of interest, even though, in practice, it may be difficult to measure.

2.2.3. Stochastic Collocation

There are many options for the selection of the collocation operator. Although a simpler formulation would also suffice for the numerical experiments below, it is advantageous to aim for efficient computation even in high-dimensional parameter spaces. With this in mind, an anisotropic Smolyak-type collocation operator defined with respect to a finite multi-index set (see, e.g., [29,30,31,32]) is taken as the model.

For the sake of accuracy, it is standard practice to choose the collocation points to be the abscissae of orthogonal polynomials; see [33]. The collocation method is formulated using Legendre polynomials, which are the optimal choice when the input random variables are uniform. In the case of, for example, Gaussian random variables, one should use Hermite polynomials instead, but otherwise the collocation method remains the same.

Let $L_p$ be the univariate Legendre polynomial of degree p. Denote by ${\{{\chi }_k^{\left(p\right)}\}}_{k=0}^p$ the zeros of $L_{p+1}$ and by ${\left\{w_k^{\left(p\right)}\right\}}_{k=0}^p$ the associated Gauss–Legendre quadrature weights. The one-dimensional Lagrange interpolation operators ${\mathcal{I}}_p^{\left(m\right)}$ are defined via

$$\left({\mathcal{I}}_p^{\left(m\right)}v\right)\left({\xi }_m\right)=\sum _{k=0}^{p}v\left({\chi }_k^{\left(p\right)}\right){\ell }_k^{\left(p\right)}\left({\xi }_m\right),$$

(12)

where ${\{{\ell }_k^{\left(p\right)}\}}_{k=0}^p$ are the related Lagrange basis polynomials of degree p.

Now, let $\mathcal{A}\subset {\mathbb{N}}_0^M$ be a finite set of multi-indices. For $\alpha ,\beta \in \mathcal{A}$ write $\alpha \le \beta$ if ${\alpha }_m\le \beta m$ for all m = 1, …, M. Further, let us assume that $\mathcal{A}$ is monotone in the following sense: $\exists \alpha \in \mathcal{A}$ such that $\beta \le \alpha \Rightarrow \beta \in \mathcal{A}.$ The sparse collocation operator is defined as

$$\begin{array}{c}\hfill {\mathcal{I}}_{\mathcal{A}}:=\sum _{\alpha \in \mathcal{A}}\underset{m=1}{\overset{M}{⨂}}\left({\mathcal{I}}_{{\alpha }_m}^{\left(m\right)}-{\mathcal{I}}_{{\alpha }_m-1}^{\left(m\right)}\right),\end{array}$$

(13)

where ${\mathcal{I}}_{-1}^{\left(m\right)}:=0$. The collocation points are of the form

$$\begin{array}{c}\hfill {\chi }_{\gamma }^{\left(\alpha \right)}=\left({\chi }_{{\gamma }_1}^{\left(\alpha 1\right)},\dots ,{\chi }_{{\gamma }_M}^{\left(\alpha M\right)}\right)\in \mathrm{\Gamma }\end{array}$$

(14)

for some $\gamma \in {\mathbb{N}}_0^M$ such that $\gamma \le \alpha \in \mathcal{A}$. Similarly, the tensorised quadrature weights $w_{\gamma }^{\left(\alpha \right)}=w_{\gamma 1}^{\left(\alpha 1\right)}\dots w_{\gamma M}^{\left(\alpha M\right)}$ are defined. Statistics, such as the expected value and variance, for the collocated solution may now be computed by applying the quadrature rule

$${\mathbb{E}}_{\mathrm{\Gamma }}\left[\underset{m=1}{\overset{M}{⨂}}{I}_{\alpha m}^{\left(m\right)}v\right]=\sum _{\gamma \le \alpha }v\left({\chi }_{\gamma }^{\left(\alpha \right)}\right)w_{\gamma }^{\left(\alpha \right)}$$

(15)

on the terms in (13).

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 $\mathcal{A}\subset {\mathbb{N}}_0^M$. For a detailed analysis, refer to [30,32]. In this paper, tensor product grids defined by

$$\begin{array}{c}\hfill \mathcal{A}=\{\alpha \in {\mathbb{N}}_0^M|\alpha m\le p,m=1,\dots ,M\},p\in {\mathbb{N}}_0,\end{array}$$

(16)

are used.

3. Free Vibration of Perforated Shells

Shells of revolution have symmetry-induced clusters of eigenmodes. With the exception of torsion modes, in free vibration, every eigenmode is a member of a cluster of size at least two. This is a simple consequence of the eigenmodes having integer-valued wave numbers in the angular direction. In a non-perforated case, this wave number increases with the rate ∼$1/\sqrt[4]{t}$ as $t\to 0$.

For regular perforation patterns, the asymptotics for the lowest eigenmode are given in terms of a conjecture.

Conjecture 1

(Asymptotics of Perforated Parabolic Shells of Revolution, [7]) Given a $g\times g$ regular perforation pattern, there exists a critical dimensionless thickness $t_c$ at which the angular wave number $k_c\approx g/2$. For thicknesses $t<t_c$, the asymptotics of the lowest eigenmode do not conform to those of the non-perforated case.

Geometrically, this means that there exists a thickness at which, within the subspace containing the lowest eigenmode, at least one pair of eigenmodes is in perfect alignment with the perforation grid so that one mode acts (has maximal amplitudes) on the holes and the other in the areas in between.

In many structures, the perforation patterns may be locally regular—for instance, in shells made out of panels with local perforation patterns. Once the angular symmetry is removed, the eigenmodes are no longer pure, but linear combinations of periodic functions in the angular direction. Removal of axial symmetry affects the local axial amplitudes only. Here, the trommel screens with precisely these properties are considered. In Figure 1, an example of the axial effect (Figure 1a,b) and the saturation with high wave numbers in response (Figure 1d) are illustrated. Currently, even a case of two adjacent but different regular perforation grids is open in the sense that the asymptotics cannot be predicted.

4. Low-Rank Approximation

In this section, the choice of the Rayleigh damping model and its effect on the low-rank approximation are motivated. Let us recall the equation of motion (8) for the system integrated at a given collocation point $\xi$. Set $\mathbf{C}=\zeta \left(\alpha \mathbf{S}\right(\xi )+\beta \mathbf{M})$, where the parameters $\alpha >0$, $\beta >0$, $\alpha +\beta =1$, and $\zeta >0$ are chosen on the basis of experimental results and previous experience, and can vary for the same structural model depending on the external conditions. These facts pose an additional requirement for simulation—in other words, the free variation of $\alpha$ and $\beta$ without having to repeat the reduction procedure. This choice of the damping is often referred to as proportional damping and it follows immediately that it preserves the eigenspace of the original undamped problem. Another option would be to use simple diagonal damping.

Given an eigenspace preserving damping model, the task is to find a suitable subspace V of (small) rank r such that the solution of the reduced system where the subscript r denotes reduced matrix, e.g., ${\mathbf{M}}_r=V^T\mathbf{M}V$,

$${\mathbf{M}}_r\ddot{\mathbf{v}}\left(\xi \right)+{\mathbf{C}}_r\dot{\mathbf{v}}\left(\xi \right)+{\mathbf{S}}_r\left(\xi \right)\mathbf{v}\left(\xi \right)=\mathbf{f},$$

(17)

is sufficiently close to the solution of the full problem.

4.1. Krylov Subspace Construction

A remarkably efficient choice for the subspace is the Krylov subspace [9]. Using standard notation from numerical linear algebra, a Krylov basis of dimension n is denoted as

$$colspan\left(V\right)={\mathcal{K}}_n(-{\mathbf{S}}^{-1}\left(\xi \right)\mathbf{M},{\mathbf{S}}^{-1}\left(\xi \right)\mathbf{f}).$$

(18)

The reduced system in Equation (17) matches the first n moments of the original full-order system [16]. Here, the usefulness of proportional damping becomes apparent; it is not needed in its explicit form, and therefore the construction is valid for any choice of $\alpha$ and $\beta$.

4.2. Deflation Preconditioning

The Krylov subspace has to be constructed for every collocation point $\xi$. In order to simplify the notation, the dependence on $\xi$ explicit is not made explicit. It is necessary to solve systems in terms of bending ${\mathbf{S}}_B$ and membrane ${\mathbf{S}}_M$ (membrane and shear summed together)

$$({\mathbf{S}}_M+t^2{\mathbf{S}}_B)v=b,$$

(19)

for every thickness t. In order to avoid factorisation of $\mathbf{S}$ at every $\xi$, it would be natural to consider iterative methods. Unfortunately, shell problems are singularly perturbed and the parameter-independent preconditioning of such systems is an open problem. However, if one considers a sequence of problems, it is possible to transform the problem as follows: Let ${\mathbf{S}}_M=\mathbf{L}{\mathbf{L}}^T$ (Cholesky decomposition); then,

$$\mathbf{L}(\mathbf{I}+t^2{\mathbf{L}}^{-1}{\mathbf{S}}_B{\mathbf{L}}^{-T}){\mathbf{L}}^{T}v=b,$$

(20)

where the subspace defined by ${\mathbf{L}}^{-1}{\mathbf{S}}_B{\mathbf{L}}^{-T}$ is invariant over all thicknesses t. The inner systems $(\mathbf{I}+t^2{\mathbf{L}}^{-1}{\mathbf{S}}_B{\mathbf{L}}^{-T})\widehat{v}=\widehat{b},$ have their spectra bounded by 1 from below, making it sufficient to collect the largest eigenvectors into subspace W, say. Once the subspace W has been constructed, the deflated conjugate gradient method can be applied [34,35]. Here, collecting Lanczos vectors is not sufficient since the construction of the Krylov subspace leads to multiple right-hand sides. This can be taken still further by constructing the subspace only once, corresponding to the point $\xi =(0,\dots ,0)$, and using this as an approximative preconditioner for every other collocation point, leading to only a small number of extra iteration steps.

Remark 1.

Since the system can be singularly perturbed, i.e.,${\mathbf{S}}_M$is not necessarily invertible, it is possible to use the stabilised version

$$(({\mathbf{S}}_M+ϵ{\mathbf{S}}_B)+(t^2-ϵ){\mathbf{S}}_B)v=b,$$

(21)

where$ϵ\in [0,t^2]$. The choice of optimal$ϵ$is problem-dependent.

Remark 2.

This kind of preconditioner requires multiple sparse triangular solves. In high-level systems such as Matlab and Mathematica, the built-in solver routines typically have superior performance, and the benefits of the approach proposed here are difficult to realise [36].

5. Numerical Experiments: Trommel Screen

By design, the perforation patterns of trommel screens are not uniform, which makes the problem more interesting since the existing asymptotic results on eigenmodes are not directly applicable here since they rely on the material uniformity assumption. In the reference case (see Figure 2), the representative parabolic profile function is $\varphi \left(x\right)=1$. In all cases, $x\in [-\pi ,\pi ]$ so that the 2D computational domain is $D=[-\pi ,\pi ]\times [0,2\pi ]\backslash \left\{\mathrm{perforations}\right\}.$ The hole coverage is very high, 55%, with one half of the shell with a regular $20\times 10$ grid, and another with two $5\times 5$ triangular pattern panels. Even though the mesh appears rather coarse, the p-version of the finite element method is used with uniform $p=4$, which is sufficient to prevent numerical locking effects from dominating the solution.

The quantity of interest is the relative accuracy of the low-rank approximation over a range of thicknesses. As a byproduct of this, one can obtain asymptotic rates for quantities of interest such as the total energy or the frequency of the dominant response as the thickness tends to zero. The size of the Krylov subspace is determined dynamically via orthogonalisation as vectors are added. Using notation of (18), ${\mathcal{K}}_4$ was found to be adequate over the whole set of experiments.

Material constants adopted for all simulations are: $E=2.069\times {10}^{11}$ MPa, $\nu =1/3$, and $\rho =7868$ kg m⁻³, unless otherwise specified. For damping, the selected weights are simply $\alpha =\beta =1/2$, and $\zeta =1/2000$, which is as in earlier work [8]. The deterministic parameter ranges of $t\in [1/100,1/1000]$ (the so-called practical range), and frequency $f\in [5,240]$ Hz (chosen experimentally) with angular frequency $\omega =2\pi f$, are covered. The loads or excitation functions act only on the transverse direction w. The first three are Fourier modes in the angular direction and constant in the axial direction, i.e., $\widehat{f}(x,y)=Ccos\left(Ky\right)$ N, where $C=1000$ is a scaling parameter, and $K\in \{5,10,20\}$ is the angular wave number. The fourth loading is a concentrated load $\widehat{f}(x,y)=C^{2}exp(-100(x^2+{(y-\pi )}^2))$ N. The concentrated load is not a point load, since, formally, the convergence of the finite element method is not established for point loads. The random Young’s modulus depends only on the axial direction and has the form $E(x,y,\xi )=E_0\left(x\right)+{\sum }_{m=1}^M{E}_m\left(x\right){\xi }_m,(x,y)\in D,$$\xi \in {[-1,1]}^M,$ where $E_0\left(x\right)=E$ (constant), and $E_m\left(x\right)=\sqrt{{\lambda }_m}Esin\left(mx\right)$ with decaying sequence of coefficients ${\lambda }_m=1/{(m+1)}^4$.

The angular wave numbers are chosen as factors of the regular grid size $=20$ in the angular direction. In particular, the results on free vibration indicate (Conjecture 1) that the performance profiles for $K=10$ should mark a transition to dominant local (boundary layer) features as the thickness $t\to 0$ [7].

5.1. Effect of Boundary Layers

For parabolic shells, there exist axial layers with a characteristic length scale ∼$\sqrt{t}$ and angular layers ∼$\sqrt[4]{t}$. In Figure 3 and Figure 4, transverse deflection profiles for one Fourier load ($K=20$) and the concentrated load are shown, respectively. In the Fourier load, the large perforations dictate the profile, except in the case of $t=1/1000$ and $f=40$. However, for the concentrated load, the angular wave structure and its dependence on the thickness is clearly visible. In both Figure 3d and Figure 4d for $t=1/1000$ and $f=40$, the local features dominate completely.

It also worth mentioning that even though different types of solution features dominate in different parameter configurations, the energies may not necessarily differ significantly. In both of the examples, the relative change is greatest within the thick ($t=1/100$) case, a fact which is practically impossible to deduce from the transverse displacement field only.

5.2. Low-Rank Approximation Concerns

As demonstrated above, for perforated shells, the interactions between loading, perforation patterns, and deterministic parameters, such as thickness and frequency in this context, are complex. In particular, in those situations where the local features become dominant, it is intuitively clear that the low-rank approximation may lead to unexpected errors. On the other hand, in the global feature range of problems, there is no reason not to employ low-rank approximations. Given the mesh of Figure 2 with uniform $p=6$, one ends up with a system of 481,140 degrees of freedom after the boundary conditions are deployed. Remarkably, the Krylov basis has a rank $=4$ after orthogonalisation for every instance tested here.

In Figure 5, the frequency responses of the cases $K=20$ and the concentrated load at $t=1/1000$ are shown. For the Fourier case, the relative errors in the total energy are fairly large, but the frequencies correct, with the exception of the secondary peak at the low frequencies. In the concentrated load case, the low-rank approximation fails and the dominant frequency is not correct.

Let us further focus on the case of Figure 5b. In Figure 6, the convergence of the low-rank approximation as the number of basis vectors grows has been added in two cases. The basis has a rank $=4$, so the sequence rank $=2,3,4$ is complete. In the case of $K=5$, $t=1/100$, already rank $=3$ is sufficient. However, in Figure 6b, strong variation within the sequence is observed. In particular, in the case of rank $=3$, the observed energy is very small $=0.001$ exactly at the frequency $f=70$, where the full version has its maximum $=0.072$. In this case, the $L^2$-norms of the different vector field components of the solution (see Appendix A.2) reveal that, in the low-rank approximation, the boundary layers in the axial direction have amplitudes that are an order of magnitude smaller than those of the full solution. These data are tabulated in Figure 6c.

An interesting question is whether stochastic modelling of the material parameters could recover or contain the right solution in the case shown in Figure 5d. In order to test this, the collocation with $M=6$ and second-order polynomials (729 collocation points in tensor quadrature) was used with pointwise low-rank approximation. It appears that the approximation error is orders of magnitude larger than anything that could result from simple perturbations of the material parameters. This is not necessarily a general phenomenon; in Figure 7, an example where this compensation actually happens is shown.

It has to be emphasised that the message here is that, indeed, there are extreme configurations where the low-rank approximations may fail. It is reasonable to expect that, in usual design situations, these cases can be identified a priori.

5.3. Asymptotics

As stated above, there are two different interesting asymptotics related to the most energetic response that follow directly from this study. These are the observed frequency and the total energy as functions of the thickness. Two features of the experimental setup are central. First, the energy contributions of the regular grid part and the triangular pattern panels are parameter-dependent; second, for a given loading, there exists a critical thickness $t_c$ at which the solution is aligned with the regular grid so that maximal amplitudes are exactly in between the perforation rows leading to the energetic response. For thicknesses $\ll t_c$, local bending occurs at the scale of the smallest perforations and thus total energy growth saturates.

In Figure 8 and 10, both sets of asymptotics are shown. Let us first concentrate on the observed frequencies. Both $K=5$ and the concentrated load have one dominant parameter-dependent response over the whole range, with respective rates of ∼1 and ∼1/2. However, for $K=10$ and $K=20$, there are two competing ones, with the transition to the higher-frequency one occurring at some critical thickness $t_c$. For $K=5$, this transition would also eventually occur, but with a critical thickness below the range considered here. Moreover, for $K=10$, this transition is sharp; for $K=20$, it is more gradual; in fact, in Figure 8c, around $t=3/1000$, the low-rank approximation picks higher frequencies already before the transition. This is explained with graphs in Figure 9.

The energy asymptotics of Figure 10 tell a similar story. For $K=10$, the sharpness of the transition leads to an almost monotone increase in the observed total energy. In the case of $K=5$, the large holes dominate first, leading to a higher initial rate of energy increase. A similar change in rate is also seen for $K=20$, but here, the reason is saturation in the higher frequencies with local features dominating.

5.4. Other Considerations

For the iterative scheme discussed above, it is the largest thickness which is the problematic one, since automatically also the condition number is the largest one. For the computationally most efficient case with only one factorisation over the whole set, the iteration counts for the Krylov basis for $t=1/100$ are $(33,29,25,23)$, i.e., constant in relative terms. This is also the worst-case result. In the current implementation, this approach is 50-times slower (!) than the built-in solver in Mathematica 13 on MacOS 12.2.

The Krylov basis has to be computed separately for every thickness, and using the fastest option took on average two seconds on a modern workstation. Naturally, solving a resulting $4\times 4$ system is orders of magnitude faster to solve once the reduction has been performed.

The collocation solutions were performed on Aalto University cluster Triton using the low-rank approximation.

6. Conclusions

Low-rank approximation techniques are highly effective in reducing the computational complexity of frequency response analysis. In parameter-dependent problems such as shell problems, where boundary layer-induced local features may actually dominate the response, it is possible that the approximation error is so large that the observed frequency response is misleading. This curious phenomenon should not be alarming since it is likely to occur only in cases where either the shell is very thin or the loading is challenging—for instance, highly concentrated.

The addition of uncertainties to models increases the need for efficient solution techniques, including efficient low-rank approximation. As always, effective implementations depend on the underlying platforms and the current state of the art in this context is still very much evolving.