Improved non-dimensional dynamic influence function method based on tow-domain method for vibration analysis of membranes

Abstract

This article introduces an improved non-dimensional dynamic influence function method using a sub-domain method for efficiently extracting the eigenvalues and mode shapes of concave membranes with arbitrary shapes. The non-dimensional dynamic influence function method (non-dimensional dynamic influence function method), which was developed by the authors in 1999, gives highly accurate eigenvalues for membranes, plates, and acoustic cavities, compared with the finite element method. However, it needs the inefficient procedure of calculating the singularity of a system matrix in the frequency range of interest for extracting eigenvalues and mode shapes. To overcome the inefficient procedure, this article proposes a practical approach to make the system matrix equation of the concave membrane of interest into a form of algebraic eigenvalue problem. It is shown by several case studies that the proposed method has a good convergence characteristics and yields very accurate eigenvalues, compared with an exact method and finite element method (ANSYS).

Introduction

The authors developed the non-dimensional dynamic influence function (NDIF) method for free vibration analysis of arbitrarily shaped membranes.^1,2 Although the NDIF method has the feature that it yields highly accurate solutions compared with the finite element method (FEM) ³ and the boundary element method,^4,5 the final system matrix equation of the NDIF method does not have a form of the algebraic eigenvalue problem ⁶ unlike the FEM and the boundary element method. As a result, the NDIF method needs the inefficient procedure of searching values of the frequency parameter that make the system matrix singular by sweeping the frequency parameter in the range of interest.

The authors introduced an improved NDIF method to eliminate the above inefficient procedure by changing the final system matrix equation in the NDIF method into a form of the algebraic eigenvalue problem.^7,8 However, the improved NDIF method has the weak point that it does not give good solutions for concave membranes with high concavity and multi-connected membranes with holes. This article employs the sub-domain method of dividing the entire domain of a membrane into two convex sub-domains. Although the proposed method is similar to the authors’ previous method ² in employing the sub-domain approach, the two methods differ in whether the final system matrix equation has a form of the algebraic eigenvalue problem.

A vast literature exists for obtaining analytical solutions of free vibration of membranes and plates having no exact solution as surveyed in the authors’ previous articles.^1,2,7 However, a survey of the literature performed by the authors reveals that only a very few articles of them dealt with concavely shaped membranes and plates.

Recently, many researchers have studied new numerical methods for more accurate solutions of concave membranes and plates than the FEM. For instance, Wu et al. ⁹ developed the local radial basis function–based quadrature method for the vibration analysis of arbitrarily shaped membranes and solved highly concave-shaped membranes without the use of any domain decomposition. Shu et al. ¹⁰ applied the two-dimensional least-square-based finite difference method for solving free vibration problems of arbitrarily shaped plates with simply supported and clamped edges. It should be noted that the above methods^9,10 still have a limitation in accuracy owing to a large amount of numerical calculation. In this article, a simple and practical approach, which is applicable to arbitrary shapes and offers a highly accurate solution, is proposed by extending the authors’ previous research.^2,7 Note that the proposed method is occasionally ill-posed when too many nodes are used. On the other hand, the proposed method does not consider the problem of spurious eigenvalues, which have been studied in many articles related to plate vibrations,^11–16 because the spurious eigenvalues do not appear in the NDIF method for fixed membranes.

NDIF method reviewed

NDIF

The NDIF primarily satisfies the governing equation of the eigen-field of interest and physically describes the displacement response of a point in an infinite domain due to a unit displacement excited at another point.^1,7 In the case of an infinite membrane (see Figure 1), the NDIF between the excitation point $P_k$ and the response point $P$ is given by the Bessel function of the first kind of order zero as follows ¹

$$\mathrm{NDIF}=J_0\left(\Lambda |\mathbf{r}-{\mathbf{r}}_k|\right)$$

(1)

Figure 1.

Infinite membrane with harmonic excitation points that are distributed along the fictitious contour (dotted line) with the same shape as the finite-sized membrane of interest.

in which $\Lambda$ represents a frequency parameter; $\mathbf{r}$ and ${\mathbf{r}}_k$ denote the position vectors for $P_k$ and $P$ , respectively; and $|\mathbf{r}-{\mathbf{r}}_k|$ denotes the distance between $P_k$ and $P$ . The NDIF satisfies the governing equation (Helmholtz equation)

$${\nabla }^{2}W(\mathbf{r})+{\Lambda }^{2}W(\mathbf{r})=0$$

(2)

where $W(\mathbf{r})$ is the transverse displacement of the finite-sized membrane depicted as the dotted line in Figure 1, $\Lambda =\omega /\sqrt{T/\rho }$ denotes the wavenumber in terms of the angular frequency $\omega$ , the uniform tension per unit length $T$ , and the mass per unit area $\rho$ . Detailed illustrations on the NDIF are given in the authors’ previous article. ¹

NDIF method

For free vibration analysis of an arbitrarily shaped membrane, of which the boundary is illustrated by the dotted line in Figure 1, N boundary points are first distributed at points $P_1,P_2,\dots ,P_N$ along the boundary depicted in an infinite membrane. Assuming that harmonic displacements of amplitudes $A_1,A_2,\dots ,A_N$ are, respectively, generated at points $P_1,P_2,\dots ,P_N$ , the total displacement response at the point $P$ may be obtained by the sum of responses (the linear combination of NDIFs given in equation (1)) that have resulted from each boundary point, that is

$$W(\mathbf{r})={\displaystyle \sum _{k=1}^N}{A}_k{J}_0\left(\Lambda |\mathbf{r}-{\mathbf{r}}_k|\right)$$

(3)

where $A_k$ may be called a contribution coefficient. Equation (3) is employed as an approximate solution for the eigen-field of the finite-sized membrane. Note that the approximate solution also satisfies the governing equation because the NDIF satisfies the governing equation.

Next, the boundary condition given for the membrane is discretized at boundary points $P_1,P_2,\dots ,P_N$ as follows

$$W({\mathbf{r}}_i)=U_i,i=1,2,\dots ,N$$

(4)

where $U_i$ denotes a boundary displacement given in boundary point $P_i$ . Then, applying the discrete boundary condition (equation (4)) to the approximate solution (equation (3)) gives

$$W({\mathbf{r}}_i)={\displaystyle \sum _{k=1}^N}{A}_k{J}_0\left(\Lambda |{\mathbf{r}}_i-{\mathbf{r}}_k|\right)=U_i,i=1,2,\dots ,N$$

(5)

Finally, equation (5) may be written in a simple matrix form

$$\mathbf{SM}(\Lambda )\mathbf{A}=\mathbf{U}$$

(6)

where elements of the $N\times N$ system matrix $\mathbf{SM}(\Lambda )$ are given by ${\mathrm{SM}}_{\mathrm{ik}}=J_0(\Lambda |{\mathbf{r}}_i-{\mathbf{r}}_k|)$ , the contribution vector $\mathbf{A}$ represents the contribution strength of the NDIFs defined at each boundary point, and the displacement vector $\mathbf{U}$ represents the discrete boundary condition. In the case of the fixed boundary condition ( $\mathbf{U}=0$ ), equation (6) leads to

$$\mathbf{SM}(\Lambda )\mathbf{A}=0$$

(7)

which is the system matrix equation of the membrane.

Since the system matrix $\mathbf{SM}(\Lambda )$ in equation (7) depends on the frequency parameter $\Lambda$ , the inefficient procedure of calculating the singularity of the system matrix in the frequency range of interest is required to obtain the eigenvalues of the membrane. In the following section, an improved theoretical formulation for the NDIF method is presented to make the system matrix independent of the frequency parameter. The improved formulation is applicable to not only convex membranes but also concave one unlike the authors’ previous method. ⁷

Improved formulation

Assumed solution and boundary conditions

A concave membrane, illustrated as the solid line in Figure 2, is divided into the two domains, $D_I$ and $D_{\mathrm{II}}$ , which have the common boundary ${\Gamma }_c$ . The domain $D_I$ is surrounded by the fixed boundary ${\Gamma }_1$ and the common boundary ${\Gamma }_c$ . The domain $D_{\mathrm{II}}$ is surrounded by the fixed boundary ${\Gamma }_2$ and the common boundary ${\Gamma }_c$ . If the NDIF method is first applied for $D_I$ , an approximate solution in $D_I$ is obtained by the linear combination of NDIFs defined at boundary points on ${\Gamma }_1$ and ${\Gamma }_c$ . As a result, the approximate solution may be in the same manner as equation (3)

$$W_I(\mathbf{r})={\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{J}_0\left(\Lambda |\mathbf{r}-{\mathbf{r}}_k^{(1)}|\right)$$

(8)

Figure 2.

Concavely shaped membrane subdivided into two domains, $D_I$ and $D_{\mathrm{II}}$ .

where $N1$ and $\mathrm{Nc}$ denote the numbers of boundary points on ${\Gamma }_1$ and ${\Gamma }_c$ , respectively; $A_k^{(1)}$ denotes contribution coefficient for domain $D_I$ ; and ${\mathbf{r}}_k^{(1)}$ denotes a position vector for the kth boundary point for $D_I$ . The fixed boundary condition at boundary points $P_1^{(1)},P_2^{(1)},\dots ,P_{N1}^{(1)}$ on ${\Gamma }_1$ may be expressed as

$$W_I({\mathbf{r}}_i^{(1)})=0,i=1,2,\dots ,N1$$

(9)

where ${\mathbf{r}}_i^{(1)}$ denotes the position vector for $P_i^{(1)}$ . Applying the boundary condition, equation (9), to the assumed solution, equation (8), gives

$${\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(I)}{J}_0(\Lambda R_{i,k}^{(1)})=0,i=1,2,\dots ,N1$$

(10)

where $R_{i,k}^{(1)}=|{\mathbf{r}}_i^{(1)}-{\mathbf{r}}_k^{(1)}|$ denotes the distance between boundary points $P_i^{(1)}$ and $P_k^{(1)}$ .

In order to decouple $\Lambda$ from the Bessel function in equation (10), it is expanded in a Taylor series ¹⁷ having the number of terms of the series $M$ as follows

$$J_0(\Lambda R_{i,k}^{(1)})\approx {\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\varphi }_j^{(1)}(i,k)$$

(11)

where function ${\varphi }_j^{(1)}(i,k)$ is

$${\varphi }_j^{(1)}(i,k)=\frac{{(-1)}^j{(R_{i,k}^{(1)}/2)}^{2j}}{{[\Gamma (j+1)]}^2}$$

(12)

Substituting equation (11) into equation (10) yields

$${\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\phi }_j^{(1)}(i,k)=0,i=1,2,\dots ,N1$$

(13)

which is rearranged as follows

$${\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{\phi }_j^{(1)}(i,k)=0,i=1,2,\dots ,N1$$

(14)

In the same manner as equation (8), an approximate solution in $D_{\mathrm{II}}$ is assumed as

$$W_{\mathrm{II}}(\mathbf{r})={\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{J}_0\left(\Lambda |\mathbf{r}-{\mathbf{r}}_k^{(2)}|\right)$$

(15)

where $N2$ denotes the number of boundary points on ${\Gamma }_2$ in Figure 2, $A_k^{(2)}$ denotes contribution coefficient for domain $D_{\mathrm{II}}$ , and ${\mathbf{r}}_k^{(2)}$ denotes a position vector for the kth boundary point for $D_{\mathrm{II}}$ . If the fixed boundary condition at boundary points $P_1^{(2)},P_2^{(2)},\dots ,P_{N2}^{(2)}$ on ${\Gamma }_2$

$$W_{\mathrm{II}}({\mathbf{r}}_i^{(2)})=0,i=1,2,\dots ,N2$$

(16)

is applied to equation (15), one can obtain

$${\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{J}_0(\Lambda R_{i,k}^{(2)})=0,i=1,2,\dots ,N2$$

(17)

where $R_{i,k}^{(2)}=|{\mathbf{r}}_i^{(2)}-{\mathbf{r}}_k^{(2)}|$ denotes the distance between boundary points $P_i^{(2)}$ and $P_k^{(2)}$ . Expanding the Bessel function in equation (17) in a Taylor series and rearranging yield

$${\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{\phi }_j^{(2)}(i,k)=0,i=1,2,\dots ,N2$$

(18)

where function ${\varphi }_j^{(2)}(i,k)$ is

$${\varphi }_j^{(2)}(i,k)=\frac{{(-1)}^j{(R_{i,k}^{(2)}/2)}^{2j}}{{[\Gamma (j+1)]}^2}$$

(19)

Next, continuity conditions in displacement and slope are considered on the common boundary ${\Gamma }_c$ . Displacements at boundary points $P_{N1+1}^{(1)},P_{N1+2}^{(1)},\dots ,P_{N1+\mathrm{Nc}}^{(1)}$ on ${\Gamma }_c$ may be obtained by equation (8) as follows

$$W_I({\mathbf{r}}_{N1+i}^{(1)})={\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{J}_0\left(\Lambda |{\mathbf{r}}_{N1+i}^{(1)}-{\mathbf{r}}_k^{(1)}|\right),i=1,2,\dots ,\mathrm{Nc}$$

(20)

Displacements at boundary points $P_{N2+1}^{(2)},P_{N2+2}^{(2)},\dots ,P_{N2+\mathrm{Nc}}^{(2)}$ on ${\Gamma }_c$ may be obtained by equation (15) as follows

$$\begin{array}{cc}\multicolumn{1}{c}{}& W_{\mathrm{II}}({\mathbf{r}}_{N2+i}^{(2)})={\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{J}_0\left(\Lambda |{\mathbf{r}}_{N2+i}^{(2)}-{\mathbf{r}}_k^{(2)}|\right),\hfill \\ \multicolumn{1}{c}{}& i=1,2,\dots ,\mathrm{Nc}\hfill \end{array}$$

(21)

Since displacements at $P_{N1+1}^{(1)},P_{N1+2}^{(1)},\dots ,P_{N1+\mathrm{Nc}}^{(1)}$ are, respectively, the same as displacements at $P_{N2+1}^{(2)},P_{N2+2}^{(2)},\dots ,P_{N2+\mathrm{Nc}}^{(2)}$ , one can obtain, by using equations (20) and (21)

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{J}_0\left(\Lambda |{\mathbf{r}}_{N1+i}^{(1)}-{\mathbf{r}}_k^{(1)}|\right)\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{J}_0\left(\Lambda |{\mathbf{r}}_{N2+i}^{(2)}-{\mathbf{r}}_k^{(2)}|\right),i=1,2,\dots ,\mathrm{Nc}\hfill \end{array}$$

(22)

Expanding the two Bessel functions in equation (22) in a Taylor series and rearranging yield

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{\phi }_j^{(1)}(N1+i,k)\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{\phi }_j^{(2)}(N2+i,k),i=1,2,\dots ,\mathrm{Nc}\hfill \end{array}$$

(23)

Since slopes at boundary points $P_{N1+1}^{(1)},P_{N1+2}^{(1)},\dots ,P_{N1+\mathrm{Nc}}^{(1)}$ are the same as slopes at boundary points $P_{N2+1}^{(2)},P_{N2+2}^{(2)},\dots ,P_{N2+\mathrm{Nc}}^{(2)}$ , respectively, one can obtain, by differentiating equation (22)

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}\frac{\partial }{\partial n_{N1+i}^{(1)}}{J}_0\left(\Lambda |{\mathbf{r}}_{N1+i}^{(1)}-{\mathbf{r}}_k^{(1)}|\right)\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}\frac{\partial }{\partial n_{N2+i}^{(2)}}{J}_0\left(\Lambda |{\mathbf{r}}_{N2+i}^{(2)}-{\mathbf{r}}_k^{(2)}|\right),i=1,2,\dots ,\mathrm{Nc}\hfill \end{array}$$

(24)

where $n_{N1+i}^{(1)}$ and $n_{N2+i}^{(2)}$ denote the normal directions from points $P_{N1+i}^{(1)}$ and $P_{N2+i}^{(2)}$ on ${\Gamma }_c$ , respectively. Equation (24) is rewritten as

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}\Lambda J_1(\Lambda R_{N1+i,k}^{(1)})\frac{\partial R_{N1+i,k}^{(1)}}{\partial n_{N1+i}^{(1)}}\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}\Lambda J_1(\Lambda R_{N2+i,k}^{(2)})\frac{\partial R_{N2+i,k}^{(2)}}{\partial n_{N2+i}^{(2)}},i=1,2,\dots ,\mathrm{Nc}\hfill \end{array}$$

(25)

The Bessel functions in equation (25) are expanded in a Taylor series as follows

$$J_1(\Lambda R_{\mathrm{Na}+i,k}^{(a)})\approx {\displaystyle \sum _{j=0}^M}\frac{{(-1)}^j{(\Lambda R_{\mathrm{Na}+i,k}^{(a)}/2)}^{1+2j}}{\Gamma (j+1)\Gamma (j+2)}$$

(26)

where $a=1\mathrm{or}2$ . Substituting equation (26) into equation (25) leads to

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\displaystyle \sum _{k=1}^{N1+\mathrm{Nc}}}{A}_k^{(1)}{\psi }_j^{(1)}(N1+i,k)\hfill \\ \multicolumn{1}{c}{}& ={\displaystyle \sum _{j=0}^M}{\Lambda }^{2j}{\displaystyle \sum _{k=1}^{N2+\mathrm{Nc}}}{A}_k^{(2)}{\psi }_j^{(2)}(N2+i,k),i=1,2,\dots ,\mathrm{Nc}\hfill \end{array}$$

(27)

where ${\psi }_j^{(1)}(N1+i,k)$ and ${\psi }_j^{(2)}(N2+i,k)$ are given by

$${\psi }_j^{(a)}(\mathrm{Na}+i,k)=\frac{{(-1)}^j{(R_{\mathrm{Na}+i,k}^{(a)}/2)}^{1+2j}}{\Gamma (j+1)\Gamma (j+2)}\frac{\partial R_{\mathrm{Na}+i,k}^{(a)}}{\partial n_{\mathrm{Na}+i}^{(a)}}$$

(28)

System matrix equation and eigenvalues

Equations (14), (18), (23), and (27) may be expressed in simple matrix equations, respectively

$$({\mathbf{\Phi }}_0^{(1)}+\lambda {\mathbf{\Phi }}_1^{(1)}+{\lambda }^2{\mathbf{\Phi }}_2^{(1)}+\cdots +{\lambda }^M{\mathbf{\Phi }}_M^{(1)}){\mathbf{A}}^{(1)}=0$$

(29)

$$({\mathbf{\Phi }}_0^{(2)}+\lambda {\mathbf{\Phi }}_1^{(2)}+{\lambda }^2{\mathbf{\Phi }}_2^{(2)}+\cdots +{\lambda }^M{\mathbf{\Phi }}_M^{(2)}){\mathbf{A}}^{(2)}=0$$

(30)

$$\begin{array}{cc}\multicolumn{1}{c}{}& ({\hat{}\mathbf{\Phi }}_0^{(1)}+\lambda {\hat{}\mathbf{\Phi }}_1^{(1)}+{\lambda }^2{\hat{}\mathbf{\Phi }}_2^{(1)}+\cdots +{\lambda }^M{\hat{}\mathbf{\Phi }}_M^{(1)}){\mathbf{A}}^{(1)}\hfill \\ \multicolumn{1}{c}{}& =({\hat{}\mathbf{\Phi }}_0^{(2)}+\lambda {\hat{}\mathbf{\Phi }}_1^{(2)}+{\lambda }^2{\hat{}\mathbf{\Phi }}_2^{(2)}+\cdots +{\lambda }^M{\hat{}\mathbf{\Phi }}_M^{(2)}){\mathbf{A}}^{(2)}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\multicolumn{1}{c}{}& ({\mathbf{\Psi }}_0^{(1)}+\lambda {\mathbf{\Psi }}_1^{(1)}+{\lambda }^2{\mathbf{\Psi }}_2^{(1)}+\cdots +{\lambda }^M{\mathbf{\Psi }}_M^{(1)}){\mathbf{A}}^{(1)}\hfill \\ \multicolumn{1}{c}{}& =({\mathbf{\Psi }}_0^{(2)}+\lambda {\mathbf{\Psi }}_1^{(2)}+{\lambda }^2{\mathbf{\Psi }}_2^{(2)}+\cdots +{\lambda }^M{\mathbf{\Psi }}_M^{(2)}){\mathbf{A}}^{(2)}\hfill \end{array}$$

(32)

In equations (29)–(32), the ith row and kth column element of matrices ${\mathbf{\Phi }}_j^{(1)}$ , ${\mathbf{\Phi }}_j^{(2)}$ , ${\hat{}\mathbf{\Phi }}_j^{(1)}$ , ${\hat{}\mathbf{\Phi }}_j^{(2)}$ , ${\mathbf{\Psi }}_j^{(1)}$ , and ${\mathbf{\Psi }}_j^{(2)}$ are given by ${\varphi }_j^{(1)}(i,k)$ , ${\varphi }_j^{(2)}(i,k)$ , ${\varphi }_j^{(1)}(N1+i,k)$ , ${\varphi }_j^{(2)}(N2+i,k)$ , ${\psi }_j^{(1)}(N1+i,k)$ , and ${\psi }_j^{(2)}(N2+i,k)$ , respectively. Equations (29)–(32) may be merged in a single matrix equation as follows

$$({\mathbf{Q}}_0+\lambda {\mathbf{Q}}_1+{\lambda }^2{\mathbf{Q}}_2+\cdots +{\lambda }^M{\mathbf{Q}}_M)\left\{\begin{array}{c}\hfill {\mathbf{A}}^{(1)}\hfill \\ \multicolumn{1}{c}{{\mathbf{A}}^{(2)}}\end{array}\right\}=0$$

(33)

where matrix ${\mathbf{Q}}_j$ is

$${\mathbf{Q}}_j=\left[\begin{array}{c}\hfill {\mathbf{\Phi }}_j^{(1)}\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{{\hat{}\mathbf{\Phi }}_j^{(1)}}\\ \multicolumn{1}{c}{{\mathbf{\Psi }}_j^{(1)}}\end{array}\begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{{\mathbf{\Phi }}_j^{(2)}}\\ \multicolumn{1}{c}{-{\hat{}\mathbf{\Phi }}_j^{(2)}}\\ \multicolumn{1}{c}{-{\mathbf{\Psi }}_j^{(2)}}\end{array}\right]$$

(34)

Note that equation (33) is called the higher order polynomial eigenvalue problem. ¹⁸ Equation (33) may be changed into a linear matrix equation as follows

$${\mathbf{SM}}_L\mathbf{B}=\lambda {\mathbf{SM}}_R\mathbf{B}$$

(35)

where the system matrices ${\mathbf{SM}}_L$ and ${\mathbf{SM}}_R$ are given, using the diagonal matrix $\mathbf{I}$ , by

$${\mathbf{SM}}_L=\left[\begin{array}{ccccc}\hfill 0\hfill & \hfill \mathbf{I}\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill \mathbf{I}\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill \mathbf{I}\hfill \\ \multicolumn{1}{c}{-{\mathbf{Q}}_0}& \hfill -{\mathbf{Q}}_1\hfill & \hfill -{\mathbf{Q}}_2\hfill & \hfill \cdots \hfill & \hfill -{\mathbf{Q}}_{M-1}\hfill \end{array}\right]$$

(36)

$${\mathbf{SM}}_R=\left[\begin{array}{ccccc}\hfill \mathbf{I}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill \mathbf{I}\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill \mathbf{I}\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill {\mathbf{Q}}_M\hfill \end{array}\right]$$

(37)

and the vector $\mathbf{B}$ is given by

$$\mathbf{B}=\left\{\begin{array}{ccccc}\hfill \mathbf{A}\hfill & \hfill \lambda \mathbf{A}\hfill & \hfill {\lambda }^2\mathbf{A}\hfill & \hfill \cdots \hfill & \hfill {\lambda }^{M-1}\mathbf{A}\hfill \end{array}\right\}$$

(38)

where $\mathbf{A}=\{\begin{array}{cc}\hfill {\mathbf{A}}^{(1)}\hfill & \hfill {\mathbf{A}}^{(2)}\hfill \end{array}\}$ .

Finally, multiplying equation (35) by the inverse matrix of ${\mathbf{SM}}_R$ yields

$${\mathbf{SM}}_R^{-1}{\mathbf{SM}}_L\mathbf{B}=\lambda \mathbf{B}$$

(39)

which may be expressed as the algebraic eigenvalue problem

$$\mathbf{SM}\mathbf{B}=\lambda \mathbf{B}$$

(40)

where the system matrix is given by

$$\mathbf{SM}={\mathbf{SM}}_R^{-1}{\mathbf{SM}}_L$$

(41)

Note that the system matrix $\mathbf{SM}$ in equation (40) is independent of the frequency parameter $\lambda$ . As a result, the eigenvalues of the membrane of interest can simply be extracted from equation (40) without the inefficient procedure required in the original NDIF method. ²

Case studies

The validity of the proposed method is shown through numerical tests of circular, rectangular, and highly concave membranes. For each case, the eigenvalues obtained by the proposed method are compared with those computed by the exact solution and FEM (ANSYS).

Circular membrane

The proposed method is first applied to a circular membrane of unit radius where the exact solution ¹⁹ is known. As shown in Figure 3, the membrane is intentionally divided into two domains $D_I$ and $D_{\mathrm{II}}$ for the proposed method. The boundary of each domain is discretized with 12 boundary point and as a result 24 boundary points are used for the membrane.

Figure 3.

Circular membrane divided into two domains using 24 boundary points.

Eigenvalues obtained by the proposed method using $M=12$ , $M=16$ , and $M=20$ are presented in Table 1, which also shows the eigenvalues given by the exact method ¹⁹ and FEM (ANSYS). In Table 1, it may be said that the eigenvalues by the proposed method in the case of $M=20$ converge accurately to those by the exact method. Furthermore, it should be noted in Table 1 that the proposed method using only 24 boundary points yields more accurate eigenvalues than FEM (ANSYS) using 3668 nodes. In addition, the ith natural frequency $f_i$ of the membrane can be calculated from ${\Lambda }_i=2\pi f_i\sqrt{\rho /T}$ where ${\Lambda }_i$ is the ith eigenvalue in Table 1, $\rho$ is the mass per unit area, and $T$ is the uniform tension per unit length.

Table 1.

Eigenvalues of the circular membrane obtained by the proposed method, the exact method, and FEM (parenthesized values denote errors (%) with respect to the exact method).

No.	Proposed method (24 points)			Exact method 19	FEM (ANSYS)
	$M=12$	$M=16$	$M=20$		3668 nodes	2378 nodes	1024 nodes
1	2.405 (0.00)	2.405 (0.00)	2.405 (0.00)	2.405	2.405 (0.00)	2.406 (0.05)	2.417 (0.50)
2	3.832 (0.00)	3.832 (0.00)	3.832 (0.00)	3.832	3.833 (0.04)	3.835 (0.07)	3.851 (0.50)
3	None	5.135 (0.02)	5.136 (0.00)	5.136	5.139 (0.05)	5.141 (0.10)	5.174 (0.74)
4	None	5.520 (0.00)	5.520 (0.00)	5.520	5.528 (0.14)	5.535 (0.27)	5.552 (0.58)
5	None	6.390 (0.16)	6.380 (0.00)	6.380	6.386 (0.10)	6.390 (0.16)	6.461 (1.27)
6	None	None	7.016 (0.00)	7.016	7.029 (0.19)	7.040 (0.35)	7.059 (0.61)

FEM: finite element method.

On the other hand, the accuracy of an eigenvalue obtained by the proposed method can be verified by plotting its mode shape, which is omitted in the article. If the plotted mode shape does not satisfy exactly the given boundary condition (the fixed boundary condition), it may be said that the eigenvalue is not accurate and larger number of nodes and series functions are required to improve its accuracy.

Rectangular membrane

As shown in Figure 4, a rectangular membrane with dimensions $1.2\mathrm{m}\times 0.9\mathrm{m}$ is divided into two domains of the boundaries discretized with 16 points, respectively. In Table 2, eigenvalues obtained by the proposed method for $M=12$ , $M=16$ , and $M=20$ are compared with eigenvalues obtained by the exact method ¹⁹ and FEM (ANSYS). In Table 2, it may be observed that only the first three eigenvalues are obtained for $M=12$ , but the first six eigenvalues are successfully obtained when M increases (for $M=16$ and $M=20$ ) and fully converge to eigenvalues given by the exact method.

Table 2.

Eigenvalues of the rectangular membrane obtained by the proposed method, the exact method, and FEM (parenthesized values denote errors (%) with respect to the exact method).

No.	Proposed method (32 points)			Exact method 19	FEM (ANSYS)
	$M=12$	$M=16$	$M=20$		2806 nodes	1813 nodes	1089 nodes
1	4.363 (0.00)	4.363 (0.00)	4.363 (0.00)	4.363	4.364 (0.03)	4.364 (0.03)	4.365 (0.05)
2	6.293 (0.00)	6.293 (0.00)	6.293 (0.00)	6.293	6.295 (0.03)	6.296 (0.06)	6.301 (0.13)
3	7.458 (0.03)	7.456 (0.00)	7.456 (0.00)	7.456	7.461 (0.07)	7.464 (0.11)	7.467 (0.15)
4	None	8.595 (0.00)	8.595 (0.00)	8.595	8.602 (0.08)	8.606 (0.13)	8.621 (0.30)
5	None	8.727 (0.00)	8.727 (0.00)	8.727	8.732 (0.06)	8.736 (0.10)	8.741 (0.16)
6	None	10.51 (0.00)	10.51 (0.00)	10.51	10.518 (0.07)	10.523 (0.13)	10.54 (0.29)

FEM: finite element method.

Figure 4.

Rectangular membrane divided two domains using 32 boundary points.

Highly concave membrane

Figure 5 shows a $1.2\mathrm{m}\times 0.9\mathrm{m}$ rectangular membrane with a V-shaped notch. The membrane is divided into two domains to apply the proposed method. Each domain is discretized with 16 boundary points as shown in Figure 5. Eigenvalues obtained by the proposed method and FEM (ANSYS) are summarized in Table 3 where it is shown that eigenvalues by the proposed method fully converge for $M=16$ . Since the current membrane has no exact solution, eigenvalues by the proposed method are compared with those by FEM. It may be observed in Table 3 that the five eigenvalues except the first eigenvalue by the proposed method have very small differences within 0.81% with respect to FEM using 1701 nodes. Interestingly, it is noticed that the first eigenvalues by the proposed method and FEM have relatively large difference. The reason why this difference happens will be revealed in a future study.

Table 3.

Eigenvalues of the highly concave membrane obtained by the proposed method and FEM (parenthesized values denote differences (%) with respect to FEM using 1701 nodes).

No.	Proposed method (32 points)			FEM (ANSYS)
	$M=12$	$M=16$	$M=20$	1701 nodes	976 nodes	451 nodes
1	5.89 (3.15)	5.89 (3.15)	5.89 (3.15)	5.71	5.72	5.74
2	6.41 (0.16)	6.41 (0.16)	6.41 (0.16)	6.42	6.43	6.44
3	8.19 (0.24)	8.19 (0.24)	8.19 (0.24)	8.17	8.18	8.21
4	None	8.88 (0.23)	8.88 (0.23)	8.89	8.90	8.92
5	None	9.95 (0.81)	9.95 (0.81)	9.87	9.89	9.94
6	None	11.23 (0.71)	11.23 (0.71)	11.31	11.34	11.43

FEM: finite element method.

Figure 5.

Highly concave membrane divided two domains using 32 boundary points.

Conclusion

A practical, improved NDIF method for the free vibration analysis of concave membranes with arbitrary shapes was proposed in the article. It was revealed that the proposed method yields much more accurate eigenvalues, which converge to exact values for membranes having exact solutions, than FEM using a large number of nodes. On the other hand, the proposed method cannot be directly applied to a multiply connected membrane because it divides the region of the membrane of interest into only two regions. In order to overcome this weak point, an extended way of dividing the region of the membrane into more than two regions will be developed in future research.