Flutter of a Plate at High Supersonic Speeds

1. Introduction

Fluid-elastic structure interactions are ubiquitous in many engineering disciplines. There are historically famous examples that demonstrate 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 which, if exceeds a certain value (flutter velocity), will cause instability. Main thrust of fluid-elastic structure interaction research is the determination of flutter velocity.

The type of problem considered in this study is especially important in aerospace structures and is usually given the name panel flutter. Uzal et al. deal with incompressible and irrotational flow in a cylindrical channel to find the flutter velocity analytically [1]. Epureanu et al. investigated vibration-based damages with the chances in material and/or stiffness properties of structures. Kapkin et al. investigate the membrane vibrations located at the stagnation point of the flow and they give the flutter velocity for the system [2]. Vedeneev considers panel flutter at low supersonic speeds by using piston theory [3]. E. H. Dowell gives flutter velocity graphs for an infinite plate lays with the same direction of fluid motion by using the potential theory for different Mach Numbers [4]. Uzal et al. give an analytical solution for a plate placed in a rigid channel which fluid flows in [5]. Durak B. investigated the plate vibrations by using potential flow theory in his PhD thesis [6]. Also, some of the researchers are also try to avoid this flutter point by applying a force or moving one boundary to control the motion of the plate. Uzal and Korbahti control resonance frequencies of a rectangular plate vibrations by applying discrete force by measuring the displacement of the plate at a point [7]. Sezgin et al give a boundary backstepping control method to stabilize the flow-induced vibrations for a membrane [8]. Tubaldi et al. investigate a periodically supported flexible plate under flowing fluid axial flow which is in a bounded channel by a rigid wall [9]. The effects of the system parameters on the stability of the plate are discussed.

2. Infinite Strip Plate

Figure 1 shows a schematic of the first problem. The flow occupies the region $z>0$ and is in $x-$direction with uniform velocity U. The region in $xy-$plane $-b/2<x<b/2$, $-\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 will cause small perturbations on flow velocity. The vibrations of the plate are governed by

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

(1)

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

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

(2)

is flexural rigidity; E is young modulus and $\nu$ is the Poisson ratio. The last term in 1 denotes the fluid pressure on the plate. Although the plate is moving, 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 lack of flow separation. In general, the linearized form of compressible potential equation can be used, but here a simpler approximation called piston theory will be adopted. Piston theory is widely used and basically states that local pressure is proportional to local slope of the plate [10].

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

(3)

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$$

(4)

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 will be 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}$$

(5)

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

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

(6)

Starred quantities are non-dimensional. Substituting in (4) and (5) and simplifying, the result is, 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$$

(7)

and the boundary conditions

$$\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}$$

(8)

here

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

(9)

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}$$

(10)

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

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

(11)

Then, (7) and (8) become

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

(12)

$$\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}$$

(13)

where

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

(14)

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

(15)

(12) and (13) is an eigenvalue problem for $\omega$. The vibrations of the plate will 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$ the non-dimensional flutter velocity. The problem defined by (12) and (13) will be solved using the collocation method. For this purpose, the approximate solution is expressed as

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

(16)

where ${\varphi }_n$ suitably chosen base functions and $C_n$ are constants to be determined. Base functions are to be chosen so as to satisfy the boundary conditions (12). Here, a family of polynomials will be 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}$$

(17)

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

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

(18)

where, for brevity, we defined

$$\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}$$

(19)

where

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

(20)

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

(21)

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,...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}$$

(22)

For 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)& ...& Q_N(x_N,\omega ,\mu ,u)\end{array}\right]=0\end{array}$$

(23)

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

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

(24)

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

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

(25)

Figure 1. Fluid and strip plate coupled system.

Figure 1. Fluid and strip plate coupled system.

3. Rectangular Plate

Figure 2 shows a schematic of the second problem considered. Now the plate also has a finite width 2d; $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$$

(26)

The plate is assumed to be simply-supported at all sides, so 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}$$

(27)

$$\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}$$

(28)

Non-dimensionalization is defined similarly

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

(29)

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,$$

(30)

$$\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}$$

(31)

$$\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}$$

(32)

with the same $\mu$ and u as before. In this case, the solution is assumed 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}$$

(33)

This is the first term of a Fourier expansion in , but the common wisdom is that the higher terms do not affect the results ([4,7]). (30)–(32) 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,$$

(34)

$$\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}$$

(35)

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

(36)

The solution of (34,35) is carried out exactly as before, only the expression $Q_n$ changes.

4. Analytical Solution

The solutions of both problems were performed by collocation method above. These problems can also be solved analytically; but the problem with 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 its zeroes is difficult. The collocation method is simpler to apply and is at similar to analytical solution in terms of performance. The solutions of both (12) and (34) are sought in the form

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

(37)

Substituting (37) and (12) gives

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

(38)

and (34) gives

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

(39)

Denoting the roots of these equations $r_1,r_2,r_3,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}$$

(40)

Applying the boundary conditions (13) and (35) 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}$$

(41)

The difference between the two solutions is due to the fourth-degree algebraic Equations (38) and (39). For non-trivial solution, the determinant of the coefficients in (41) 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}$$

(42)

This gives omega as a function of the other parameters in the problem.

5. Results and Discussion

Collocation solutions were performed for $N=10,12$ and 14 and sufficient convergence was observed for $N=10$. Table 1 shows convergence of the results for flutter velocity $u_f$ while increasing the number of collocation points. The collocation points were chosen to be equally spaced between $x=-1/2,..+1/2$ in all cases.

The solution of (38) (analytical solution) was carried out by an iterative shooting method [9,10,11]. Table 2 shows comparisons between collocation, and analytical solutions.

The collocation solution agrees with the analytical solution. As was mentioned, collocation method gives results quickly; analytical solution here is meant to check the correctness of the collocation solution. As another validation of the results presented here, comparison with Dowell 1966 shows similar results, bearing in mind that in the mentioned work, full potential theory was used, and the solution was carried out for a plate infinite in the direction of flow The results obtained in this study for small d agree with Dowell.

5.1. Infinite Strip

Figure 3 shows the non-dimensional flutter velocity as a function of mass ratio. As expected, flutter velocity decreases with increasing mass ratio. Since the curve in Figure 3 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 and the result is

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

(43)

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

5.2. Rectangular Plate

For rectangular plate, there are two parameters that the flutter velocity depends on; mass ratio and plate width ratio. Figure 4 shows the flutter velocity as a function of plate width ratio for various values of mass ratio, for $d>1$. The results for $d<1$ are shown in a separate Figure 4. It is observed that the flutter velocity decreases and asymptotically converges to the value for the infinite strip as the plate width increases.

To generalize (43) to rectangular plate, it was found to be useful to look at the graph of $d^2{u}_f$ (Figure 5). Assuming the relation 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}$$

(44)

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

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

(45)

(44) gives the flutter velocity in Figure 5 within 1 percent for $d>1$ and $\mu <25$ To give an example, for Aluminum (with density $\rho =2720$ $kg/m^3$, Elasticity Modulus $E=70$ GPa and Poisson ratio $\nu =0.3$ plate of thickness 5 mm for $b=1m$ and $d=3m$ Equation (44) gives the non-dimensional flutter velocity as 91.7847 (91.166737 with collocation method) and the actual flutter velocity is 704.5209 m/s (699.7770 m/s with collocation method).

6. Conclusions

To summarize, the flutter velocity of a plate structure, under very varied geometric conditions can be estimated with simple algebraic expressions given in this study. These expressions were derived by curve fitting to data obtained by collocation method. To validate the collocation results, analytical solution was also carried out and seen to give the same results.