Modeling and Simulation of a Fluttering Cantilever in Channel Flow

TL;DR

This paper develops a coupled fluid-structure model for a fluttering cantilever in channel flow, validating it against simulations and exploring its accuracy across various parameters, with implications for energy harvesting and biomedical applications.

Contribution

The paper introduces a new coupled fluid-structure model for viscous channel flow that accurately predicts flutter boundaries over a wide parameter range, extending beyond previous models.

Findings

Model accurately predicts flutter boundaries within valid parameter ranges.

Model predictions converge to inviscid results at high Reynolds numbers.

Validation against DNS confirms model's applicability for small channel heights.

Abstract

Characterizing the dynamics of a cantilever in channel flow is relevant to applications ranging from snoring to energy harvesting. Aeroelastic flutter induces large oscillating amplitudes and sharp changes with frequency that impact the operation of these systems. The fluid-structure mechanisms that drive flutter can vary as the system parameters change, with the stability boundary becoming especially sensitive to the channel height and Reynolds number, especially when either or both are small. In this paper, we develop a coupled fluid-structure model for viscous, two-dimensional channel flow of arbitrary shape. Its flutter boundary is then compared to results of two-dimensional direct numerical simulations to explore the model's validity. Provided the non-dimensional channel height remains small, the analysis shows that the model is not only able to replicate DNS results within the parametric limits that ensure the underlying assumptions are met, but also over a wider range of Reynolds numbers and fluid-structure mass ratios. Model predictions also converge toward an inviscid model for the same geometry as Reynolds number increases.

Tables

Table 1. Table 1 : Table of fluid-structure dimensional parameters.

Variable	Description	Dimension
$\delta$	beam displacement	$l$
$x$	beam length coordinate	$l$
$t$	time	$t$
$U_c$	characteristic velocity	$l\ast t^{-1}$
$L$	beam length	$l$
$h_b$	beam thickness	$l$
$b$	beam width	$l$
$\overline{h}$	throat height	$l$
${\rho }_f$	fluid density	$m\ast l^3$
${\mu }_f$	fluid viscosity	$m\ast l^{-1}\ast t^{-1}$
${\rho }_s$	beam density	$m\ast l^{-3}$
$E$	Young’s modulus	$m\ast l^{-1}\ast t^{-2}$
$\nu$	Poisson’s ratio	ND

Table 2. Table 2 : Table of solutions to the characteristic equation 28 .

${\beta }_{1}L$	${\beta }_{2}L$	${\beta }_{3}L$	${\beta }_{4}L$	${\beta }_{5}L$	${\beta }_{6}L$
1.8751	4.6941	7.8548	10.9955	14.1372	17.2788

Table 3. Table 3 : Table of clamped-free fluid-structure nondimensional parameters.

Variable	Expression	Description
$M^{\ast }$	$\frac{{\rho }_{f}L}{{\rho }_s{h}_b}$	mass ratio
$U^{\ast }$	$\frac{q_{x0}L}{\overline{h}}{\left(\frac{{\rho }_s{h}_b}{E{I}_{2D}}\right)}^{\frac{1}{2}}$	stiffness ratio
$\hat{h}$	$\frac{\overline{h}}{L}$	gap or throat ratio
${\hat{h}}^{2}R{e}_L$	${\left(\frac{\overline{h}}{L}\right)}^2\left(\frac{{\rho }_f{q}_{x0}L}{\overline{h}{\mu }_f}\right)$	viscous parameter
$St$	$\frac{\text{Im}\left[\lambda \right]}{2\pi }=\frac{f_{s}L\overline{h}}{q_{x0}}$	Strouhal number
$f_s^{\ast }$	$U^{\ast }St$	nondimensional frequency

Table 4. Table 4 : Table of cases for constant channel flow simulations with U ∗ superscript 𝑈 U^{*} as the bifurcation parameter and varying M ∗ superscript 𝑀 M^{*} and h ^ 2 ​ R ​ e L superscript ^ ℎ 2 𝑅 subscript 𝑒 𝐿 \hat{h}^{2}Re_{L} .

Case $\#$	$\hat{h}$	${\hat{h}}^{2}R{e}_L$	$M^{\ast }$
1	0.025	0.500	[$0.01-0.50$]
2	0.050	0.500	[$0.01-1.00$]
3	0.050	1.250	[$0.01-1.00$]
4	0.050	2.500	[$0.01-1.00$]
5	0.125	0.500	[$0.01-1.00$]
6	0.125	1.250	[$0.01-1.00$]
7	0.050	[0.10-4.50]	0.01
8	0.125	[0.10-9.50]	0.01

Equations

$${\rho }_s{h}_b\frac{{\partial }^2}{\partial t^2}\delta (x,t)+\frac{{\partial }^2}{\partial x^2}\left(E{I}_{2D}\frac{{\partial }^2}{\partial x^2}\delta (x,t)\right)=P^{\mathrm{bot}}(x,t)-P^{\mathrm{top}}(x,t),$$

$$I_{2D}=\frac{h_b^3}{12(1-{\nu }^2)}.$$

$$\delta (0,t)=0,\frac{\partial }{\partial x}\delta (0,t)=0,\frac{{\partial }^2}{\partial x^2}\delta (L,t)=0,\frac{{\partial }^3}{\partial x^3}\delta (L,t)=0.$$

$$\delta (x,t)={\delta }^{\mathrm{top}}(x,t)=-{\delta }^{\mathrm{bot}}(x,t).$$

$$-\frac{\partial \delta }{\partial t}+\frac{\partial Q_x}{\partial x}=0,$$

$$Q_x={\int }_{\delta }^{h_0}udy,$$

$$\begin{array}{r}\frac{\partial }{\partial t}\left({\int }_{\delta }^{h_0}\left[\begin{array}{c}u\\ v\end{array}\right]dy\right)+\frac{\partial }{\partial x}\left({\int }_{\delta }^{h_0}\left[\begin{array}{c}{u}^2\\ uv\end{array}\right]dy\right)=-\frac{1}{{\rho }_f}\left\{\frac{\partial }{\partial x}\left({\int }_{\delta }^{h_0}\left[\begin{array}{c}P\\ 0\end{array}\right]dy\right)-\left[\begin{array}{c}{h}_0^{\prime }\\ -1\end{array}\right]P{|}_{y=h_0}+\left[\begin{array}{c}{\delta }^{\prime }\\ -1\end{array}\right]P{|}_{y=\delta }-{\mathbf{F}}_{\mathrm{visc}}\right\},\end{array}$$

$$\frac{\partial Q_x}{\partial t}+\frac{\partial {\mathcal{N}}_x}{\partial x}=-\frac{1}{{\rho }_f}\left\{\frac{\partial }{\partial x}\left({\int }_{\delta }^{h_0}Pdy\right)-h_0^{\prime }P{|}_{y=h_0}+\frac{\partial \delta }{\partial x}P{|}_{y=\delta }-F_{\mathrm{visc},x}\right\},$$

$${\mathcal{N}}_x={\int }_{\delta }^{h_0}{u}^{2}dy,$$

$${\mathcal{N}}_x=\frac{6}{5}\frac{Q_x^2}{h_0-\delta }\text{, and }{F}_{\mathrm{visc},\mathrm{x}}=-12{\mu }_f\frac{Q_x}{{\left(h_0-\delta \right)}^2},$$

$${\mathcal{N}}_x$$

$$F_{\mathrm{visc},\mathrm{x}}$$

$$f=\{\begin{array}{ll}48R{e}_h^{-1}& \text{if }R{e}_h<1000\\ 0.26R{e}_h^{-0.24}& \text{if }R{e}_h\ge 1000\end{array},$$

$${\xi }_x=\{\begin{array}{ll}6/5& \text{if }R{e}_h<1000\\ 1& \text{if }R{e}_h\ge 1000\end{array},$$

$$\frac{\partial Q_x}{\partial t}+{\xi }_x\frac{\partial }{\partial x}\left(\frac{Q_x^2}{h_0-\delta }\right)=-\frac{h_0-\delta }{{\rho }_f}\frac{\partial P}{\partial x}-\frac{f}{4}\frac{Q_x^2}{(h_0-\delta {)}^2}.$$

$$P(t){|}_{x=0}=P_{\mathrm{in}}-\frac{{\zeta }_{\mathrm{in}}}{2}{\rho }_f{\left[{\left(\frac{Q_x}{h_0-\delta }\right)}^2\right]}_{x=0},P(t){|}_{x=L}=P_{\mathrm{out}}+\frac{{\zeta }_{\mathrm{out}}}{2}{\rho }_f{\left[{\left(\frac{Q_x}{h_0-\delta }\right)}^2\right]}_{x=L}.$$

$$\delta (x,t)$$

$$P(x,t)$$

$$Q_x(x,t)$$

$$f(Q_x)$$

$$\approx f_0+\epsilon \eta q_{x1}(x,t)+\dots ,$$

$$E{I}_{2D}\frac{{\mathrm{d}}^4}{\mathrm{d}{x}^4}{\delta }_0(x)=p_0^{\mathrm{bot}}-p_0^{\mathrm{top}},$$

$$p_0(x)$$

$$q_{x0}$$

$${\rho }_s{h}_b\frac{{\partial }^2{\delta }_1}{\partial t^2}+E{I}_{2D}\frac{{\partial }^4{\delta }_1}{\partial x^4}=p_1^{\mathrm{bot}}-p_1^{\mathrm{top}},$$

$$\frac{\partial q_{x1}}{\partial t}=\frac{\partial {\delta }_1}{\partial x},$$

$$\begin{array}{r}\frac{\partial p_1}{\partial x}=\frac{{\rho }_f}{h_e}\left\{\left[\frac{{\xi }_x{q}_{x0}^2}{h_e^2}\frac{\partial }{\partial x}-\frac{3p_0^{\prime }}{{\rho }_f}\right]{\delta }_1-\left[\frac{\partial }{\partial t}+\frac{2{\xi }_x{q}_{x0}}{h_e}\left(\frac{\partial }{\partial x}-\frac{h_e^{\prime }}{h_e}\right)+\frac{q_{x0}}{2{h_e}^2}\left(f_0+\frac{\eta q_{x0}}{2}\right)\right]q_{x1}-\frac{\partial q_{x1}}{\partial t}\right\}.\end{array}$$

$$p_1(0,t)=\frac{2\left(P_{\mathrm{in}}-p_0\left(0\right)\right)}{h_e\left(0\right)}{\delta }_1\left(0,t\right)-{\zeta }_{\mathrm{in}}\frac{{\rho }_f{q}_0}{{h_e\left(0\right)}^2}{q}_{x1}\left(0,t\right)$$

$$p_1(L,t)=\frac{2\left(P_{\mathrm{out}}-p_0\left(L\right)\right)}{h_e\left(L\right)}{\delta }_1\left(L,t\right)+{\zeta }_{\mathrm{out}}\frac{{\rho }_f{q}_0}{{h_e\left(L\right)}^2}{q}_{x1}\left(L,t\right),$$

$${\delta }_1(x,t)=\sum _{i=0}^n{a}_i(t)g_i(x)$$

$$g_i(x)=\{\begin{array}{ll}0& \text{for }i=0\\ {\varphi }_i(x)& \text{for }i=[1,n]\end{array},$$

$$\frac{{\mathrm{d}}^4\varphi }{\mathrm{d}{x}^4}+{\beta }^4\varphi =0.$$

$$\phi_{k}(x)=\cosh\!\left(\beta_{k}\,x\right)-\cos\!\left(\beta_{k}\,x\right)+\left[\frac{\cos\!\left(\beta_{k}\,{L}\right)+\cosh\!\left(\beta_{k}\,{L}\right)}{\sin\!\left(\beta_{k}\,{L}\right)+\sinh\!\left(\beta_{k}\,{L}\right)}\right]\Big{(}\sin\!\left(\beta_{k}\,x\right)\,-\sinh\!\left(\beta_{k}\,x\right)\,\Big{)},$$

Full text

††journal: Journal of Fluids and Structures

1 Introduction

The stability of an elastic member within a channel, or as part of the channel, has been studied for many decades (Johansson, 1959, Miller, 1960, Inada & Hayama, 1988, 1990, Nagakura & Kaneko, 1991). Applications include wind instruments (Sommerfeldt & Strong, 1988, Backus, 1963), human snoring (Balint & Lucey, 2005, Tetlow & Lucey, 2009), vocalization (Tian et al., 2014), and enhanced heat transfer (Shoele & Mittal, 2016, Hidalgo et al., 2015). Recently, this geometry has also been used for flow-energy harvesting, with devices specifically targeting power generation for remote sensor networks (Sherrit et al., 2014, 2015, Lee et al., 2015, 2016). For most of these applications, the flutter instability boundary is the essential result sought, as the functional requirements of engineering designs (i.e. flow-energy harvester, heat management systems) or the manifestation of sound in natural systems (i.e instruments, snoring, voice) depend on it. A particular challenge with characterizing the flutter onset is its dependence on the specific flow regimes, which vary between and within applications noted. In this paper, we target the prediction of the flutter instability over laminar and turbulent regimes, and show that with specific restrictions to fluid-to-structure mass and throat-to-length ratios, our model agrees well with both viscous numerical simulations and inviscid modeling results.

With recent advancements, it has become possible to directly simulate the fluid-structure interaction (FSI) numerically by solving the Navier-Stokes equations coupled to a model of the structure. Two dimensional FSI algorithms have been used to study the stability of an elastic member within channel flow (Tetlow & Lucey, 2009, Shoele & Mittal, 2016), and more recently to assess the effect of Reynolds number on the ensuring flutter boundary (Cisonni et al., 2017). Yet, one of the challenges with fluid-structure systems is the large number of parameters necessary to describe the subset of possible system regimes; fully coupled computational approaches require considerable computing time to capture the flutter instability, often only being able to span a fraction of this parameter space. A more tractable but less accurate (or versatile) approach involves modeling of the structure displacement and velocities, with fluid forces approximated via simplified equations of motion. The early works of Miller (Miller, 1960) and Johansson (Johansson, 1959) address the divergence instability in this light specifically targeting channels within nuclear reactor cooling systems. More recent work by Guo & Paidoussis (Guo & Paidoussis, 2000) takes an inviscid approach to understanding the onset of flutter in a symmetric channel. Alben (Alben, 2015) solves for the inviscid flutter boundary using a vortex sheet model, and Shoele & Mittal (Shoele & Mittal, 2016) extend the plane wake vortex sheet method by Alben (Alben, 2008, 2015) to asymmetric channel flow.

As recently evident from (Cisonni et al., 2017), the inviscid treatment of channel flow for small throat-to-length ratios is unable to predict stability boundaries that more generally depend on the Reynolds number. Large deviations are observed between their results and the inviscid models, particularly when both Reynolds number and throat-to-length ratio are small. Moreover, Cisonni et al. (Cisonni et al., 2017) are unable to address the bridge between the viscous and inviscid stability regimes due to the Reynolds number restriction in their numerical simulations.

Apart from fully-coupled numerical simulations that include viscous effects, it is possible to include viscosity in an approximate way in reduced-order models. Nagakura & Kaneko (Nagakura & Kaneko, 1991) and Wu & Kaneko (Wu & Kaneko, 2005) employ models with a viscous formulation of fluid forces for elastic beams in channel flow. This form for fluid forces was originally proposed by Inada & Hayama (Inada & Hayama, 1988) for rigid plates in converging or diverging channels, with Fujita & Shintani further expanding the analysis to flexible, cylindrical constant channel flows (Fujita & Shintani, 1999, 2001, 2007). Yet, the parametric bounds for validity of this model remain largely untested. This includes not only viscosity effects, but also geometrical parameters (i.e. throat-to-length ratio), and fluid-to-structure mass and stiffness ratios.

This paper addresses these parameter bounds first by deriving a fully coupled fluid-structure model that considers flow confinement in arbitrarily shaped channels through a set of assumptions grounded in parametric limits. We then compare the model predicted flutter stability boundaries to those from fluid-structure direct numerical simulations (DNS), as well as to results of the inviscid model from Shoele & Mittal (Shoele & Mittal, 2016). Though derived differently, our model accounts for fluid forces in a similar fashion to (Nagakura & Kaneko, 1991, Wu & Kaneko, 2005) when the constant channel geometry is considered, and to Inada & Hayama (Inada & Hayama, 1988) when a rigid beam is considered. The difference comes as a factor in the advection term, and careful attention to the parametric validity established through the control volume analysis and subsequent approximations of the geometry and flow field.

Our DNS algorithm uses an immersed-boundary projection formulation developed by Taira & Colonius in (Taira & Colonius, 2007, Colonius & Taira, 2008). The complementary FSI algorithm is a strongly-coupled formulation between the flow and the structure developed by Goza & Colonius (Goza & Colonius, 2017). Both modeling and simulations only consider two-dimensional flow, with the modeling further simplifying the problem into a quasi-1 dimensional framework. An extension to three-dimensional flow is possible (Tosi, 2018), and will the subject of forthcoming publications.

The paper is organized as follows. We first derive the quasi-1 dimensional fluid-structure model, which captures the linearized incompressible Navier-Stokes and solid equations of a cantilever beam in an arbitrarily shaped channel configuration. Next, the fluid-structure direct numerical simulation algorithm is described, along with the dynamic mode decomposition employed to extract the frequency, growth rate, and beam shapes of unstable modes. Lastly, modeling and simulation results for a cantilever beam in constant channel flow are compared over a broad range of parameters. Our model results are then compared to those from the inviscid model as Reynolds number is varied.

2 Quasi-1 Dimensional Fluid-Structure Model

The geometry illustrated in figure 1(a) and dimensional parameters in table 1 are inspired by the flow energy harvester configurations in (Sherrit et al., 2014, 2015).

We consider the beam displacement as the output of a system defined by a characteristic velocity (or flow rate), geometrical parameters, and material properties, for a total of 10 possible nondimensional groups that determine its dynamics. This large number of parameters makes a thorough numerical or experimental investigation difficult. The purpose of this section is to provide a simple model that allows us to analytically identify the most important properties affecting the stability of the system.

2.1 Structure equations of motion

Considering the cantilever in figure 1(a) in transverse vibration, we apply the undamped, two-dimensional Euler-Bernoulli beam equation (Inman, 2008),

[TABLE]

where the two-dimensional area moment of inertia is

[TABLE]

$P^{\mathrm{bot}}(x,t)$ and $P^{\mathrm{top}}(x,t)$ are the pressures acting on the bottom and top of the beam, respectively (i.e. per superscript). The beam is clamped at its leading edge and free at its trailing edge, so the boundary conditions are

[TABLE]

In considering the flow separately in the top and bottom channels, we write the geometrical constraint,

[TABLE]

2.2 Fluid Equations of Motion

To develop a relation between pressure, beam displacement and its derivatives, we consider the control volume defined in $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ illustrated in figure 1(b), which corresponds to a small section of one of the channels in figure 1(a). The channel is formed between the moving cantilever boundary at a variable height $\delta(x,t)$ and the upper surface at a specified location, $h_{0}(x)$.

We apply mass and momentum conservation to this control volume under the simplifying assumptions of constant fluid density and a gradually-varying channel in the streamwise direction, $h_{0}^{\prime 2}\ll 1$ and $\delta^{\prime 2}\ll 1$, such that $\sqrt{1+h_{0}^{\prime 2}}\approx 1$ and $\sqrt{1+\delta^{\prime 2}}\approx 1$ for $x\in[0,L]$. From mass conservation, we obtain

[TABLE]

where

[TABLE]

is the volume flow rate. From momentum conservation,

[TABLE]

where the $\mathbf{\hat{x}}$ component, after substituting $Q_{x}$,

[TABLE]

[TABLE]

and $F_{\mathrm{visc},x}$ is the net viscous stress acting on the walls.

To make further progress, we need a closure relation for $\mathcal{N}_{x}$ and $F_{\mathrm{visc,x}}$ in term of $Q_{x}$ and $\delta$, as well as a relation between the integrated and evaluated quantities of $P$ in equation 8. Such a closure can be rigorously determined for a narrow gap where both $\frac{\bar{h}}{L}=\hat{h}\ll 1$ and $\hat{h}^{2}Re_{L}\ll 1$ so that inertial terms and streamwise diffusion of momentum can be locally neglected, i.e. the lubrication limit (Kundu et al., 2012, p. 319). This lubrication scaling of the N-S equations leads to the pressure being constant across the gap, $P=P(x,t)$, and the streamwise velocity becoming parabolic. Using the definition in equation 6 and integrating a parabolic $u$ leads to expressions for $Q_{x}$ in terms of $\delta$ and $P$. $\mathcal{N}_{x}$ and $F_{\mathrm{visc,x}}$ can then be evaluated as

[TABLE]

with the latter taking the form for a Newtonian fluid.

Even when $\hat{h}^{2}Re_{L}$ is not small, a closure can be determined under the weaker assumption that the profile shape is unchanging with $x$ and $t$. For a particular velocity profile $u$, we define a profile factor, $\xi_{x}$ and Fanning friction factor, $f$, so that the appropriate relations become

[TABLE]

respectively. Following (Shimoyama & Yamada, 1957, Inada & Hayama, 1988, Nagakura & Kaneko, 1991), we take

[TABLE]

whereas we model the profile factor as

[TABLE]

where the laminar value ($Re_{h}<1000$) coincides with the aforementioned lubrication result, and the turbulent case follows from the blunted mean velocity profile in the outer region and neglects the thin inner region.

With this closure and $P=P(x,t)$, we need not consider the $\mathbf{\hat{y}}$ component of equation 7. Equation 14 can be simplified as

[TABLE]

In order to solve equations 5 and 14 uniquely, two pressure boundary conditions are required. A simple and common treatment is to use the extended Bernoulli equation with empirical loss coefficients associated with the specific geometry and flow conditions (including Reynolds number) near the inlet and outlet. We adopt here the approach used to treat leakage-flow instabilities (Nagakura & Kaneko, 1991, Inada & Hayama, 1988, 1990),

[TABLE]

where $\zeta_{\mathrm{in}}\geq 1$ and $\zeta_{\mathrm{out}}\geq 0$ are loss coefficients, and the departure from equality represents non-isentropic processes. $P_{\mathrm{out}}$ and $P_{\mathrm{out}}$ are constants.

To summarize, equations 5 and 14, together with the boundary conditions 15 relate the pressure, flow rate, and deflection in either channel. Both sets are closed with the geometrical constraint, equation 4, the common boundary conditions ($P_{\mathrm{in}}$ and $P_{\mathrm{out}}$), and the common equation of motion for the beam, equation 1. The top and bottom channels correspond to distinct channel shapes $h_{0}^{\mathrm{top}}$ and $h_{0}^{\mathrm{bot}}$, respectively. These constitute the closed quasi-1D fluid-structure model.

Before proceeding, we review the assumptions underpinning the model. In addition to flow two-dimensionality, we require that the gap is thin $\hat{h}\ll 1$ and that axial variations in the channel gap are small: $h_{0}^{\prime 2}\ll 1$ and $\delta^{\prime 2}\ll 1$. In addition, we require the velocity profile to be slowly changing with both $x$ and $t$, so that $\xi_{x}$ can be regarded as a constant, as would strictly be the case when $\hat{h}^{2}Re_{L}\ll 1$. As we will show, the model also performs well in circumstances where $\hat{h}^{2}Re_{L}$ is not small. Our goal, after linearizing the model, is to test the validity of these bounds by contrasting them with results from numerical simulations.

2.3 Linearized model

The primary goal of this paper is to predict the linear stability (flutter boundary) of an equilibrium beam shape $\delta_{0}(x)$, as a function of parameters on table 1. We begin this process by expanding the dependent variables about their respective equilibrium values in a small parameter, $\varepsilon$, representing the amplitude of the beam displacement. That is, we take

[TABLE]

as well as the linearized friction factor

[TABLE]

determined from laminar and turbulent relations in equation 11 with the $Re_{h}=\hat{h}Re_{L}$.

At zeroth order, we obtain a differential equation describing the equilibrium beam shape

[TABLE]

together with the homogeneous beam boundary conditions of equation 3. Once again, the superscripts top and bot refer to parameters associated with $h_{0}^{\mathrm{top}}$ and $h_{0}^{\mathrm{bot}}$ as the channel shapes above and below the beam, respectively. The pressure distribution and flow rate in either channel (i.e. top and bottom) are given by

[TABLE]

where we have set $h_{e}(x)=h_{0}(x)-\delta_{0}(x)$, the equilibrium channel height. We note that at equilibrium, the flow rate in each channel is a constant independent of $x$. As an example, $p_{0}^{\mathrm{top}}$ in equation 16 refers to equations 17 and 18 where $h_{e}^{\mathrm{top}}$ is substituted for $h_{e}$, while $h_{e}^{\mathrm{bot}}$ substituted into $p_{0}^{\mathrm{bot}}$. The same is true for the subsequent higher order terms in $\varepsilon$ (i.e. $p_{1}$ and $q_{x1}$ with top and bot superscripts).

The steady equilibrium equations here are similar to those obtained by Inada & Hayama (Inada & Hayama, 1988), assuming $\xi_{x}=1$, and $\delta_{0}=0$. Equation 17 has three distinct terms in the parenthesis, the first with $f_{0}$ as a factor represents the pressure drop due to viscous effects as a function of the integral equilibrium gap shape over $x$; the second integral, multiplied by $\xi_{x}$, comes from the inertia term and is only a function of the initial ($x=0$) and the current (at $x$) gap size, representing the pressure change due to the gap expansion or contraction; and the third is solely due to system inlet conditions. If $\zeta_{\mathrm{in}}=1$ (isentropic), then the only pressure drop comes from accelerating the flow to the average inlet velocity.

Next we collect and equate coefficients to linear order in $\varepsilon$. For the beam, we obtain

[TABLE]

together with its homogeneous boundary conditions. For the pressure and flow rate, we have

[TABLE]

[TABLE]

To obtain a closed form for $p_{1}$, equation 20 is first integrated in $x$, solved for $q_{x1}$ in terms of $\delta_{1}$ and $q_{x1}(0,t)$, and substituted into equation 21. The result is, once again, integrated in $x$, and together with the boundary conditions

[TABLE]

and

[TABLE]

and with the beam equation, can be solved solely for $\delta_{1}$ and $q_{x1}$.

2.4 Numerical solution of the perturbation equations

To numerically solve the linear system of PDEs given by equations 19 to 21, we expand the first-order beam displacement in a series of basis functions

[TABLE]

where

[TABLE]

and $\phi_{i}(x)$ are solutions of the homogeneous (unforced) beam equation,

[TABLE]

in the domain $x\in[0,L]$. Its solution, $\phi_{k}$, $k\in\mathbb{Z}:[1,\infty]$, when subject to the clamped-free boundary conditions is

[TABLE]

with characteristic equation

[TABLE]

The first six values are listed in table 2.

Equation 24 is substituted into the coupled beam-pressure equation (equation 19 and solution of 21), the result which is then projected onto the same basis functions (Galerkin projection). By construction, the equation for the time evolution of $q_{x1}(0,t)$ is not a function of $x$, but rather of the boundary points at $x=0$ and $x=L$. The final total number of unknowns is $2n+2$. These include the expansion coefficients $a_{i}(t)$, their time derivatives, ${\dot{a}}_{i}(t)$, and the remaining two unknowns are the entrance flow rate $q_{z1}(0,t)$ for the top and bottom channels. We write the solution vector as

[TABLE]

the resulting system of ODE are

[TABLE]

The entries of ${\bf A}$ are given in the appendix. The eigenvalues and eigenvectors of ${\bf A}$ are computed in section 3.1 to determine the flutter boundary for the coupled FSI system. We note here that the model, once nondimensionalized, involves four independent nondimensional groups given in table 3: $U^{*}$, $M^{*}$, $\hat{h}$, and $\hat{h}^{2}Re_{L}$. The table also contains the beam frequency response parameters, including the Strouhal number $St$ as a function of the imaginary part of the eigenvalue $\lambda$ or the dimensional beam frequency response $f_{s}$.

3 Immersed-Boundary Direct Numerical Simulation and Data Processing

In order to validate aspects of the quasi-1D model developed in section 2, we employ a two-dimensional fluid-structure algorithm (Goza & Colonius, 2017) that utilizes the immersed boundary (IB) projection method (Taira & Colonius, 2007, Colonius & Taira, 2008) along with Newton-Raphson approach to solve the strongly-coupled fluid-structure system. Strong-coupling ensures that the nonlinear constraint between the fluid and the structure is enforced at each time step, and is necessary for accurate computation of large structural deformations.

The current implementation uses a co-rotational formulation of the structural equations for the beam, where strains are assumed small within constituent material equations in the beam-local frame (Criesfield, 1991). The beam-fluid nondimensional parameters that govern the dynamics are $Re_{L}$, $M^{*}$, and $U^{*}$ and the same as those defined in the quasi-1 dimensional model in table 3.

The internal flow passage is created by immersed boundaries, including one spanning the inlet where the streamwise velocity is set to a given profile and the normal velocity is set to zero.

The initial discrete delta function used in the IB method is constructed as follows. The 3-point kernel of Yang et al. (Yang et al., 2009) is smoothed three times using their recursive formula such that its total support is 6 points. This kernel was chosen as a compromise between computational efficiency and accuracy and smoothness of the local stresses (Goza et al., 2016).

The algorithm has been verified extensively in (Goza, 2018) for external flows, where regimes of the standard and inverted flags were explored and compared to results from other strongly-coupled fluid-structure solvers. It has also been verified (Tosi, 2018) for internal flows using the suggested benchmark of an elastic member in an internal, incompressible, laminar flow and compared to (Turek & Hron, 2006, Tian et al., 2014, Bhardwaj & Mittal, 2012, Shoele & Mittal, 2014). Results including frequency and amplitude response, along with the drag coefficient, were consistent with the other compared schemes for the two cases at $Re_{L}=350$ and 700.

In order to evaluate the linear stability of the FSI system, we employ an empirical approach based on the Dynamic Mode Decomposition (DMD) to calculate growth/decay rates, frequencies, and mode shapes directly from the (nonlinear) simulation data. DMD is a data decomposition technique that approximates the eigenmodes of the linear operator that best describes the dynamics of the system, i.e. from one time instance to the next (Schmid, 2010, Chen et al., 2012, Tu et al., 2013, Goza & Colonius, 2018). We take the exact-DMD algorithm and augment our data matrix, comprising of beam $\mathbf{\hat{y}}$ positions along its Lagrangian IB coordinate, with the eigensystem realization algorithm (ERA) to ensure robustness of our eigenvalue calculations (Tu et al., 2013). DMD eigenvalues and modes for the time series near the beam equilibrum are used to compare with quasi-1D model eigenvalues and modes.

3.1 Cantilever in Constant Channel

Although the quasi-1D model is formulated for arbitrarily shaped channels, flow in a constant, symmetric channel presents a simple geometry that can be simulated and compared to existing inviscid models (Alben, 2015, Shoele & Mittal, 2016), and provides a relevant problem configuration to many of the aforementioned applications.

The FSI DNS computational domain is illustrated in figure 2. The coordinate system has its origin at the beam leading edge, and the channel is defined at a constant half width, $\bar{h}$. The initial beam position $\delta(s,0)=0$, and its velocity is $\frac{\partial}{\partial t}\delta(s,0)=0$, where $s$ is the Lagrangian coordinate that describes the beam parameterized by its arc length.

The boundary conditions are clamped and free for leading and trailing edges of the beam, respectively. A uniform velocity profile is specified at the inlet as $U_{\mathrm{in}}$ and used as the reference parameter for non-dimensionalization of other quantities. The flow is impulsively started and the beam is perturbed by a small body force at the initial time step in order to break symmetry.

Figure 2 is contrasted with the domain in figure 1(a), where quasi-1D model parameters are illustrated ($U_{c}=\frac{q_{x0}}{\bar{h}}$). First, the steady flow rate, $q_{x0}$, specified for the quasi-1D model, constitutes half the integrated $\mathbf{\hat{x}}$ flow velocity (in the channel) over $\mathbf{\hat{y}}$ from $-\bar{h}$ to $\bar{h}$ immediately downstream of the inlet. The beam boundary conditions for the quasi-1D model are the same as the DNS. Results also do not depend on inlet and outlet distances except through loss coefficients $\zeta_{\mathrm{in}}$ and $\zeta_{\mathrm{out}}$ defined in equations 15. We assume that stagnation pressure loss at either end is negligible and take the no loss coefficients values ($\zeta_{\mathrm{in}}=1$ and $\zeta_{\mathrm{out}}=0$).

Lastly, we are interested in the dynamics near the DNS initial conditions, which, apart from the small body force to deflect the beam, represents a (potentially unstable) equilibrium for both the DNS and the quasi-1D model. In this limit the coordinate $x$, shown in both figures 1(a) and 2, coincides with the Lagrangian coordinate $s$ in figure 2. For small perturbations, the Lagrangian beam shape then becomes the beam $\mathbf{\hat{y}}$ displacement, $\delta(x,t)$.

We will use DMD to find the least-stable mode of the DNS results using the procedure discussed. These results are directly comparable with the eigenvalues and eigenmodes of the quasi-1D linear operator $\mathbf{A}$ in equation 30. In order to compare to DMD growth rate, the dynamically significant eigenvalues of the quasi-1D will be shown as well (the least, and sometimes second least, stable eigenvalues for the quasi-1D model). If this eigenvalue has an imaginary part (and complex conjugate pair since the data matrix is real), we will track the positive frequency counterpart. Given the parameters in table 3, the model eigenvalues and DMD spectrum, represented as $\lambda$, are scaled with the inverse of nondimensional convective time units.

3.2 FSI DNS Discretization

The FSI DNS Eulerian mesh is uniform in $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ with grid spacing $\Delta x^{*}=\frac{\Delta x}{L}$. Because the parametric study spans a wide range of $Re_{L}$ and $\hat{h}$ values, with over 4000 simulations carried out, the $\Delta x^{*}$ is automatically determined by the most restrictive of three conditions: the grid Reynolds number $Re_{\Delta x}\leq 2$; the minimum number of grid elements in $\hat{h}$ is 20; the minimum number of elements on the beam surface is 160,

[TABLE]

The Lagrangian grid spacing is always $\Delta s^{*}=2\Delta x^{*}$, as suggested in (Goza et al., 2016). The time step size $\Delta t^{*}$ is determined by holding the $CFL=\frac{\Delta t^{*}}{\Delta x^{*}}=0.2$ for the $\Delta x^{*}$ that satisfies the criteria 31. These conditions were determined by trial and error to capture at least 200 time steps per beam oscillating cycle for all results. The grid Reynolds number chosen captures fluid advection and diffusion terms well. The resulting $\Delta x^{*}$, $\Delta t^{*}$ combination produces results within acceptable wall-time for the number of simulations ran in this study. We explore the effect of grid refinement from criteria 31 in our results next.

3.3 Grid Convergence and Effective Beam Thickness

The convergence of the FSI DNS beam dynamics as a function of the spatial discretization is explored through a grid refinement study. The parameters are $\hat{h}=0.05$, $\hat{h}^{2}Re_{L}=0.5,\ M^{*}=0.01$, and the corresponding $\Delta x=0.0025$. This parameter set is chosen specifically because it insures that the closure approximations introduced in section 2 are satisfied; thus the quasi-1D model is expected to accurately represent the fully-coupled system. $U^{*}$ is varied as the bifurcation parameter, while $M^{*}$, $\hat{h}$, and $\hat{h}^{2}Re_{L}$ are held constant. DMD eigenvalues are calculated for $\Delta x^{*}$, $\frac{\Delta x^{*}}{2}$, and $\frac{\Delta x^{*}}{3}$, where $\Delta x^{*}$ satisfies criteria 31. The DMD procedure is applied to the beam displacement results from simulations for each refined $\Delta x^{*}$ value.

Figure 3(a) shows the real and imaginary parts of the least-damped eigenvalue for all $U^{*}$ and $\Delta x^{*}$ grid values. The leading quasi-1D model eigenvalues are also shown. As the FSI DNS grid is refined, the DMD spectrum appears to be converging to the model eigenvalues: the real part of the DMD spectrum are monotonically moving toward the real part of the model eigenvalues; yet most notably, the imaginary part of the DMD spectrum is moving down for points where $U^{*}<5.27$ and up for $U^{*}>5.27$, approximating the quasi-1D curve shape.

To better understand the slow $\lambda$ convergence, we consider the beam thickness in light of the immersed boundary projection method. If the Eulerian mesh does not resolve the physical thickness of the beam, the immersed boundary produces an “effective” beam thickness that is proportional to $\Delta x^{*}$. This phenomenon is caused by the discrete kernel of the delta function being always positive, in addition to the unidirectional flow conditions on both sides of the beam (i.e. moving in the $+\mathbf{\hat{x}}$ direction). The IB projection method ensures that the no-slip condition is exactly enforced at each Lagrangian IB point. Once the discrete delta function smears it onto the flow, the Navier-Stokes equations are altered by the wall forcing for Eulerian points within the delta function kernel support. The flow outside its support, however, behaves as if the no-slip condition had been applied at the IB point locations. This effective thickness is not necessarily a physical quantity and cannot be systematically measured. An upper bound for this effective thickness is the number of Eulerian support points of delta function kernel (6 in our case) and expected to vary with $\Delta x^{*}$. By refining $\Delta x^{*}$, we are also decreasing the beam’s effective thickness, and slightly increasing the channel size. This effect would be more pronounced as $\hat{h}$ decreases, when the dynamics become a stronger function of $\hat{h}$ (Guo & Paidoussis, 2000, Shoele & Mittal, 2014, 2016, Cisonni et al., 2017).

Running the full parametric study with the finest grid in figure 3(a) presents issues due to computational resource restrictions and time constraints111$\Delta x^{*}$ grid results require approximately 32 hours of wall-time computation, while $\frac{\Delta x^{*}}{3}$ grid requires approximately 13 days.. These problems would be significantly amplified should the actual beam thickness need to be resolved by the Eulerian mesh.

Hence, we explore instead altering the channel thickness within the quasi-1D model to understand whether the slight change in channel width could explain the slow convergence. Figure 3(b) shows the comparison of the DMD spectrum for $\Delta x^{*}$ and the quasi-1D model with a corrected channel width,

[TABLE]

which gives an effective beam thickness slightly smaller than 3.5 Eulerian cells (or slightly larger than 1.5 Lagrangian cells). This result replicates the DMD spectrum extremely well, both in shape and in the instability boundary at $U^{*}_{\mathrm{cr}}=4.90$; it also confirms the sensitivity of the dynamics to $\hat{h}$, and lends credibility to this convergence hypothesis. Hence, we apply the grid criteria 31 in subsequent results presented in section 4.1. We assume that the effective beam thickness is captured by equation 32 when comparing between FSI DNS DMD spetra and quasi-1D eigenvalues for all grids in this study.

4 Results

4.1 Comparison of FSI DNS and Quasi-1D Model Results

We begin by assessing the validity of the quasi-1D model for parameter values compatible with the closure assumptions enumerated in section 2. We consider the specific cases listed in table  4. Similar to section 3.2, FSI DNS simulations are run for each set of parameters, with $U^{*}$ treated as the bifurcation parameter. For each parameter trio $\left[\hat{h},\ \hat{h}^{2}Re_{L},\ M^{*}\right]$, we find the pair $U^{*}_{\mathrm{cr}}$ and Im$\left[\lambda\right]_{\mathrm{cr}}$ as the critical values for the flutter boundary by linearly interpolating between the two nearest parametric mesh points between stable and unstable results. All cases in table 4 initially explored a $U^{*}$ range that spanned at least two orders of magnitude, with refined cases near the bifurcation point that yielded appropriate results for linear interpolation to find $U^{*}_{\mathrm{cr}}$ and $f_{s\mathrm{cr}}^{*}$ values. For all cases, $Re_{L}$ values were chosen such that $\hat{h}^{2}Re_{L}$ remains constant across the different $\hat{h}$ values tested, consistent with the fact that $Re_{L}$ only appears in the model equations through the group $\hat{h}^{2}Re_{L}$. All cases in table 4 fall within the laminar flow regime (i.e. parabolic flow profile), consistent with the definition in equations 12 and 13 for the linearized $f_{0}$, where $\hat{h}Re_{L}<1000$.

We first focus on cases 1 to 6, where $U^{*}_{\mathrm{cr}}$ is tracked as the mass ratio is finely incremented. The results are presented in figures 4 to 9. In addition to the flutter boundary, each figure shows the corresponding frequency of the unstable mode and beam mode shapes for selected mass ratios. Figure 4 shows the narrowest channel at $\hat{h}=0.025$. The quasi-1D model predicts the flutter boundary exceptionally well for the range of $M^{*}$ simulated. The corresponding beam mode shapes are given in figures 4(c) and 4(d). For the modes, the corresponding $U^{*}$ values were chosen from the nearest supercritical value from the available results for both the quasi-1D and FSI DNS computations. The beam shapes are qualitatively similar in the DMD and quasi-1D results. We attribute small oscillations superimposed onto the primary mode shape to the DMD data matrix having components from the impulsive start and body force perturbation at $t=0$. Mode switching is evident as $M^{*}$ increases from 0.01 to 0.3: not only is there an abrupt jump in frequency, but an additional effective node appears in both the real and imaginary parts of the modes shown. Compared to the orthogonal beam modes in vacuum, described in equation 27 we see a strong resemblance to beam mode two at $M^{*}=0.01$ and the third mode at $M^{*}=0.3$. The mode numbers refer to the eigenvalue index on table 2.

Similar results are obtained when the channel width is doubled to $\hat{h}=0.05$ holding $\hat{h}^{2}Re_{L}$ constant (case 2, figure 5). In particular, we confirm that quasi-1D model replicates the flutter boundary well for all $M^{*}$ values simulated. Quasi-1D results show multiple bifurcations at a single value of $M^{*}$. This occurs when, at a given $M^{*}$, a second eigenvalue crosses the stability boundary in the stable-to unstable direction as $U^{*}$ increases.

Figure 6 shows similar results for case 3 where $\hat{h}^{2}Re_{L}$ is raised to 1.5, holding $\hat{h}=0.05$. The flutter boundary has moved to higher values of $U^{*}$, indicating stabilization with increasing $Re_{L}$. This trend is reversed, however, when Reynolds numbers are high enough and inertial forces begin to dominate the dynamics. Quasi-1D and FSI DNS modes still mirror each other, but the mode switching inflection point has moved to a higher $M^{*}$ (relative to cases 1 and 2), so that $M^{*}=0.01$ and $M^{*}=0.3$ are primarily composed of beam mode two.

A further increase $\hat{h}^{2}Re_{L}=2.5$ yields results in figure 7 (case 4). The quasi-1D model marginally under-predicts the FSI DNS boundary, with the discrepancy increasing with $M^{*}$. Consequently, the mode switching $M^{*}$ point and higher mode boundary are under-predicted, yet $f^{*}_{s\mathrm{cr}}$ remains in close agreement. The mode shapes at $M^{*}=0.01$ also show some disagreement, with the DMD results appearing to have a larger contribution resembling orthogonal beam mode three.

Further increasing the channel width to $\hat{h}=0.125$ causes an increasing discrepancy between the DNS and quasi-1D model prediction, particularly at higher $M^{*}$ values. For case 5 (figure 8) the model accurately predicts flutter properties for $M^{*}<0.2$, but underestimates the higher mode boundary and critical frequencies for $M^{*}>0.4$. This is also evident in the mode shapes, as the $M^{*}=0.01$ modes are in agreement, but the $M^{*}=0.3$ modes are not. For case 6 (figure 9), $\hat{h}^{2}Re_{L}$ is increased to 1.25, where we also begin to see larger differences in the lower $M^{*}$ range. The quasi-1D model over-estimates $U^{*}_{\mathrm{cr}}$ relative to the FSI DNS simulations for $M^{*}<0.2$ but under-estimates it for $M^{*}>0.4$, also missing the mode switching $M^{*}$ point. $f^{*}_{s\mathrm{cr}}$ values remain well predicted through all values of $M^{*}$, however, as long as the system is in the correct branch.

Overall, the level of agreement between the quasi-1D model and the simulation results is consistent with the assumptions underpinning the model. When $\hat{h}$ and $\hat{h}^{2}Re_{L}\leq 0.5$ are sufficiently low, the quasi-1D model predicts critical flutter values well even for higher $M^{*}$ values, where condition $\delta^{\prime 2}\ll 1$ is less apt as the beam is oscillating in higher modes. Results where $\hat{h}^{2}Re_{L}=0.5$ is kept constant and $\hat{h}$ increases from 0.025 to 0.125 illustrate that, indeed, as $\hat{h}$ approaches $O(1)$, the quasi-1D boundary predictions worsen. Their results miss the mode switching $M^{*}$ along with the critical properties for $M^{*}\geq 0.4$. Yet even at $\hat{h}=0.125$, critical values for $M^{*}<0.2$ remain well approximated by the model. This indicates that as long as $\delta^{\prime 2}\ll 1$ holds, with only the lowest mode is considered, the model remains accurate. Also of note is that at lower $M^{*}$ values, inertial terms associated with the beam, i.e. moving channel walls, tend to dominate over the fluid inertia. The former is well captured with our closure, while the latter may not be due to the dependence of $\xi_{x}$ in $t$ and $x$. Cases where we hold $\hat{h}=0.05$ constant and increase $\hat{h}^{2}Re_{L}$ from 0.5 to 2.5 ($Re_{L}$ from 200 to 1000), show the quasi-1D model is not restricted to $\hat{h}^{2}Re_{L}\ll 1$ for accurate predictions. This is true even as the beam is excited at higher modes through increasing $M^{*}$. However, considering case 6, we see the $\hat{h}$ effect as $\delta^{\prime 2}$ increases through mode switching (i.e. increasing $M^{*}$). In summary, condition $\hat{h}\ll 1$ appears to be the most stringent restriction: as long as $\hat{h}$ remains small (i.e. $\hat{h}\leq 0.05$), the quasi-1D model predicts the flutter boundary for a wide range of $M^{*}$ and $\hat{h}^{2}Re_{L}$.

Cases 7 and 8 explore a larger range of ${\hat{h}}^{2}Re_{L}$ for a low mass ratio $M^{*}=0.01$ (heavy beam) where only the lowest beam mode is excited. Figure 10 presents results for case 7, where $\hat{h}=0.05$, with the flutter boundary now depicted as a function of $\hat{h}^{2}Re_{L}$, and representative mode shapes depicted for $\hat{h}^{2}Re_{L}=3$. At the upper limit of $\hat{h}^{2}Re_{L}$, the equivalent $Re_{L}$ reaches 1800. The quasi-1D model closely replicates the simulation results, with a small discrepancy in $U^{*}_{\mathrm{cr}}$ apparent when $\hat{h}^{2}Re_{L}>3$. Representative modes shapes also agree. In case 8 (figure 11) $\hat{h}$ is increased to 0.125 and $\hat{h}^{2}Re_{L}$ ranges from 0.1 to 9.5, with an equivalent $Re_{L}=608$ at the upper limit. Similar to case 7, quasi-1D results agree reasonably well over the range of $\hat{h}^{2}Re_{L}$, with a slight deviation near $\hat{h}^{2}Re_{L}=2$. Most notable, however, is the agreement at $\hat{h}^{2}Re_{L}>5$, as both results appear to asymptote to $U^{*}_{\mathrm{cr}}\sim\left(\hat{h}^{2}Re_{L}\right)^{1/2}$. Modes are shown in figure 11(c), with FSI DNS mode again showing some a larger component of the third orthogonal beam mode than the quasi-1D model predicts.

In summary, results show that the quasi-1D model well predicts the flutter boundary up to $\hat{h}\leq 0.125$ and $\hat{h}^{2}Re_{L}\leq 10$ for the heavy beam. The large values of $\hat{h}^{2}Re_{L}$ suggest that $\hat{h}^{2}Re_{L}\ll 1$ need not hold strongly in the limit $M^{*}\ll 1$. The impact of the viscous term on the flutter boundary though non-negligible, is well captured by the parabolic profile description of $\mathcal{N}_{x}$ and $F_{\mathrm{visc}x}$. This suggests that for low $M^{*}$ values, the inertia associated with the walls largely dominates over that associated with changing the velocity profile shape, i.e. $u$ largely remains locally parabolic as $\delta$ evolves in $t$ and $x$.

4.2 Quasi-1D Flutter Boundary Comparison to Inviscid Model

The ability of the quasi-1D model to predict the flutter stability behavior even when $\hat{h}^{2}Re_{L}$ is not small appears promising. We would like to compare quasi-1D predictions in cases where $\hat{h}^{2}Re_{L}\gg 1$ while $\hat{h}\ll 1$ to an inviscid flow solution. The inviscid model in (Shoele & Mittal, 2016) presents the flutter stability boundary as a function of $\hat{h}$, with $\hat{h}=0.05$ as the lowest provided value. Thus, taking the case for $\hat{h}=0.05$ and $\hat{h}^{2}Re_{L}=1.25$ as the starting point, we compare $U^{*}_{\mathrm{cr}}$ as a function of $M^{*}$ for $\hat{h}^{2}Re_{L}=[1.25-50]$ ($Re_{L}=[500-2\times 10^{4}]$) for the lowest frequency mode branch. Results are shown in figure 12. Notable trends appear as $\hat{h}^{2}Re_{L}$ increases: first, as $\hat{h}^{2}Re_{L}$ increases from 1.25 to 2.50, the system is stabilized as the stability boundary shifts upwards. Yet as $\hat{h}^{2}Re_{L}$ is further increased to 12.50 and thereafter the system is destabilized for $M^{*}>0.03$, with the boundary eventually disappearing for $0.03<M^{*}<0.2$ at $\hat{h}^{2}Re_{L}=50$. This means that no matter the stiffness of the system, the first mode is unstable if the beam is heavy enough, but not too heavy. The original stabilization trend for increasing $\hat{h}^{2}Re_{L}$ remains true, however, for $M^{*}<0.03$. Furthermore, as $\hat{h}^{2}Re_{L}$ becomes large, the quasi-1D flutter boundary appears to near the inviscid results acquired from (Shoele & Mittal, 2016) for $0.1<M^{*}<1$, with the mode switching inflection nearly matching over all $\hat{h}^{2}Re_{L}$ boundaries shown.

Shoele & Mittal (Shoele & Mittal, 2016) conjectured based on earlier DNS studies (Shoele & Mittal, 2014) that $Re_{L}\approx 200$ was sufficient to observe inviscid behavior for the $\hat{h}$ values in their study. Though this may be true for $\hat{h}>0.125$, beyond which the quasi-1D model’s predictions cannot be trusted, the inviscid behavior boundary for $\hat{h}=0.05$ is predicted at about $Re_{L}\approx 2\times 10^{4}$ (figure 12). Predictions of the quasi-1D model indicate that the choice of $Re_{L}$ for inviscid treatment is a strong function of $\hat{h}$.

In light of these results and those in section 4.1, we utilize the quasi-1D model to produce the complete flutter boundary for $M^{*}=0.01$ in figure 13 and $M^{*}=0.1$ in figure 14. Both figures only show the flutter boundary for the lowest frequency mode branch. Three trends become apparent from these plots: first, as $M^{*}$ increases, the lowest frequency mode is destabilized significantly; as $\hat{h}^{2}Re_{L}$ increases, the mode is stabilized. This stabilization is accelerated at higher $\hat{h}$, for values for $\hat{h}^{2}Re_{L}<10$. Figure 12 shows that opposite is true as $\hat{h}^{2}Re_{L}$ is increased further. Lastly, as $\hat{h}$ increases, the mode is stabilized, as has been pointed out previously (Shoele & Mittal, 2014).

Though $f^{*}_{s\mathrm{cr}}$ remains within a narrow range for the lowest frequency mode, an interesting pattern arises as $\hat{h}^{2}Re_{L}$ and $\hat{h}$ are varied. Bands of lower frequency appear alternating with higher frequency states in both $M^{*}$ values. This indicates that these parameters have an effect on the frequency response, but it is much less pronounced than their effects on the stability boundary as judged by $U^{*}_{\mathrm{cr}}$, where it should be noted that over the entire range of $\hat{h}$ and $\hat{h}^{2}Re_{L}$ considered, there is only about 10% variation in $f^{*}_{s\mathrm{cr}}$.

5 Conclusion

In this study, we derived a model for the fully coupled fluid-structure dynamics of a cantilever in channel flow, and predicted the flutter stability boundary in a constant, symmetric channel. The model takes a quasi-one-dimensional approach based on a slowly varying, thin-gap approximation ($h_{0}^{\prime 2}\ll 1$, $\delta^{\prime 2}\ll 1$, and $\frac{h}{L}=\hat{h}\ll 1$), and is closed by assuming that the velocity profile across the gap is unchanging, as would strictly be true in the lubrication limit where ${\hat{h}}^{2}Re_{L}\ll 1$. We examined the range of validity of the model by performing full numerical simulations of the two-dimensional incompressible Navier-Stokes equations coupled to the large-deflection Euler-Bernoulli beam.

Our results show that the model validity bounds are flexible: by enforcing $\hat{h}\ll 1$ strictly, the model is able to predict flutter properties at both large $\hat{h}^{2}Re_{L}$ and $M^{*}$ values. The agreement with FSI DNS values at high $\hat{h}^{2}Re_{L}$ is surprising; yet it is corroborated by the comparison between the inviscid model flutter boundary by Shoele & Mittal (Shoele & Mittal, 2016) as $\hat{h}^{2}Re_{L}\gg 1$. This suggests that as the viscous term diminishes, the modeled inertial term gives rise to similar balances between the quasi-1D model presented and the inviscid models at small $\hat{h}$ values. Notable also is that increasing $\hat{h}^{2}Re_{L}$ at small $\hat{h}$ values presents first a stabilizing effect for both heavy and light cantilevers, followed by a destabilizing effect for light cantilevers. Heavy cantilevers appear to always become more stable as $\hat{h}^{2}Re_{L}$ increases.

Furthermore, the model is shown as a reliable alternative to expensive fluid-structure numerical simulations, with the potential of handling various geometries of the flow passage, including asymmetric ones. It may serve as the first instance to explore the parametric space for device designs or understanding fluid-structure resonance phenomena in internal flows.

Acknowledgement

The authors would like to thank Andres Goza for his invaluable insights on fluid-structure physics and numerical methods that largely comprise the code base for the DNS simulations in this work. We would also like to acknowledge Bosch Energy Research Network (BERN) grant 13.01.CC17 and the NASA Jet Propulsion Laboratory for their support of this research.

Appendix A Two Dimensional Model Coefficients

The linear fluid-structure coupled operator in equation 30 is,

[TABLE]

where,

[TABLE]

exist in $\mathbb{R}^{n+1\ \times\ n+1}$, and

[TABLE]

exist in $\mathbb{R}^{n+1\ \times\ 1}$. The coefficients are defined in terms of the basis expansion as follows. Let

[TABLE]

where $g_{\alpha}$ is the basis function of index $\alpha$ defined in equation 25, $h_{e}$ is the equilibrium channel gap, $l$, $m$, and $n$ are their exponents, and $\mathcal{L}_{i}\left[\cdot\right]$ is the linear operator of order $i\in\mathbb{Z}:>-2$, defined as

[TABLE]

We also define the ratio,

[TABLE]

Given equations 36, 37, and 38, the following are the coefficients from equation 34: added fluid mass,

[TABLE]

fluid damping,

[TABLE]

fluid stiffness,

[TABLE]

and boundary flow rate forcing,

[TABLE]

The structure coefficients are

[TABLE]

The boundary flow-rate forcing components of matrices in 33 are,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For laminar flow, substituting $f_{0}$ from equation 12, we can show that

[TABLE]

and that viscous coefficients scales as

[TABLE]