The flow-induced fluttering motion of a flexible reed inside a heated channel is modeled numerically and used to investigate the relationship between the aeroelastic vibration of the reed and heat-transfer enhancement. An immersed boundary method is developed to solve the coupled flow-structure-thermal problem, and the simulations show that the vibrating reed significantly increases the mean heat flux through the channel, as well as the thermal performance, quantified in terms of the thermal enhancement factor. The effect of reed material properties on vibratory dynamics and heat transfer is studied. Changes in material properties produce a rich variety of vibratory behavior, and the thermal performance is found to depend more strongly on the reed inertia than its bending stiffness. The effects of both the Reynolds number and channel confinement are examined and it is found that the thermal performance is maximized when the reed creates large modulations in the boundary layer of the channel, while at the same time avoiding the creation of strong vortices.

## I. INTRODUCTION

An important challenge in the design of high-performance electronic systems, such as high-speed microprocessors and power converters is the demand for high-power-density heat dissipation to ensure proper operation and to prevent thermal damage. The current air-cooled heat sinks mostly operate at low characteristic channel Reynolds numbers *Re* < 3000 to minimize flow losses and cooling power. This, however, comes with a nearly laminar and therefore less effective rate of heat transfer. The efficiency and heat transfer rate of these systems can either be enhanced by unsteady fluctuations and mixing through increased Reynolds number^{1} or by introduction of secondary vortices through vortex generators.^{2–4} These and similar approaches, however, result in much higher flow losses and a reduced thermal enhancement factor (*TEF*) which is the ratio of the mean heat transfer of a channel with a reed to the mean heat transfer of a channel without a reed at the same pumping power.

It has been observed recently that direct actuation of small-scale motions within heat sink channels leads to significant improvement in their thermal performance. This improvement also continues at significantly lower flow rates for which the conventional channel type heat sinks stop being effective.^{5,6} Recently Hidalgo *et al.*^{7} showed that by inducing time-periodic small-scale vortices using the vibration of thin, planar piezoelectrially driven reeds, it is possible to increase the channel’s performance without affecting its characteristic Reynolds number. In particular, they reported that for a matched thermal dissipation conditions, the ratio of thermal power dissipated to the mechanical power invested in driving the flow known as the coefficient of performance (*CoP*) increases by 140%. Despite this superior performance compared to conventional air-side heat exchangers, the actuated reeds require power and wiring which add complexity of the system. An alternative is the use of aeroelastically driven, self-oscillating flexible reeds.^{8} By installing these flexible systems in a heat sink channel, a portion of the kinetic energy from the core flow within the channel is transformed to aeroelastic fluttering motion of reeds which then results in small-scale fluctuations near the heated surfaces and increase of heat transfer through the boundaries. While flexible reeds frequently are employed in different instruments such as clarinets, oboes, and saxophones,^{9–11} there are few studies exploring their potential as an alternative heat transfer enhancement technique.

The aeroelastic vibration of a reed in a channel also has similarities to the classical problem of the fluttering or flapping motion of a flag in a parallel, unconfined flow that has been the subject of a large number of experimental,^{12,13} numerical,^{14–16} and analytical studies.^{17–20} The fluttering instability of the flag results in a periodic or chaotic motion of the system. The temporal variation of the motion of the flag is from the competing effects of the fluid dynamic pressure forces and structural mass and flexural rigidity. However, all the above analyses have been for flags in unconfined flow and the effect of wall confinement on the dynamics of flag flutter remains mostly unexplored. In addition to flag flutter, the flow-induced vibration of plates/membranes also presents in many engineering and biomedical applications^{21,22} such as snoring,^{23} wing flutter,^{24} propulsion,^{25–28} and energy-harvesting.^{29} Moreover, hydroelastic motion of plates is important in some industrial cooling systems.^{30–32}

The purpose of the present study is to conduct a numerical investigation of the flow physics underlying the aeroelastic vibration of a reed in a channel and associated enhancement in convective heat transfer. We investigate the performance of this coupled flow-structure-thermal (CFST) system for a range of reed and flow parameters and explore the effects of oscillatory modes on the heat transfer enhancement. We show that proper tuning of these parameters can lead to significant enhancement in net heat transfer rates as well as the efficiency of heat transfer as assessed via the *TEF*.

In Secs. II–IV, we describe the two dimensional (2D) computational model used in this study. Three dimensional effects might ultimately found to be important in this configuration. However, given the tremendous computational expense of 3D simulations as well as the expansion in parameter space associated with 3D configurations, a necessary precursor is 2D models that allow us to explore a larger parameter space and provide a baseline for interpreting the results of the 3D simulations in the future.

The immersed boundary method (IBM) based approach that has been developed to solve the underlying CFST problem is described in detail. The computational modeling approach is validated with three benchmark problems, the details of which are included in the Appendix. Using this approach, we examine the effect of reed material properties (mass and bending rigidities) on the reed motion, the vortex dynamics and the thermal performance of the system. The roles of key flow parameters as well as channel aspect-ratio on the reed and flow are also examined and described in detail.

## II. PROBLEM FORMULATION

### A. Problem description

A modular setup of a flexible reed inside a rectangular channel is shown in Fig. 1(a). In this study, we assume that the spanwise width of the channel (*W*) is much larger than the channel height *H* and the reed length *L*. It is also assumed that the side-wall effects on the reed are negligible and the thermal performance of the system can be investigated by studying the fluid-structure-thermal interactions via a two dimensional computational model.

A schematic of the computational configuration is shown in Fig. 1(b). The reed is assumed to be fixed at the leading-edge and free at the trailing edge. A Lagrangian coordinate is chosen to be the arc-length along the reed starting at the leading-edge and ending at the trailing edge. The flow into the left boundary of the channel is assumed to be parabolic with a mean velocity of *U*. The reed is modeled as a two dimensional body with zero-thickness and assumed to be elastic and inextensible with its dynamics governed by

where *s* is a Lagrangian coordinate along the reed, $X= x , y $ describes the material position of the reed, $t= \u2202 X \u2202 s $ is the unit tangent vector along the coordinate *s*, and **n** is the unit normal vector. Furthermore, *m _{s}* is the excess mass per unit length, per unit width of the reed defined as

*m*=

_{s}*m*−

*ρh*, where

_{s}*m*is the mass per unit length, per unit width of the reed,

*ρ*is the density of the fluid, and

*h*denotes the cross sectional thickness of the reed. On the right-hand side of the equation,

_{s}*σ*is the tension and

*q*is the transverse stress. Finally,

**F**is the force density applied by the reed on the surrounding fluid and

**F**

_{c}denotes the force on the reed due to contact with the channel walls.

The inextensibility constraint of the reed is expressed as $ d d t \u2202 X \u2202 s =0$,^{33,34} which is reformulated as an equation for the tension *σ* as

Using the Euler-Bernoulli assumption^{35} for a slender inextensible reed (*M* = *k _{b}κ*, with

*M*the elastic moment along reed,

*k*the bending rigidity, and $\kappa =\u2212n\u22c5 \u2202 2 X \u2202 s 2 $ the local curvature along the reed), the transverse stress

_{b}*q*becomes

To account for the support mechanism at the leading edge, we have also assumed that the reed is fixed at its first 6% of its length; and the free boundary condition at the trailing edge is enforced by imposing ∂^{2}**X**/∂*s*^{2} = 0 and ∂^{3}**X**/∂*s*^{3} = 0.

Following previous studies on flexible filaments,^{36–38} we use the following nondimensional parameters to represent the inertia and bending rigidity of the reed:

where *ρ* is the density of the fluid. Physically, *M*^{*} represents the relative importance of the structural inertial force and *U*^{*} can be thought of as the ratio between time scale related to the reed’s natural oscillation in vacuum and the convective time scale of the fluid passing the reed.

The governing equations for incompressible flow and forced thermal convection in the entire fluid domain are expressed in Eulerian variables as

where in Eqs. (5) and (6), **x** is the location vector in the fluid domain, $u x , t $ is the fluid velocity, *p* is the dynamic pressure, *μ* is the dynamic viscosity, and **f** is the momentum forcing term to enforce the boundary conditions on the surface of the reed. In Eq. (7), $T x , t $ denotes the temperature, *c _{p}* the heat capacity,

*k*the coefficient of thermal conductivity, and

**q**the external heat source from the immersed structures.

Assuming *U* denotes the characteristic velocity, *L* the characteristic length, *T _{s}* the temperature of boundaries, and

*T*

_{0}the baseline temperature of the incoming flow temperature, Eqs. (5) and (6) can be written in their non-dimensional form. For convenience, in the following, we use same notations for the dimensionless quantities as their dimensional counterparts. Equations (5) and (6), then become

Other non-dimensional coefficients are defined in Table I along with their range of the values in this study.

Related equations . | Name . | Notation . | Range . |
---|---|---|---|

Fluid and thermal parameters Eqs. (5)–(7) | Reynolds number | $Re= \rho L U \mu $ | 50 − 800 |

Prandtl number | $Pr= c p \mu k $ | 1 | |

Reed parameters Eqs. (1)–(3) | Mass ratio | $ M * = \rho L m s $ | 10^{−1.5} − 10^{1.5} |

Reduced flow velocity | $ U * =UL m s / k b $ | 3 − 18 | |

Non-dimensional channel height | H^{*} = H/L | 0.25 − 2 |

Related equations . | Name . | Notation . | Range . |
---|---|---|---|

Fluid and thermal parameters Eqs. (5)–(7) | Reynolds number | $Re= \rho L U \mu $ | 50 − 800 |

Prandtl number | $Pr= c p \mu k $ | 1 | |

Reed parameters Eqs. (1)–(3) | Mass ratio | $ M * = \rho L m s $ | 10^{−1.5} − 10^{1.5} |

Reduced flow velocity | $ U * =UL m s / k b $ | 3 − 18 | |

Non-dimensional channel height | H^{*} = H/L | 0.25 − 2 |

Unless otherwise stated, the Reynolds number is *Re* = 400 and without lost generality, the Prandtl number (*Pr*) is chosen to be 1. Our sensitivity analysis of Prandtl number shows that the results are not appreciably different as long as *Pr* = *O*(1). These values are typical of the applications in electronic cooling that are the primary motivation for the current study. We also neglect natural convection effects and focus only on the forced convection regime which is more relevant for the applications of interest in electronic cooling. In addition, given the typical high temperature difference between *T _{s}* and

*T*

_{0}(on the order of O(10

^{2})) in air-cooled electronic systems, and the domination of heat convection and heat conduction, the increase in the thermal energy due to viscous dissipation is neglected in the thermal equation.

### B. Immersed boundary formulation

The Lagrangian (reed) and Eulerian (momentum and temperature) variables are related together by the Dirac delta function based immersed boundary formulation. For example, the Lagrangian fluid velocity **U** and Lagrangian fluid temperature Γ at the position of the beam are obtained by

where **Ω**_{F} is the fluid domain. A similar relation is employed to relate the Lagrangian and Eulerian descriptions of force densities and heat sources, i.e.,

where **Ω**_{S} is the domain of the immersed boundary (which in the case of the immersed reed is the contour along the reed).

The interaction between fluid variables and structural variables is calculated using a feedback law^{27,34,39} as

where *c*_{1} = *ραU*^{2}/*L* and *c*_{2} = *ρβU* with *α* and *β* are large penalty parameters to enforce no-slip and no-flux conditions along the reed. Furthermore, $V s , t $ and $Y s , t $ are the velocity and temperature of the reed, respectively.

Finally, in the present study, it is assumed that because of the fluid lubrication, the reed does not actually contact but rather interacts repulsively with the walls via the intervening fluid when the two are in close proximity.^{40} This short-range repulsive force can be formulated using the Dirac delta function as

where *S _{w}* is contour line along the walls and

*k*is the repulsive coefficient. This approach has successfully been used in previous studies such as collective dynamics of immersed particles,

_{c}^{41–43}the interaction between several flexible filaments,

^{34,44}and bacterial flagella bundling and tumbling in a viscous fluid.

^{45}

### C. Numerical implementation

A finite difference method is employed to solve the governing Eqs. (1)–(16). The Navier-Stokes and thermal convection Eqs. (5) and (6) are discretized on a non-staggered grid with **u**, *p*, and *T* are all defined at the cell centers.^{46,47} In addition, the face-center velocities, $ u \xaf $, are computed and updated separately to guarantee strong coupling between velocity and pressure, and to impose the incompressibility constraint (Eq. (5)) in each cell. Since there is a high degree of similarity between the numerical methods that are employed for solving thermal convection Eq. (7) and momentum Eq. (5), hereafter, we primarily explain the main steps in time-marching the momentum equation along with the key difference for the thermal equation.

A Crank-Nicolson scheme is used for implicitly advancing the diffusion terms, while all other terms, including convective terms, marched with the second-order Adams-Bashforth scheme. To summarize, the resulting formula of momentum equation is

where *H _{i}* represents the

*i*th component of convective terms,

*G*is the

_{i}*i*th component of the negative gradient operator divided by the fluid density, and

*D*is the diffusion operator divided by the density. To discretize the convective term

*H*in the momentum equation, the QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme is employed in which the upwind direction is determined by the magnitude of $ u \xaf i $.

_{i}^{48}The convective terms in the thermal convection equation is discretized with the SHARP (Simple High Accuracy Resolution Program) scheme.

^{49}The monotonic behavior of the SHARP scheme improves the stability of solution in the cases where there is a sharp temperature gradient in the flow domain.

^{49}The central difference is applied for the diffusion terms and a four-step fractional-step time advancement scheme proposed by Choi and Moin

^{50}is used to solve Eq. (18).

Staggered grids are used on the reed along *s*; the reed tension is defined on the centroid of the grid $ ( i \u02c6 ) $ nodes and other variables are defined on the node vertices (*i*). Equations (1) and (2) are then solved in the following sequence:

where *D _{s}* is the difference operator along

*s*,

**X**

^{*}is the predicted position of the reed at time

*n*+ 1 and is used for calculation of all other intermediate variables with superscripts “*”. The tension $ \sigma i \u02c6 $ is calculated at the intermediate step and is used to update the position of the reed,

**X**

^{n+1}. The interaction force between the reed and the fluid,

**F**

^{n}, is obtained explicitly from the Eq. (15) using the fluid and structural velocities up to the time step

*n*

where **V**^{n} is the velocity of the beam and **U**^{n} is the fluid velocity at the structural grid points along the reed, both calculated at time step *n*. The fluid velocity along the reed is evaluated by the smoothed delta function

with Δ*x* and Δ*y* are the grid sizes around the reed in *x* and *y* directions, respectively. *δ _{h}* is chosen to be a four-point smoothed delta function

^{51}

with

The other variables are transferred between Eulerian and Lagrangian descriptions similarly. The smoothened delta function shown in Eq. (17) is also used for regularizing the contact force.

The method has been tested and validated for several benchmark problems involving fluid-structure interaction as well as forced convection with stationary and moving boundaries, and the results of these benchmark tests are presented in the Appendix.

### D. Boundary conditions

Figure 2 shows the computational domain as well as the specified boundary conditions. The channel length is chosen to be 20*L*, which, as will be shown later, is sufficient for the effects of the reed to decay to nearly ambient conditions for all the cases simulated here. At the channel walls, we assume *T* = *T _{w}* with Δ

*T*=

_{w}*T*−

_{w}*T*

_{0}denoting the increase with respect to the temperature of the incoming flow. The origin of the coordinate system is placed at the leading-edge of the reed. The parabolic flow profile is imposed at the inflow boundary condition at

*x*= − 5

*L*and Neumann boundary condition is applied at the outflow boundary at

*x*= 15

*L*. A Neumann boundary condition is imposed for the pressure at all boundaries.

### E. Output variables

The amplitude of the reed oscillation is denoted by *A*/*L* and the reed vibration Strouhal number with *St _{A}* =

*f*/U in which

_{A}L*f*denotes the dominant oscillating frequency of the reed. Unless otherwise noted, all output fluxes are calculated at a cross section

_{A}*X*= 10

_{o}*L*downstream of the leading-edge of the reed. In particular, the net change in the thermal energy of the fluid is

where the integration is performed at *x* = *X _{o}*, as indicated in Fig. 2. The time-averaged value of the

*Q*is indicated with $ Q \xaf $, the magnitude of its temporal variation with $\Delta Q=s.d. Q $ (the standard deviation of

*Q*), and a thermal Strouhal number with

*St*=

_{Q}*f*/U where

_{Q}L*f*is the dominant frequency of the heat flux at

_{Q}*x*=

*X*. In addition, the net energy loss due to pressure drop between the inflow and output plane is calculated as

_{o} where *p*_{∞} is the ambient pressure.

The *TEF*, defined as the ratio of the mean heat transfer through the output plane of a channel with the reed $ ( Q \xaf ) $ to the mean heat transfer through the output plane of a channel without the reed $ ( Q \xaf 0 ) $ at an identical pumping power,^{52–55} is used to compare the performance of the reed enhanced channel to a baseline channel with no reed. This *TEF* is given by

where $ E \u0307 \xaf 0 $ is the mean energy loss between the inflow and output plane in the channel without a reed and a *TEF* value greater than 1 indicates higher heat transfer for the same input mechanical power. Further details about the derivation of *TEF* coefficient can be found in previous studies.^{52–55} In addition, the convective heat transfer on the channel walls is quantified via the Nusselts number (*Nu*) defined as

with *n* is the normal direction to the walls pointing into the fluid region.

## III. RESULTS

After describing the results of grid refinement studies, we examine the performance of the reed in a channel within the parametric space of *M*^{*} and *U*^{*} for the baseline configurations for which, *Re* = 400, *H*/*L* = 1, and *Pr* = 1.0. This is followed by a detailed investigation of the effects of mode switching on the oscillation and thermal performance of the reed, concentrating on the correlation between the characteristics of flow field and deformation of the reed. We then discuss the effects of Reynolds numbers and investigate how the height of the channel influences the thermal performance of the system.

### A. Grid refinement and IB parameters

The baseline grid has a 830 point non-uniform grid in the x-direction, with a higher grid density in the region around the reed (extending from −0.5*L* to 3*L* as shown with a shaded area in Fig. 2). The y-direction is discretized with a uniform grid and the number of grid points in the vertical direction depends on the channel height. Near the reed, the grid sizes are Δ*x* = Δ*y* = 0.0125*L*, and a uniform grid with 96 points is used to discretize the reed. Grid refinement studies have been performed to ensure that the simulation results are independent of the grid sizes for the flow as well as the Lagrangian grid along the reed. For example, by checking grid sizes of Δ*x* = Δ*y* = 0.0167*L*, 0.0125*L*, and 0.01*L* for a typical case (with *M*^{*} = 1 and *U*^{*} = 12), we predict the non-dimensional mean heat flux, $ Q \xaf / \rho c p U H \Delta T w $ computed at *x* = 10*L* to be 0.516, 0.484, and 0.476, respectively. Thus the change in this key quantity is less than 1.6%, as the grid size is reduced from 0.0125*L* to 0.01*L*. Other flow and structural parameters show similar convergence and their maximum error is always less than ∼2% as the fluid grid size is reduced below 0.0125*L*. By choosing the size of the grid on the reed to be smaller than the nearby fluid grid sizes, our convergence study shows that the results are insensitive to the resolution of the reed grid. Similar convergence behavior has also been observed previously.^{27,28}

A temporal sensitivity analysis was also carried out for selected cases where Δ*t* and IB (immersed boundary) parameters $ \alpha , \beta $ were varied until the results were found to be relatively independent of these parameters. The final time-step employed in all the simulations was Δ*t* = 0.0002*L*/*U* and we also chose *α* = − 3 × 10^{5} and *β* = − 2 × 10^{2}. The time-step chosen here is limited by the explicit treatment of the feedback law in Eqs. (15) and (16). This chosen time-step size provides at least 5000 time-steps per vibration cycles for the range of cases simulated here. The repulsive coefficient *k _{c}* is specified as 1000. This choice has been tested numerically to ensure that it is sufficiently large to enforce a reasonable contact constraint and small enough so as to maintain numerical stability. With this value, the reed does not cross the channel walls, even for the extreme cases (i.e., small

*M*

^{*}and large

*U*

^{*}), and the closest approach of the reed to the wall is about 0.01

*L*, which is of the same order as the grid size near the wall. We have also examined the effect of further increase in

*k*and have found negligible effects on thermal characteristics as well as the dynamic responses of the reed.

_{c}The inflow velocity amplitude and the wall temperature are initially increased exponentially as $ 1 \u2212 exp \u2212 t / t 0 $ to reduce the transient effects with *t*_{0} chosen to be 0.04*L*/*U*. The sensitivity of the numerical results to changes in *t*_{0} has also been investigated, and we find that a stationary state attained by the system beyond about *t* = 20*L*/*U* is independent of the startup transient. All statistics presented here have been computed between *t* = 25*L*/*U* and *t* = 40*L*/*U*.

### B. Reed dynamics and thermal performance in *M*^{*} − *U*^{*} space

Figures 3(a) and 3(b) show the peak-to-peak amplitude of the motion of the trailing edge (*A*/*L*) and the Strouhal number $ S t A $ of the reed within a range of 10^{−1.5} ≤ *M*^{*} ≤ 10^{1.5} (with Δlog_{10}*M*^{*} = 0.2) and 3 ≤ *U*^{*} ≤ 16 (with Δ*U*^{*} = 1). The width of microchannels (*H*) used in electronic or power cooling applications that are of interest here vary from a few millimeters to tens of centimeters with flow velocities up to 10 m/s.^{8} The reed thickness would be in the range of $O 10 \u2212 100 \mu m$ with lengths in the range of $O 10 \u2212 100 mm$. Therefore, *M*^{*} is expected to vary between 0.01 and 1.0 for a copper reed, 0.04 and 4.0 for an aluminum reed, and 0.1 and 10.0 for a mylar reed. Similarly, *U*^{*} is anticipated to be larger than 0.5, 0.3, and 1.3 for copper, aluminum, and mylar reeds, respectively. The range of *M*^{*} and *U*^{*} for the current study are therefore chosen to be inline with the applications of interest to us.

Figure 3 summarizes the reed response over the range of *M*^{*} and *U*^{*} studied here. Figure 3(a) indicates that when the reed inertia is very large (small *M*^{*}) or very small (large *M*^{*}), then irrespective of *U*^{*}, the reed is very stable and is not actuated by the fluid dynamic forces. As *M*^{*} becomes larger than a threshold value (which can be approximated by a line from $ log 10 M * , U * \u223c \u2212 1 . 1 , 16 $ to $ 0 , 5 $), a sharp increase is observed in *A*, which rapidly reaches its limit of *A* = 0.5*L* (in this case, this also equals 0.5*H*). In this regime, which extends almost to log_{10}*M*^{*} ∼ − 0.1, the reed undergoes large oscillations, but the presence of the channel walls on both sides restricts the motion to be at most 0.5*H*. Further increase in *M*^{*}, leads to an increase in the oscillation amplitude. This decrease is more rapid for smaller *U*^{*} where the bending rigidity is larger than for larger *U*^{*} where the reed is more flexible.

For intermediate values of *U*^{*}, *A*/*L* exhibits a double peak with variation in *M*^{*}. This phenomenon is related to mode switching of the reed oscillation and is more noticeable in Fig. 3(b) where the frequency of the reed is found to increase beyond log_{10}*M*^{*} ∼ 0.25.

For better illustration of this mode switching behavior, we plot in Fig. 4(a), the observed flapping mode shapes at four points marked in Fig. 3(b). By changing $ log 10 M * , U * $ from $ 0 . 1 , 12 $ (marked in Fig. 3(b) with ◊) to $ 0 . 3 , 13 $ (⧫), which are two neighboring points in the surface plot, there is a clear transition in the observed mode shapes. In the case of ⧫, not only does the reed show more deflection, its deformation is also more concentrated near the trailing edge. A similar transition is also observed between $ log 10 M * , U * = 0 . 7 , 14 ( \u2218 ) $ and $ 0 . 9 , 15 ( \u2022 ) $ in which the deformation of the reed changes from a two-node shape to a three-node shape. In Fig. 4(b), we have plotted the eigenmodes of the reed in vacuum and while there is some superficial similarity between the modes observed for the reed in the channel and the eigenmodes, the significant differences between the two indicate that the effect of the fluid on the reed clearly cannot be ignored.

Figures 5(a)–5(d) show the mean heat flux per unit length of the height of the channel at the output plane $ ( Q \xaf ) $, the mean flow power dissipated $ ( E \u0307 \xaf ) $, and the *TEF*. In the stable cases where the reed remains fixed, the presence of the reed only has a negligible effect on the heat flux (with $ Q \xaf / \rho c p U H \Delta T w =0.33$ compared to 0.32 in the channel without the reed). However, in the finite range of log_{10}*M*^{*} and *U*^{*} (almost coincident with high amplitude region in Fig. 3(a)), the channel with the flexible reed shows much higher heat flux at the output plane with the maximum $ Q \xaf / \rho c p U H \Delta T w =0.56$ occurring at $ log 10 M * , U * \u223c \u2212 0 . 3 , 8 $.

It is also observed that in a large area of the *M*^{*} − *U*^{*} space, especially for 10^{−0.9} ≤ *M*^{*} ≤ 10^{0.1}, the system delivers large $ Q \xaf $ with the value $ Q \xaf / \rho c p U H \Delta T w >0.55$, which is very close to the maximum value observed in the *M*^{*} − *U*^{*} parametric space. Except near the stability boundaries in the *M*^{*} − *U*^{*} space, the resultant heat flux from the system is not very sensitive to changes in *M*^{*} and *U*^{*}. In particular, $ Q \xaf $ shows very similar values for *M*^{*} ≤ 10^{0.1} when the reed is in the oscillatory state, and there is only a sharp change in the value of $ Q \xaf $ near the stability boundaries.

As *M*^{*} increases beyond 10^{0.1}, $ Q \xaf $ decreases toward the value corresponding to the channel with a stationary reed. Overall, there is a very good correlation between $ Q \xaf $ and the oscillation amplitude of the reed *A*, with the exception near $ log 10 M * , U * \u223c \u2212 0 . 3 , 9 $ where, despite the fact that *A* has very similar values to surrounding cases, $ Q \xaf $ is slightly higher.

Although the reed shows high $ Q \xaf $ in almost all the regions where *A* is near its peak value, *TEF* attains its optimal value near *M*^{*} = 1. The maximum of *TEF* is 1.16 compared to ∼0.9 for the stationary reed and this occurs at $ log 10 M * , U * \u223c \u2212 0 . 1 , 10 $. In fact, there is a range around *M*^{*} = 1 where the reed continues to deliver near-optimal *TEF* for a wide range of values of *U*^{*}. This suggests that the reed inertia is a factor that is more important than bending stiffness for achieving high *TEF*.

### C. Detailed analysis of effect of *U*^{*} and *M*^{*}

We now focus on the effect of *U*^{*} on the thermal performance of the channel with a reed. From the discussion in Part A, two representative *M*^{*} values have been chosen: *M*^{*} = 0.1, where the reed reaches its maximum flapping amplitude, and *M*^{*} = 1 where we observed the highest *TEF*. In the case of *M*^{*} = 1 (black line in Fig. 6(a)), when *U*^{*} is less than 4, the reed remains stationary. This is associated with the large bending stiffness of the reed in this regime, which makes the reed very stable. As *U*^{*} increases from 4 to ∼6, there is a sudden increase in the amplitude of the motion from 0 to ∼0.3*L*. A similar trend is also present for *M*^{*} = 0.1 (blue-diamond line in Fig. 6(a)), although in this case, the transition initiates at a larger *U*^{*} value (*U*^{*} ∼ 9) endures until *U*^{*} ∼ 13.5 and it shows a higher amplitude (*A* ∼ 0.5*L*). A sudden transition from a stationary state to an oscillatory mode has been observed previously in flapping flags inside unbounded as well as bounded domains.^{56,57}

In the case of *M*^{*} = 1, for a large interval of *U*^{*} from ∼6 to ∼12, the reed continues to deliver a near-constant amplitude of motion with a shallow maximum (*A*/*L* = 0.34) around *U*^{*} ∼ 9. The amplitude plot also indicates that the system undergoes a mode transition around *U*^{*} ∼ 13 which results in a slight increase in *A*.

The onset of the mode transition is clearer in Fig. 6(b), where there are rapid changes in the vibrational Strouhal numbers at the transition points. In the case of *M*^{*} = 0.1, the system displays only a single jump in the amplitude from *A*/*L* ∼0 to ∼0.5 and subsequent changes in the oscillatory modes are absent, at least up to *U*^{*} ∼ 18. Also in this case, the dominant frequency of oscillation is much smaller than what is present in *M*^{*} = 1, mostly associated with the increase of the inertia of the reed. For *M*^{*} = 1, the largest transition is around *U*^{*} ∼ 13 where the reed Strouhal number jumps from about 0.5 to 0.8. There is also a second transition around *U*^{*} ∼ 16 where the Strouhal number jumps to about 0.9. We note that there is no significant change in amplitude with this latter transition.

Figure 7(a) displays the dependence of $ Q \xaf $ on *U*^{*} and it is noted that the *M*^{*} = 1 case has high $ Q \xaf $ in a large range of *U*^{*} (*U*^{*} ≥ 7) studied here. On the other hand, the case with *M*^{*} = 0.1 shows enhancement in heat transfer only for *U*^{*} ≥ 15. The overall trend of $ Q \xaf $ is very similar to the curve of the oscillation amplitude, *A* (Fig. 7(a)), with the main difference being that here the transition happens at a slightly larger *U*^{*}.

The variation of *TEF* for *M*^{*} = 1 (Fig. 7(c)) shows that the range of *U*^{*} that exhibit high *TEF* is coincident with the range with the highest $ Q \xaf $ and the range in which the mechanical energy deficit still has not undergone a rapid increase (Fig. 7(b)). Moreover, the *TEF* in this range is quite insensitive to the changes in *U*^{*} (or the bending stiffness) and the reed continues to provide excellent thermal performance even if its bending stiffness deviates from the optimal value. This relative insensitivity to structural flexibility is a very useful property since it implies a higher tolerance in the manufacturing of these reeds. The case of *M*^{*} = 0.1, on the other hand, provides only a slight improvement in *TEF* and it is, therefore, not the optimal case for enhancing the heat transfer.

In Fig. 7(d), we plot the variation of the Strouhal number corresponding to the variation of heat flux at *x* = 10*L* versus *U*^{*} for the two cases. Superposed on these are the Strouhal numbers corresponding to the reed oscillation. It is interesting to note that for *M*^{*} = 1, the reed and the downstream flow lock-on to the same frequencies for oscillatory parts of the *U*^{*} range. In contrast, for the *M*^{*} = 0.1 case, the reed and flow are “detuned” in some ranges of *U*^{*}.

To examine this detuning further, in Figs. 8(a) and 8(b), we compare the power spectra of tail motion and heat flux for two selected cases chosen from Fig. 7(d). In the case with *M*^{*} = 1.0 and *U*^{*} = 17 (Fig. 8(a)), the power spectra of both the tail motion and the heat flux have only one peak, and the dominant heat flux frequency is equal to twice of the dominant vibrational frequency. For the heavier case with *M*^{*} = 0.1 and *U*^{*} = 17 (Fig. 8(b)) which corresponds to a “detuned” case, we can still see one dominant vibrational frequency, but the power spectrum of the heat flux consists of several peaks at smaller and larger frequencies, with the largest peak occurring at a frequency lower than the vibrational frequency. The origin of the richer spectra for this latter case can be traced to a more complex motion of the reed which is apparent in the phase-plane plot of the reed trailing-edge as shown in Figs. 8(b) and 8(c) for the two cases, respectively. The motion of the reed tip for the *M*^{*} = 0.1 and *U*^{*} = 17 case is highly stochastic due to interaction with vortices shed from upstream portions of the reed as well as interaction with the wall. This produces a complex vortex shedding with a relatively large range of vortex scales. These vortices interact with each other as well and also decay at different rates as they convect downstream leading to the generation of a more complex spectrum. It is interesting to note that the dominant frequency in the heat flux is 1/3 of the dominant vibrational frequency suggesting the presence of a subharmonic mechanism such as vortex merging.

In Fig. 9(a), we plot the distribution of the mean Nusselts number $ N u \xaf $ along the channel wall. This quantity is calculated by averaging the top and the bottom walls at the corresponding longitudinal positions along the channel, and for comparison, the Nusselts number of a channel without the reed is also included. There is a small increase in $ N u \xaf $ for *U*^{*} = 4 (for which the reed is nearly stationary) compared to the baseline channel, especially in the area where the reed is located. For larger values of *U*^{*} (*U*^{*} = 8, 12, and 16), the local $ N u \xaf $ shows the highest value in the near-wake region of the reed. For *U*^{*} = 16, the maximum occurs around ∼1.5*L* whereas for *U*^{*} = 8 and 12, the maximum of $ N u \xaf $ occurs at ∼2.5*L*. The transition in the position of maximum *Nu*^{1} is related to the mode switching of the reed to a higher frequency of oscillation.

The increase in the heat transfer from the channel walls however results in larger rate of increase in the mean cross-sectional temperature of the channel, especially for the cases with larger $ N u \xaf $ enhancement at the walls (Fig. 9(b)). The plot shows that the inclusion of the reed causes an immediate decrease of mean channel temperature after the leading edge of the reed reaching its minimum value at the trailing edge of the reed. This is followed by a rapid increase of the mean cross-sectional temperature after the trailing edge. As a consequence of higher mean temperature of flow at farther downstream points in the channel, the $ N u \xaf $ starts decreasing from its maximum value, eventually approaching zero as the temperature of the fluid equalizes with that of the walls. This has an important implication in selecting the length of channel in order to obtain the best heat transfer performance in the presence of the reed. In fact, by studying the dependence of *TEF* on *x*/*L* for the most efficient case with *U*^{*} = 8, it is observed that the maximum *TEF* is achieved if *x* ∼ 10*L*. This is the minimum length for obtaining the highest achievable *TEF* and can be considered as the optimal length of the channel for enhancing the heat transfer with flexible reeds for design purposes. Beyond this channel length, the temperature of the fluid inside the channel rises approximately to the temperature of the walls thereby reducing heat flux from the wall, but the pressure drop along the channel continues to increase with the increase of the channel length. Therefore, *TEF* does not improve with a further increase in the channel length.

To relate the heat enhancement capacity to the characteristics of the flow, in Fig. 10, we plot vorticity contour lines and isotherms in the near wake of the reed with *M*^{*} = 1 at one time-instance. The temperature of the flow is shown in grayscale, and snapshots of the corresponding reed deformation are also included.

For *U*^{*} = 4, the reed stays stationary, and the flow exhibits features of a steady wake (Fig. 10(a)). For *U*^{*} = 8 and 12, the reed oscillates with a large amplitude, and the amplitude of motion in this case grows rapidly from the leading edge to the trailing edge. These cases generate a well-organized vortex street where vortices of opposite rotation entrain hot fluid away from the wall. In both these cases, the counterclockwise (CCW) vortices are concentrated mostly near the top wall and the clockwise (CW) close to the bottom wall. As the reed becomes more flexible (*U*^{*} = 16), it vibrates in a new mode shape with a higher frequency (*St _{A}* = 0.75 compares to 0.5 for

*U*

^{*}= 12) and its deflection occurs predominantly near the trailing edge (Fig. 10(d)). There is also a clear transition in the vortex structures in the wake with stronger vortices that can penetrate into the flow in the middle of the channel. Creation of these large vortex structures with large scale vertical velocity perturbations causes larger energy deficit and lowers the

*TEF*in this case.

This effect becomes more apparent by examining the distribution of the viscous dissipation inside the channel. In Fig. 11, we show the instantaneous and mean contour plots of mean viscous dissipation (calculated as **τ**:∇**v** with **τ** being the stress tensor) in the channel for the cases with *U*^{*} = 12 and *U*^{*} = 16. Comparison of the contours of instantaneous viscous dissipation in Figs. 11(a) and 11(c) shows that while in the case of *U*^{*} = 12, the regions with high dissipation coincides with the locations of large heat flux, in the case of *U*^{*} = 16, there is significant increase in the area of large **τ**:∇**v** in the middle of the channel and between successive vortices resulting in a larger energy deficit in this case. Comparison of the mean viscous dissipation for these cases indicates that there is a large increase in **τ**:∇**v** (dark contour color in Figs. 11(c) and 11(d)) in both the near wake of the reed, and far downstream for *U*^{*} = 16 compared to *U*^{*} = 12.

The above analysis indicates that the more flexible reed creates stronger vortices, but these in turn create significant losses, thereby significantly increasing $ E \u0307 $. The current results therefore suggest that effective heat enhancement mechanism should center on modulating the boundary layer of the channel in a constructive way rather than imposing strong distinct flow features that bring a disproportionate rise in the mechanical loss of the system.

### D. Effect of Reynolds number and channel height

The simulations described in Secs. III A, III B, and III C have all been conducted at a fixed Reynolds number (*Re* = 400) and channel height *H*/*L* = 1. In this section, we describe the effect of these two parameters on the reed dynamics and thermal performance. For this study, we choose three representative cases: Case 1: a suboptimal case with *M*^{*} = 1 and *U*^{*} = 8; Case 2: a case with the highest observed *TEF* with *M*^{*} = 1 and *U*^{*} = 12; and Case 3: the case which delivers the maximum heat transfer flux at the output plane with *M*^{*} = 1 and *U*^{*} = 14. In each of Secs. III D 1 and III D 2, we alter only one parameter (either the Reynolds number or the normalized channel height) while all other parameters are assumed to be similar to their original values.

#### 1. Effect of Reynolds number

Changes in the size of the channel and/or the flow rate have a direct influence on the Reynolds number at a fixed Prandtl number of *Pr* = 1, and we explore the effect of this parameter on the dynamics and heat transfer of the system. In particular, Reynolds number is expected to have a significant effect on the vortex shedding and dynamics in the wake of the reed and this in turn should impact the heat transfer.

The amplitude of the motion, *A*, of the reed initially increases with the increase in *Re* but reached a nearly constant value beyond *Re* ∼ 200 (Fig. 12). This behavior is very similar for all *U*^{*} cases. Fig. 13 shows the deformation patterns for the *U*^{*} = 12 case for various Reynolds numbers, and we note that between a *Re* of 100 and 600, the deformation pattern is essentially unchanged. At *Re* = 50, the pattern is also the same, albeit with a smaller amplitude. At *Re* = 800, the reed exhibits a transition to a different mode shape with much smaller amplitude. Thus, mode transitions are not only driven by changes in reed properties but also by flow parameters.

In Figs.14(a)-14(c), we plot $ Q \xaf $, $ E \u0307 \xaf $, and *TEF* as a function of *Re*. In addition, for comparison, the variation of $ Q \xaf $ and $ E \u0307 \xaf $ for the channel without a reed is also presented (Figs.14(a) and 14(b)). In all three cases with the reed as well as the case without the reed, $ Q \xaf $ reduces monotonically with increasing *Re*. This decrease is attributed to the reduction in the rate of heat diffusion at the walls with increasing *Re*, which reduces the effective heat flux from the walls. Figure 14(c) shows the variation of *TEF* for these cases, and this plot shows a rapid initial increase in *TEF* for all cases. The rate of this increase however decreases at larger *Re* values with the *TEF* approaching an asymptotic value. More flexible reeds (higher *U*^{*}) reach their maximum *TEF* values at the smaller *Re* numbers. The *U*^{*} = 14 has the lowest *TEF* among all the cases, while the *TEF* for *U*^{*} = 8 and *U*^{*} = 12 shows a very similar trend with *Re* up to *Re* = 600. The maximum gain in thermal performance as quantified by *TEF* reaches up to about 25% for the cases studied here and these occur at the higher Reynolds numbers.

The snapshots of the flow near the reed for the cases with *M*^{*} = 1 and *U*^{*} = 12, but different *Re* numbers are shown in Fig. 15. At *Re* = 50, the changes in the flow due to the reed oscillation dissipate rapidly, and the flow returns back to the baseline channel flow. Since diffusion plays an important role in this case, the temperature of the flow at the channel-center rises rapidly from the inflow, thereby reducing the effective temperature difference and the heat transfer between the near-wall flow and the wall. At higher Reynolds numbers, it can be seen that vorticity produced by the reed dominated the flow field in the wake of the reed, and the heat transfer is mostly concentrated in near-wall flow entrained by these vortices. Also, there is noticeable vortex separation at the leading edge of the reed at *Re* = 800 case with *M*^{*} = 1 and *U*^{*} = 12 in which the reed shows a large oscillation. The leading edge vortex separation seems to reduce the strength of the main vortexes at the trailing edge, and therefore weakens the heat transfer at the walls, resulting in small decrease in *TEF* is this case (Fig. 14(c)). These phenomenological features explain to some extent, the thermal performance of the reed at different Reynolds numbers.

#### 2. Effects of channel height

In this subsection, we describe how the height of the channel relative to the length of the reed affects the reed motion and the heat transfer performance of the system. Figure 16 shows the variation in normalized reed amplitude with the channel height, and the plots show a different trend for *U*^{*} = 8 as compared to *U*^{*} = 12 and 14. This difference is due to the fact that in the case of *U*^{*} = 8, the reed bending stiffness is large enough to prevent it from contacting the walls for all cases except *H*/*L* = 0.2. On the other hand, in the cases of *U*^{*} = 12 and 14, when *H*/*L* ≤ 0.75, the reed contacts with the walls and the contact force limits its deflection. As *H*/*L* increases, *A* initially shows a small decrease and then rises again gradually toward its asymptotic value at higher *H*/*L*.

In Fig.17(a), we plot the dependence of the mean heat flux per unit height of the channel at output plane as a function of *H*/*L*. In general, $ Q \xaf / \rho c p U H \Delta T w $ reduces monotonically with increasing *H*/*L*. The different values of *U*^{*} show very similar trends in $ Q \xaf / \rho c p U H \Delta T w $, with the main difference occurring at the intermediate values of *H*/*L*. In this case, more flexible reeds (*U*^{*} = 12 and 14) deliver slightly higher values. When we compare $ Q \xaf / \rho c p U H \Delta T w $ between the channel with and without the reed, we can see that at *H*/*L* between 0.75 and 1, the channel with the reed substantially increases the heat flux compared to the channel without the reed. For example, this increase for *U*^{*} = 12 at *H*/*L* = 1 is ∼51% (changing from $ Q \xaf / \rho c p U H \Delta T w =0.32$ for the channel without the reed to 0.49 for the one with the flexible reed). As *H*/*L* increases beyond *H*/*L* = 1, the mean heat transfer per unit height of the channel with a reed monotonically decreases similarly to the channel without reed. However, these cases still deliver higher mean heat flux.

The *TEF* generally shows a similar trend for various *U*^{*} for *H*/*L* < 0.75 and *H*/*L* > 1.25 where it rises slowly with *H*/*L* and asymptotes to a value between 1.1 and 1.2. For 0.75 < *H*/*L* < 1.25, the *U*^{*} = 8 and 12 cases show a maximum whereas *U*^{*} = 14 shows a minimum. The main reason behind the smaller *TEF* for *U*^{*} = 14 at *H*/*L* ∼ 1 is the appearance of a higher frequency mode shape in this case.

Figure 18 shows the comparison between the flow fields of the cases with *M*^{*} = 1 and *U*^{*} = 12 and with two different channel heights, *H*/*L* = 0.5 and 1.5. When *H* is small (here the case with *H*/*L* = 0.5, Fig.18(a)), the wall interferes with the dynamics of the reed causing the vortical structures of the reed to be less coherent. In addition, the flow perturbations by the reed in this case interact rapidly with the channel walls and disappear shortly after the trailing edge in the wake, resulting in a very small enhancement in the heat flux. The presence of the reed in this case, however, imposes subsequent pressure loss and therefore larger viscous dissipation along the channel. The resultant *TEF* in this case is lower than 1.0 showing that the system performance is in fact worse than the baseline channel without a reed.

In the cases with larger *H* (*H*/*L* = 1.5 in Fig. 18(b)), the reed oscillates with large amplitude and initiates large vortical structures which reorganize themselves in a reverse Karman vortex street. The resultant flow in the channel exhibits strong entrainment of hot fluid from the channel wall and increases the mean heat flux and *TEF* of the channel.

## IV. CONCLUSIONS

Computational modeling studies have been conducted to examine the dynamics of a flexible reed in a heated channel and the effect of flow-induced reed oscillations on heat transfer. In the two-dimensional model employed here, the reed is modeled as an inextensible plate with fixed boundary condition at the leading edge and free boundary condition at the trailing edge. An immersed boundary approach on collocated grids has been employed to model the fluid. The model has been verified through testing with several benchmark problems.

A parametric study of the effect of reed inertia (represented by *M*^{*}) and bending stiffness (represented by *U*^{*}) is presented, and the effect of these parameters on vibration amplitude, frequency, and mode transition is described. In addition, regimes with large increase in the heat flux and higher thermal performance (characterized by *TEF*) are distinguished and correlated with the reed vibration and resulting vortex dynamics. Compared with the channel without the reed, significant increase in mean heat flux is observed in a large area of the *M*^{*} − *U*^{*} parametric space with a very good correlation with the amplitude of the vibration of the reed. The *TEF* of the system shows a peak at *M*^{*} ∼ 1, and it is found that *TEF* is more sensitive to the inertia of the reed than its bending stiffness. This suggests that as long as the mass properties of the reed are adjusted to maintain *M*^{*} ∼ 1, the system continues performing with near maximum *TEF* over a large range of structural bending stiffness.

The simulations indicate that the thermal performance of the system is maximized if the wake of the reed interacts with the channel walls in a way that creates large modulations in the boundary layer. On the other hand, conditions that produce vortex streets comprising of nearly circular vortices that occupy a significant proportion of the channel are not optimal for heat transfer. The simulations also show that the beneficial effects of the vibrating reed on heat transfer only extend a finite distance downstream of the reed. In terms of the design of heat sinks, this would indicate an optimal channel length for a given reed. The effect of Reynolds number has also been studied and it is shown that thermal performance improves with Reynolds number, at least in the range of *Re* that has been simulated here. Finally, the comparison between the performance of the system with different channel heights demonstrates that the maximum $ Q \xaf $ increase compared to the channel without the reed happens when *H* is approximately equal to *L*. This also coincides with the situation in which the *TEF* shows its maximum value. This along with other parametric studies performed here provides potential configurations of a channel with a self-oscillatory reed that can perform optimally.

Although the two-dimensional rendition used in the current study reveals important key mechanisms that are contributing to the heat transfer enhancement of the system, intrinsic three dimensionality of the flow, as well as extrinsic 3D effects associated with the channel sidewalls, could very well alter the system performance in ways that is not revealed by the current study. These effects are being explored currently and will be presented in the future.

## Acknowledgments

This study was supported by the Air Force Office of Scientific Research under Grant No. FA9550-13-0053 monitored by Dr. Douglas Smith. We would also like to acknowledge discussion with Dr. Ari Glezer (Georgia Tech.) and Dr. Silas Alben (U. Michigan) on various aspects of this problem.

### APPENDIX A: VALIDATION FOR FLOW PAST A CIRCULAR CYLINDER WITH HEAT TRANSFER

Flow pass a stationary cylinder with forced convective heat transfer is simulated. The Reynolds number is $Re= \rho D U \mu =20$, Prandtl number is chosen to be $Pr= c p \mu k =1$, and Grashof number is $Gr= g \beta T s \u2212 T 0 D 3 \rho 3 \mu 2 =0$. This problem has been solved previously by Chang and Finlayson,^{59} Gan *et al.*,^{60} Wang *et al.*,^{61} and by Turek and Hron^{62} using different numerical approaches. The size of computation domain is chosen to be 28*D* × 20*D*, with *D* being the diameter of the cylinder. The computational domain and boundary conditions are illustrated in Fig.19(a). Uniform flow with zero temperature is imposed at the inflow and zero derivative Neumann boundary condition is applied at the outflow. In all boundaries in y direction, the far-field boundary condition is applied through Neumann boundary condition. At the surface of the cylinder, a uniform temperature is applied. The simulation time step is Δ*t* = 0.0002*D*/*U*; the grid size near the cylinder is uniform with $\Delta x=\Delta y= 1 96 D$; and the number of Lagrangian points along the cylinder is *N _{s}* = 323. Finally, the penalty parameters are

*α*= − 10

^{5}and

*β*= − 10

^{2}.

### APPENDIX B: FLOW-INDUCED VIBRATION OF AN ELASTIC BEAM BEHIND A STATIONARY CYLINDER

We validate the fluid structure coupling model with flow-induced vibration of a flexible beam behind a fixed cylinder in a confined channel, originally proposed by Dennis *et al.*^{58} as a benchmark validation for two-dimensional fluid structure interaction models. The model setup is shown in Fig. 20(a).

The cylinder is placed at 2*D* from both the inlet and bottom wall in a rectangular domain with longitudinal and transversal sizes of, respectively, *L _{x}* = 11

*D*and

*L*= 4.1

_{y}*D*with

*D*being the diameter of the cylinder. The parabolic velocity profile with mean velocity $ U \xaf $ is imposed at the inlet boundary condition, no-slip wall boundary conditions at the top and bottom boundaries, and Neumann boundary condition at the outlet. The computational domain is discretized by 320 × 129 nonuniform grids with Δ

*x*

_{min}= Δ

*y*

_{min}= 0.018

*D*, wherein high resolution is provided in the region where the immersed structures are expected to present. A 248-point uniform Lagrangian grid is used along the beam, and the time step is set to $\Delta t=0.005D/ U \xaf $.

The beam is inextensible with bending stiffness $ k b =0.933\rho D 3 U \xaf 2 $ which is equivalent to assuming Young modules $E=1.4\xd71 0 \u2212 3 \rho U \xaf 2 $ for a beam with thickness *h* = 0.2*D* originally used by Dennis *et al.*^{58} In the current test, We consider two cases with (1) *ρ _{s}*/

*ρ*= 10 and $Re=\rho U \xaf D/\mu =100$, and (2)

*ρ*/

_{s}*ρ*= 1 and

*Re*= 200. The time history of the tip deflection of the beam of case-1 is illustrated in Fig. 20(b) along with the previous results by Dennis

*et al.*

^{58}and Bhardwaj and Mittal,

^{63}and Table II shows quantitative comparisons of key quantities. Overall, our predictions agree very well with others, both qualitatively and quantitatively, thereby providing confidence in the fidelity of the solver.

Cases . | Sources . | A/_{m}D
. | St
. | C
. _{D} |
---|---|---|---|---|

ρ/_{s}ρ = 10, Re = 100 | Present result | 0.89 | 0.19 | 4.19 |

Dennis et al.^{58} | 0.83 | 0.19 | 4.13 | |

Bhardwaj and Mittal^{63} | 0.92 | 0.19 | 3.56 | |

Tian et al.^{64} | 0.78 | 0.19 | 4.11 | |

ρ/_{s}ρ = 1, Re = 200 | Present result | 0.44 | 0.27 | 2.48 |

Dennis et al.^{58} | 0.36 | 0.26 | 2.30 | |

Bhardwaj and Mittal^{63} | 0.41 | 0.28 | 2.20 | |

Tian et al.^{64} | 0.32 | 0.29 | 2.16 |

Cases . | Sources . | A/_{m}D
. | St
. | C
. _{D} |
---|---|---|---|---|

ρ/_{s}ρ = 10, Re = 100 | Present result | 0.89 | 0.19 | 4.19 |

Dennis et al.^{58} | 0.83 | 0.19 | 4.13 | |

Bhardwaj and Mittal^{63} | 0.92 | 0.19 | 3.56 | |

Tian et al.^{64} | 0.78 | 0.19 | 4.11 | |

ρ/_{s}ρ = 1, Re = 200 | Present result | 0.44 | 0.27 | 2.48 |

Dennis et al.^{58} | 0.36 | 0.26 | 2.30 | |

Bhardwaj and Mittal^{63} | 0.41 | 0.28 | 2.20 | |

Tian et al.^{64} | 0.32 | 0.29 | 2.16 |