## Abstract

Eigenvalues of fourth-order elliptic operators feature prominently in stability analysis of elastic structures. This paper considers out-of-plane modes of vibration of a thin elastic plate perforated by a collection of small clamped patches. As the radius of each patch shrinks to zero, a *point constraint eigenvalue problem* is derived in which each patch is replaced by a homogeneous Dirichlet condition at its centre. The limiting problem is consequently not the eigenvalue problem with no patches, but a new type of spectral problem. The discrepancy between the eigenvalues of the patch-free and point constraint problems is calculated. The dependence of the point constraint eigenvalues on the location(s) of clamping is studied numerically using techniques from numerical algebraic geometry. The vibrational frequencies are found to depend very sensitively on the number and centre(s) of the clamped patches. For a range of number of punctures, we find spatial clamping patterns that correspond to local maxima of the base vibrational frequency of the plate.

## 1. Introduction

Determining the spectra of fourth-order elliptic operators is a ubiquitous problem arising in structural analysis of rigid objects. Two central problems (cf. [1]) are the determination of modes of vibration *u* of a thin plate _{k} and the corresponding eigenfunctions *u*_{k}(**x**) satisfying (1.1). The governing equations (1.1) are supplemented with boundary conditions that vary with particular applications. A common assumption is the clamped condition

The perforation structure, denoted *Ω*_{ε}, can be formally defined as a collection of small patches of common radius *ε* and written explicitly as
*Ω*_{k}. The plate itself then occupies the domain *Ω*∖*Ω*_{ε} (cf. figure 1).

The problem considered in this work is the determination of the free modes of vibration (solutions to (1.1a)) in the presence of the perforation structure (1.2). This leads to consideration of the fourth-order partial differential equation (PDE)
*Ω*∖*Ω*_{ε} is the extent of the plate and *Ω*_{ε} is the collection of patches. Clamped boundary conditions (1.3b) are applied to the periphery of the domain and on each perturbing perforation.

The focus of this work is to study the dependence of the eigenvalues λ_{ε} on *Ω*_{ε}, and more specifically on the number *N* and locations *ε*→0.

Before detailing the main results of this study, it is useful to review established key results. A good deal of mathematical interest in the bi-Laplacian eigenvalue problem has focused on the qualitatively different behaviours of (1.3) compared with its well understood second-order elliptic counterpart. A key property of problem (1.3) is that the fundamental eigenfunction, i.e. that associated with the lowest eigenvalue, is not guaranteed to be of one sign [10–15]. This contrasts sharply with classical results for the eigenvalue problem for the Laplacian [16,17]. In the case of the annular region *ε*<*r*<1, detailed properties of Bessel functions [11] determine that for *ε*^{−1}>762.36, the fundamental eigenvalue is not simple and the corresponding eigenspace is two-dimensional. In domains with a corner, the first eigenfunction may possess an infinite number of nodal lines [13]. The numerical treatment of fourth-order eigenvalue problems, with and without perforations, requires careful attention to account for these complexities and has led to many specialized methodologies [18–22].

There are additional remarkable phenomena associated with the bi-Laplacian eigenvalue problem in domains with small subdomains removed. One may expect that as *ε*→0, (*u*_{ε},λ_{ε}) would approach (*u*^{★},λ^{★}), the problem with no patches satisfying
*ε* and centred at **x**_{0} is removed, the limiting behaviour is *u*_{0} is a distinct *point constraint eigenvalue problem*
*ε*→0 through a point constraint *u*_{0}(**x**_{0})=0 on the limiting eigenfunction. Note that the condition on the gradient of the eigenfunction is not retained. For the case of a single patch, the limiting behaviour of the eigenvalues λ_{ε} of (1.3) was determined [23,24] to be
*ε*→0, provided |∇*u*_{0}(**x**_{0})|≠0. Higher-order corrections to the behaviour (1.6) were also calculated in [23,24] for both the case |∇*u*_{0}(**x**_{0})|≠0 and the degenerate scenario |∇*u*_{0}(**x**_{0})|=0. In the exactly solvable case of the annular region *ε*<*r*<1, this singular limiting behaviour is manifested in a plot (cf. figure 2) of λ_{ε} against *ε*→0 as a jump discontinuity at *ε*=0. In the case of a rectangular domain, numerical studies (cf. [25]) of (1.5) have found that a single puncture has a significant localizing effect on the eigenfunctions. Specifically, the presence of a puncture may subdivide *Ω* into distinct regions of high and low energy.

We remark that the limiting behaviour of eigenvalues to (1.3) as *ε*→0 is qualitatively very different from the well-established behaviour of the eigenvalues of the Laplacian (cf. [26–30]) in the presence of small perturbing holes. In such problems, the limiting behaviour is to the patch-free problem. This dichotomy between limiting behaviours between the second- and fourth-order eigenvalue problems has recently been studied in a mixed-order problem (cf. [31]).

In this paper, we consider problem (1.3) in the presence of *N* punctures and in the limit *ε*→0. In §2, we study a simplified model problem in the annulus *ε*<*r*<1 with clamped boundary conditions. This serves to illustrate the discontinuous behaviour of (1.3) as *ε*→0 in an exactly solvable situation. These reduced problems also serve as a guide for the analysis of the general case (1.3) where exact solutions are not available.

In §3, we analyse the limiting behaviour of (1.3) as *ε*→0 and assuming the non-degeneracy condition *ε*→0, where *μ*=*O* (*ν*^{M}) for any integer *M*. The limiting point constraint problem for (*u*_{0},λ_{0}) is
_{k} in (3.13). The one term case is given by
*ε*→0. In §3a, an explicit expression for λ^{★}−λ_{0}, the magnitude of the jump discontinuity (cf. figure 2) in λ_{ε} as *ε*→0 is derived. This deviation depends on the number and locations

In §4, we develop a numerical methodology to solve for the limiting eigenfunctions and eigenvalues and investigate the dependence of the eigenvalues on the configuration of centres *N*. An analytical treatment of this question is challenging owing to a lack of an explicit Green's function with which to express closed form solutions of (1.8). To overcome issues relating to high condition numbers, notorious in discretizations of fourth-order partial differential equations, we employ the adaptive precision solver Bertini (cf. [32,33]). In particular, our numerical methodology allows us to determine particular arrangements of punctures _{0}. For the range *N*=1,…,9, we determine (locally) optimal configurations of patch centres. In §5, we discuss future avenues of investigation which arise from this work.

## 2. A model problem with a point constraint

In this section, we consider the reduced model problem for the deflection of the annular plate *ε*<*r*<1 under a load *f*(**x**). The purpose of this example is to clarify the presence and consequences of a point constraint in the limiting problem in the simplest possible and exactly solvable setting. This also serves to highlight the phenomena that will appear in (1.3) and gives explicit information on the limiting form of solutions. The problem considered is the clamped plate under load
*f*. This problem is simpler than the corresponding eigenvalue problem for annulus as there is no unknown λ to be determined. The eigenvalue problem for the annulus has been considered in [11], and some key results are reviewed later in §4.

### (a) Case *f*=1

In the case *ε*=0 in which there is no perturbing patch, the solution of (2.1) is denoted *u*^{★}, and the radially symmetric solution is calculated to be
*ε*>0 and there is a circular patch of radius *ε* centred at the origin, the general solution of (2.1) satisfying the conditions at *r*=1 is given by
*A*_{ε},*B*_{ε}} is determined by applying the boundary conditions at *r*=*ε* and the limiting behaviour as *ε*→0 is
*u*_{0}(*r*) satisfying
*u*^{★}(*r*) of (2.2) shows that the limiting problem does not converge to the patch-free problem, and, in particular, *u*_{0}(*r*) features the point constraints *u*_{0}(0)=∂_{r}*u*_{0}(0)=0. We also remark that the limiting function *u*_{0}(*r*) is not smooth and satisfies *r*→0. The limiting form of (2.4) reveals that the expansion for *u*_{ε} of (2.1) in the case *f*=1 is of form
_{r}*u*_{0}(0)=0. To investigate the non-degenerate case ∂_{r}*u*_{0}(0)≠0, we consider an asymmetric forcing *f*.

### (b) Case f = cos θ

In the absence of a patch (*ε*=0), the smooth solution of problem (2.1) is
*u*^{★}(0)=0. As before, the annulus *ε*<*r*<1 is considered, and the limit *ε*→0 is investigated. The problem (2.1) for *ε*>0 has a general solution spanned by *r*=1 gives
*A*_{ε},*B*_{ε}} determined by enforcing *u*_{ε}=∂_{r}*u*_{ε}=0 on *r*=*ε*. The limiting behaviour of {*A*_{ε},*B*_{ε}} as *ε*→0 is
*ν*, a so-called *logarithmic series*. The series can be expanded in *ν* to give
*u*_{0}(0)=0 and

These model problems indicate two qualitative solution regimes that we can expect to encounter in the eigenvalue problem (1.3). First, in the case of a single clamping point **x**_{0}, if *u*^{★}(**x**_{0})=0, then the limiting problem *u*_{0} will coincide with the unperturbed problem. This result is encapsulated later in theorem 3.2. The second is that we can expect two different types of expansions, depending wherever the gradient of the function at the clamping point ∇*u*_{0}(**x**_{0}) is zero or not. In the non-degenerate scenario ∇*u*_{0}(**x**_{0})≠0, the expansion is a logarithmic series, whereas if ∇*u*_{0}(**x**_{0})=0, then the expansion is of form (2.6).

## 3. Derivation of the limiting behaviour

In this section, we consider the eigenvalue problem (1.3) and determine the limiting behaviour of a simple (multiplicity one) eigenpair (*u*_{ε},λ_{ε}) of (1.3) in the limit as *ε*→0. As motivated by the simplified example of the previous section, we will assume that the corresponding eigenfunction of the limiting point constraint problem satisfies a certain non-degeneracy condition.

To construct solutions of (1.3) in the limit as *ε*→0, we first note that *ε*→0, which implies that the patch centred at **x**_{k} must be replaced by an equivalent local condition on the limiting eigenfunction as **x**→**x**_{k}. To establish this condition at the *k*th patch, we consider the rescaled problem in the stretched variables
*v*_{k}(**y**), satisfying
*Ω*_{k} and can be identified explicitly in a few simple cases. When *Ω*_{k} is the unit disc, the exact solution of (3.2) is
*I* is the 2×2 identity matrix. In the example of an ellipse with semi-major axis *a*, semi-minor axis *b* for *a*>*b* and where the semi-minor axis is inclined at an angle *α* to the horizontal coordinate, the matrix entries of *Ω*_{k},

The logarithmic terms in the far-field behaviour (3.2b) of the canonical function *v*_{k} indicate that the inner solution *v* near **x**=**x**_{k} should have an infinite series expansion in powers of **a**_{m,k} are vector-valued constants whose value will be determined. Upon using the far-field behaviour (3.2b) in this infinite sum, and writing the resulting expression in the variable (**x**−**x**_{k}) through (3.1), we obtain the matching condition
**x**→**x**_{k}. The form of (3.3) suggests that the limiting eigenfunction for (1.3) should be expanded as
*M*. Note that the form of this asymptotic expansion corresponds to the model example of §2b in the non-degenerate case. Upon substituting this expansion into (1.3), we equate terms at each order of *ν* to obtain a sequence of problems. At leading order, the pair (λ_{0},*u*_{0}) satisfies
**x**=**x**_{k} then provides the conditions that
*u*_{0}(**x**) can be represented as
*G*(**x**;** ξ**,λ) is the bi-Laplacian Helmholtz Green's function satisfying

*R*is a regular function at

**x**=

**. The constants**

*ξ***={**

*α**α*

_{1},…,

*α*

_{N}} represent the strengths of singularity associated with each point constraint. The (

*N*+1) unknowns (

**,λ**

*α*_{0}) are determined by the (

*N*+1) conditions

_{m},

*u*

_{m}) to this leading-order problem. The equation satisfied by

*u*

_{m}at order

*k*=1,…,

*N*. The sequence of problems (3.8) recursively determines the unknown vectors

**a**

_{m,k}for

*m*≥1 and

*k*=1,…,

*N*. In the case

*m*=1, for example, the strength of singularity is prescribed as

**a**

_{0,k}=∇

*u*

_{0}|

_{x=xk}which uniquely determines

*u*

_{1}. Thereafter, matching the regular part of

*u*

_{1}as

**x**→

**x**

_{k}to the behaviour (3.8b) provides a condition for

**a**

_{1,k}. The recursive process continues, where matching the regular part at

_{m}for

*m*≥1, we first need to establish the following lemma:

### Lemma 3.1

*Let* (*u*_{0},λ_{0}) *be an eigenpair of* (3.4) *with multiplicity one. Then, a necessary condition for the problem*
*and*
*to have a solution is that* λ_{m} *satisfies*
*Here*, **a**_{0,k}=∇_{x}*u*_{0}|_{x=xk}, *and we have defined the inner product*

### Proof.

To begin the proof, we consider the integral identity
*u*_{0} and integrating over *Ω*∖*B*_{σ}, where *σ*→0
**x**=**x**_{k}, where *r*=|**x**−**x**_{k}| and _{n}=−∂_{r} on |**x**−**x**_{k}|=*σ*. In terms of these quantities, we use the local behaviour (3.9b) for *u*_{m} and from (3.5) *r*→0 to calculate for *r*→0 that
*σ*→0:

By applying this lemma to (3.8), the coefficients λ_{m} in the expansion for the eigenvalue λ_{ε} are obtained as follows:

### Principal Result 2

Let (*u*_{0},λ_{0}) be an eigenpair of (3.4) with multiplicity one, and assume _{ε} of the perturbed problem (1.3) admits the expansion
*M* and where λ_{m} for *m*≥1 are defined by
**a**_{m−1,k} for *m*≥2 and *u*_{m} for *m*≥1 are determined from (3.8).

The key feature of the analysis is that the difference λ_{ε}−λ_{0} depends on the gradient of the base eigenfunction *u*_{0} at each point of clamping. Crucially, this base eigenfunction *u*_{0} is not determined by the problem in the absence of perturbing patches, but by the following problem with *N* additional point constraints:
_{0}, and in terms of the solution, we calculate **a**_{0,k}≡∇*u*_{0}|_{x=xk} for *k*=1,…,*N* from the regular part *R*_{k}(**x**) of the limiting eigenfunction *u*_{0}. In the case where a clamping location **x**_{k} is a point at which *u*_{0} has zero gradient, i.e. ∇*u*_{0}(**x**_{k})=0, then that location does not contribute to λ_{1} in (3.12).

### (a) Calculation of jump in eigenvalue

As the limiting problem does not correspond to a patch-free eigenvalue problem for (*u*^{★},λ^{★}) satisfying
_{0}−λ^{★}, we establish the following result:

### Theorem 3.2

*Consider a simple eigenpair (u*_{0}*,λ*_{0}*) of the limiting problem (3.14) with local behaviour (3.14b) as* **x***→***x**_{k} *for k=1,…,N and a simple eigenpair (u*^{★}*,λ*^{★}*) of the patch-free problem (3.15). In the case where 〈u*_{0}*,u*^{★}*〉≠0, then
*

### Proof.

To incorporate the singularity structure (3.14b) into the calculation of the jump in the eigenvalue, we integrate over the region *B*(**x**_{k},*σ*) is a ball of radius *σ* centred at **x**_{k}∈*Ω*. Repeated integration by parts yields the following identity
*r*=|**x**−**x**_{k}| and _{n} of ∂*B*(**x**_{k},*σ*) satisfies ∂_{n}=−∂_{r}. From direct calculation, the local behaviour (3.14b) yields the expressions
**c**_{k}=∇*R*_{k}|_{x=xk}. After substituting this local behaviour into (3.17), together with ∂_{n}=−∂_{r}, the limit as *σ*→0 is taken. Several terms tend to zero in this limit leaving the final result
*u*_{0},*u*^{★}〉≠0, the result (3.16) is recovered. ▪

The formula (3.18) can be interpreted as a generalization of the standard eigenfunction orthogonality relationship 〈*u*_{i},*u*_{j}〉=*δ*_{ij} to the point constraint eigenvalue problem (3.4). This result also confirms the intuitive expectation that owing to the common boundary conditions *u*_{0}|_{∂Ω}=*u*^{★}|_{∂Ω}=0, the contribution owing to a particular puncture **x**_{k} vanishes as **x**_{k}→∂*Ω*. Moreover, the jump vanishes whenever the configuration *u*^{★}(**x**).

A casual inspection of (3.18) might suggest that the puncture locations _{0}−λ^{★}) correspond to extrema of the unperturbed eigenfunction *u*^{★}. However, the strengths of each singularity *α*_{k}=*α*_{k}(**x**_{1},…,**x**_{N}) also carry a dependence on the puncture locations which is not known explicitly.

The lack of an explicit solution to the limiting problem (3.14) is a significant hindrance in understanding the dependence of λ_{0} on the perforation pattern _{0} from (3.7).

## 4. Numerical investigation of the limiting problem

In this section, we numerically study the dependence of the principal eigenvalue λ_{0} of the point constraint eigenvalue problem (3.14) on the number and location(s) of the patch centres

### (a) Description of the numerical algorithm

The first step in our method exploits the linearity of the problem to separate the regular and singular portion of *u*_{0} with the decomposition *u*_{0}=*u*_{S}+*u*_{R}, where
*u*_{R} satisfies
*N*+1) unknowns (*α*_{1},…,*α*_{N},λ_{0}) are found by enforcing the point constraints *u*_{0}(**x**_{k})=0 for *k*=1,…,*N* together with the normalization condition *α*={*α*_{1},…,*α*_{N}} and eigenvalue λ_{0}(**x**_{1},…,**x**_{N}). The normalization condition adopted in (4.3) is for numerical convenience only. Once a solution is obtained, it can be rescaled to satisfy 〈*u*_{0},*u*_{0}〉=1, or any other normalization condition. Our numerical experiments are performed on the unit disc domain *r*,*θ*) with a uniform polar mesh consisting of *N*_{r}×*N*_{θ} points.

To track the first eigenvalue as the isolated points *N*_{r},*N*_{θ}) increases. To maintain convergence of the numerical scheme and accurately track the paths, we have used adaptive higher precision arithmetic [39,40]. In table 1, we demonstrate second-order convergence of our numerical scheme on an exactly solvable test case, once adaptive precision is used.

At this point, we briefly review closed form solutions of the patch-free (1.4) and punctured (3.14) problems. The factorization Δ^{2}−*μ*^{4}=(Δ−*μ*^{2})(Δ+*μ*^{2})=0 indicates that for the unit disc scenario, the eigenfunctions take the form
*m*=0,±1,±2,… and *J*_{m}(*z*),*Y* _{m}(*z*),*I*_{m}(*z*) and *K*_{m}(*z*) are the standard Bessel functions. In the patch-free problem (1.4), the smooth eigenfunctions satisfying *u*^{★}=∂_{r}*u*^{★}=0 on *r*=1 are

In the punctured case, we can develop a closed form solution of (3.14) for *N*=1 and **x**_{0}=(0,0). The radially symmetric eigenfunctions (*m*=0) of (3.14) are spanned by (4.4) and those satisfying *u*_{0}(0)=0 and *u*_{0}(1)=∂_{r}*u*_{0}(1)=0 are
*A* is a constant set by the normalization condition 〈*u*_{0},*u*_{0}〉=1. The smallest positive root of (4.7b) is

### (b) Results

Numerical determination of the configurations _{0}(**x**_{1},…,**x**_{N}) is challenging on account of several factors. Principally, for any given configuration, λ_{0} is determined by solving a nonlinear system of equations, each evaluation of which requires solving the fourth-order PDE (4.2). Numerical optimization of this objective function is therefore computationally intensive, particularly as the number of unknowns increases. As an initial step in analysing the optimizing locations *r*, the centres are given explicitly by
*r* over which to optimize.

In the case *N*=1, we examine the dependence of the first two eigenvalues on *r*. In figure 4*a*, the first two eigenvalues are seen to converge to *r*→1, in agreement with (3.18). The eigenfunctions associated with the lowest eigenvalue makes a connection between the second (mode 1) eigenfunction *r*→1. In an opposite fashion, the eigenfunctions associated with the second eigenvalue are continuously deformed from the point constraint problem (4.7) to the second (mode 1) eigenfunction

In figure 5*a*, the principal eigenvalue λ_{0}(*r*) is plotted against *r* for the case *N*=2, and a unique global maximum is observed at (*r*_{c},λ_{c}). The open circles represent the eigenvalues of the unperturbed problem (1.4), and in agreement with (3.18), the curve λ_{0}(*r*) tends to this value as *r*→1. As *r*→0, we see that λ_{0}(*r*)→λ_{0,0,} where λ_{0,0} is the exact value obtained in (4.8).

Table 3 displays the numerically computed maximizing ring radius along with the associated maximum principal eigenvalue.

For a single ring with *N* punctures, we observe

In table 3, we observe that a ring with a puncture at the centre generates a larger maximum eigenvalue over a single ring case for *N*≥5. As with the single ring case (4.9), the maximum eigenvalue saturates as *N* increases with limiting behaviour (*r*_{c},λ_{c})∼(0.50,3838.4) as _{c}. The family of two ring puncture patterns is given by
*r*_{1}>*r*_{2} and *N*_{1}≥*N*_{2}. The system for λ_{0} is solved over all 0<*r*_{2}<*r*_{1}<1 and the maximum value recorded (cf. figure 6). In table 5, the maximum first eigenvalue is displayed for a range of two ring configurations.

For *N*=1,…,9, we display in figure 7 the configurations of patch centres which generate the local maximum first eigenvalue over the patterning schemes with one ring (4.9), one ring with a centre puncture (4.10) and two rings (4.11).

## 5. Conclusion

In this paper, we have studied the problem for the modes of vibration of a thin elastic plate with a collection of clamped patches. As the patch radius shrinks to zero, the limiting problem is *not* the patch-free problem, but a point constraint eigenvalue problem. The contribution of this work is threefold. First, we have obtained detailed information on this limiting behaviour in the form of an asymptotic expansion in the limit of small patch radius. Second, exact expressions for the jump between the point constraint eigenvalues (1.8) and patch-free eigenvalues (1.4) were obtained in (3.18). The third contribution of this work is a numerical study on the dependence of the vibrational modes with respect to the number and location of clamping points.

The numerical study of §4 indicates that the vibrational characteristics of the plate depend sensitively on the clamping configuration *N*=1,…,9, we have numerically calculated (locally) optimal configurations for disc geometry.

There are many avenues of future investigation which arise from this study. Prominent among them is the question of determining the patch centres *N* and on non-radially symmetric domains. Analytical progress towards this goal can be made by establishing detailed knowledge of the bi-Laplacian Helmoltz Green's function
*Ω* is a disc and **x**_{0}=(0,0), we know of no closed form exact solutions to (3.6). Such a solution would be very useful in determining optimizing configurations of clamping points

In scenarios where a numerical solution of the full PDE (4.2) is required for each solution of the eigenvalue problem (4.3), it is highly desirable to have an accurate and efficient solution technique which is adaptable to a variety of geometries. An integral equation approach seems well suited to accomplish this goal (cf. [22]).

This paper has largely focused on the dependence of the fundamental eigenfunction and eigenvalue on the perforation configuration. We also expect a significant effect on localization in higher modes, as seen in [25] for the case of a rectangular domain and *N*=1. An interesting problem would study this phenomenon and determine a theory for the placement of punctures to achieve a desired localization pattern.

## Authors' contributions

A.E.L. conceived the problem, performed analysis and wrote manuscript. W.H. and A.J.S. performed numerical simulations with Bertini software package [32,33].

## Competing interests

We have no competing interests.

## Funding

A.E.L. acknowledges support from the National Science Foundation under grant no. DMS-1516753. W.H. has been supported by the Mathematical Biosciences Institute and the National Science Foundation under grant no. DMS-0931642. A.J.S. acknowledges support from the National Science Foundation under grant no. ACI-1440607.

## Acknowledgements

We acknowledge the assistance of the Notre Dame Center for Research Computing (CRC).

- Received July 10, 2015.
- Accepted November 20, 2015.

- © 2015 The Author(s)