Simulation of the unsteady vortical flow of freely falling plates

Abstract

An inviscid vortex shedding model is numerically extended to simulate falling flat plates. The body and vortices separated from the edge of the body are described by vortex sheets. The vortex shedding model has computational limitations when the angle of incidence is small and the free vortex sheet approaches the body closely. These problems are overcome by using numerical procedures such as a method for a near-singular integral and the suppression of vortex shedding at the plate edge. The model is applied to a falling plate of flow regimes of various Froude numbers. For $\text {Fr}=0.5$, the plate develops large-scale side-to-side oscillations. In the case of $\text {Fr}=1$, the plate motion is a combination of side-to-side oscillations and tumbling and is identified as a chaotic type. For $\text {Fr}=1.5$, the plate develops to autorotating motion. Comparisons with previous experimental results show good agreement for the falling pattern. The dependence of change in the vortex structure on the Froude number and its relation with the plate motion is also examined.

Graphical abstract

1 Introduction

Solid objects that fall freely in fluids such as falling paper, coins, leaves, and seeds are commonly observed in nature. Maxwell [1] provided insightful observations and explanations on the turnover of falling papers over 150 years ago. Since then, this subject has attracted the attention of researchers [2]. Free-body motion has been widely investigated in the context of unsteady aerodynamics [3, 4], biomechanics [5], and sedimentology [6].

Among various shapes of falling bodies [2, 7], rectangular plates and circular disks have been studied extensively. Smith [8] carried out the first systematic experiment on a freely falling rectangular plate and classified it into three falling patterns: steady falling, fluttering (periodic oscillation), and tumbling (autorotation). Tanabe and Kaneko [9] suggested a chaotic trajectory that transits between fluttering and tumbling motions from their model of a dynamical system. Field et al. [10] performed comprehensive measurements on circular disks and found chaotic motion in a range of the dimensionless moment of inertia.

Belmonte et al. [11] conducted an experiment for rectangular flat plates and identified a transition from fluttering to tumbling at a critical Froude number. (The Froude number plays essentially the same role as the dimensionless moment of inertia.) Subsequently, Mahadevan et al. [12] considered a free tumbling card with a high aspect ratio and used dimensional arguments to provide a scaling for the average rate of descent and the tumbling frequency. Many investigators conducted numerical simulations of the two-dimensional Navier–Stokes equations for cylinders with rounded edges [13], elliptical cylinders [3], flat plates [14], and flexible plates [15].

An inviscid vortex shedding model was used to simulate falling flat plates by Jones and Shelley [16]. In this model, both the solid body and separated vortices are described as vortex sheets [17]. The vortex shedding model was developed for a moving rigid body by Jones [18]. It has been extended to various problems of vortex-body interactions, including falling flexible sheets [19], flapping plates and flags [20, 21], hovering and insect flights [22,23,24,25], and fish-like swimming [26,27,28]. A different type of vortex shedding model based on discrete point vortices was also used for falling plates [29].

Although the vortex shedding model has been successfully applied to various problems, there are limitations in the numerical simulations. If the free vortex sheet approaches the body closely, the computation becomes unstable and stops. Moreover, when the angle of incidence is small, the condition of vortex shedding that moves away from the plate edges is violated, and this results in a breakdown of the computation. Due to these problems, computations of the model for falling flat plates [16] stopped relatively early, demonstrating only the starting phase of fluttering motion. The plate may change its falling pattern from fluttering to tumbling or chaotic motion, or is reversed at later times. In this study, we overcome the limitations of the vortex shedding model by adopting several numerical procedures and apply them to conduct stable and accurate computations of falling plates.

Recently, Nitsche [30] reported that when a vortex sheet is located closely to the body, the integral of the velocity field of the model is a near-singular type, and demonstrated that a standard numerical integration for the velocity field is inaccurate. The author presented a numerical method for a near-singular integral for accurate calculation. We adopt this method for a near-singular integral in the falling-plate problem. In addition, if the shedding condition is violated, we suppress vortex shedding at the plate edge, detaching the existing free vortex sheet from the edge. The numerical method of the model is improved by using these procedures, and the computation proceeds for a significantly longer time than in the previous study [16]. This enables to investigate the unsteady dynamics of the vortex-body interaction of the falling plate such as the falling pattern of the plate and detailed vorticity field, for various Froude numbers.

In Sect. 2, we describe the inviscid vortex shedding model for a falling plate. In Sect. 3, the numerical method of the model is presented. The numerical results for falling plates for the flow regimes of various Froude numbers are presented in Sect. 4. Section 5 gives conclusions.

2 Model description

We consider a flat plate with chord length 2L, thickness T, width W, and density $\rho _b$. It moves under a gravitational field of strength g in an incompressible fluid of density $\rho _f$. We focus on a thin plate of a high aspect ratio, $T \ll L \ll W$, and thereby the motion of the plate can be regarded as approximately two-dimensional.

The flow is nondimensionalized by using the half-chord length L, fluid density $\rho _f$ and undetermined characteristic velocity U. The balance of the gravitational force on the plate, $2LTW (\rho _b - \rho _f)g$, and the magnitude of the aerodynamic force on the plate, of the order $2LW \rho _f U^2$, provides an estimate of the characteristic velocity U. By identifying $U^2 = Tg(\rho _b - \rho _f)/\rho _f$, the only remaining non-dimensional parameters are the Reynolds number and Froude number. We define the Froude number as

$$\begin{aligned} \textrm{Fr} = \frac{\rho _b U^2}{(\rho _b - \rho _f)gL} = \frac{m_b}{2 \rho _f L^2 W } = \frac{\rho _b T}{\rho _f L}, \end{aligned}$$

(1)

where $m_b$ denotes the mass of the plate. As we consider an effectively inviscid fluid, the Reynolds number is assumed to be large.

The Froude number has various physical meanings, as indicated by Eq. (1). In the context of a falling plate, it would be natural to interpret it as the ratio of the mass of the plate over the mass of the fluid column immediately surrounding the plate [16]. The Froude number thus characterizes the relative importance of plate inertia: a relatively heavy plate ($\textrm{Fr} \gg 1$) is less affected by the motion of the fluid, while a relatively light plate ($\textrm{Fr} \ll 1$) is highly dependent on the surrounding fluid motion.

In this section, we describe the equations for a flow with two free vortex sheets, i.e., one leading and one trailing edge vortex. The model is extended to a flow with multiple free vortex sheets at Sect. 3.3.

2.1 The body motion

The plate is parameterized by the arc-length s, $-L \le s \le L$, and expressed as

$$\begin{aligned} \zeta (s,t) = c(t) + s e^{\text {i} \theta }, \end{aligned}$$

(2)

where c(t) denotes the location of the center of mass of the plate and $\theta (t)$ denotes the orientation of the plate with the horizon in the laboratory frame.

We assume that the free vortex sheets are separated from the plate edges. The plate is also regarded as a bound vortex sheet. Recall that a vortex sheet is a surface across which the tangential velocity is discontinuous. The free vortex sheets are denoted by $\zeta _{\pm } (\Gamma , t)$ where $\Gamma$ represents the circulation as a Lagrangian parameter, which is explained shortly. $\Gamma _+(t)$ and $\Gamma _-(t)$ denote the total circulations of the free vortex sheets emanating from the (initial) right and left edges of the plate, respectively.

We introduce the body frame of reference which has its origin at the center of mass of the body. See Fig. 1 for a schematic. The plate in the laboratory frame $\zeta (s,t)$ is related to the body frame parametrization $s = (\zeta (s,t) - c(t)) e^{-\text {i} \theta }$. We denote $z_{\pm } (\Gamma ,t)$ as the position of the free vortex sheets in the body frame, where $z_{\pm } (\Gamma , t)= (\zeta _{\pm } (s,t) - c(t)) e^{-\text {i} \theta }$.

Fig. 1

Schematic of the body and free vortex sheets in a the laboratory frame and b body frame centered at c(t)

The location of the center of mass c(t) and angle $\theta (t)$ of the plate are governed by the linear and angular momentum balance equations, in non-dimensional form:

$$\begin{aligned} \textrm{Fr} \, \ddot{c}(t)&= \frac{1}{2} \, {\mathbb {F}}(t) - \textrm{i}, \end{aligned}$$

(3)

$$\begin{aligned} \frac{2}{3} \, \textrm{Fr} \, \ddot{\theta }(t)&= {\mathbb {T}} (t), \end{aligned}$$

(4)

where ${\mathbb {F}}(t)$ and ${\mathbb {T}} (t)$ denote the force and torque acting on the body, respectively. The relevant moment of inertia of the body is taken to be $2 \rho _b L^3 TW/3$. The dot represents the time derivative.

2.2 Equations for free vortex sheets

The evolution of the free vortex sheets is described by the Birkhoff-Rott equation [31]:

$$\begin{aligned} \frac{\partial z_{\pm }}{\partial t} (\Gamma ,t) = \bar{w} (z_{\pm } (\Gamma ,t), t) - \dot{c}(t) e^{-\text {i} \theta } - \textrm{i} \dot{\theta }(t) z_{\pm } (\Gamma ,t), \end{aligned}$$

(5)

where $\bar{w} (z,t)$ represents the velocity field in the body frame of reference, and the bar represents the complex conjugate. We define the vortex sheet strength $\gamma (s,t)$ as the jump in the tangential velocity across the bound vortex sheet. The circulation $\Gamma (s,t)$ is the integral of the vortex sheet strength over the free vortex sheet. According to the Helmholtz law, the circulation $\Gamma$ at material points on the free vortex sheet is conserved and does not change with time. We thus can label material points on the free vortex sheet by $\Gamma$. The velocity w(z, t) is then expressed as the boundary integral

$$\begin{aligned} w (z,t) = \frac{1}{2\pi \text {i}} \int _{-1}^1 \frac{\gamma (\lambda ,t) }{\lambda - z} \textrm{d} \lambda - \frac{1}{2\pi \text {i}} \left[ \int _0^{\Gamma _{\pm } (t)} K_{\delta } ( z_{\pm } (\Lambda ,t) - z ) \textrm{d} \Lambda \right] _{-}^{+}, \end{aligned}$$

(6)

where $K_{\delta } (z)$ is a regularized kernel for the singular one $K(z) = 1/z$, which is explained in detail in Sect. 3, and $\Lambda$ is the integration variable for the circulation. The square bracket denotes a difference, $[\alpha _{\pm }]_-^+ = \alpha _+ - \alpha _-$. The first integral in Eq. (6) represents the contribution of the bound vortex sheet, and the second integral represents the contribution of the free vortex sheets. Using the circulation to label points on the free vortex sheet prevents a separate evolution equation for $\gamma (s,t)$ on the free sheet. The bound vortex sheet strength $\gamma (s,t)$ and edge circulations $\Gamma _{\pm } (t)$ are determined as part of the solution.

2.3 Equations for the edge circulations

To find the bound vortex sheet strength and edge circulations, we apply the kinematic condition and unsteady Kutta condition. The unsteady Kutta condition is imposed to ensure that the fluid velocity at the edges of the bound vortex sheet remains bounded [18, 26].

The kinematic condition means that fluid does not penetrate the plate on either side. That is to say, the normal component of the fluid velocity on either side of the body should be the same as the body’s normal velocity. From Jones [26], the kinematic condition gives the following equation:

$$\begin{aligned} \frac{1}{\pi } \int _{-1}^1 \frac{\gamma (\lambda ,t) }{\lambda - s} \textrm{d} \lambda = f(s,t), \qquad \text{ for }~~ -1 \le s \le 1, \end{aligned}$$

(7)

where the right hand side takes the form

$$\begin{aligned} f(s,t) = 2 \nu (s,t) + \frac{1}{\pi } \left[ \int _0^{\Gamma _{\pm } (t)} \text{ Re } \left\{ K_\delta ( z_{\pm } (\Lambda ,t) - s ) \right\} \textrm{d} \Lambda \right] _{-}^{+}. \end{aligned}$$

(8)

Here $\nu (s,t)$ is the normal component of the velocity of the plate and is given by

$$\begin{aligned} \nu (s,t) = \text{ Im } \{ \dot{c}(t) e^{-\text {i} \theta } \} + \dot{\theta }(t) s. \end{aligned}$$

(9)

The integral Eq. (7) can be solved analytically [32]; however, the solution is not unique. By imposing the Kutta condition,

$$\begin{aligned} \gamma (\pm 1,t) = 0, \end{aligned}$$

(10)

we obtain the bounded solution as

$$\begin{aligned} \gamma (s,t) = -\frac{1}{\pi } \int _{-1}^1 \sqrt{\frac{1-s^2}{1-\lambda ^2}} \frac{f(\lambda , t)}{\lambda -s} \textrm{d} \lambda , \end{aligned}$$

(11)

with a constraint

$$\begin{aligned} \int _{-1}^1 \frac{f(\lambda , t)}{\sqrt{1-\lambda ^2}} \textrm{d} \lambda = 0. \end{aligned}$$

(12)

The constraint (12) is rewritten as

$$\begin{aligned} \nu (0,t) = -\frac{1}{2 \pi ^2} \int _{-1}^{1} \frac{1}{\sqrt{1 - \lambda ^2} } \left[ \int _0^{\Gamma _{\pm } (t)} \text{ Re } \left\{ K_\delta ( z_{\pm } (\Lambda ,t) - \lambda ) \right\} \textrm{d} \Lambda \right] _{-}^{+} \textrm{d} \lambda . \end{aligned}$$

(13)

This means that the normal velocity of the plate’s center of mass should be the same as the integral of the contribution of free vortex sheets on the mean normal fluid velocity over the plate. These equations are combined with the conservation of the total circulation of the flow,

$$\begin{aligned} \int _{-1}^1 \gamma (s,t) \textrm{d}s = \Gamma _+ (t) - \Gamma _- (t). \end{aligned}$$

(14)

Then, Eqs. (12) and (14) give a system of equations for the unknowns $\Gamma _{\pm } (t)$. The bound vortex sheet strength $\gamma (s,t)$ is obtained from Eq. (11), once $\Gamma _{\pm } (t)$ are determined.

2.4 Interaction of the body and fluid

To close the equations of the body motion (3) and (4), we need to calculate the force and torque. The only force exerted on the body is the normal force that acts with the pressure difference across the two sides of the body, since the fluid is inviscid. (There is an additional force when vortex shedding is suppressed at the plate edge, which is explained in Sect. 3.3.) Thus, we have

$$\begin{aligned} \textrm{Fr} \, \ddot{c}(t)&= - \frac{1}{2} \text {i} e^{\text {i} \theta (t)} \int _{-1}^1 [p (\lambda ,t)] \, \textrm{d} \lambda - \text {i} , \end{aligned}$$

(15)

$$\begin{aligned} \frac{2}{3} \textrm{Fr} \, \ddot{\theta }(t)&= - \int _{-1}^1 \lambda \, [p (\lambda ,t)] \textrm{d} \lambda . \end{aligned}$$

(16)

The pressure difference across the two sides of the bound vortex sheet is given by

$$\begin{aligned} {[}p (s,t)] = - \frac{\textrm{d} \Gamma }{\textrm{d}t} (s,t) - \gamma (s,t) (\mu (s,t) - \tau (t)). \end{aligned}$$

(17)

Here $\mu (s,t)$ denotes the tangential component of the average of the fluid velocities on either side of the body, which has the form

$$\begin{aligned} \mu (s,t) = - \frac{1}{2 \pi } \left[ \int _0^{\Gamma _{\pm } (t)} \text{ Im } \left\{ K_\delta ( z_{\pm } (\Lambda ,t) - s ) \right\} \textrm{d} \Lambda \right] _{-}^{+}, \end{aligned}$$

(18)

and $\tau (t) = \text{ Re } \{ \dot{c}(t) e^{-\text {i} \theta } \}$ is the tangential component of the body velocity. The circulation on the body is expressed as

$$\begin{aligned} \Gamma (s,t) = \Gamma _{-} (t) + \int _{-1}^s \gamma (\lambda ,t) \textrm{d} \lambda . \end{aligned}$$

(19)

Differentiating of this equation with respect to time gives

$$\begin{aligned} \frac{\textrm{d} \Gamma }{\textrm{d}t} (s,t) = \dot{\Gamma }_{-} (t) + \int _{-1}^s \dot{\gamma } (\lambda ,t) \textrm{d} \lambda . \end{aligned}$$

(20)

To determine $\ddot{c}(t)$ and $\ddot{\theta }(t)$, we have to find the time derivatives of the edge circulations and sheet strength on the body. Several approaches can be applied to find these quantities [16, 23, 27]. Here, we use the solution (7) and constraint (12). Differentiating Eqs. (7) and (12) with respect to time, we have

$$\begin{aligned} \gamma _t (s,t) =&-\frac{1}{\pi } \int _{-1}^1 \sqrt{\frac{1-s^2}{1-\lambda ^2}} \frac{f_t (\lambda , t)}{\lambda -s} \textrm{d} \lambda , \end{aligned}$$

(21)

$$\begin{aligned}&\int _{-1}^1 \frac{f_t (\lambda , t)}{\sqrt{1-\lambda ^2}} \textrm{d} \lambda = 0. \end{aligned}$$

(22)

In Eqs. (21) and (22), $f_t (s,t)$ is obtained as

$$\begin{aligned} f_t(s,t) =&~ 2 \nu _t (s,t) - \frac{1}{\pi } \left[ \int _0^{\Gamma _{\pm } (t)} \text{ Re } \left\{ L_\delta ( z_{\pm } (\Lambda ,t) - s ) \right\} \textrm{d} \Lambda \right] _{-}^{+} ~~~ \nonumber \\&+ \frac{1}{\pi } \left[ \dot{\Gamma }_{\pm } (t) \, \text{ Re } \left\{ K_\delta ( z_{\pm } (\Gamma _{\pm },t) - s ) \right\} \right] _{-}^{+}, \end{aligned}$$

(23)

by applying the Leibniz integral rule, where $L_\delta (z)$ is the time derivative of the regularized kernel $K_\delta (z)$. The differentiation of the circulation conservation gives

$$\begin{aligned} \int _{-1}^1 \gamma _t (s,t) \textrm{d}s = \dot{\Gamma }_+ (t) - \dot{\Gamma }_- (t). \end{aligned}$$

(24)

Then, Eqs. (15), (16), (21), (22), and (24) provide a system of equations for the unknowns $\ddot{c}(t)$, $\ddot{\theta }(t)$, $\dot{\Gamma }_{\pm } (t)$, and $\dot{\gamma } (s,t)$.

3 Numerical method

In this section, we present the numerical procedures required for stable and accurate computation of a falling plate. The discretization and time-integration of the model are described in the Appendix.

3.1 $\delta$-regularization

The vortex sheet evolution suffers from the Kelvin-Helmholtz instability for all disturbance wavenumbers [31]. It is well known that in a free shear flow, the vortex sheet develops a singularity in finite time [34]. To circumvent the ill-posedness, the vortex blob method in which the singular kernel $K(z)=1/z$ is replaced by a smoothed kernel [35, 36] is commonly applied. A widely used blob-regularization is to give a constant parameter $\delta =\delta _0$ in the kernel,

$$\begin{aligned} K_{\delta } (z) = \frac{\bar{z}}{z \bar{z} + \delta ^2}. \end{aligned}$$

(25)

Note that the blob-regularization should not be applied to the bound vortex sheet, because such regularization renders the equation ill-posed. Thus, there is a jump in the regularization parameter at the edges of the bound vortex sheet when the uniform regularization is applied. The jump of the regularization parameter often yields oscillations in the solution at an early time. To prevent oscillations, a small value of $\delta _0$ should be used; however, the computation becomes expensive because of the high resolution of spiral cores.

Alben [27] proposed a tapered smoothing using nonuniform $\delta (s)$, to reconcile the regularization of the free vortex sheet with the absence of regularization on the bound vortex sheet. In this approach, $\delta (s)$ gradually decreases as the arc-length s approaches the edges of the bound vortex sheet. For the free vortex sheet separated from the edge of the body, $\delta (s)$ is given by

$$\begin{aligned} \delta (s) = \delta _0 [1 - (1-\tau ) e^{- \tilde{s}^2 (\Gamma , t)/\epsilon ^2} ], \end{aligned}$$

(26)

where $\tilde{s} (\Gamma , t)$ is the arc-length of the point $z_{\pm } (\Gamma ,t)$ from the plate edge. The parameter $\tau$ is introduced to prevent vanishing of $\delta$ in the limit, and $\delta (s)$ decreases to $\delta _0 \tau$ as $\tilde{s} \rightarrow 0$. With the nonuniform smoothing, $L_\delta (z)$ in Eq. (23) takes the form

$$\begin{aligned} L_\delta (z)&= \frac{ \dot{\bar{z}} - 2 K_\delta (z) [ \text{ Re } \{ \bar{z} \dot{z} \} + \delta \dot{\delta } ] }{ |z|^2 + \delta ^2 }. \end{aligned}$$

(27)

Here, $\delta = \delta (s,t)$, since $\delta$ in Eq. (26) is also a function of t. In our simulations, the parameters are set to $\epsilon = 2 \delta _0$ and $\tau =0.2$, with $\delta _0 = 0.2$.

3.2 Near-singular integral

The vortex sheet model encounters a computational problem when the free vortex sheet is located near a body. Jones [18] commented on this problem, demonstrating the breakdown of the computation for an example of the close interaction between a body and free vortex sheets. Several authors also reported difficulties in near-boundary simulations [24, 37, 38]. In the falling plate problem, the near-boundary situation occurs frequently and should therefore be overcome to conduct long-time simulations.

The problem in the near-boundary simulation arises from the boundary integral over a body to evaluate the velocity field. That is to say, the first integral in Eq. (6), becomes near-singular and loses accuracy. Recently, Nitsche [30] proposed a numerical method to accurately evaluate near-singular integrals. This method is based on the Taylor-series approximations of the integrands that capture the near-singular behavior and are integrated in closed forms. We apply this method to a falling plate and summarize it briefly.

Using the transformation $s = \cos \alpha$, $0 \le \alpha \le \pi$, the boundary integral over the plate is converted to

$$\begin{aligned} \int _{-1}^1 \frac{\gamma (\lambda ,t)}{\lambda - z} \textrm{d} \lambda = \int _0^{\pi } \frac{ \lambda (\alpha ) - \bar{z} }{ | \lambda (\alpha ) - z |^2} \gamma (\alpha ,t) \sin \alpha \, \textrm{d} \alpha . \end{aligned}$$

(28)

We rewrite this integral in the vector form,

$$\begin{aligned} \int _0^\pi \frac{ \textbf{x} ( \alpha ) - \bar{\textbf{x}}_0 }{ | \textbf{x} ( \alpha ) - \textbf{x}_0 |^2} \omega (\alpha ) \, \textrm{d} \alpha , \end{aligned}$$

(29)

by denoting $\textbf{x} (\alpha ) = (\lambda (\alpha ), 0)$ and $\omega (\alpha ) = \gamma (\alpha ,t) \sin \alpha$ and changing $z \rightarrow \textbf{x}_0$. The flat plate is discretized uniformly with respect to $\alpha$ by the points $\alpha _j = jh$, where h is the mesh size. The distance from $\textbf{x}_0$ to the plate is denoted by d. The main idea is to approximate the integrand to leading-order for $d \ll h$ by a function $H (\alpha )$ which captures the near-singularity and can be integrated exactly. Denoting the integrand of (29) by G, we approximate the integral of G as follows:

$$\begin{aligned} \int G&= \int (G - H + H) \end{aligned}$$

(30a)

$$\begin{aligned}&\approx T [G-H] + \int H \end{aligned}$$

(30b)

$$\begin{aligned}&= T [G] + \left( \int H - T [H] \right) \end{aligned}$$

(30c)

$$\begin{aligned}&= T [G] + E [H]. \end{aligned}$$

(30d)

By adding and subtracting the approximation from G, the smoother difference $G-H$ is approximated by a trapezoid rule T, while integrating H exactly. The inaccurate evaluation of the trapezoidal rule for G is improved by applying the same approximation to H. The rearrangement of the terms yields the final expression (30d), which adds the correction $E[H] = \int H - T [H]$ to the trapezoidal approximation of the uncorrected integral.

For T, we use the second-order trapezoidal rule. The approximation H consists of a Taylor-series approximation of $G(\alpha )$ about $\alpha _p$, where $\textbf{x}_p = \textbf{x} (\alpha _p)$ is the projection of $\textbf{x}_0$ onto the plate, i.e., $\textbf{x} (\alpha _p) = ( \lambda (\alpha _p), 0)$. The expansion of the numerator and denominator of G and its approximation H are given by

$$\begin{aligned} G(\alpha )&= \frac{\textbf{x} ( \alpha ) - \bar{\textbf{x}}_0}{| \textbf{x} (\alpha ) - \textbf{x}_0 |^2} \omega (\alpha ) = \frac{c_0 + c_1 \hat{\alpha } + c_2 \hat{\alpha }^2 + O (\hat{\alpha }^3)}{d + c^2 \hat{\alpha }^2 + e_3 \hat{\alpha }^3 + O (\hat{\alpha }^4)} \end{aligned}$$

(31a)

$$\begin{aligned} \approx H(\alpha )&= \frac{c_0 + c_1 \hat{\alpha } + c_2 \hat{\alpha }^2}{d + c^2 \hat{\alpha }^2} - e_3 \frac{c_0 \hat{\alpha }^3 + c_1 \hat{\alpha }^4}{(d + c^2 \hat{\alpha }^2)^2}, \end{aligned}$$

(31b)

where $\hat{\alpha } = \alpha - \alpha _p$. The coefficients of the denominator in Eq. (31a) are

$$\begin{aligned} c^2&= | \textbf{x}_p^{\prime } |^2 + \textbf{x}_p^{\prime \prime } \cdot ( \textbf{x}_p - \textbf{x}_0 ), \end{aligned}$$

(32a)

$$\begin{aligned} e_3&= \textbf{x}_p^{\prime } \cdot \textbf{x}_p^{\prime \prime } + \frac{1}{3} ( \textbf{x}_p - \textbf{x}_0 ) \cdot \textbf{x}_p^{\prime \prime \prime }, \end{aligned}$$

(32b)

where the prime denotes the derivative with respect to $\alpha$. The second terms on the right hand sides of Eq. (32) are zero in the case of a flat plate, owing to orthogonality. Further applying $\lambda = \cos \alpha$, the coefficients are expressed as $c^2 = \sin ^2 \alpha _p$ and $e_3 = \sin \alpha _p \cos \alpha _p$. The functions of $H (\alpha )$ can be integrated analytically in closed forms, and the correction E[H] then can be computed. The approximation $G \approx H$ yields a method of order $O(h^2)$, with the property that the error vanishes as $d \rightarrow 0$.

It is necessary to evaluate the error E[H] only over an interval centered on $\alpha _p$,

$$\begin{aligned} I = [ \alpha _{j_0 - n_w}, \alpha _{j_0 + n_w} ], \end{aligned}$$

(33)

where $\alpha _{j_0}$ is the grid point closest to $\alpha _p$. The number of intervals for the method is set to $n_w = 20$, in our computations. The numerical implementation of this method is as follows:

Step 1: Compute T[G].

Step 2: If $\textbf{x}_0$ is too close to the plate (here, $d < 4 \Delta s$),

• (i) find $\alpha _p, d, \textbf{x}_p, \textbf{x}_p^{\prime }, \textbf{x}_p^{\prime \prime }$ and $\omega ' (\alpha _p)$.

• (ii) find the coefficients $c_0, c_1, c_2, c$, and $e_3$.

• (iii) compute $E[H]= \int H - T [H]$ over the interval $I = [ \alpha _{j_0 - n_w}, \alpha _{j_0 + n_w} ]$.

• (iv) correct $\int G \approx T[G] + E[H]$.

Here, $\Delta s$ is the spacing between the two grid points closest to $\textbf{x}_p$. As $c_0 = O(d)$, this method requires only the first-order derivative of $\omega (\alpha )$. For $\omega '(\alpha _p)$, $\gamma '(\alpha _p)$ is needed and can be provided by interpolating $\gamma (\alpha _j)$ on the body.

3.3 Suppression of vortex shedding

In addition to the near-singular integral, the model has another limitation. The numerical procedure of vortex shedding is valid if new vortex particles are separated tangentially from the body edge. In other words,

$$\begin{aligned} \mu _+ (t) - \tau (t) > 0 \qquad \text{ and } \qquad \mu _- (t) - \tau (t) < 0, \end{aligned}$$

(34)

where $\mu _{\pm } (t)$ denotes the value of $\mu (t)$ at the plate edges. This condition does not hold in certain situations such as small angles of attack. To handle this, we do not shed vortices at the corresponding edge if one of the inequalities in (34) is not satisfied. Therefore, the free vortex sheet can shed from one, both, or neither of the plate edges during the computation. At the edges that do not shed vortices, the Kutta condition is not applied and a flow singularity is allowed. In this case, the solution and constraint of the integral equation for $\gamma (s)$, Eqs. (11) and (12), are given in different forms. We refer to Muskhelishvili [32] or Sohn [28] for the detailed expressions.

When vortex shedding is suppressed at the edge, the flow velocity and pressure diverge as the inverse square root of distance from the edge. The divergent pressure produces a finite suction force on the body, which is the integral of the pressure over an infinitesimal small circle at the plate edge [33]. The force in the linear momentum balance Eq. (3) has an additional term from the suction force and is given by

$$\begin{aligned} {\mathbb {F}} (t) = {\mathbb {F}}_{\text {n}} (t) \text {i} e^{\text {i} \theta (t)} \pm \frac{\pi }{8} \sigma ^2 (\pm 1, t) e^{\text {i} \theta (t)}, \end{aligned}$$

(35)

where $\sigma (s,t)$ denotes the regular part of the sheet strength, $\gamma (s,t) = \sigma (s,t) / \sqrt{1-s^2}$, and ${\mathbb {F}}_{\text {n}} (t)$ represents the normal force on the plate.

For the edge with no vortex shedding, a free vortex sheet is newly shed when the condition (34) is recovered to satisfy. Therefore, several free vortex sheets may exist in the flow, at late times. The model is easily extended to the flow with multiple free vortex sheets; two, one, or none free vortex sheets are attached to the bound vortex sheet, and the others are detached from the bound vortex sheet. In Eqs. (6), (8), (18), and (23), the integration on the free vortex sheet should be taken for each one. In Eq. (14) of the circulation conservation, the edge circulations are replaced by adding the circulations of all free vortex sheets separated from the corresponding edges, i.e., $\Gamma _+ \rightarrow \sum _k \Gamma _{+, k}$ and $\Gamma _- \rightarrow \sum _{k} \Gamma _{-, k}$. In Eqs. (23) and (24), the time derivatives of $\Gamma _+$ and $\Gamma _+$ are taken only on the edge circulations of the free vortex sheets attached to the bound vortex sheet, because the circulations of the detached free vortex sheets do not change with time.

4 Numerical results

We present numerical results for the falling plate by dividing the flow regime into $\text{ Fr }=1, \text{ Fr } < 1$, and $\text{ Fr } > 1$. The plate is released from rest with the initial angle $\theta = \pi /256$ for all results in this section. The plate is initially discretized by $N=83$.

4.1 $\text{ Fr } = 1$

We first consider the case of $\text{ Fr } = 1$. The falling plate and free vortex sheets for $\text{ Fr } = 1$ up to $t=15$ are shown in Fig. 2. Two small starting vortices are shed from the plate edges early, and the plate becomes unstable and its oscillations increase. The result of Fig. 2 is similar to that in Jones and Shelley [16], but the y-location of the plate in Fig. 2 is marginally lower than their result. This difference is due to the nonuniform regularization for the Birkhoff-Rott equation in our simulation, while they used the uniform regularization.

Fig. 2

Falling plate and free vortex sheets for $\text{ Fr } = 1$ up to $t=15$

The asymmetry of the initial condition immediately produces the normal force on the plate, which accelerates the plate to move to left. The right vortex is located slightly higher than the left vortex. Owing to the torque of this asymmetry, the plate rotates back to the horizontal position and passes the horizon by inertia. Then, the plate decelerates and the right vortex is located lower than the left vortex. This process of the vortex-body interaction continues in the opposite direction and is repeated.

The second phase of the falling plate and free vortex sheets for $\text{ Fr } = 1$ are shown in Fig. 3. The flow topology changes during this phase. The shedding condition (34) is violated at the leading edge, immediately after the plate starts to move to right at $t=15.06$. Thus, vortex shedding is suppressed at the leading edge, and the free vortex sheet is detached. In Fig. 3, the detached free vortex sheet is colored red, to distinguish it from the separating vortex sheet. The end of the detached vortex sheet rolls up on the plate and is washed down to the large structure of the trailing edge vortex. During the second phase, the plate swings to right and rotates rapidly after passing the local minimum. As the plate steepens, the free vortex sheet is newly shed from the leading edge, at $t=19.05$, recovering the shedding condition (34) to satisfy. The plate turns over at $t=19.26$, and immediately after that, the trailing edge vortex is detached. The new free vortex sheet, of short length, is colored green in Fig. 3, and the circle denotes the initial left edge of the plate.

Fig. 3

Falling plate and free vortex sheets for $\text{ Fr } = 1$ at $16 \le t \le 20$. The detached vortex sheet is colored red. At $t=20$, the newly shed vortex sheet at the leading edge appears and is colored green. The circle denotes the initial left edge of the plate (color figure online)

Fig. 4

Falling plate and free1 vortex sheets for $\text{ Fr } = 1$ at $21 \le t \le 29$. The circle denotes the initial left edge of the plate

Figure 4 shows the third phase of the falling plate for $\text{ Fr } = 1$. After the plate overturns at $t=19.26$, the rotation decelerates and the plate swings to right, where the roles of the trailing and leading edges are switched. The newly shed vortex (green) at the trailing edge is nearly attached to the primary vortices (red) at $t=21$. The plate turns over again at about $t=23.4$ and then swings back to left. A relatively strong trailing edge vortex is formed and the bound vortex sheet crosses over the free vortex sheets while traveling to left.

Fig. 5

Trajectory of a falling plate for $\text{ Fr } = 1$ from $t=0$ to $t=29$ in the increment of 0.5

The trajectory of the falling plate for $\text{ Fr } = 1$ from $t=0$ to $t=29$ in the increment of 0.5 is plotted in Fig. 5. The plate falls with side-to-side oscillations at the early stage. After the angle of the plate changes largely with respect to the left edge, the plate rotates end-over-end and extends to right. At the end of the swing, the plate turns over again and moves to left. The plate motion of $\text{ Fr } = 1$ is a combination of side-to-side oscillations and autorotation, and is therefore identified as a chaotic type.

Figure 6 shows the temporal variations of the x- and y-locations, and the angle of the plate and the total edge circulations for $\text{ Fr } = 1$. $c_x (t)$ and $c_y (t)$ denote the x- and y-components of the center of mass of the plate c(t), respectively. Interestingly, the x-location of the plate has a similarity with the angle $\theta (t)$. We observe that the total circulations of the initial right edge $\Gamma _+$ and the initial left edge $\Gamma _-$ are stationary at $15.47< t < 19.05$, and $t > 19.31$, respectively. The edge circulation $\Gamma _-$ increases significantly at $16< t < 17.5$, when the plate swings to right after detaching the leading edge vortex.

The similarity of the x-location of the plate with the plate angle would be explained by a simple pendulum. The early-time side-to-side oscillations behave like a stable pendulum, neglecting the y-location of the plate. Even when the plate rotates end-over-end, or rotates $180^\circ$ and swings back, the x-location of the plate increases or decreases, respectively, with the plate angle, which is analogous to an unstable pendulum. Kuznetsov [39] showed from a dynamical system approach that the dynamics of a falling body, under certain assumptions, can be expressed as the pendulum equation.

Fig. 6

x- and y-locations, and angle of the plate and total edge circulations with time for $\text{ Fr } = 1$

We have checked that the plate motion is insensitive to the initial condition. The direction of rotation depends on the initial angle, and the early-time transient of side-to-side oscillations is decreased as the initial angle increases. However, the falling style depends little on the initial plate angle. (This also holds for the other Froude number regimes.) This result of the independence on the initial condition from our model is consistent with previous studies of the falling plate [4, 9, 14].

4.2 $\text{ Fr } < 1$

We choose $\text{ Fr } = 0.5$ as a representative case for $\text{ Fr } < 1$. The flow regime $\text{ Fr } < 1$ was not considered in Jones and Shelley [16], because of the close interaction between vortices and the plate.

Figure 7 shows the falling plate and free vortex sheets for $\text{ Fr } = 0.5$ up to $t=11$. The plate falls similar to the case of $\text{ Fr } = 1$ at early times. The oscillations grow rapidly after the early stage and soon becomes unstable. The shedding condition is violated at the leading edge, and the free vortex sheet is detached at $t=11.37$, as the plate begins to move to right. The plate then becomes steep, while rotating clockwise with respect to the left edge. Subsequently, the plate swings largely to right, as shown in Fig. 8. The vortices fall down and become larger. They remain near $x=0$ and are located away from the plate.

Fig. 7

Falling plate and free vortex sheets for $\text{ Fr } = 0.5$ up to $t=11$

Fig. 8

Falling plate and free vortex sheets for $\text{ Fr } = 0.5$ at $12 \le t \le 17.5$

As the plate approaches the right end, it decelerates and rotates slightly, and then falls freely to left, as shown in Fig. 9. (See also Fig. 10.) We describe more details of the flow at $17.5< t < 19$. The plate begins to fall down at $t=17.89$. Immediately thereafter, a new free vortex sheet is shed at the trailing edge, which is colored green in Fig. 9, and the leading edge vortex is detached. No leading edge vortex is formed in the free-fall phase.

Fig. 9

Falling plate and free vortex sheets for $\text{ Fr } = 0.5$ at $19 \le t \le 23$

Fig. 10

Trajectory of a falling plate for $\text{ Fr } = 0.5$ from $t=0$ to $t=23$ in the increment of 0.5. Circles are not drawn as the plate is steepened at the right end, for a clear view

The trajectory of the falling plate for $\text{ Fr } = 0.5$ from $t=0$ to $t=23$ in the increment of 0.5 is plotted in Fig. 10. The plate for $\text{ Fr } = 0.5$ oscillates with a side-to-side motion. However, as the oscillation is amplified, it swings in a large scale. After reaching the right end, the plate accelerates and falls freely. Figure 10 shows the remarkable differences of the motion of the plate of $\text{ Fr }=0.5$ from that of $\text{ Fr } =1$. The plate motion for $\text{ Fr } = 0.5$ until reaching the left end is qualitatively similar to that for $\text{ Fr } =1$. However, after the left end, it persists side-to-side motion, whereas the plate for $\text{ Fr } =1$ is turned over. Figure 11 plots the x- and y-locations, and the angle of the plate and the total edge circulations for $\text{ Fr } = 0.5$. The circulation of the left edge $\Gamma _-$ increases only slightly after the detachment of the leading edge vortex at $t=11.37$, compared with the case of $\text{ Fr } =1$ in Fig. 6.

Fig. 11

x- and y-locations, and angle of the plate and total edge circulations with time for $\text{ Fr } = 0.5$

4.3 $\text{ Fr } > 1$

We take the two Froude numbers, $\text{ Fr } = 1.5$ and 2, as representative cases of $\text{ Fr } > 1$. The motion of the plate for $\text{ Fr } = 1.5$ at $t < 24$ is similar to that of $\text{ Fr } = 1$, differing only in time scale. Thus, we show the flow of $\text{ Fr } = 1.5$ only at late times, in Fig. 12. The plate of $\text{ Fr } = 1.5$ has already turned over at $t<24$, although not shown here. A key difference between the plate motions of the two Froude numbers is that the plate of $\text{ Fr } = 1.5$ continues to rotate and travels to right, whereas the plate of $\text{ Fr } = 1$ swings back to left at late times. Therefore, the plate for $\text{ Fr } = 1.5$ demonstrates tumbling motion.

In Fig. 12, we observe more small-scale vortices in the flow of $\text{ Fr } = 1.5$, compared to the flow of $\text{ Fr } = 1$. More importantly, no rolled-up trailing edge vortex appears at $t \ge 28$ for $\text{ Fr } = 1.5$, whereas a relatively strong trailing edge vortex evolves for $\text{ Fr } = 1$. This difference in the vortex structure may be responsible for $360^\circ$ rotation of the plate of $\text{ Fr } = 1.5$, and $180^\circ$ rotation and swinging back in the case of $\text{ Fr } = 1$.

Fig. 12

Falling plate and free vortex sheets for $\text{ Fr } = 1.5$ at late times. The circle denotes the initial left edge of the plate

Figure 13 shows the falling plate and free vortex sheets for $\text{ Fr } = 2$ up to $t=20$. The plate falls with side-to-side oscillations. The computation of $\text{ Fr } = 2$ in the previous study by Jones and Shelley [16] stopped at $t=21.19$. However, our computation proceeds until $t=25.16$, satisfying the shedding condition, as shown in Fig. 13. This is owing to the adoption of the nonuniform regularization and the method for a near-singular integral. In Fig. 13, we observe that the free vortex sheets are in close proximity of the plate, at $22 \le t \le 25$. This indicates that the near-singular integral method is successfully applied. The subsequent motions are plotted in Fig. 14. We see that the trailing edge vortex is detached. At $23 \le t \le 25$, the plate rotates with respect to the leading edge. (See also Fig. 15). Subsequently, the plate falls steeply and then changes its angle sharply, thereby making a circular turn. No leading edge vortex is formed during the free fall and circular turn.

Fig. 13

Falling plate and free vortex sheets for $\text{ Fr } = 2$ up to $t=25$

Fig. 14

Falling plate and free vortex sheets for $\text{ Fr } = 2$ at $26 \le t \le 32$

The trajectories of the falling plates of $\text{ Fr } = 1.5$ and 2 are shown in Fig. 15. Tumbling motion of the plate of $\text{ Fr } = 1.5$ is clearly observed. The plate makes $180^\circ$ turn first and then $360^\circ$ turn. For $\text{ Fr } = 2$, the side-to-side oscillations persist longer. The plate falls nearly vertically during the transition between the side-to-side oscillations and rotation. However, our computation is not long enough to determine the falling type of $\text{ Fr } = 2$. Note that the computational time of the vortex sheet model increases dramatically at late times, because of $O(N^2)$ computation cost of the Birkhoff-Rott equation. Figure 16 plots the x- and y-locations, and the angle of the plate and the total edge circulations for $\text{ Fr } = 2$. The x-location of the plate of $\text{ Fr } = 2$ also has a similarity to the angle $\theta$, as explained previously.

Fig. 15

Trajectory of falling plates for $\text{ Fr } = 1.5$ and 2 in the increment of 0.5. Times are $0 \le t \le 29$ for $\text{ Fr } = 1.5$, and $0 \le t \le 32$ for $\text{ Fr } = 2$

Fig. 16

x- and y-locations, and angle of the plate and total edge circulations for $\text{ Fr } = 2$

4.4 Comparison with experiments

For validation, we compare our results with previous experiment studies. From an experiment for thin rectangular flat plates, Smith [8] showed that a transition from fluttering to tumbling occurs at the dimensionless moment of inertia, $I^* \approx 0.2-0.3$, where $I^*$ is defined as

$$\begin{aligned} I^* = \frac{32 I}{\pi \rho _f D^4} = \frac{8 \rho _b}{3 \pi \rho _f} \left[ \frac{T}{D} + \left( \frac{T}{D} \right) ^3 \right] , \end{aligned}$$

(36)

and the chord length $D = 2L$. Chaotic motion was not identified in his study. Neglecting the third-order term for a thin plate, the dimensionless moment of inertia is expressed as $I^* \approx 4 \, \textrm{Fr} / (3 \pi )$. We have checked that in our model, the motion of the plate transits from fluttering to overturning at the critical Froude number, $\textrm{Fr}_c = 0.6$. This corresponds to $I^* \approx 0.26$ and is in good agreement with Smith [8].

From a quasi-two-dimensional experiment on rectangular plates, Belmonte et al. [11] reported the transition from fluttering to tumbling at the critical Froude number, $\textrm{Fr}_c = 0.9$, in our definition of the Froude number. (Chaotic motion was also not found in their study.) The critical value $\textrm{Fr}_c = 0.9$ is rather larger than that from our model. This difference may arise from their experimental setup where the flow tank is narrow and the plates touch the sides of the container to constrain their motion to be two-dimensional.

In addition, Anderson et al. [3] conducted essentially a two-dimensional experiment for rectangular plates, with different thickness-to-length ratios. They showed the trajectories of the plate to be fluttering for $I^* = 0.16$, chaotic for $I^* = 0.39$, and tumbling for $I^* = 0.48$. The Froude numbers in our simulations, $\textrm{Fr} = 0.5, 1$, and 1.5 correspond to $I^* \approx 0.21, 0.42$, and 0.64, and the plate motions appear to be fluttering, chaotic, and tumbling, respectively. Therefore, the results of our model are consistent with their experiment.

5 Conclusions

We have improved the numerical method of the inviscid vortex shedding model to investigate the unsteady motion of falling plates, which extends the previous study [16]. The main contribution of this work is to identify the falling pattern of the plate and the detailed vorticity field for various regimes of the Froude number. The limitations of the model are overcome by adopting a method for a near-singular integral and suppression of vortex shedding, and a nonuniform regularization. The numerical computation of the model successfully simulates the motion of falling plates at late times. In particular, the close interactions of the plate and vortices are calculated accurately by using the near-singular integral method.

Our results show remarkable differences of the falling pattern for varying the Froude number. The plate develops large-scale side-to-side oscillations for a small Froude number, whereas transverse oscillations of the plate decrease significantly for a large Froude number. We find that the plate of $\text{ Fr } = 0.5$ corresponds to fluttering motion, and $\text{ Fr } = 1$ is apparently chaotic motion, while the case of $\text{ Fr } = 1.5$ develops into autorotating motion. The transition from fluttering to overturning occurs at the critical Froude number, $\textrm{Fr}_c = 0.6$.

The vortex structure shed from the plate varies with the Froude number. Topological changes of the vortex structure with the Froude number are summarized as followings.

• (i) $\text{ Fr } = 0.5, 1, 1.5$: LEV & TEV $\Rightarrow$ LEV detached $\Rightarrow$ new separation & detachment,

• (ii) $\text{ Fr } = 2$: LEV & TEV $\Rightarrow$ TEV detached,

where the LEV and TEV denote the leading edge vortex and trailing edge vortex, respectively. In all the cases, the first phase of the coexistence of the LEV and TEV corresponds to side-to-side oscillations. For $\text{ Fr } = 0.5, 1$, and 1.5, the plate swings from one side to the other after the LEV is detached, whereas for $\text{ Fr } = 2$, the plates falls and rotates, after the TEV is detached. For the cases of $\text{ Fr } = 0.5, 1$, and 1.5, new separation and detachment of vortices occur within a short duration when the plate rotates, following a new phase of swing.

Although topological changes of the vortex structure for $\text{ Fr } = 0.5, 1$, and 1.5 are similar, the interactions of vortices with the plate are markedly different in details. For $\text{ Fr } = 0.5$, the primary vortices are located far from the plate at late times, and thus interact little with the plate. In the cases of $\text{ Fr } = 1$ and 1.5, vortices strongly interact with the plate. At late times, a relatively strong trailing edge vortex is developed for $\text{ Fr } = 1$, whereas no rolled-up trailing edge vortex is shed for $\text{ Fr } = 1.5$. This difference of the vortex structure between $\text{ Fr } = 1$ and 1.5 may determine the falling patterns. These observations indicate that the change of the vortex structure plays a significant role in the motion of the falling plate.

The three Froude number regimes exhibit different trends in the circulation. For $\text{ Fr }=0.5$, the edge circulations increase only slightly at late times, while for $\text{ Fr } \gtrsim 1$, the edge circulations grow significantly larger. For $\text{ Fr }=1$, both edge circulations are stationary for a while at later times, but for $\text{ Fr }=2$, one or both of the edge circulations grow continuously. These differences indicate that the variation of the circulation is closely related to the evolution of the flow field.

Despite our extension of the model, the computational time for the falling plates is not sufficient to observe the steady-state behavior. There is a limitation in the long-time simulation of the model, due to $O(N^2)$ computation for the Birkhoff-Rott integral, as commented before. An efficient computation of the integration of the model can be carried out by applying fast summation techniques such as a tree code algorithm [40] or fast multipole method [23]. On the other hand, several recent studies [38, 41, 42] have developed an approach of representing free vortex sheets into point vortices with regularized cores. This approach reduce the computation cost significantly. Furthermore, the problem of close interactions between the plate and vortices would be lessened to some extent. We leave these issues for a future study.

The model and numerical method in this study can be extended to cambered or flexible plates. In these cases, the kinematic condition includes the geometric effect of the body and can be written as a singular Fredholm integral equation of the first kind for the bound vortex sheet strength [26, 28]. There have been few studies on the free falling motion of cambered and flexible plates, and the exploration on these bodies will be interesting.

Discretization and time-integration

We discretize the free and bound vortex sheets by N lagrangian point vortices. We here assume that the free vortex sheets are separated from the leading and trailing edges and there is no free vortex sheet detached to the bound vortex sheet. The numerical method can be easily extended to multiple free vortex sheets.

Using the circulation $\Gamma$ as a Lagrangian variable, we denote the location of the free vortex sheet at current time $t^n$ by $z_j^n = z_{\pm } (\Gamma _j,t^n)$ for $1 \le j \le P-1$ and $Q+1 \le j \le N$. The Lagrangian indices of the leading and trailing edges of the bound vortex sheet are denoted by the integer P and Q, respectively. The bound vortex sheet is discretized by using Gauss-Lobatto collocation nodes [18], which are given by

$$\begin{aligned} s_j = \frac{1}{2} \left[ 1 - \cos \left( \frac{\pi (j-P)}{Q-P} \right) \right] , \qquad P \le j \le Q. \end{aligned}$$

(A1)

The edge circulation is expressed as $\Gamma _P^n = \Gamma _+ (t^n)$ and $\Gamma _Q^n = \Gamma _- (t^n)$.

We assume that $z_j^{n}$, $c^{n}$, $\theta ^{n}$, $\dot{c}^{n}$ and $\dot{\theta }^{n}$ are given at time $t^n$. At each time step, we first determine $\Gamma _P^n$, $\Gamma _Q^n$, and $\gamma _j^n$ by solving the linear Eq. (12) and (14), and Eq. (11), and calculate the velocity of the free vortex sheet from Eq. (5). In Eq. (6), the integration over the free vortex sheet is approximated by using the trapezoidal rule, while the integration for the bound vortex sheet is obtained by using the trapezoidal rule, with the numerical method of a near-singular integral if point vortices in the free vortex sheet are close to the plate.

From Eq. (23), $\dot{f}_j^n$ are the function of $\ddot{c}^n$, $\ddot{\theta }^n$, $\dot{\Gamma }_P^n$ and $\dot{\Gamma }_Q^n$, and in turn, $\dot{\gamma }_j^n$ and $\dot{\Gamma }_j^n$ are also functions of these quantities, from Eqs. (21) and (20). Therefore, Eqs. (15), (16), (22) and (24) are a system of linear equations for the unknowns $\ddot{c}^n$, $\ddot{\theta }^n$, $\dot{\Gamma }_P^n$ and $\dot{\Gamma }_Q^n$. Finding these quantities from the linear equation, we are ready for time-integration. We employ the classical fourth-order Runge–Kutta (RK) method for time-integration to obtain $z_j^{n+1}$, $c^{n+1}$, $\theta ^{n+1}$, $\dot{c}^{n+1}$ and $\dot{\theta }^{n+1}$.

To find the time derivatives of $\Gamma _P^n$ and $\Gamma _Q^n$ is the key step in the numerical method for this type of vortex-body interaction. Our method of finding them differs from previous studies that used the Chebyshev series [16], or an implicit scheme for solving a nonlinear system of equations [23, 27]. We use the formulation of the time derivative of the kinematic condition and determine $\dot{\Gamma }_P^n$ and $\dot{\Gamma }_Q^n$ by solving a system of linear equations.

At each time step, two new point vortices are released from both edges of the body. The free vortex sheet lacks resolution at late times, due to the nonuniform distribution of point vortices. To handle this problem, an adaptive point insertion procedure is applied to maintain the resolution of free vortex sheets [17, 18]. The third-order local polynomial interpolation is used to insert points whenever the distance between two consecutive points exceeds a given threshold.

The solution procedure for the equations requires the initial conditions at a time $t_1 > 0$. We use the self-similar solution of the free vortex sheet separated from a plate for the initial conditions [16]. Two starting point vortices with zero circulation $\Gamma _1 = \Gamma _N = 0$ are placed initially and their locations $z_1$ and $z_N$ are given by the small-time asymptotic expansion of the self-similar solution.

The numerical method is summarized as the following time-integration algorithm:

For given $z_j^{n}$, $c^{n}$, $\theta ^{n}$, $\dot{c}^{n}$ and $\dot{\theta }^{n}$,

• (1) At the first step of the RK method,

  • (i) determine $\Gamma _P^n$, $\Gamma _Q^n$, and $\gamma _j^n$.

  • (ii) compute the velocity of the free vortex sheet, $\dot{z}_j^n$.

  • (iii) find $\ddot{c}^n$, $\ddot{\theta }^n$, $\dot{\Gamma }_P^n$ and $\dot{\Gamma }_Q^n$.

  • (iv) compute $z_j^{n+1/2}$, $c^{n+1/2}$, $\theta ^{n+1/2}$, $\dot{c}^{n+1/2}$ and $\dot{\theta }^{n+1/2}$.

• (2) Repeat steps (i)-(iv) in (1) similarly to the other steps of the RK method.

• (3) Find $z_j^{n+1}$, $c^{n+1}$, $\theta ^{n+1}$, $\dot{c}^{n+1}$ and $\dot{\theta }^{n+1}$ by combining the intermediate values.

• (4) Find the location of new point vortices released from the plate edges.

• (5) A point insertion procedure is performed, when necessary.

To validate the numerical method of the present model, we applied it to the falling plate with the uniform $\delta$-regularization with $\delta =0.2$ and checked that the results of the present model agree well with Jones and Shelley [16].