1. Introduction
The flapping motion of a flexible flag or a flexible plate immersed in a parallel flow is a fundamental problem in flow-induced vibrations (Shelley & Zhang Reference Shelley and Zhang2011), with numerous applications such as propulsion (Lighthill Reference Lighthill1975; Michelin & Llewellyn Smith Reference Michelin and Llewellyn Smith2009; Shoele & Zhu Reference Shoele and Zhu2012, Reference Shoele and Zhu2013), thermal management (Miller Reference Miller1960; Kim & Davis Reference Kim and Davis1995; Guo & Païdoussis Reference Guo and Païdoussis2000a ), airfoil flutter (Tang, Yamamoto & Dowell Reference Tang, Yamamoto and Dowell2003), snoring (Balint & Lucey Reference Balint and Lucey2005; Tetlow & Lucey Reference Tetlow and Lucey2009) and musical instruments such as clarinets, oboes, saxophones, etc. (Sommerfeldt & Strong Reference Sommerfeldt and Strong1988; Backus Reference Backus2005; Dalmont & Frappe Reference Dalmont and Frappe2007). In almost all of these dynamic systems (such as human snoring, propulsion, etc.), the fluttering response of a flexible plate in a confined or unconfined flow field plays a major role, and identification of its dominant oscillatory frequency and most unstable mode shape is essential. The stability threshold and the dominant oscillation frequency of such systems are also the pivotal design features of new passive aeroelastic heat transfer enhancement technology in microchannels (Herrault et al. Reference Herrault, Hidalgo, Ji, Glezer and Allen2012; Shoele & Mittal Reference Shoele and Mittal2014).
While experiments and numerical simulations can be used to study this problem, theoretical analysis of the corresponding inviscid problem offers the attractive possibility of conducting a comprehensive parameter survey, and to develop deep insight into the dynamics of this problem. Inviscid analysis is particularly relevant since the initial stages of the development of the flutter instability are expected to be mostly inviscid in nature. A variety of theoretical approaches have been proposed by Theodorsen (Reference Theodorsen1935), Kornecki, Dowell & O’Brien (Reference Kornecki, Dowell and O’Brien1976), Argentina & Mahadevan (Reference Argentina and Mahadevan2005), Eloy, Souilliez & Schouveiler (Reference Eloy, Souilliez and Schouveiler2007), Alben (Reference Alben2008a ) and Alben & Shelley (Reference Alben and Shelley2008) to study the stability of a single plate in unbounded flow, but the introduction of wall confinement complicates the analysis significantly.
The stability of a flexible plate inside a two-dimensional channel has been studied numerically with fully viscous simulations at low and moderate Reynolds number (Tetlow & Lucey Reference Tetlow and Lucey2009; Shoele & Mittal Reference Shoele and Mittal2014), or using inviscid simulations with a boundary-element method (Howell et al. Reference Howell, Lucey, Carpenter and Pitman2009). These studies report a strong correlation between the distance of the plate to the channel walls and the dynamics of the plate. The three-dimensional effects of confinement in the spanwise direction have been studied experimentally (Doaré, Mano & Ludena Reference Doaré, Mano and Ludena2011a ) and also analytically (Doaré, Sauzade & Eloy Reference Doaré, Sauzade and Eloy2011b ) and it has been concluded that the reduction in the gap between the plate and the walls in the spanwise direction has a destabilizing effect on the dynamics of the system, especially for plates with larger inertia. Finally, Auregan & Depollier (Reference Auregan and Depollier1995) used a quasi-static frictionless fluid model to study a similar problem of aeroelastic instability of the human soft palate and provided an explanation for the occurrence of snoring.
One study of direct relevance to the current work is the recent work of Alben (Reference Alben2015) who used a vortex-sheet method coupled with a large-deflection model for the plate to simulate the large deflection dynamics of a flag in a channel. The study indicated that for heavier flags, greater channel confinement increases the instability of the flag whereas for lighter flags, confinement has minimal influence. The advantage of this method is that is allows for the study of large-deflection dynamics, however, the non-trivial expense of these simulations make them less suitable for the large parameter survey of the stability boundaries of this configuration.
Another study of particular relevance to the current work is the analysis of the stability of a plate inside a two-dimensional channel by Guo & Païdoussis (Reference Guo and Païdoussis2000b ). This study however is limited to the situation in which the plate is placed in the middle of the channel. Moreover, this study employs an upstream wake model as a way to deal with the leading-edge singularity. For the problem of stability of a single plate in unbounded flow, Eloy et al. (Reference Eloy, Lagrange, Souilliez and Schouveiler2008) and Michelin & Smith (Reference Michelin and Smith2009) showed that the upstream wake model could lead to a good prediction of the critical velocity of instability for light plates, but deviates from the analytical solution for heavier plates. The reason behind this discrepancy is that the characteristic spatial scales of deformation in heavy plates are of the order of the plate length and therefore, the dynamics are more sensitive to the description of the wake. This problem is more severe when the plate is placed inside a confined channel, as the inclusion of an auxiliary non-physical upstream wake, while simplifying the mathematical analysis, alters the dynamical properties of the system and changes its stability behaviour in ways that are not fully understood. To eliminate the non-physical effects of the upstream wake, Howell et al. (Reference Howell, Lucey, Carpenter and Pitman2009) employed a regularized boundary-element method to discretize both the plate and the channel walls and studied the transient dynamics of the plate. The computational expense of this approach however increases with time as a result of the continuous creation of free vortices at the trailing edge and their convective propagation to the wake, and thus is not suited for the prediction of the long-time system response.
In the current paper, an analytical approach is proposed that overcomes the shortcomings of the method of Guo & Païdoussis (Reference Guo and Païdoussis2000b ) for calculation of the stability of a flexible plate confined inside an inviscid two-dimensional channel flow. The current method employs a more physically realistic representation of the plate wake within the context of an inviscid flow, and provides a theoretical prediction of the stability boundaries and vibrational characteristics of a plate confined in a channel. In particular, a novel Green’s function for the vortex sheet of the plate wake inside the channel, as well as its rapidly convergent series approximation, are proposed. The stability of the system is then investigated by solving a nonlinear eigenvlaue problem for a large set of system parameters. A similar method for a plate in unbounded flow has previously been proposed by Alben (Reference Alben2008a ); thus the current work may be considered an extension of this previous work to a plate with wall confinement.
The roles of the plate inertia and bending stiffness on its flutter stability for different channel confinements are examined and the effect of asymmetric placement of the plate inside the channel on its stability is also quantified. The scope of the current paper is to analytically predict the stability boundary of the problem as opposed to numerically solving for the nonlinear dynamics of the flag in a confined channel, as in Alben (Reference Alben2015). In fact, the analytical methodology in current study allows a more refined portrait of the stability boundary of the system compared to what has been previously reported using more expensive numerical simulations; the current method also provides an efficient tool to optimize the system in presence of many free parameters, and this is one of the key motivations for conducting this study.
2. Methodology
In the current analysis, we focus on a two-dimensional problem in which a thin, flexible plate of length
$L$
, mass per unit length
$m_{s}$
and bending rigidity
$k_{b}$
is clamped at its leading edge (figure 1
a), and placed at a distance
${\it\eta}$
from the lower wall of a channel with height
$H$
. The equations are non-dimensionalized using the fluid density
${\it\rho}$
, incoming flow velocity
$U_{0}$
and
$L$
, resulting in four characteristic non-dimensional parameters:
$M^{\ast }={\it\rho}L/m_{s}$
(the mass ratio);
$U^{\ast }=U_{0}L\sqrt{m_{s}/k_{b}}$
(the reduced velocity);
$h=H/L$
(the channel confinement ratio); and
$d={\it\eta}/H$
(the plate placement ratio) (Guo & Païdoussis Reference Guo and Païdoussis2000b
; Eloy et al.
Reference Eloy, Lagrange, Souilliez and Schouveiler2008).
Figure 1. (a) Schematic view of a cantilever beam in a channel subjected to axial flow and (b) image systems of bounded and free vortex sheets in the
$y$
-direction used to impose no-penetration boundary conditions on the walls. Blue lines have similar sign to the actual vortex sheet and red lines have opposite sign.
The flow is assumed to be potential; the span of the plate is infinite; and the effect of viscosity is confined to the thin two-dimensional dipole bound vortex sheet along the plate and the deformable free-vortex sheet associated with the wake of the plate. Furthermore, the displacement of the plate from rest is assumed to be small with the non-dimensional vertical displacement of the plate being of the order of
${\it\epsilon}$
. The equation of motion of the plate to
$O({\it\epsilon})$
can be expressed as the Euler–Bernoulli beam equation,
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \partial _{t}^{2}{\it\xi}=-\frac{1}{{U^{\ast }}^{2}}\partial _{x}^{4}{\it\xi}-M^{\ast }[p],\\ \displaystyle {\it\xi}(0,t)=\partial _{x}{\it\xi}(0,t)=\partial _{x}^{2}{\it\xi}(1,t)=\partial _{x}^{3}{\it\xi}(1,t)=0,\end{array}\right\}\end{eqnarray}$$
where
$t$
is non-dimensional time,
$x$
is the longitudinal coordinate measured from the leading edge of the plate which to
$O({\it\epsilon})$
also coincides with the Lagrangian coordinate along the plate (
$0\leqslant x\leqslant 1$
).
$[p]$
is the pressure difference between the upper and lower surfaces of the plate and
${\it\xi}$
is the displacement of the plate.
${\it\xi}$
and
$[p]$
are related to the vortex sheet strength
${\it\gamma}$
through the Biot–Savart kernel and the unsteady Bernoulli equation (Saffman Reference Saffman1992) which, to
$O({\it\epsilon})$
, are
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (\partial _{t}+\partial _{x}){\it\xi}=\frac{1}{2{\rm\pi}}\,\unicode[STIX]{x2A0D}_{0}^{1}{\it\gamma}(x^{\prime },t)G(x,x^{\prime })\,\text{d}x^{\prime }+w(x,t;\infty ),\\ \displaystyle \partial _{x}[p]=(\partial _{t}+\partial _{x}){\it\gamma},\end{array}\right\}\end{eqnarray}$$
respectively, where
$G(x,x^{\prime })$
is the kernel function for the linearized problem satisfying the no-penetration boundary conditions on the top and bottom walls, and
$w(x,t;x_{M})=(1/2{\rm\pi})\int _{1}^{x_{M}}{\it\gamma}(x^{\prime },t)G(x,x^{\prime })\,\text{d}x^{\prime }$
, is the vertical velocity induced at
$x$
by the free-vortex sheet extending from 1 to
$x_{M}$
.
$G(x,x^{\prime })$
for a vortex sheet between two walls can be computed by the ‘method of images’ by adding an infinite image system in the
$y$
direction as shown in figure 1(b) (Greengard Reference Greengard1990). The image system consists of vortex sheets of the same strength placed at
$(x,\text{d}h+2jh)$
and vortex sheets with opposite strength placed at
$(x,-\text{d}h+2jh)$
with
$j=(-\infty ,\ldots ,\infty )$
, as shown in figure 1(b). After some algebra, as described in appendix A, the
$O({\it\epsilon})$
approximation of
$G(x,x^{\prime })$
can be expressed as
$$\begin{eqnarray}\displaystyle G(x,x^{\prime })=\frac{{\it\alpha}}{2}\coth \left(\frac{{\it\alpha}}{2}r\right)-\frac{\displaystyle \frac{{\it\alpha}}{2}\sinh ({\it\alpha}r)}{\cosh ({\it\alpha}r)-\cos (2{\rm\pi}\,d)},\quad \text{with }r=x-x^{\prime },\text{and }{\it\alpha}=\frac{{\rm\pi}}{h}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The dependence of
$G(x,x^{\prime })$
on
$r=x-x^{\prime }$
for different values of
$h$
and
$d$
is shown in figure 2. It appears that the reduction of either
$h$
or
$d$
results in more rapid decrease of
$G(x,x^{\prime })$
with
$r=x-x^{\prime }$
, and consequently from (2.2), the Biot–Savart kernel is mostly affected by the near-field vortices and less by the free-vortex sheet in the wake of the plate.
In addition, the pressure jump satisfies
$[p](1,t)=0$
, and
${\it\gamma}$
is finite at
$x=1$
due to the Kutta condition (Saffman Reference Saffman1992). It is known that the solution of
${\it\gamma}(x,t)$
has inverse-square-root singularities at
$x=0,1$
(Alben Reference Alben2008b
; Muskhelishvili & Radok Reference Muskhelishvili and Radok2008). Therefore, by denoting the finite part of
${\it\gamma}$
by
${\it\nu}$
, defined as
${\it\gamma}={\it\nu}/\sqrt{x(1-x)}$
, the Kutta condition can be written as
${\it\nu}(1,t)=0$
.
Finally, the total circulation in the free-vortex sheet beyond
$x$
is defined as
${\it\lambda}(x,t)=\int _{x}^{\infty }{\it\gamma}(x^{\prime },t)\,\text{d}x^{\prime }$
. Since the circulation associated with fluid particles is conserved, known as Kelvin’s theorem (Saffman Reference Saffman1992), and the free-vortex sheet originates at the trailing edge, the circulation of other points in the wake can be determined from
${\it\lambda}(1,t)$
through,
$$\begin{eqnarray}{\it\lambda}(x,t)={\it\lambda}(1,t_{0}),\end{eqnarray}$$
where
$t_{0}$
is the time when the current fluid particle departed from the trailing edge.
Figure 2. The dependence of
$O({\it\epsilon})$
approximation of the kernel function
$G(x,x^{\prime })$
(2.3), on
$r=x-x^{\prime }$
for (a) different
$h$
with
$d=0.5$
, and (b) different
$d$
with
$h=1.0$
.
To study the stability of the system, we assume a time-periodic solution for the vertical displacement and pressure jump and decompose them into their normal modes represented by
$Y$
and
$[P]$
respectively. Thus for example, the vertical position of the plate is defined as
${\it\xi}=\text{Re}(Y(x)\text{e}^{\text{i}{\it\omega}t})$
, where
${\it\omega}$
is the complex eigenfrequency of the non-dimensional system of equations. In this description,
$\text{Re}({\it\omega})$
is the angular frequency and
$-\text{Im}({\it\omega})$
the growth rate of the normal mode. An unstable mode corresponds to
$\text{Im}({\it\omega})<0$
and vice-versa. The resultant governing system of equations is:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle -{\it\omega}^{2}Y=-\frac{1}{{U^{\ast }}^{2}}\partial _{x}^{4}Y-M^{\ast }[P],\\ \displaystyle (\text{i}{\it\omega}+\partial _{x})Y=\frac{1}{2{\rm\pi}}\,\unicode[STIX]{x2A0D}_{0}^{1}{\it\Gamma}(x^{\prime })G(x,x^{\prime })\,\text{d}x^{\prime }+W(x,{\it\omega};\infty ),\\ \displaystyle \partial _{x}P=(\text{i}{\it\omega}+\partial _{x}){\it\Gamma},\quad \text{with }P(1)=0,\\ \displaystyle {\it\Upsilon}(x=1)=0,\quad \text{with }{\it\Gamma}(x)=\frac{{\it\Upsilon}(x)}{\sqrt{x(1-x)}},\end{array}\right\}\end{eqnarray}$$
where
${\it\Upsilon}$
and
${\it\Gamma}$
are the normal modes of
${\it\nu}$
and
${\it\gamma}$
, respectively. The effect of the free-vortex sheet
$W(x,{\it\omega};\infty )$
takes a much simpler form due to the conservation properties of circulation in the free-vortex wake. In addition, the total circulation in the plate is also periodic,
${\it\lambda}(x=1,t)=\text{Re}({\it\Lambda}_{0}\text{e}^{\text{i}{\it\omega}t})$
where
${\it\Lambda}_{0}$
is the amplitude of
${\it\lambda}$
at
$x=1$
. Given that any fluid particle at
$(x,t)$
was located at the trailing edge (
$x=1$
) at time
$t-(x-1)$
, its circulation can be expressed as
${\it\lambda}(x,t)={\it\Lambda}_{0}\text{e}^{\text{i}{\it\omega}(t-(x-1))}$
and consequently we have,
$$\begin{eqnarray}W(x,{\it\omega};\infty )=-\frac{{\it\omega}{\it\Lambda}_{0}}{2{\rm\pi}}\int _{1}^{\infty }\text{e}^{-\text{i}{\it\omega}(x^{\prime }-1)}G(x,x^{\prime })\,\text{d}x^{\prime }.\end{eqnarray}$$
Because of the nonlinearity of
${\it\omega}$
in (2.6), a nonlinear eigenvalue problem based on (2.5) and (2.6) mush be solved to find
${\it\omega}$
and its associated eigenvectors.
While it is typical in such analysis to truncate the integral in (2.6) to a large but finite limit, this can cause errors for the neutrally stable case. This is because in this state, the system sustains its vibration at a constant amplitude and the strength of the wake does not diminish with downstream distance. To further improve the predictions of the eigenvalue and eigenvectors of the neutrally stable case, we decompose
$W(x,{\it\omega};\infty )$
into a sum of a finite and an infinite integral with the infinite integral computed by the shifting technique (Greengard Reference Greengard1990). In particular, an arbitrary point
$x_{0}>1$
is chosen in the wake, and the integral is decomposed as
$W(x,{\it\omega};\infty )=W(x,{\it\omega};x+x_{0})+\mathscr{M}(x,x_{0})$
where
$W(x,{\it\omega};x+x_{0})$
can be calculated using numerical integration and
$\mathscr{M}(x,x_{0})$
, the contribution of the wake beyond
$x+x_{0}$
on the velocity of the point
$x$
, can be estimated based on the following convergent series:
$$\begin{eqnarray}\mathscr{M}(x,x_{0})=-\frac{{\it\omega}{\it\Lambda}_{0}}{h}\text{e}^{-\text{i}{\it\omega}(x+x_{0}-1)}\mathop{\sum }_{j=1}^{J_{M}}\frac{\text{e}^{-{\it\alpha}jx_{0}}}{{\it\alpha}j+{\it\omega}\text{i}}\sin ^{2}({\rm\pi}\,dj)\quad \text{with }J_{M}\geqslant \frac{h}{{\rm\pi}x_{0}}\sqrt{\frac{{\it\kappa}}{{\it\sigma}}},\end{eqnarray}$$
where
${\it\kappa}=|{\it\omega}{\it\Lambda}_{0}/{\rm\pi}|$
, and
$J_{M}$
is the minimum number of terms needed to achieve a relative precision of
${\it\sigma}$
, and is calculated
a priori
to the solution based on the initial estimation of
${\it\Lambda}_{0}$
as explained in § 3. Details of the derivations of
$\mathscr{M}(x,x_{0})$
and
$J_{M}$
are presented in appendices B and C.
Based on convergence studies, 60 Chebyshev–Lobatto nodes along the plate are chosen to discretize the governing equations. Equations (2.5) and (2.6) along with the series expansion based on (C 3) with
${\it\sigma}=10^{-8}$
are solved iteratively using a Newton–Kantorovitch algorithm (Boyd Reference Boyd2001) (similar to the method used previously by Alben & Shelley (Reference Alben and Shelley2008) and Michelin & Llewellyn Smith (Reference Michelin and Llewellyn Smith2009) for studying the stability of the flag in an unbounded flow), to investigate the stability of the system in the
$(M^{\ast },U^{\ast })$
-plane. In particular, for each channel height and plate placement we initially consider five points with
$M^{\ast }=(-2.0,-1.0,0,1.0,2.0)$
and identify the corresponding
$U^{\ast }$
values of the neutrally stable cases in unbounded flow using the analytical method similar to Kornecki et al. (Reference Kornecki, Dowell and O’Brien1976). The first 20 eigenfrequencies are identified and used to calculate the range of anticipated frequencies of the current system; this range is then divided in to more than 2000 intervals (with 100 steps between successive frequencies) and used for initial estimations of the anticipated
${\it\omega}_{0}$
s. We use a two-stage predictive–corrective approach to solve the nonlinear eigenvalue problem. In the predictive phase, we only consider the effects of the wake for
$x\leqslant 11$
(e.g.
$W(x,{\it\omega};\infty )\simeq W(x,{\it\omega};11L)$
), which has been previously shown to be sufficient for estimation of the wake effects (Alben Reference Alben2008a
), and discretize the wake into 1000 equal distance grid intervals. For each
${\it\omega}_{0}$
value, the equations are then separated based on their power of
${\it\omega}$
and written as a quadratic eigenvalue problem of the form,
$$\begin{eqnarray}[{\it\omega}^{2}\boldsymbol{A}_{2}+{\it\omega}\boldsymbol{A}_{1}+\boldsymbol{A}_{0}({\it\omega}_{0})]\,\boldsymbol{X}=\mathbf{0},\end{eqnarray}$$
where the vector
$\boldsymbol{X}$
is the eigenvector consisting of the plate curvature and hydrodynamic pressure of Chebyshev–Lobatto nodes as well as the scalar
${\it\Lambda}_{0}$
. This equation is solved iteratively and in each iteration the nonlinear term,
$\boldsymbol{A}_{0}({\it\omega}_{0})$
, is kept constant to its value from the most recent calculation of
${\it\omega}$
. The calculated eigenpairs of
${\it\omega}$
and
$\boldsymbol{X}$
are used as initial conditions for more efficient Newton–Kantorovitch calculations of other
$U^{\ast }$
values of each particular
$M^{\ast }$
assuming that the wake is truncated to
$x\leqslant 11$
. For each
$M^{\ast }$
, the
$U^{\ast }$
case with
$\text{Im}({\it\omega})\simeq 0$
is further refined in order to reach to the first estimate of
$U_{c}^{\ast }$
(
$U_{c}^{\ast }$
is the critical reduced velocity that separates stable and unstable regimes).
In the corrective step, the series representation of the wake for
$x\geqslant 11$
is incorporated in order to eliminate the potential effects of wake truncation on the calculation of the stability boundary, and several more Newton–Kantorovitch calculations are performed to find the
$U_{c}^{\ast }$
with less than
$10^{-4}$
error. At the final stage, smaller values of
$U^{\ast }$
are checked to ensure that the calculated
$U_{c}^{\ast }$
is indeed the highest
$U^{\ast }$
for which the plate stays neutrally stable.
Figure 3. The critical velocity
$U_{c}^{\ast }$
as a function of mass ratio
$M^{\ast }$
for a plate in an unbounded flow. The symbols are for the relevant experimental results. Data from several other numerical and analytical studies have also been included for comparison with the current method.
Subsequently, the most recent calculation of
$U_{c}^{\ast }$
for the particular
$M^{\ast }$
is used as an initial prediction of
$U_{c}^{\ast }$
at
$M^{\ast }+{\rm\Delta}M^{\ast }$
with
$\log _{10}({\rm\Delta}M^{\ast })=0.005$
and the above two-stage solving method repeated for this new condition. Each case takes, on-average, approximately 5 s to solve on a personal computer.
3. Results
Figure 3 shows the prediction of the stability of the plate placed in the unbounded uniform flow in
$(M^{\ast },U^{\ast })$
-plane by the current technique and its comparison with previous experimental, analytical and numerical results (e.g. Abderrahmane et al.
Reference Abderrahmane, Païdoussis, Fayed and Dick Ng2012; Huang Reference Huang1995; Michelin, Llewellyn Smith & Glover Reference Michelin, Llewellyn Smith and Glover2008). The unbounded domain is approximated here by choosing a large channel height (
$H=10L$
). The predicted stability curve shows a reasonably good agreement with the previous studies especially those which are based on the similar formulation for unbounded flow, analytical results by Kornecki et al. (Reference Kornecki, Dowell and O’Brien1976) and numerical eigenvalue perditions by Alben (Reference Alben2008a
). The minor mismatch between the current prediction of the stability curve and the analytical prediction by Kornecki et al. (Reference Kornecki, Dowell and O’Brien1976) at small
$M^{\ast }$
cases is primarily due to the truncation error that arises in nonlinear eigenvalue analyses and partially because of the domain truncation at
$H=10L$
. In general, the current method agrees very well with other methods in capturing the peaks of the stability curves, but slightly differs with the other methods at small
$M^{\ast }$
. While there are vertical differences between predicted
$U_{c}^{\ast }$
from the current analysis and others at small
$M^{\ast }$
, the associated growth rates only change to a very small degree with
$U^{\ast }$
, indicating that the frequency estimations of all methods are quite similar.
The presentation of the confined plate results begins by showing in figure 4 a comparison between the current prediction of the stability boundary of the plate placed in the middle of a channel with
$H/L=1$
, with the fully nonlinear fluid–structure interaction simulations by Shoele & Mittal (Reference Shoele and Mittal2014), the previous analytical method by Guo & Païdoussis (Reference Guo and Païdoussis2000b
) with the upstream wake assumption and the stability boundary predicted by Alben (Reference Alben2015) with the regularized large-deflection vortex method. For the fully viscous simulation, it is known from Shoele & Mittal (Reference Shoele and Mittal2014) that in the case of the symmetric placement of the plate inside a channel, the plate dynamics are almost independent of the Reynolds number (based on the mean flow velocity in the channel) as long as
$Re>200$
. We therefore selected
$Re=400$
and a series of 250 Navier–Stokes simulations were carried out to cover a range of
$M^{\ast }$
from 0.03 to 30 and
$U^{\ast }$
from 3 to 16. Each simulation was integrated in time till approximately
$t=100L/U$
in order to eliminate transient effects and obtain representative asymptotic deflection amplitudes. The stability boundary from the current method shows a good agreement with the Navier–Stokes simulations both for the lower range, as well as the higher range, of
$M^{\ast }$
values. In contrast, the predictions from Guo & Païdoussis (Reference Guo and Païdoussis2000b
) compare well with the Navier–Stokes simulations and the current method only for
$M^{\ast }>1$
but depart significantly from them for lower
$M\ast$
values. As previously discussed, this mismatch for
$M^{\ast }<1$
is due to the upstream wake assumption used in Guo & Païdoussis (Reference Guo and Païdoussis2000b
). The current method is also in general agreement with the method of Alben (Reference Alben2015) over the entire range of
$M^{\ast }$
. For
$0.01<M^{\ast }<1$
, the method of Alben (Reference Alben2015) does however predict a instability boundary that is somewhat lower than that from the Navier–Stokes simulations and the current method, and the reason for this is not clear.
Figure 4. Comparison of the stability curves of the confined plate with the confinement ratios
$H/L=1$
and
${\it\eta}/H=0.5$
between the present method (solid line), analytical prediction by Guo & Païdoussis (Reference Guo and Païdoussis2000b
) with the upstream wake assumption (dashed line), numerical prediction using the large deflection vortex method by Alben (Reference Alben2015) and the fully nonlinear fluid–structure interaction simulations at
$Re=400$
from Shoele & Mittal (Reference Shoele and Mittal2014) (contour plot).
$A$
is the peak-to-peak amplitude of the motion at the trailing edge of the plate.
In figure 5 we compare the stability boundaries of the plate for different confinements ranging from
$H/L=10$
down to
$1/10$
. The cases with
$H/L>0.5$
have very similar stability curves, with the main differences appearing at small
$M^{\ast }$
values. When
$M^{\ast }$
reduces below 1,
$U_{c}^{\ast }$
first reduces rapidly and then increases with further decrease in
$M^{\ast }$
. Furthermore, the onset of this increase in
$U_{c}^{\ast }$
happens at smaller
$M^{\ast }$
values for increased confinements. For heavier plates (smaller
$M^{\ast }$
), the confinement effects by the channel walls reduce the critical velocity
$U_{c}^{\ast }$
, indicating that for a similar inflow speed and plate mass ratio, confinement increases the instability of heavy flag. A similar behaviour has been found in our previous viscous flow simulations (Shoele & Mittal Reference Shoele and Mittal2014) as well as the studies of Howell et al. (Reference Howell, Lucey, Carpenter and Pitman2009) and Alben (Reference Alben2015). It is also noted that for small
$H/L$
ratios, here less than
$1/4$
as shown in figure 5, the stability curves are modified significantly; the peaks of the
$M^{\ast }$
–
$U_{c}^{\ast }$
curves move toward smaller
$M^{\ast }$
values, and the corresponding
$U_{c}^{\ast }$
values for these peaks are altered and even new peaks appear in the stability curve. For example, while the
$U_{c}^{\ast }$
curve of
$H/L=1$
(shown in figure 5 with black colour) for the range of
$M^{\ast }$
presented in the figure has four peaks, each corresponding to a particular mode transition, the
$U_{c}^{\ast }$
curve of
$H/L=1/10$
(shown with light green colour in figure 5) exhibits six peaks.
Figure 5. Stability curves for the confined plate for different confinement ratios
$H/L$
. The data in this plot are for
${\it\eta}/H=0.5$
.
The availability of Navier–Stokes model (Shoele & Mittal Reference Shoele and Mittal2014) for this configuration allows us to go beyond a qualitative description of behaviour to a quantitative evaluation of the prediction of confinement on the stability of the flag. Figure 6 shows the variation of the peak-to-peak vertical amplitude of the motion at the trailing edge of the plate with
$U^{\ast }$
at one selected value of
$M^{\ast }=0.2$
, obtained from the fully viscous simulations. As before, for these comparisons, we choose
$Re=UL/{\it\nu}=400$
in the viscous simulations, with
${\it\nu}$
being the kinetic viscosity of the fluid. The plots also identify for each case, the critical
$U^{\ast }$
as predicted from the current method. The Navier–Stokes simulations exhibit, for each case, a steady rise in the asymptotic flutter amplitude (of approximately 2–2.5 orders of magnitude) with
$U^{\ast }$
until a saturation state is reached. As
$H/L$
increases (decreasing confinement), the asymptotic amplitude curves shift to the right, indicating increasing stability for decreasing confinement. The
$U_{c}^{\ast }$
predicted from the current inviscid method falls within the nonlinear growth phase for each case and these values also shift to the right as confinement is decreased. As would be expected of any linear stability analysis, in no case does the
$U_{c}^{\ast }$
predicted from the current method exceed the value of
$U^{\ast }$
, where amplitude saturation is observed to occur. The inset figure compares the
$U_{c}^{\ast }$
predicted from the current analysis with the
$U^{\ast }$
at which the amplitude saturates for the Navier–Stokes model. The comparison indicates that the inviscid linear stability analysis employed here provides a reasonably good lower bound for the
$U_{c}^{\ast }$
in realistic cases.
Figure 6. Comparison of stability boundaries between the current model and fully viscous Navier–Stokes fluid–structure interaction (NSFSI) numerical calculations by Shoele & Mittal (Reference Shoele and Mittal2014) (solid lines) for a heavy flag (
$M^{\ast }=0.2$
). Solid lines with symbols show the non-dimensional amplitude
$A/H$
from the NSFSI simulations for
$H/L=1$
,
$4/5$
,
$1/2$
,
$1/3$
,
$1/4$
and
${\it\eta}/H=1/2$
. Vertical dashed arrows indicate the critical
$U^{\ast }$
predicted by the current analysis for each case. The NSFSI simulations correspond to cases with
$Re=UL/{\it\nu}=400$
, with
${\it\nu}$
being the kinematic viscosity of the fluid. The inset figure compares
$U_{c}^{\ast }$
for the current method with the lowest
$U^{\ast }$
from NSFSI simulation that
$A/H$
reaches its saturated value.
In figure 7, we present the oscillation frequency and the mode shape of the neutral mode shape at both the peaks and troughs of the curves form
$H/L=1$
to
$1/10$
. At each peak, the neutral mode shape of the plate undergoes a gradual mode transition between the corresponding modes of instability of its surrounding troughs. For example, the most unstable mode at the first peak is the second mode with
${\it\omega}=3.11$
and 18.2 for
$H/L=1$
and
$1/10$
, respectively. Also the most unstable mode shapes of the peaks show larger deformation around the leading half of the plate when
$H/L$
ratio is small. Reduction in
$H/L$
also causes mode transitions to appear at smaller
$M^{\ast }$
intervals and increases the oscillation frequency of the neutral mode shape along the stability curve, primarily as a result of the reduction in
$U_{c}^{\ast }$
and higher bending stiffness of the plate.
Figure 7. The neutrally stable mode shapes of the plate at the peaks and troughs of the
$U_{c}^{\ast }$
curve for the cases of (a)
$H/L=1$
, (b)
$H/L=1/2$
, (c)
$H/L=1/4$
, and (d)
$H/L=1/10$
. The real parts of the neutral mode shapes are shown with thick red dashed lines, and the imaginary parts of the neutral mode shapes are shown with thick blue dotted lines. All cases shown here correspond to
${\it\eta}/H=0.5$
. Corresponding symbols for each cases are marked in figure 5.
Figure 8. Stability curves for the confined plate with
$H/L=1.0$
for (a)
${\it\eta}/H=1/5$
, (b)
${\it\eta}/H=1/10$
, (c)
${\it\eta}/H=1/20$
and (d)
${\it\eta}/H=1/40$
. The corresponding curve for
${\it\eta}/H=1/2$
is shown with dashed line for comparison. Snapshots of the mode shapes of the neutral mode at the peaks and troughs of the
$U_{c}^{\ast }$
curve are included with the real part of the neutral mode shape shown with thick red dashed lines, and the imaginary part of the neutral mode shape shown with thick blue dotted lines.
The current model also allows for an examination of the effect of asymmetric placement of the plate in the channel on its stability, and figure 8 shows the effect of plate placement as denoted by
${\it\eta}/H$
on flutter instability for the cases with
$H/L=1$
. The stability curve does not change if the plate moves slightly away from the centre of the channel. However, as the plate approaches one wall, the peaks of
$U_{c}^{\ast }$
start shifting toward smaller
$M^{\ast }$
and the magnitude of
$U_{c}^{\ast }$
at the peaks decreases more for smaller
$M^{\ast }$
than higher
$M^{\ast }$
. It is found that as the plate is placed ever nearer to the wall, new peaks emerge in the stability curve, and the frequency of the neutral mode at these peaks increases more rapidly with the increase of
$M^{\ast }$
. For example, by comparing the frequency of neutral modes of
${\it\eta}/H=1/40$
and
$1/5$
(shown in figure 8
a,d), it can been seen that
${\it\omega}$
of the first three crests of the stability curve are 2.76, 4.41 and 5.76 for
${\it\eta}/H=1/40$
whereas they are at
${\it\omega}=2.32$
, 3.67 and 5.00 for
${\it\eta}/H=1/5$
. In addition, heavy plates with small
$M^{\ast }$
become unstable at smaller
$U^{\ast }$
when they are placed in close proximity of the wall. The mode shapes of the different
${\it\eta}/H$
at the crests and troughs of the stability curves are similar and only show differences in the first half of the plate at higher
$M^{\ast }$
. We note that the neutrally stable mode shapes calculated by the present stability analysis are symmetric, which is a direct consequence of the small deformation assumption in the present formulation. However, it is anticipated that as the plate undergoes large motion, it will interact strongly with the closer wall, leading to asymmetric deflection modes (Alben Reference Alben2015).
It is noted that in the limit of very close placement of the plate to one wall in a viscous channel flow, the viscous force on both sides of the plate would be quite different and the flag experiences much larger hydraulic resistance on the side that faces the closer wall. In this situation, viscosity could play a major role in the dynamic response of the plate and a fully viscous model (e.g. Shoele & Mittal Reference Shoele and Mittal2014) is required to investigate the stability of the system. Interestingly however, previous numerical calculations by Tetlow & Lucey (Reference Tetlow and Lucey2009) also predicted that offsetting the plate toward one wall has a destabilization effect on the dynamics of the plate due to the enhancement of the energy transfer between the fluid and the structure, similar to what has been observed in the current study. In the recent numerical and experimental study by Dessi & Mazzocconi (Reference Dessi and Mazzocconi2015), it is shown that the reduction of the gap between the plate to a solid boundary results in a decrease of the critical fluttering response,
$U_{c}^{\ast }$
, which is consistent to what has been found in this study. The similarities between these studies and current results provide some evidence that the instability mechanism for even highly asymmetric cases at high Reynolds number flow is inviscid in origin, but this hypothesis needs to be fully evaluated.
The above results regarding the destabilizing effect of confinement and wall proximity on heavy flags are interesting and even perhaps counterintuitive, and the physical mechanisms that drive this behaviour are worth discussing here. Following Guo & Païdoussis (Reference Guo and Païdoussis2000b
) and based on the local behaviour of the kernel function
$G(x,x^{\prime })$
in a confined domain ((2.3) and figure 2), the dynamic pressure on the plate in a confined channel can be considered to consist of three main components: the inertial or added-mass effect,
$m_{1}(x){\it\xi}_{tt}$
; the Coriolis force effect,
$2m_{2}(x)U^{\ast }{\it\xi}_{xt}$
; and the centrifugal force effect,
${U^{\ast }}^{2}m_{3}(x){\it\xi}_{xx}$
, where
$m_{1}(x)$
,
$m_{2}(x)$
and
$m_{3}(x)$
depend on the fluttering mode shape of the plate. At a small confinement ratio, it can be shown that
$m_{1}(x)=m_{2}(x)=m_{3}(x)\simeq m$
where
$m$
is a positive constant and linearly proportional to
$M^{\ast }$
. The rate of the work done by the plate on the surrounding fluid can then be well approximated as
$$\begin{eqnarray}{\dot{W}}\simeq m\int _{L}({\it\xi}_{tt}{\it\xi}_{t}+2U^{\ast }{\it\xi}_{xt}{\it\xi}_{t}+{U^{\ast }}^{2}{\it\xi}_{xx}{\it\xi}_{t})\,\text{d}x.\end{eqnarray}$$
If
${\dot{W}}<0$
, the plate receives energy from the fluid and starts oscillating in its most unstable mode shape. The plate in a highly confined channel resembles an oscillating flat plate near a wall, and since
$m$
is real at small
$H/L$
ratio, there is no dissipation of energy from the added-mass term (i.e. the first term on the right-hand side of (3.1)) (Jaiman, Parmar & Gurugubelli Reference Jaiman, Parmar and Gurugubelli2014). This suggests that unlike a plate in unbounded flow in which the energy lost from the added-mass term due to pressure radiation to the far field (physically in terms of acoustic waves), in a very confined channel, the flow acts as a virtual mass and exchanges energy with the plate in a non-dissipative manner. For a typical travelling wave deflection of the plate, it is possible to show that at small
$H/L$
ratio, the second term on the right-hand side of (3.1) (the Coriolis term proportional to
$U^{\ast }$
) always stabilizes the motion of the plate and the third term (the centrifugal term proportional to
${U^{\ast }}^{2}$
) always destabilizes the oscillation of the plate. Therefore as
$H/L$
decreases, the stabilized added-mass term diminishes causing a smaller
$U^{\ast }$
to be required to balance the Coriolis term and the centrifugal term such that the plate can receive energy from the fluid for its oscillation. It equivalently means that a larger bending stiffness at the given flow velocity compared to the unbounded flow is required to break the phase-lock between the pressure force and the plate vertical motion, and to prevent a single-mode fluttering response of the plate. The increase in the required bending stiffness for the flutter stability manifests itself in the higher fluttering frequency for heavy plates. To fully understand and dissect all important mechanisms that prompt fluttering response, future studies in this direction are needed.
4. Conclusions
In summary, we have proposed and implemented a novel analytical approach to study the stability of a plate confined between two parallel walls using a nonlinear eigenvalue analysis. Using this approach, we have computed the critical reduced velocity of the plate as a function of the mass-ratio of the plate for different channel confinements and placements. It is found that increased confinement reduces the critical velocity of heavy plates, shifts the stability curve toward smaller mass-ratios and causes the emergence of new mode shapes along the stability curve. New mode shapes and a decrease in the critical reduced velocity at low mass-ratios (i.e. for heavy plates) is also observed if the plate is asymmetrically placed closer to one wall. The frequency of the unstable mode in these cases also depends strongly on the distance of the plate to the nearby wall. Comparisons of the results of the current model with the orders-of-magnitude more expensive fully coupled Navier–Stokes fluid–structure interaction models indicate that the current inviscid model quite accurately predicts the stability boundaries for a range of confinement ratios.
Finally, the current study suggests that the strong correlation between the vibrational frequency of the plate and its placement in the channel can be utilized to modify, improve or control new heat enhancement technologies that rely on the aeroelastic instabilities of such configurations (Herrault et al. Reference Herrault, Hidalgo, Ji, Glezer and Allen2012) and could also provide insights into the physics of snoring as well as the functioning and design of reed instruments. The current work also provides a theoretical framework that could be utilized to study similar problems in the future.
Acknowledgements
The authors thank Dr S. Alben of the University of Michigan for helpful discussions and for providing data for validation. We also acknowledge helpful discussions with Dr A. Glezer from the Georgia Institute of Technology and support from Air Force Office of Scientific Research (Grant FA9550-13-1-0053) as well as NSF (Grant CBET-1357819) for this work.
Appendix A. Derivation of
$G(x,x^{\prime })$
Assuming that the flow is inviscid and the effect of viscosity is only confined to a thin sheet along the plate and a continuous deformable surface behind the plate with vortex sheet strength
${\it\gamma}(s,t)$
, the flow velocity in the fluid domain can be computed using the Biot–Savart kernel (Saffman Reference Saffman1992).
${\it\gamma}$
can also be considered as a jump in flow tangential velocity across the plate. The resultant equation in its dimensionless form is,
$$\begin{eqnarray}\overline{u}(z_{0})=1+\frac{1}{2\text{i}{\rm\pi}}\int _{C_{b}+C_{f}}\frac{{\it\gamma}(s,t)}{z_{0}-{\it\zeta}(s,t)}\,\text{d}s+\mathscr{I}(z_{0}),\end{eqnarray}$$
where
$z_{0}$
is the complex representation of the position of the field point defined as
$z_{0}=x_{0}+\text{i}y_{0}$
(the bottom and top walls, respectively, are at
$y_{0}=0$
and
$y_{0}=h$
),
$C_{b}$
is the line contour along the plate and
$C_{f}$
is the free-vortex sheet in the wake starting from the trailing edge.
$\overline{u}(z_{0})=u_{x}(z_{0})-\text{i}u_{y}(z_{0})$
is the complex conjugate of the fluid velocity at
$z_{0}$
. In addition,
$\mathscr{I}$
function is added to satisfy the zero normal flow condition at the bottom and top walls.
To find the function
$\mathscr{I}$
, an infinite image system is employed along the
$y$
direction (Greengard Reference Greengard1990). In particular, the image system consists of infinite vortex sheets with positive strength placed vertically at
${\it\xi}(s,t)+2jh\text{i}$
and infinite vortex sheets with negative strength placed vertically at
$\overline{{\it\xi}}(s,t)+2jh\text{i}$
with with
$j=(-\infty ,\ldots ,\infty )$
as shown in figure 1(b).
$\mathscr{I}(z_{0})$
can be explicitly written as,
$$\begin{eqnarray}\mathscr{I}(z_{0})=\frac{1}{2\text{i}{\rm\pi}}\mathop{\sum }_{\substack{ j=-\infty \\ j\neq 0}}^{\infty }\int _{C_{b}+C_{f}}\frac{{\it\gamma}(s,t)}{z_{0}-{\it\zeta}(s,t)+2jh\text{i}}\,\text{d}s-\frac{1}{2\text{i}{\rm\pi}}\mathop{\sum }_{j=-\infty }^{\infty }\int _{C_{b}+C_{f}}\frac{{\it\gamma}(s,t)}{z_{0}-\overline{{\it\zeta}}(s,t)+2jh\text{i}}\,\text{d}s,\end{eqnarray}$$
which allows (A 1) to be expressed as,
$$\begin{eqnarray}\displaystyle \overline{u}(z_{0}) & = & \displaystyle 1+\frac{1}{2\text{i}{\rm\pi}}\int _{C_{b}+C_{f}}\left\{\frac{{\it\gamma}(s,t)}{z_{0}-{\it\zeta}(s,t)}+\mathop{\sum }_{j=1}^{\infty }\frac{2{\it\gamma}(s,t)(z_{0}-{\it\zeta}(s,t))}{(z_{0}-{\it\zeta}(s,t))^{2}+4h^{2}j^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\frac{{\it\gamma}(s,t)}{z_{0}-\overline{{\it\zeta}}(s,t)}-\mathop{\sum }_{j=1}^{\infty }\frac{2{\it\gamma}(s,t)(z_{0}-\overline{{\it\zeta}}(s,t))}{(z_{0}-\overline{{\it\zeta}}(s,t))^{2}+4h^{2}j^{2}}\right\}\text{d}s.\end{eqnarray}$$
Knowing that (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014),
$$\begin{eqnarray}\coth ({\rm\pi}z)=\frac{1}{{\rm\pi}}\left(\frac{1}{z}+2z\mathop{\sum }_{j=1}^{\infty }\frac{1}{z^{2}+j^{2}}\right),\end{eqnarray}$$
we can combine (A 3) and (A 4) and express
$\overline{u}(z_{0})$
as,
$$\begin{eqnarray}\overline{u}(z_{0})=1+\frac{1}{2\text{i}{\rm\pi}}\int _{C_{b}+C_{f}}{\it\gamma}(s,t)\underbrace{\displaystyle \frac{{\rm\pi}}{2h}\left(\coth \left[\displaystyle \frac{{\rm\pi}}{2h}(z_{0}-{\it\zeta}(s,t))\right]-\coth \left[\frac{{\rm\pi}}{2h}(z_{0}-\overline{{\it\zeta}}(s,t))\right]\right)}_{\mathscr{G}(z_{0},{\it\zeta})}\text{d}s.\end{eqnarray}$$
Taking the limit of (A 5) as
$z$
approaches a point on
$C_{b}$
, the flow velocity can be written as,
$$\begin{eqnarray}\overline{u}(s,t)=1\!+\frac{1}{2\text{i}{\rm\pi}}\,\unicode[STIX]{x2A0D}_{0}^{1}{\it\gamma}(s^{\prime },t)\mathscr{G}({\it\zeta}(s,t),{\it\zeta}(s^{\prime },t))\,\text{d}s^{\prime }\!+\underbrace{\displaystyle \frac{1}{2\text{i}{\rm\pi}}\int _{1}^{\infty }{\it\gamma}(s^{\prime },t)\mathscr{G}({\it\zeta}(s,t),{\it\zeta}(s^{\prime },t))\,\text{d}s^{\prime }}_{w(s,t;\infty )},\end{eqnarray}$$
with
$\unicode[STIX]{x2A0D}$
being the principal Cauchy integral. Assuming the displacement of the plate and its wake is small of the order of
${\it\epsilon}$
, the position of the plate and its wake to
$O({\it\epsilon})$
is,
$$\begin{eqnarray}{\it\zeta}(s,t)=x(s)+\text{i}\,dh,\end{eqnarray}$$
and
$\mathscr{G}({\it\zeta}(s,t),{\it\zeta}(s^{\prime },t))$
to order
$O(1)$
can be written as,
$$\begin{eqnarray}\displaystyle \mathscr{G}(x,x^{\prime }) & = & \displaystyle \frac{{\rm\pi}}{2h}\left(\coth \left[\frac{{\rm\pi}}{2h}(x-x^{\prime })\right]-\coth \left[\frac{{\rm\pi}}{2h}(x-x^{\prime }+2\text{i}\,dh)\right]\right)\nonumber\\ \displaystyle & = & \displaystyle \frac{{\rm\pi}}{2h}\coth \left(\frac{{\rm\pi}}{2h}(x-x^{\prime })\right)-\frac{\displaystyle \frac{{\rm\pi}}{2h}\sinh \left(\displaystyle \frac{{\rm\pi}}{h}(x-x^{\prime })\right)}{\cosh \left(\displaystyle \frac{{\rm\pi}}{h}(x-x^{\prime })\right)-\cos (2{\rm\pi}\,d)}\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}\frac{{\rm\pi}}{4h}\frac{\sin (2{\rm\pi}d)}{\sinh ^{2}\left(\displaystyle \frac{{\rm\pi}}{2h}(x-x^{\prime })\right)+\sin ^{2}({\rm\pi}\,d)}.\end{eqnarray}$$
Defining
$G(x,x^{\prime })=\mathbb{R}(\mathscr{G}(x,x^{\prime }))$
as the real part of
$\mathscr{G}$
, we find (2.3),
$$\begin{eqnarray}G(x,x^{\prime })=\frac{{\it\alpha}}{2}\coth \left(\frac{{\it\alpha}}{2}r\right)-\frac{\displaystyle \frac{{\it\alpha}}{2}\sinh ({\it\alpha}r)}{\cosh ({\it\alpha}r)-\cos (2{\rm\pi}\,d)},\quad \text{with }r=x-x^{\prime },\text{and }{\it\alpha}=\frac{{\rm\pi}}{h}.\end{eqnarray}$$
Appendix B. Power expansion of
$G(x,x^{\prime })$
Knowing that the point on the plate is always upstream of the points in the wake and choosing
${\it\sigma}$
as a real positive number, for any pairs of complex points on the plate (
$z$
) and in the wake (
$z^{\prime }$
) we can write (Greengard Reference Greengard1990),
$$\begin{eqnarray}\displaystyle \coth ({\it\sigma}(z-z^{\prime })) & = & \displaystyle \frac{\text{e}^{{\it\sigma}(z-z^{\prime })}+\text{e}^{-{\it\sigma}(z-z^{\prime })}}{\text{e}^{{\it\sigma}(z-z^{\prime })}-\text{e}^{-{\it\sigma}(z-z^{\prime })}}=1-\frac{2}{1-\text{e}^{2{\it\sigma}(z-z^{\prime })}}\nonumber\\ \displaystyle & = & \displaystyle -1-2\mathop{\sum }_{j=1}^{\infty }\text{e}^{2j{\it\sigma}(z-z^{\prime })},\end{eqnarray}$$
provided that
$\text{e}^{2{\it\sigma}(z-z^{\prime })}<1$
, which is satisfied in the current case as the wake points are always downstream of the plate points and
${\it\sigma}=({\rm\pi}/2h)>0$
. Substituting (B 1) into
$G$
expression in (A 8) we obtain,
$$\begin{eqnarray}\mathscr{G}(z,z^{\prime })=\frac{{\rm\pi}}{h}\mathop{\sum }_{j=1}^{\infty }(\text{e}^{-({\rm\pi}/h)z^{\prime }j}-\text{e}^{-({\rm\pi}/h)\overline{z^{\prime }}j})\text{e}^{({\rm\pi}/h)zj},\end{eqnarray}$$
and,
$$\begin{eqnarray}G(x,x^{\prime })=\mathbb{R}(\mathscr{G}(x,x^{\prime }))=\frac{2{\rm\pi}}{h}\mathop{\sum }_{j=1}^{\infty }\text{e}^{({\rm\pi}/h)j(x-x^{\prime })}\sin ^{2}({\rm\pi}j\,d).\end{eqnarray}$$
By truncating
$G(x,x^{\prime })$
after
$J_{M}$
terms in the series, it is possible to show that the error in the expansion is bounded by,
$$\begin{eqnarray}\left|G(x,x^{\prime })-\frac{2{\rm\pi}}{h}\mathop{\sum }_{j=1}^{J_{M}}\text{e}^{({\rm\pi}/h)j(x-x^{\prime })}\sin ^{2}({\rm\pi}j\,d)\right|\leqslant \frac{2{\rm\pi}}{h}\frac{{\it\eta}^{J_{M}+1}}{1-{\it\eta}},\end{eqnarray}$$
where
${\it\eta}=\text{e}^{(-{\rm\pi}/h)(x^{\prime }-x)}$
.
Appendix C. Power expansion of the wake contribution
$(W(x,{\it\omega};\infty ))$
in the neutral stable case
We define
$\mathscr{F}(x,x^{\prime })$
as
$$\begin{eqnarray}W(x,{\it\omega};\infty )=-\frac{{\it\omega}{\it\Lambda}_{0}}{2{\rm\pi}}\int _{1}^{\infty }\underbrace{\text{e}^{-\text{i}{\it\omega}(x^{\prime }-1)}G(x,x^{\prime })}_{\mathscr{F}(x,x^{\prime })}\,\text{d}x^{\prime }.\end{eqnarray}$$
We can further decompose
$W(x,{\it\omega};\infty )$
into,
$$\begin{eqnarray}W(x,{\it\omega};\infty )=\underbrace{-\displaystyle \frac{{\it\omega}{\it\Lambda}_{0}}{2{\rm\pi}}\int _{1}^{x_{0}+x}\mathscr{F}(x,x^{\prime })\,\text{d}x^{\prime }}_{W(x,{\it\omega};x+x_{0})}\underbrace{\displaystyle -\frac{{\it\omega}{\it\Lambda}_{0}}{2{\rm\pi}}\int _{x_{0}+x}^{\infty }\mathscr{F}(x,x^{\prime })\,\text{d}x^{\prime }}_{\mathscr{M}(x,x_{0})},\end{eqnarray}$$
where
$x_{0}>1$
is an arbitrary splitting point in the wake. The first integral in the right-hand side is finite and can be numerically calculated without any difficulty. The second integral
$\mathscr{M}(x,x_{0})$
however is infinite but still can be solved analytically using the shifting technique (Greengard Reference Greengard1990). In particular, by incorporating (B 3) in the
$\mathscr{M}(x,x_{0})$
and changing the limits of integral to 0 and
$\infty$
, we obtain the series representation of
$\mathscr{M}(x,x_{0})$
as,
$$\begin{eqnarray}\mathscr{M}(x,x_{0})=-\frac{{\it\omega}{\it\Lambda}_{0}}{h}\text{e}^{-\text{i}{\it\omega}(x+x_{0}-1)}\mathop{\sum }_{j=1}^{\infty }\frac{\text{e}^{-({\rm\pi}/h)x_{0}j}}{\displaystyle \frac{{\rm\pi}}{h}j+{\it\omega}\text{i}}\,\sin ({\rm\pi}\,dj)^{2}.\end{eqnarray}$$
Similar to (B 4),
$\mathscr{M}(x,x_{0})$
is also truncated after
$J_{M}$
terms in the series,
$$\begin{eqnarray}\mathscr{M}(x,x_{0})=-\frac{{\it\omega}{\it\Lambda}_{0}}{h}\text{e}^{-\text{i}{\it\omega}(x+x_{0}-1)}\mathop{\sum }_{j=1}^{J_{M}}\frac{\text{e}^{-({\rm\pi}/h)x_{0}j}}{\displaystyle \frac{{\rm\pi}}{h}j+{\it\omega}\text{i}}\,\sin ({\rm\pi}\,dj)^{2}+{\it\delta}\mathscr{M}(x,x_{0}),\end{eqnarray}$$
and using the truncation error in the power expansion of
$G(x,x^{\prime })$
from (B 4), the minimum number of terms needed such that
${\it\delta}\mathscr{M}(x,x_{0})\leqslant {\it\sigma}$
for
$0\leqslant x\leqslant 1$
can be determined
a priori
as follows,
$$\begin{eqnarray}\displaystyle |{\it\delta}\mathscr{M}(x,x_{0})| & {\leqslant} & \displaystyle \left|\frac{{\it\omega}{\it\Lambda}_{0}}{h}\text{e}^{-\text{i}{\it\omega}(x+x_{0}-1)}\text{e}^{-({\rm\pi}/h)(J_{M}+1)x_{0}}\int _{0}^{\infty }\frac{\text{e}^{-(\text{i}{\it\omega}+(J_{M}+1)({\rm\pi}/h))x^{\prime }}}{1-\text{e}^{-({\rm\pi}/h)(x_{o}+x^{\prime })}}\,\text{d}x^{\prime }\right|\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \left|\frac{{\it\omega}{\it\Lambda}_{0}}{h}\right|\left|\frac{\text{e}^{-({\rm\pi}/h)x_{0}}}{1-\text{e}^{-({\rm\pi}/h)x_{0}}}\right|\frac{1}{\displaystyle (J_{M}+1)\frac{{\rm\pi}}{h}}\text{e}^{-({\rm\pi}/h)J_{M}x_{0}}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \underbrace{\left|\displaystyle \frac{{\it\omega}{\it\Lambda}_{0}}{{\rm\pi}}\right|}_{{\it\kappa}}\frac{1}{\left(\displaystyle \frac{{\rm\pi}}{h}J_{M}x_{0}\right)^{2}}\end{eqnarray}$$
and therefore,
$$\begin{eqnarray}J_{M}\geqslant \frac{h}{{\rm\pi}x_{0}}\sqrt{\frac{{\it\kappa}}{{\it\sigma}}}.\end{eqnarray}$$
