Flutter of a Plate at High Supersonic Speeds

Abstract

The vibrations of plate structures placed in a supersonic flow was considered. The undisturbed fluid flow was parallel to the plate. This type of problem is especially important in the aerospace industry, where it is named panel flutter. It has been noticed for a long time that panel flutter may be problematic at high speeds. In this article, two specific problems were treated: in the first one, the plate was in the form of an infinite strip and the flow was in the direction of its finite length. Rigid walls indefinitely extended from the sides of the plate. In the second problem, the plate was a finite rectangle and the flow was parallel to one of its sides. The rest of the plane of the rectangle was again rigid. The first problem was a limiting case of the second problem. The flow was modeled by piston theory, which assumes that the fluid pressure on the plate is proportional to its local slope. This approximation is widely used at high speeds (supersonic speeds in the range of M > 1), and reduces the interaction between the fluid flow and the vibrations of the plate to an additional term in the vibration equation. The resulting problem can be solved by assumed mode methods. In this study, the solution was also found by using the collocation method. The contribution of this study is the correlation between the flutter velocity and the other parameters of the plate. The main result is the flutter velocity of the free fluid flow under which the plate vibrations become unstable. Finally, simple expressions are proposed between the various non-dimensional parameters that allows for the quick estimation of the flutter velocity. These simple expressions were deduced by least squares fits to the computed flutter velocities.

1. Introduction

Fluid–elastic structure interactions are ubiquitous in many engineering disciplines. There are historically famous examples that demonstrate that the interaction between fluid flow and the vibrations of the structure can cause catastrophic failure. It is observed that the character of vibrations typically depend on a characteristic flow velocity that, if it exceeds a certain value (the flutter velocity), will cause instability. The main thrust of fluid–elastic structure interaction research is the determination of the flutter velocity. This study’s innovation lies in the correlation between the flutter velocity and other model properties. The literature defines supersonic speeds as M > 1.2, and thus, the speed range pertains to Mach numbers greater than 1.2.

The type of problem considered in this study is especially important in aerospace structures and is usually given the name panel flutter. Panel flutter has been an active area of research, starting with the development of high-speed projectiles and rockets around the 1950s. A well-known publication that summarizes the theory and the general methods of solution is [1]. A more recent text is [2]. The AGARD (Advisory Group for Aerospace Research and Development—NATO) publication is an engineering manual dealing with the detection and alleviation of aeroelasticity problems, including experimental methods and various collected tables and graphs [3]. The literature on panel and plate flutter is immense; some of the recent work relevant to the problems considered in this study is summarized here: Uzal et al. dealt with incompressible and irrotational flow in a cylindrical channel to analytically find the flutter velocity [4]. Epureanu et al. investigated vibration-based damage with their chances in material and/or stiffness properties of structures. Kapkin et al. investigated the membrane vibrations located at the stagnation point of the flow and they give the flutter velocity for the system [5]. Vedeneev considered panel flutter at low supersonic speeds by using piston theory [6]. E. H. Dowell gave flutter velocity graphs for an infinite plate lying in the same direction as the fluid motion by using the potential theory for different Mach numbers [7]. Durak B. investigated plate vibrations by using potential flow theory in his PhD thesis [8]. Kapkin et al. tried to find a similar equation to predict the flutter velocity for a membrane with a stagnation flow [9]. Uzal et al. give an analytical solution for a plate placed in a rigid channel that fluid flows in [10]. Also, some of the researchers are also trying to avoid this flutter point by applying a force or moving one boundary to control the motion of the plate. Uzal and Korbahti controlled the resonance frequencies of a rectangular plate’s vibrations by applying a discrete force by measuring the displacement of the plate at a point [11]. Sezgin et al. give a boundary backstepping control method to stabilize the flow-induced vibrations for a membrane [12]. Tubaldi et al. investigated a periodically supported flexible plate under flowing fluid axial flow that is in a channel bounded by a rigid wall [13]. Lin et al. investigated how to control the oscillations of a two-dimensional ribbed plate with free boundary conditions by using a point-wise control method [14]. Li used the Rayleigh–Ritz method to determine the modal characteristics of a rectangular plate with elastic boundary supports [15]. The effects of the system parameters on the stability of the plate are discussed.

In this study, we considered the vibrations of infinite strip plates and rectangular finite plates under the effect of a high-speed flow parallel to the plate. The plate vibrations were modeled by the well-known linear (infinite strip) and two-dimensional (finite rectangle) plate vibration equations. The infinite strip was infinite in the direction perpendicular to the flow and the infinite sides were assumed to be simply supported. The finite rectangular plate was simply supported at all its sides. In both cases, the plates were made up of homogeneous, isotropic, linear material. The flow was taken to be in the high supersonic regime, and was modeled by employing piston theory, according to which, the local pressure was proportional to the local slope of the plate. This allowed for including the coupling between the plate vibrations and the fluid flow simply as an additional term in the plate vibration equation. The bare vibration equations were fourth order in space and second order in time; therefore, their solutions were in the form of undamped vibrations, neither growing nor decreasing in time. In contrast, the vibration equations including the piston theory terms contained first-order derivatives in time and space. This allowed for the solutions to decrease or blow-up in time. We were interested in the boundary between these two regimes, which gives the stability boundary, and the flutter velocity, which was the main quantity of interest.

This paper is organized as follows: In Section 2, we present general formulations for the fluid–structure systems under consideration. And a useful approximate theory of aerodynamics, the piston theory approach, is given in Section 3. In Section 4, the governing equations for the infinite strip case are given and non-dimensionalized in a suitable fashion, and the solution using the point collocation method is presented. In Section 5, the case of the finite rectangular plate is also similarly described. Both cases can also be solved analytically; these are presented in Section 6. Thus, analytical and approximate methods of solution (in the present case, the collocation method) could be compared. Section 7 presents the analytical and numerical results (the non-dimensional flutter velocity against a non-dimensional mass ratio) for both problems; also, the convergence of the collocation results to the analytical solution is demonstrated. Finally, inspired by the structure of the flutter curves, simple equations relating the non-dimensional flutter velocity to the mass ratio are developed for both the strip and rectangle cases. The rectangle case also includes another non-dimensional parameter: the ratio of the sides of the rectangle.

2. General Formulation

Many fluid–solid coupled vibration problems have the general form shown in Figure 1.

The fluid flows (in the undisturbed case) parallel to a rigid plane (which is taken as the xy-plane), with an elastic plate occupying the region bounded by a closed curve B. The flow direction is taken as the x-axis. The flow causes vibrations on the plate, which, in turn, slightly modifies (disturbs) the flow. Under certain conditions, more mechanical energy is transferred from the fluid to the plate, causing the magnitude of the vibrations of the plate to increase, usually exponentially. In the actual nonlinear case, the growing vibrations may be bounded by nonlinear interactions (or they may lead to structural failure), but this is still an unwanted situation from an engineering standpoint. Therefore, investigating the linearized fluid–elastic body coupling phenomena is of great engineering and scientific interest. The most general case may be formulated as follows: The flow is modeled by the compressible potential equation:

$${\nabla }^2\varphi =\frac{1}{a^2}\left[\frac{{\partial }^2\varphi }{\partial t^2}+\frac{\partial }{\partial t}\left(\nabla \varphi \dot{\nabla }\varphi \right)+\nabla \varphi ·\left(\nabla \varphi ·\nabla \nabla \varphi \right)\right]$$

(1)

where $\varphi$ is the velocity potential

$$\mathbf{V}=\nabla \varphi$$

(2)

The local speed of sound a is given by

$$a^2={\left(\frac{\partial p}{\partial \rho }\right)}_S$$

(3)

This is because the predominant force that causes the plate to vibrate is the fluid pressure; viscous stresses would normally have little effect since they are mainly in the flow direction. Therefore, unless there are large-scale flow separations, which is not the case here, the viscous effects and other sources of vorticity can be ignored and the flow can be assumed to be potential. The differentiation of the pressure with respect to density is done at constant entropy (isentropic). In the potential equation, the speed of sound is to be eliminated by using the compressible Bernoulli equation:

$$\frac{\partial \varphi }{\partial t}+\frac{1}{2}\nabla \varphi ·\nabla \varphi +\int a^2\frac{dp}{\rho }=constant$$

(4)

The plate vibrations, in the case of an elastic, linear, and isotropic plate material, are governed by

$$\frac{{\partial }^{2}w}{\partial t^2}+\frac{D}{{\rho }_p}{\nabla }^{4}w=0$$

(5)

where $w(x,y,t)$ denotes the plate displacement, D is the plate bending rigidity, and ${\rho }_p$ is the density of the plate ($\mathrm{kg}/{\mathrm{m}}^2$). P is the fluid pressure on the plate surface.

$$p\left(x,y,z\right)=P\left(x,y,z,t\right){\mid }_{z=w\left(x,y,t\right)}$$

(6)

We assumed that the plate boundaries are simply supported by the rigid extensions, which means that the plate displacement and the bending moment are zero at the boundaries. Finally, at the plate surface, the normal components of the plate and fluid velocities must be the same:

$$w\left(x,y,t\right){\mathbf{e}}_{\mathbf{z}}.\mathbf{n}=\mathbf{V}\left(x,y,z,t\right).\mathbf{n}\mathrm{for}z=w\left(x,y,t\right)$$

(7)

Here, ${\mathbf{e}}_{\mathbf{z}}$ is the unit vector in the z-direction, n is the vector normal to the plate surface (while it is vibrating), and V is the fluid velocity. Note that the conditions between the fluid and the solid (pressure expression and the velocity condition) makes the problem nonlinear, in addition to the velocity potential equation. To arrive at a tractable (i.e., linear) but still physical model, we assumed that the plate vibration amplitude was small enough so that the potential equation could be linearized around the uniform flow in the x-direction:

$$\left(1-\frac{U^2}{a_{\infty }^2}\right)\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}=0$$

(8)

where $a_{\infty }$ is the sound velocity in free stream conditions. As a natural extension of this assumption, the pressure on the plate was assumed to be the pressure at $z=0$:

$$p\left(x,y,t\right)=P\left(x,y,z,t\right){\mid }_{z=0}$$

(9)

and the velocity condition on the plate and fluid was also linearized:

$$\frac{\partial \varphi }{\partial y}{|}_{y=0}=\left\{\begin{array}{cc}\frac{\partial w}{\partial t}+U\frac{\partial w}{\partial x}\hfill & \mathrm{on}\mathrm{the}\mathrm{plate}\hfill \\ 0\hfill & \mathrm{off}\mathrm{the}\mathrm{plate}\hfill \end{array}\right.$$

(10)

The problem expressed by these equations can be attacked by first expressing the velocity potential in terms of the plate displacement by means of a two-dimensional Fourier transform in x and y, computing the pressure from Bernoulli’s equation, and then substituting its value at z = 0 into the plate vibration equation. Since the expression for the velocity potential involves double integrals, including the plate displacement, the final form of the plate vibration equation is a partial integro-differential eigenvalue problem.

3. Piston Theory Approximation

The integro-differential equation mentioned in the last section is more or less the most general formulation of the plate flutter problem. But it is quite a difficult mathematical problem. This type of plate vibration problem is usually termed as panel flutter in the aerospace industry, and it is well known that panel flutter is important at high supersonic speeds [16]. A useful approximate aerodynamic theory exists in this regime known as piston theory [17]. Developed in 1956, piston theory asserts that at high supersonic speeds, the local pressure is proportional to the slope of the boundary ((13) below). With the piston theory approximation, the potential flow equation is discarded and the vibrations are governed by a modified form of the plate vibration equation ((14) below). Note that the usual plate vibration equation contains only second-order derivatives; therefore, it can only exhibit steady time-harmonic vibrations. The modified equation above also has first-order derivatives, and it can, in principle, model attenuated or amplified vibrations; the latter case corresponds to panel flutter.

4. Infinite Strip Plate

Figure 2 shows a schematic of the first problem. The flow occupies the region $z>0$ and is in the x-direction with uniform velocity U. The region in the $xy$-plane with $-b/2<x<b/2$ and $-\infty <y<\infty$ is an elastic plate (infinite strip) of width b; the rest of the $xy$-plane is rigid. The vibrations of the plate cause small perturbations in the flow velocity. Some different equations can be found in the literature to define the vibrations of the plates [18,19,20,21,22,23,24]. The vibrations of the plate are governed by [18,19,20]

$$D\frac{{\partial }^{4}w}{\partial x^4}+{\rho }_p{h}_p\frac{{\partial }^{2}w}{\partial t^2}{+p|}_{z=0}=0$$

(11)

where $w=w(x,t)$ is the displacement, $h_p$ is the thickness, ${\rho }_p$ is the density of the plate, and

$$D=\frac{E{h}_p^3}{12(1-{\nu }^2)}$$

(12)

is the flexural rigidity, where E is Young’s modulus and $\nu$ is Poisson’s ratio. The last term in (11) denotes the fluid pressure on the plate. Although the plate is moving, the fluid pressure can be assumed to have its value at $z=0$ within the linear theory. The fluid flow is assumed to be inviscid since viscous effects are negligible due to the lack of flow separation. In general, the linearized form of the compressible potential equation can be used, but here a simpler approximation called piston theory is adopted. Piston theory is widely used and basically states that the local pressure is proportional to the local slope of the plate [16,25,26].

$${p|}_{z=0}=\frac{{\rho }_{f}U}{M}\left(\frac{\partial w}{\partial t}+U\frac{\partial w}{\partial x}\right)$$

(13)

Thus, the coupled fluid–plate vibrations obey

$$D\frac{{\partial }^{4}w}{\partial x^4}+{\rho }_p{h}_p\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\rho }_{f}U}{M}\left(\frac{\partial w}{\partial t}+U\frac{\partial w}{\partial x}\right)=0$$

(14)

where ${\rho }_f$ is the density of the fluid, U is the velocity of the fluid, and M is the Mach number $M=\frac{U}{c}$. The boundary conditions on the plate were taken as

$$\begin{array}{c}\hfill w(-b/2)=w(b/2)=0\\ \hfill \frac{{\partial }^{2}w}{\partial x^2}(-b/2)=\frac{{\partial }^{2}w}{\partial x^2}(b/2)=0\end{array}$$

(15)

which state that the ends of the plate are simply supported. The problem was non-dimensionalized as follows:

$$\begin{array}{c}\hfill x^{\ast }=\frac{x}{b}, w^{\ast }=\frac{w}{b}, z^{\ast }=\frac{z}{b}, t^{\ast }=\frac{t}{b^2\sqrt{\frac{{\rho }_p{h}_p}{D}}}, p^{\ast }=\frac{p}{\frac{D}{b^3}}, u=\frac{U}{\frac{1}{b}\sqrt{\frac{D}{{\rho }_p{h}_p}}}, c_s^{\ast }=\frac{u}{M}\end{array}$$

(16)

Starred quantities are non-dimensional. Substituting in (14) and (15) and simplifying, the result is, after getting rid of the stars since dimensional quantities will not be needed,

$$\frac{{\partial }^{4}w}{\partial x^4}+\frac{{\partial }^{2}w}{\partial t^2}+\mu \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}\right)=0$$

(17)

and the boundary conditions are

$$\begin{array}{c}\hfill w(-1/2)=w(1/2)=\frac{{\partial }^{2}w}{\partial x^2}(-1/2)=\frac{{\partial }^{2}w}{\partial x^2}(1/2)=0\end{array}$$

(18)

where

$$\begin{array}{c}\hfill \mu =\frac{{\rho }_{f}b}{{\rho }_p{h}_p}{c}_s\end{array}$$

(19)

is a non-dimensional number, usually called the mass ratio, and

$$\begin{array}{c}\hfill c_s=\frac{c}{\frac{1}{b}\sqrt{\frac{D}{{\rho }_p{h}_p}}}\end{array}$$

(20)

is the dimensionless sound velocity. To investigate the stability, the plate displacement was assumed to be

$$\begin{array}{c}\hfill w(x,t)=v\left(x\right)e^{i\omega t}\end{array}$$

(21)

Then, (17) and (18) become

$$\frac{{\partial }^{4}v}{\partial x^4}+\mu u\frac{\partial v}{\partial x}+K_{1}v=0$$

(22)

$$\begin{array}{c}\hfill v(-1/2)=v(1/2)=\frac{{\partial }^{2}v}{\partial x^2}(-1/2)=\frac{{\partial }^{2}v}{\partial x^2}(1/2)=0\end{array}$$

(23)

where

$$\begin{array}{c}\hfill K_1=-{\omega }^2+\mathrm{i}\omega \mu \end{array}$$

(24)

$$\begin{array}{c}\hfill u=\frac{U}{\frac{1}{b}\sqrt{\frac{D}{{\rho }_p{h}_p}}}\end{array}$$

(25)

Equations (22) and (23) give an eigenvalue problem for $\omega$. The vibrations of the plate do not grow as long as the imaginary part of $\omega$ is positive; the stability boundary is Im$\left(\omega \right)=0$. Since $\omega$ depends on u (as well as $\mu$), the condition Im$\left(\omega \right)=0$ gives the non-dimensional flutter velocity. The problem defined by (22) and (23) were solved using the collocation method. For this purpose, the approximate solution was expressed as

$$\begin{array}{c}\hfill v\left(x\right)=\sum _1^N{C}_n{\varphi }_n\left(x\right)\end{array}$$

(26)

where ${\varphi }_n$ are suitably chosen base functions and $C_n$ are constants to be determined. Base functions were chosen so as to satisfy the boundary conditions (22). Here, a family of polynomials was chosen:

$$\begin{array}{c}\hfill {\varphi }_n\left(x\right)=x^{n+3}-\frac{2n+1}{4n-2}{x}^{n+1}+\frac{2n+3}{16(2n-1)}{x}^{n-1}\end{array}$$

(27)

${\varphi }_n\left(x\right)$ satisfies all the boundary conditions (23). Substituting the approximate solution (26) into the governing (22) will not satisfy it, but will result in a “residual”:

$$\begin{array}{c}\hfill R=R(C_1,C_2,\dots ,C_N,x)=\sum _{n=1}^N{Q}_n(x,\omega ,\mu ,u)C_n\end{array}$$

(28)

where, for brevity, we define

$$\begin{array}{cc}\hfill Q_n(x,\omega ,\mu ,u)=& \sum _{n=1}^N{C}_n(n+3)(n+2)(n+1)n{x}^{n-1}-\sum _{n=1}^N{C}_n(n+1)(n-2)(n-1)n{A}_n{x}^{n-1}\hfill \\ \hfill & +\sum _{n=1}^N(n-4)(n-3)(n-2)(n-1)B_n{x}^{n-5}+K_1\sum _{n=1}^N{C}_n(x^{n+3}-A_n{x}^{n+1}+B_n{x}^{n-1})\hfill \\ \hfill & +\mu u\sum _{n=1}^N{C}_n((n+3)x^{n+2}-(n+1)A_n{x}^n+(n-1)+B_n{x}^{n-2})\hfill \end{array}$$

(29)

where

$$\begin{array}{c}\hfill A_n=\frac{2n+1}{4n-2}\end{array}$$

(30)

$$\begin{array}{c}\hfill B_n=\frac{2n+3}{16(2n-1)}\end{array}$$

(31)

In the collocation method, the free parameters $C_n$ are determined by equating the residual R to zero at N collocation points $x_1,x_2,x_3,\dots x_N$, which gives a linear homogeneous system of algebraic equations:

$$\begin{array}{c}\hfill \sum _{n=1}^N{Q}_n(x,\omega ,\mu ,u)=0\end{array}$$

(32)

For a non-trivial solution, the determinant of the coefficients should be zero:

$$\begin{array}{c}\hfill \left[\begin{array}{cccc}{Q}_1(x_1,\omega ,\mu ,u)& Q_1(x_2,\omega ,\mu ,u)& \cdots & Q_N(x_N,\omega ,\mu ,u)\\ Q_2(x_1,\omega ,\mu ,u)& Q_2(x_2,\omega ,\mu ,u)& \cdots & Q_N(x_N,\omega ,\mu ,u)\\ ⋮& ⋮& \ddots & ⋮\\ Q_2(x_1,\omega ,\mu ,u)& Q_2(x_2,\omega ,\mu ,u)& \dots & Q_N(x_N,\omega ,\mu ,u)\end{array}\right]=0\end{array}$$

(33)

The eigenvalue $\omega$ is determined from this equation in the form

$$\begin{array}{c}\hfill \omega =\omega (\mu ,u)\end{array}$$

(34)

and the flutter velocity is found as a function of the mass ratio $\mu$ from

$$\begin{array}{c}\hfill Im\left[\omega \right(\mu ,u\left)\right]=0.\end{array}$$

(35)

5. Rectangular Plate

Figure 3 shows a schematic of the second problem considered. Now, the plate also has a finite width 2d; the $d=\infty$ limit of this problem gives the first problem. Again using piston theory, the vibration equation takes the form

$$D\left(\frac{{\partial }^{4}w}{\partial x^4}+2\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+\frac{{\partial }^{4}w}{\partial y^4}\right)+{\rho }_p{h}_p\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\rho }_{f}U}{M}\left(\frac{\partial w}{\partial t}+U\frac{\partial w}{\partial x}\right)=0$$

(36)

The plate was assumed to be simply supported at all sides, and thus, the boundary conditions are

$$\begin{array}{c}\hfill w(-b/2,y,t)=w(b/2,y,t)=w(x,-d,t)=w(x,d,t)=0,\end{array}$$

(37)

$$\begin{array}{c}\hfill \frac{{\partial }^{2}w}{\partial x^2}(-b/2,y,t)=\frac{{\partial }^{2}w}{\partial x^2}(b/2,y,t)=\frac{{\partial }^{2}w}{\partial x^2}(x,-d,t)=\frac{{\partial }^{2}w}{\partial x^2}(x,d,t)=0\end{array}$$

(38)

Non-dimensionalization is similarly defined:

$$\begin{array}{c}\hfill x^{\ast }=\frac{x}{b}, y^{\ast }=\frac{y}{b}, w^{\ast }=\frac{w}{b}, z^{\ast }=\frac{z}{b}, t^{\ast }=\frac{t}{b/U}, d^{\ast }=\frac{d}{b},\end{array}$$

(39)

and again, getting rid of stars, the non-dimensional problem is

$$\frac{{\partial }^{4}w}{\partial x^4}+2\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+\frac{{\partial }^{4}w}{\partial y^4}+\frac{{\partial }^{2}w}{\partial t^2}+\mu \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}\right)=0,$$

(40)

$$\begin{array}{c}\hfill w(-1/2,y,t)=w(1/2,y,t)=w(x,-1,t)=w(x,1,t)=0,\end{array}$$

(41)

$$\begin{array}{c}\hfill \frac{{\partial }^{2}w}{\partial x^2}(-1/2,y,t)=\frac{{\partial }^{2}w}{\partial x^2}(1/2,y,t)=\frac{{\partial }^{2}w}{\partial x^2}(x,-1,t)=\frac{{\partial }^{2}w}{\partial x^2}(x,1,t)=0\end{array}$$

(42)

with the same $\mu$ and u as before. In this case, the solution was assumed to be in the following form:

$$\begin{array}{c}\hfill w(x,y,t)=v\left(x\right)sin\frac{\pi y}{d}.e^{i\omega t}\end{array}$$

(43)

This is the first term of a Fourier expansion in $\frac{y}{d}$, but the common wisdom is that the higher terms do not affect the results [7,11]. The reason for this can be explained as follows: Like any other linear vibration problem, the flutter of a flat plate or strip will have mode shapes that depend on the mode numbers (basically, the numbering of the eigenvalues starting from the simplest mode). In the simplest case, the plate will vibrate like a single sinus wave, without changing the sign between its boundaries. The higher modes will have sign changes and nodes within the rectangular region. These higher modes will also have higher vibration frequencies. From a physical point of view, it is reasonable to expect that the plate will more easily vibrate in the shape of a simple sinus wave rather than a complicated shape with nodes and sign changes. Therefore, the leading sinus wave will exhibit unstable vibrations before any of the higher modes. Therefore, using (43), (40)–(42) become

$$\frac{{\partial }^{4}v}{\partial x^4}-2{\left(\frac{\pi }{d}\right)}^2\frac{d^{2}v}{d{x}^2}+\mu u\frac{dv}{dx}+K_{2}v=0,$$

(44)

$$\begin{array}{c}\hfill v(-1/2)=v(1/2)=\frac{d^{2}v}{d{x}^2}(-1/2)=\frac{d^{2}v}{d{x}^2}(1/2)=0,\end{array}$$

(45)

$$\begin{array}{c}\hfill K_2={\left(\frac{\pi }{d}\right)}^4-{\omega }^2+i\omega \mu \end{array}$$

(46)

The solution of (44) and (45) is carried out exactly as before, where only the expression $Q_n$ changes.

6. Analytical Solution

The solutions of both problems were found by the collocation method above. These problems can also be solved analytically, but the problem with the analytical solution is that the solution procedure involves finding the roots of a quartic equation and the resulting determinant involves transcendent functions and searching for their zeroes is difficult. The collocation method is simpler to apply and is similar to the analytical solution in terms of performance. The solutions of both (22) and (44) were sought in the form

$$\begin{array}{c}\hfill v=C{e}^{rx}.\end{array}$$

(47)

Substituting (47) and (22) gives

$$\begin{array}{c}\hfill r^4+\mu ur+K_1=0,\end{array}$$

(48)

and (44) gives

$$\begin{array}{c}\hfill r^4-2{\left(\frac{\pi }{d}\right)}^2{r}^2+\mu ur+K_2=0.\end{array}$$

(49)

Denoting the roots of these equations as $r_1$, $r_2$, $r_3$, and $r_4$, the solutions can be written as

$$\begin{array}{c}\hfill v=C_1{e}^{r_{1}x}+C_2{e}^{r_{2}x}+C_3{e}^{r_{3}x}+C_4{e}^{r_{4}x}.\end{array}$$

(50)

Applying the boundary conditions (23) and (45) both give

$$\begin{array}{c}\hfill \left[\begin{array}{cccc}{e}^{\frac{1}{2}{r}_1}& e^{\frac{1}{2}{r}_2}& e^{\frac{1}{2}{r}_3}& e^{\frac{1}{2}{r}_4}\\ e^{-\frac{1}{2}{r}_1}& e^{-\frac{1}{2}{r}_2}& e^{-\frac{1}{2}{r}_3}& e^{-\frac{1}{2}{r}_4}\\ r_1^2{e}^{\frac{1}{2}{r}_1}& r_2^2{e}^{\frac{1}{2}{r}_2}& r_3^2{e}^{\frac{1}{2}{r}_3}& r_4^2{e}^{\frac{1}{2}{r}_4}\\ r_1^2{e}^{-\frac{1}{2}{r}_1}& r_2^2{e}^{-\frac{1}{2}{r}_2}& r_3^2{e}^{-\frac{1}{2}{r}_3}& r_4^2{e}^{-\frac{1}{2}{r}_4}\end{array}\right]\left[\begin{array}{c}{C}_1\\ C_2\\ C_3\\ C_4\end{array}\right]=0\end{array}$$

(51)

The difference between the two solutions is due to the fourth-degree algebraic Equations (48) and (49). For a non-trivial solution, the determinant of the coefficients in (51) should be zero.

$$\begin{array}{c}\hfill \left|\begin{array}{cccc}{e}^{\frac{1}{2}{r}_1}& e^{\frac{1}{2}{r}_2}& e^{\frac{1}{2}{r}_3}& e^{\frac{1}{2}{r}_4}\\ e^{-\frac{1}{2}{r}_1}& e^{-\frac{1}{2}{r}_2}& e^{-\frac{1}{2}{r}_3}& e^{-\frac{1}{2}{r}_4}\\ r_1^2{e}^{\frac{1}{2}{r}_1}& r_2^2{e}^{\frac{1}{2}{r}_2}& r_3^2{e}^{\frac{1}{2}{r}_3}& r_4^2{e}^{\frac{1}{2}{r}_4}\\ r_1^2{e}^{-\frac{1}{2}{r}_1}& r_2^2{e}^{-\frac{1}{2}{r}_2}& r_3^2{e}^{-\frac{1}{2}{r}_3}& r_4^2{e}^{-\frac{1}{2}{r}_4}\end{array}\right|=0\end{array}$$

(52)

This gives $\omega$ as a function of the other parameters in the problem. Although this is the analytical solution, the evaluation of the flutter velocity using this is quite difficult compared with the numerical collocation method. First, the evaluation of determinants is known to be numerically problematic, needing higher-precision computations. The second difficulty is the evaluation of the roots of the quartic equation. Although there is an exact formula for quartic equations, using numerical methods is advantageous here also.

7. Results and Discussion

Collocation solutions were found for $N=10,12$, and 14 and sufficient convergence was observed for $N=10$.

Table 1. The convergence of the results while increasing the number of collocation points.

Number of Collocation Points	${\mathit{u}}_{\mathit{f}}$ (Non-Dimensional Flutter Velocity)
8	90.9697280425894
10	91.1667370035824
12	91.1563146276537
14	91.1564970192325

shows the convergence of the results for the flutter velocity $u_f$ while increasing the number of collocation points. The collocation points were chosen to be equally spaced for $x=-1/2,...,+1/2$ in all cases.

The solution of (48) (analytical solution) was found by an iterative shooting method [27,28,29].

Table 2. Numerical results for a randomly chosen $d=3$ value of the two methods.

	Collocation Method		Analytical Results
$\mathit{\mu }$	Before Flutter	After Flutter	Before Flutter	After Flutter
1.006255	358.920571	359.571969	359.467520	359.571969
1.242290	291.097486	291.372886	291.315280	291.372886
1.572273	230.265262	230.461233	230.377137	230.461233
2.053581	176.682416	176.762617	176.671744	176.762617
2.795152	130.248788	130.299136	130.23840	130.299136
4.025020	90.6746705	91.1957893	91.156960	91.195789
6.289093	59.2604927	59.6428185	59.604512	59.642818
11.18061	36.0023423	36.1616447	36.111424	36.161644

shows the comparisons between the collocation and analytical solutions.

The collocation solution agreed with the analytical solution. To derive the solutions analytically, we used another method, the shooting method [27,29], which agreed with the collocation method solutions, as shown in

Table 3. Comparison of some results with NACA report ([30]).

Material	Experimental Result	Theory
Steel	254.7695 (panels with tension, clamped front and rear)	263.5333
Steel	254.7695 (buckled panels, clamped front and rear)	263.5333
Brass	No flutter (panels with tension, clamped front and rear)	919.6203
Aluminum	No flutter (buckled panels, clamped on four edges)	397.2558

. As mentioned, the collocation method gives results quickly; the analytical solution here is meant to check the correctness of the collocation solution. Another validation of the results is presented here by way of comparison with Dowell 1966, which shows similar results, bearing in mind that in the mentioned work, full potential theory was used, and the solution was carried out for an infinite plate in the direction of flow. The results obtained in this study for small d agreed with Dowell [7,25]. Although this is a theoretical study of flutter, a simple comparison with experiment was carried out. The type of experiment for the solution attempted here seems to be rare. The NACA Technical Note 3914 (1957) ([30]) lists some comparable measurements on thin panels. In that experimental study, some of the boundary conditions were also different. Furthermore, it was extremely difficult to experimentally imitate the plate boundary conditions. But, the magnitude of the numbers should give some idea about the validity of the theoretical solution. The notation and the system of units used were also different. But, after an adjustment, the following values for the experiment and the present computations were obtained.

Considering the differences in the setup and the difficulty in experimentally realizing the boundary conditions, it can be safely stated that the computations were in agreement with the experiment. For the cases of brass and aluminum, the computed flutter velocities were beyond the capabilities of the experiment. Therefore, it is natural that no flutter was observed.

7.1. Infinite Strip

Figure 4 shows the non-dimensional flutter velocity as a function of the mass ratio. As expected, the flutter velocity decreased with increasing mass ratio. Since the curve in Figure 4 seems to have a simple structure, a mathematical expression between $u_f$ and $\mu$ could be developed by using curve fitting. This was done by using Matlab 2023b and the result was

$$\begin{array}{c}\hfill u_f=\frac{333.7+3.249\mu +0.1461{\mu }^2}{\mu -0.02}\end{array}$$

(53)

This equation gives the same points as in Figure 4 within an error of 1 percent.

7.2. Rectangular Plate

For a rectangular plate, there are two parameters that the flutter velocity depends on: the mass ratio and plate width ratio. Figure 5 shows the flutter velocity as a function of the plate width ratio for various values of the mass ratio for $d>1$. The results for $d<1$ are shown in a separate Figure 5. It was observed that the flutter velocity decreased and asymptotically converged to the value for the infinite strip as the plate width increased.

To generalize (53) to a rectangular plate, it was found to be useful to look at the graph of $d^2{u}_f$ (Figure 6). We assumed the relationship between the flutter velocities for finite and infinite width plates to be

$$\begin{array}{c}\hfill u_f=u_{\infty }+\frac{\alpha }{d^2}.\end{array}$$

(54)

Here, $u_{\infty }$ is the expression (53), and by curve fitting, $\alpha$ was found to be

$$\begin{array}{c}\hfill \alpha =1.125+\frac{159.4}{\mu }+\frac{12.02}{{\mu }^2}.\end{array}$$

(55)

Equation (54) gives the flutter velocity in Figure 6 within 1 percent for $d>1$ and $\mu <25$. To give an example, for aluminum (with a density $\rho =2720$ $\mathrm{kg}/{\mathrm{m}}^3$, elasticity modulus $E=70$ GPa, Poisson’s ratio $\nu =0.3$) plate of thickness 5 mm for $b=1$ m and $d=3$ m, (54) gave the non-dimensional flutter velocity as 91.7847 (91.166737 with the collocation method) and the actual flutter velocity was 704.5209 m/s (699.7770 m/s with the collocation method). With the use of (54), design engineers can easily and rapidly predict the flutter velocity. This is the main novelty of this paper, i.e., to predict the flutter velocity of a plate in different situations.

8. Conclusions

To summarize, the flutter velocity of a plate structure under very varied geometric conditions can be estimated with the simple algebraic expressions given in this study. These expressions were derived by curve fitting to data obtained by the collocation method. To validate the collocation results, the analytical solution was also found and seen to give the same results. The novelty of the solution is the correlation between the flutter velocity and the other model parameters. In practice, these correlation equations should be very useful for design engineers in order to predict the flutter velocity. The range of the expressions given cover a very wide range of mass ratios (from 1 to 100) and a corresponding range of non-dimensional flutter velocities (from very low velocities up to about 6500).