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.

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 number1 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 applications21,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. IIIV, 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.

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.

FIG. 1.

A schematic of the flexible reed model with the length L, and 2D computational model of the system chosen for study.

FIG. 1.

A schematic of the flexible reed model with the length L, and 2D computational model of the system chosen for study.

Close modal

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

m s 2 X t 2 = s σ t + q n F + F c ,
(1)

where s is a Lagrangian coordinate along the reed, X = x , y describes the material position of the reed, t = X s is the unit tangent vector along the coordinate s, and n is the unit normal vector. Furthermore, ms is the excess mass per unit length, per unit width of the reed defined as ms = mρhs, where m is the mass per unit length, per unit width of the reed, ρ is the density of the fluid, and hs denotes the cross sectional thickness of the reed. On the right-hand side of the equation, σ is the tension and q is the transverse stress. Finally, F is the force density applied by the reed on the surrounding fluid and Fc 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 X s = 0 ,33,34 which is reformulated as an equation for the tension σ as

2 s 2 σ t t = m s 2 2 t 2 t t m s t t t t s q n s F + F c .
(2)

Using the Euler-Bernoulli assumption35 for a slender inextensible reed (M = kbκ, with M the elastic moment along reed, kb the bending rigidity, and κ = n 2 X s 2 the local curvature along the reed), the transverse stress q becomes

q = M s = s k b 2 X s 2 n .
(3)

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 ∂2X/∂s2 = 0 and ∂3X/∂s3 = 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:

M * = ρ L m s , U * = U L m s k b ,
(4)

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

ρ u t + uu = p + μ 2 u + f ,
(5)
u = 0 ,
(6)
ρ c p T t + u T = k 2 T + q ,
(7)

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, cp 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, Ts the temperature of boundaries, and T0 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

u t + uu = p + 1 R e 2 u + f ,
(8)
u = 0 ,
(9)
T t + u T = 1 R e P r 2 T + q .
(10)

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

TABLE I.

The key independent nondimensional parameters in the thermal, fluid, and structural equations as well as their range of values in this study.

Related equations Name Notation Range
Fluid and thermal parameters Eqs. (5)(7)  Reynolds number  R e = ρ L U μ   50 − 800 
  Prandtl number  P r = c p μ k  
Reed parameters Eqs. (1)(3)  Mass ratio  M * = ρ L m s   10−1.5 − 101.5 
  Reduced flow velocity  U * = U L 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  R e = ρ L U μ   50 − 800 
  Prandtl number  P r = c p μ k  
Reed parameters Eqs. (1)(3)  Mass ratio  M * = ρ L m s   10−1.5 − 101.5 
  Reduced flow velocity  U * = U L 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 Ts and T0 (on the order of O(102)) 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.

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

U s , t = Ω F u x , t δ x X s , t d x ,
(11)
Γ s , t = Ω F T x , t δ x X s , t d x ,
(12)

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

f x , t = Ω S F s , t δ x X s , t d s ,
(13)
q x , t = Ω S Q s , t δ x X s , t d s ,
(14)

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 law27,34,39 as

F s , t = c 1 0 t U s , t V s , t d τ + c 2 U s , t V s , t ,
(15)
Q s , t = c 1 0 t Γ s , t Y s , t d τ + c 2 Γ s , t Y s , t ,
(16)

where c1 = ραU2/L and c2 = ρβ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

F c s , t = k c S w δ x ( s ) X ( s , t ) x ( s ) X ( s , t ) | x ( s ) X ( s , t ) | d s ,
(17)

where Sw is contour line along the walls and kc is the repulsive coefficient. This approach has successfully been used in previous studies such as collective dynamics of immersed particles,41–43 the interaction between several flexible filaments,34,44 and bacterial flagella bundling and tumbling in a viscous fluid.45 

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

u i n + 1 u i n Δ t = 3 2 H i n 1 2 H i n 1 + G i p n + 1 + 1 2 D u i n + 1 u i n + f i n ,
(18)

where Hi represents the ith component of convective terms, Gi is the ith 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 Hi 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 ¯ 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 Moin50 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 ˆ ) nodes and other variables are defined on the node vertices (i). Equations (1) and (2) are then solved in the following sequence:

X * = 2 X n X n 1 ,
(19)
D s D s σ i ˆ t * t i ˆ * = m s 2 1 2 t t i ˆ n + t t i ˆ n 1 Δ t 2 m s D s V D s V i ˆ n , D s D s q n * F n + F c * t i ˆ *
(20)
m s X i n + 1 X i * Δ t 2 D s σ i ˆ D X i n + 1 + q n n + 1 i = F n + F c * ,
(21)

where Ds 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 σ i ˆ is calculated at the intermediate step and is used to update the position of the reed, Xn+1. The interaction force between the reed and the fluid, Fn, is obtained explicitly from the Eq. (15) using the fluid and structural velocities up to the time step n

F n = c 1 k = 1 n U k V k Δ t k + c 2 U n V n ,
(22)

where Vn is the velocity of the beam and Un 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

U n = x u n x δ h x X n Δ x Δ y
(23)

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 function51 

δ h x = 1 Δ x Δ y φ x Δ x φ y Δ y
(24)

with

φ r = 1 / 8 3 2 r + 1 + 4 r 4 r 2 0 r < 1 1 / 8 5 2 r + 7 + 12 r 4 r 2 1 r < 2 0 2 r .
(25)

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.

Figure 2 shows the computational domain as well as the specified boundary conditions. The channel length is chosen to be 20L, 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 = Tw with ΔTw = TwT0 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 = − 5L and Neumann boundary condition is applied at the outflow boundary at x = 15L. A Neumann boundary condition is imposed for the pressure at all boundaries.

FIG. 2.

The computational domain and assumed boundary conditions (not to scale).

FIG. 2.

The computational domain and assumed boundary conditions (not to scale).

Close modal

The amplitude of the reed oscillation is denoted by A/L and the reed vibration Strouhal number with StA = fAL/U in which fA denotes the dominant oscillating frequency of the reed. Unless otherwise noted, all output fluxes are calculated at a cross section Xo = 10L downstream of the leading-edge of the reed. In particular, the net change in the thermal energy of the fluid is

Q = 0 H ρ C p u T T 0 d y ,
(26)

where the integration is performed at x = Xo, as indicated in Fig. 2. The time-averaged value of the Q is indicated with Q ¯ , the magnitude of its temporal variation with Δ Q = s.d. Q (the standard deviation of Q), and a thermal Strouhal number with StQ = fQL/U where fQ is the dominant frequency of the heat flux at x = Xo. In addition, the net energy loss due to pressure drop between the inflow and output plane is calculated as

E ̇ = 0 H ( p p ) u d y x = 0 0 H ( p p ) u d y X 0 ,
(27)

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 ¯ ) to the mean heat transfer through the output plane of a channel without the reed ( Q ¯ 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

T E F = Q ¯ Q ¯ 0 × E ̇ ¯ 0 E ̇ ¯ 1 / 3 ,
(28)

where E ̇ ¯ 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

N u = T w T / n s Δ T w / L
(29)

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

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.

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.5L to 3L 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.0125L, 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.0167L, 0.0125L, and 0.01L for a typical case (with M* = 1 and U* = 12), we predict the non-dimensional mean heat flux, Q ¯ / ρ c p U H Δ T w computed at x = 10L 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.0125L to 0.01L. 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.0125L. 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 α , β 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.0002L/U and we also chose α = − 3 × 105 and β = − 2 × 102. 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 kc 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.01L, which is of the same order as the grid size near the wall. We have also examined the effect of further increase in kc and have found negligible effects on thermal characteristics as well as the dynamic responses of the reed.

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

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.5M* ≤ 101.5 (with Δlog10M* = 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 100 μ m with lengths in the range of O 10 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.

FIG. 3.

Effect of M* and U* on the dynamic response of the reed. (a) Variation of maximum trailing edge (TE) peak-to-peak amplitude and (b) vibrational Strouhal associated with the TE frequency.

FIG. 3.

Effect of M* and U* on the dynamic response of the reed. (a) Variation of maximum trailing edge (TE) peak-to-peak amplitude and (b) vibrational Strouhal associated with the TE frequency.

Close modal

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 * 1 . 1 , 16 to 0 , 5 ), a sharp increase is observed in A, which rapidly reaches its limit of A = 0.5L (in this case, this also equals 0.5H). In this regime, which extends almost to log10M* ∼ − 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.5H. 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 log10M* ∼ 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 ( ) and 0 . 9 , 15 ( ) 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.

FIG. 4.

(a) The flapping mode observed at the four points marked in Fig. 3(b) near the mode switching regions and (b) the first 4 mode shapes of a linearized reed in the vacuum and their corresponding natural frequencies.

FIG. 4.

(a) The flapping mode observed at the four points marked in Fig. 3(b) near the mode switching regions and (b) the first 4 mode shapes of a linearized reed in the vacuum and their corresponding natural frequencies.

Close modal

Figures 5(a)5(d) show the mean heat flux per unit length of the height of the channel at the output plane ( Q ¯ ) , the mean flow power dissipated ( E ̇ ¯ ) , 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 ¯ / ρ c p U H Δ T w = 0 . 33 compared to 0.32 in the channel without the reed). However, in the finite range of log10M* 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 ¯ / ρ c p U H Δ T w = 0 . 56 occurring at log 10 M * , U * 0 . 3 , 8 .

FIG. 5.

Effect of M* and U* on the thermal performance. (a) The non-dimensional mean heat flux per unit height, (b) mean mechanical energy deficit, (c) thermal enhancement factor. All coefficient are calculated at 10L downstream of the reed leading edge.

FIG. 5.

Effect of M* and U* on the thermal performance. (a) The non-dimensional mean heat flux per unit height, (b) mean mechanical energy deficit, (c) thermal enhancement factor. All coefficient are calculated at 10L downstream of the reed leading edge.

Close modal

It is also observed that in a large area of the M*U* space, especially for 10−0.9M* ≤ 100.1, the system delivers large Q ¯ with the value Q ¯ / ρ c p U H Δ 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 ¯ shows very similar values for M* ≤ 100.1 when the reed is in the oscillatory state, and there is only a sharp change in the value of Q ¯ near the stability boundaries.

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

Although the reed shows high Q ¯ 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 * 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.

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.3L. 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.5L). 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

FIG. 6.

Dependence of (a) peak-to-peak amplitude of the trailing edge and (b) vibrational Strouhal number with U* for the reed with M* = 1 and 0.1.

FIG. 6.

Dependence of (a) peak-to-peak amplitude of the trailing edge and (b) vibrational Strouhal number with U* for the reed with M* = 1 and 0.1.

Close modal

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 ¯ on U* and it is noted that the M* = 1 case has high Q ¯ 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 ¯ 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*.

FIG. 7.

Dependence of (a) the mean heat flux, (b) mechanical energy deficit, (c) TEF, and (d) thermal Strouhal number StQ (solid symbols) along with the vibrational Strouhal number StA (hollow symbols) and on U* for the reed with M* = 1 and 0.1.

FIG. 7.

Dependence of (a) the mean heat flux, (b) mechanical energy deficit, (c) TEF, and (d) thermal Strouhal number StQ (solid symbols) along with the vibrational Strouhal number StA (hollow symbols) and on U* for the reed with M* = 1 and 0.1.

Close modal

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 ¯ 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 = 10L 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.

FIG. 8.

Power spectrum of the tip vertical position as a function of the vibrational frequency fAL/U (red dashed line, PA) and power spectrum of the heat flux as a function of the heat flux frequency fQL/U (blue solid line, PQ) for (a) U* = 17 and M* = 1, and (b) U* = 17 and M* = 0.1 chosen from Fig. 7(d). Phase-plane plots for (c) U* = 17 and M* = 1, and (d) U* = 17 and M* = 0.1.

FIG. 8.

Power spectrum of the tip vertical position as a function of the vibrational frequency fAL/U (red dashed line, PA) and power spectrum of the heat flux as a function of the heat flux frequency fQL/U (blue solid line, PQ) for (a) U* = 17 and M* = 1, and (b) U* = 17 and M* = 0.1 chosen from Fig. 7(d). Phase-plane plots for (c) U* = 17 and M* = 1, and (d) U* = 17 and M* = 0.1.

Close modal

In Fig. 9(a), we plot the distribution of the mean Nusselts number N u ¯ 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 ¯ 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 ¯ shows the highest value in the near-wake region of the reed. For U* = 16, the maximum occurs around ∼1.5L whereas for U* = 8 and 12, the maximum of N u ¯ occurs at ∼2.5L. The transition in the position of maximum Nu1 is related to the mode switching of the reed to a higher frequency of oscillation.

FIG. 9.

Histories of (a) mean Nusselt number at the channel walls, and (b) mean flow temperature versus x/L. In these plots, x/L = 0 coincides with the leading edge of the reed. M* = 1. Purple solid line shows the longitudinal placement of the reed.

FIG. 9.

Histories of (a) mean Nusselt number at the channel walls, and (b) mean flow temperature versus x/L. In these plots, x/L = 0 coincides with the leading edge of the reed. M* = 1. Purple solid line shows the longitudinal placement of the reed.

Close modal

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 ¯ 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 ¯ 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 ∼ 10L. 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.

FIG. 10.

Vorticity (line) and temperature contours (colored) of a reed with M* = 1 and (a) U* = 4, (b) U* = 8, (c) U* = 12, and (d) U* = 16. The plot shows 37 contour levels of vorticity ranging from −10 to 10. The corresponding snapshots of the reed deformation are shown on the right of each figure.

FIG. 10.

Vorticity (line) and temperature contours (colored) of a reed with M* = 1 and (a) U* = 4, (b) U* = 8, (c) U* = 12, and (d) U* = 16. The plot shows 37 contour levels of vorticity ranging from −10 to 10. The corresponding snapshots of the reed deformation are shown on the right of each figure.

Close modal

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 (StA = 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.

FIG. 11.

Instantaneous ((a) and (c)) and mean ((b) and (d)) viscous dissipation (τ:∇v) inside the channel with the reed with M* = 1 and ((a) and (b)) U* = 12, ((c) and (d)) U* = 16. Vorticity contours are from −10 to 10 with 12 intervals. In the mean plot, four snapshot of the reed position in time have been shown.

FIG. 11.

Instantaneous ((a) and (c)) and mean ((b) and (d)) viscous dissipation (τ:∇v) inside the channel with the reed with M* = 1 and ((a) and (b)) U* = 12, ((c) and (d)) U* = 16. Vorticity contours are from −10 to 10 with 12 intervals. In the mean plot, four snapshot of the reed position in time have been shown.

Close modal

The above analysis indicates that the more flexible reed creates stronger vortices, but these in turn create significant losses, thereby significantly increasing E ̇ . 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.

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.

FIG. 12.

Dependence of the peak-to-peak amplitude of oscillation, A/L with Re for a reed with M* = 1.

FIG. 12.

Dependence of the peak-to-peak amplitude of oscillation, A/L with Re for a reed with M* = 1.

Close modal
FIG. 13.

Observed mode shapes at different Re for a reed with M* = 1 and U* = 12 (black solid line in Fig. 12).

FIG. 13.

Observed mode shapes at different Re for a reed with M* = 1 and U* = 12 (black solid line in Fig. 12).

Close modal

In Figs.14(a)-14(c), we plot Q ¯ , E ̇ ¯ , and TEF as a function of Re. In addition, for comparison, the variation of Q ¯ and E ̇ ¯ 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 ¯ 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.

FIG. 14.

Dependence of (a) Q ¯ , (b) E ̇ ¯ , and (c) TEF on the value of Re for a reed with M* = 1. The green lines are for a channel without the reed.

FIG. 14.

Dependence of (a) Q ¯ , (b) E ̇ ¯ , and (c) TEF on the value of Re for a reed with M* = 1. The green lines are for a channel without the reed.

Close modal

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.

FIG. 15.

Vorticity (line) and temperature contours (colored) of a reed at different Re for a reed with M* = 1 and U* = 12 (black solid line in Fig. 12). Contour levels are similar to Fig. 10.

FIG. 15.

Vorticity (line) and temperature contours (colored) of a reed at different Re for a reed with M* = 1 and U* = 12 (black solid line in Fig. 12). Contour levels are similar to Fig. 10.

Close modal

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.

FIG. 16.

Dependence of the peak-to-peak amplitude of reed oscillation, A/L on the relative channel height H/L for a reed with M* = 1.

FIG. 16.

Dependence of the peak-to-peak amplitude of reed oscillation, A/L on the relative channel height H/L for a reed with M* = 1.

Close modal

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 ¯ / ρ c p U H Δ T w reduces monotonically with increasing H/L. The different values of U* show very similar trends in Q ¯ / ρ c p U H Δ 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 ¯ / ρ c p U H Δ 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 ¯ / ρ c p U H Δ 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.

FIG. 17.

Dependence of (a) average heat flux per unit height of the channel, Q ¯ ; (b) average energy dissipation per unit height, E ̇ ¯ ; and (c) TEF with normalized channel height H/L for a reed with M* = 1.

FIG. 17.

Dependence of (a) average heat flux per unit height of the channel, Q ¯ ; (b) average energy dissipation per unit height, E ̇ ¯ ; and (c) TEF with normalized channel height H/L for a reed with M* = 1.

Close modal

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.

FIG. 18.

Vorticity (line) and temperature contours (colored) of a reed at (a) H/L = 0.5 and (b) H/L = 1.5 for a reed with M* = 1 and U* = 12 (black solid line in Fig. 16). Contour levels are similar to Fig. 10.

FIG. 18.

Vorticity (line) and temperature contours (colored) of a reed at (a) H/L = 0.5 and (b) H/L = 1.5 for a reed with M* = 1 and U* = 12 (black solid line in Fig. 16). Contour levels are similar to Fig. 10.

Close modal

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.

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 ¯ 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.

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.

Flow pass a stationary cylinder with forced convective heat transfer is simulated. The Reynolds number is R e = ρ D U μ = 20 , Prandtl number is chosen to be P r = c p μ k = 1 , and Grashof number is G r = g β T s T 0 D 3 ρ 3 μ 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 Hron62 using different numerical approaches. The size of computation domain is chosen to be 28D × 20D, 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.0002D/U; the grid size near the cylinder is uniform with Δ x = Δ y = 1 96 D ; and the number of Lagrangian points along the cylinder is Ns = 323. Finally, the penalty parameters are α = − 105 and β = − 102.

FIG. 19.

(a) The computational configuration and boundary conditions for the case of the flow past a fixed heated cylinder. (b) Local Nusselt number on the cylinder boundary measured from the incoming flow direction assuming Ts = 1.

FIG. 19.

(a) The computational configuration and boundary conditions for the case of the flow past a fixed heated cylinder. (b) Local Nusselt number on the cylinder boundary measured from the incoming flow direction assuming Ts = 1.

Close modal

In Fig.19(b), we compare of the local Nusselt number on the surface of the cylinder with the results of Chang and Finlayson,59 Gan et al.,60 Wang et al.,61 and Turek and Hron62 and a very good agreement is observed.

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).

FIG. 20.

(a) Schematic of the problem of flow-induced vibration of a flexible beam attached to a rigid cylinder. (b) Stationary-state time histories of X and Y displacement of the tip of the beam. Red (dashed) lines is form Bhardwaj and Mittal,63 blue (dashed-dotted) line is form Dennis et al.58 and black (solid line) is for the present method. ρs/ρ = 10 and R e = ρ U ¯ D / μ = 100 .

FIG. 20.

(a) Schematic of the problem of flow-induced vibration of a flexible beam attached to a rigid cylinder. (b) Stationary-state time histories of X and Y displacement of the tip of the beam. Red (dashed) lines is form Bhardwaj and Mittal,63 blue (dashed-dotted) line is form Dennis et al.58 and black (solid line) is for the present method. ρs/ρ = 10 and R e = ρ U ¯ D / μ = 100 .

Close modal

The cylinder is placed at 2D from both the inlet and bottom wall in a rectangular domain with longitudinal and transversal sizes of, respectively, Lx = 11D and Ly = 4.1D with D being the diameter of the cylinder. The parabolic velocity profile with mean velocity U ¯ 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 Δxmin = Δymin = 0.018D, 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 Δ t = 0 . 005 D / U ¯ .

The beam is inextensible with bending stiffness k b = 0 . 933 ρ D 3 U ¯ 2 which is equivalent to assuming Young modules E = 1 . 4 × 1 0 3 ρ U ¯ 2 for a beam with thickness h = 0.2D originally used by Dennis et al.58 In the current test, We consider two cases with (1)ρs/ρ = 10 and R e = ρ U ¯ D / μ = 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.

TABLE II.

Amplitude of y motion at the tip of the beam A m , the Strouhal number S t , and the drag coefficient C D for the flow induced vibration of the beam behind a fixed cylinder.

Cases Sources Am/D St CD
ρs/ρ = 10, Re = 100  Present result  0.89  0.19  4.19 
  Dennis et al.58   0.83  0.19  4.13 
  Bhardwaj and Mittal63   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 Mittal63   0.41  0.28  2.20 
  Tian et al.64   0.32  0.29  2.16 
Cases Sources Am/D St CD
ρs/ρ = 10, Re = 100  Present result  0.89  0.19  4.19 
  Dennis et al.58   0.83  0.19  4.13 
  Bhardwaj and Mittal63   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 Mittal63   0.41  0.28  2.20 
  Tian et al.64   0.32  0.29  2.16 
1.
F. P.
Incropera
,
A. S.
Lavine
, and
D. P.
DeWitt
,
Fundamentals of Heat and Mass Transfer
(
John Wiley & Sons
,
2011
).
2.
S.
Tiggelbeck
,
N.
Mitra
, and
M.
Fiebig
, “
Flow structure and heat transfer in a channel with multiple longitudinal vortex generators
,”
Exp. Therm. Fluid Sci.
5
,
425
436
(
1992
).
3.
G.
Biswas
,
K.
Torii
,
D.
Fujii
, and
K.
Nishino
, “
Numerical and experimental determination of flow structure and heat transfer effects of longitudinal vortices in a channel flow
,”
Int. J. Heat Mass Transfer
39
,
3441
3451
(
1996
).
4.
T.
Tanaka
,
M.
Itoh
,
T.
Hatada
, and
H.
Matsushima
, “
Influence of inclination angle, attack angle, and arrangement of rectangular vortex generators on heat transfer performance
,”
Heat Transfer—Asian Res.
32
,
253
267
(
2003
).
5.
R.
Mahalingam
and
A.
Glezer
, “
Design and thermal characteristics of a synthetic jet ejector heat sink
,”
J. Electron. Packag.
127
,
172
177
(
2005
).
6.
D. W.
Gerlach
,
D.
Gerty
,
R.
Mahalingam
,
Y. K.
Joshi
, and
A.
Glezer
, “
A modular stackable concept for heat removal from 3-D stacked chip electronics by interleaved solid spreaders and synthetic jets
,”
IEEE Trans. Adv. Packag.
32
,
431
439
(
2009
).
7.
P.
Hidalgo
,
F.
Herrault
,
A.
Glezer
,
M.
Allen
,
S.
Kaslusky
, and
B. S.
Rock
, “
Heat transfer enhancement in high-power heat sinks using active reed technology
,” in
Thermal Investigations of ICs and Systems
,
Proceedings of the 16th International Workshop On Thermal Investigations of ICs and Systems (THERMINIC)
,
Barcelona, Spain
,
2010
(
IEEE
,
2010
), pp.
1
6
.
8.
F.
Herrault
,
P.
Hidalgo
,
C.
Ji
,
A.
Glezer
, and
M.
Allen
, “
Cooling performance of micromachined self-oscillating reed actuators in heat transfer channels with integrated diagnostics
,” in
Micro Electro Mechanical Systems
,
Proceedings of the IEEE 25th International Conference On Micro Electro Mechanical Systems (MEMS)
,
Paris, France
,
2012
(
IEEE
,
2010
), pp.
1217
1220
.
9.
J.
Backus
, “
Small-vibration theory of the clarinet
,”
J. Acoust. Soc. Am.
35
,
305
313
(
2005
).
10.
S. D.
Sommerfeldt
and
W. J.
Strong
, “
Simulation of a player–clarinet system
,”
J. Acoust. Soc. Am.
83
,
1908
1918
(
1988
).
11.
J.
Dalmont
and
C.
Frappe
, “
Oscillation and extinction thresholds of the clarinet: Comparison of analytical results and experiments
,”
J. Acoust. Soc. Am.
122
,
1173
1179
(
2007
).
12.
J.
Zhang
,
S.
Childress
,
A.
Libchaber
, and
M.
Shelley
, “
Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind
,”
Nature
408
,
835
839
(
2000
).
13.
M.
Shelley
,
N.
Vandenberghe
, and
J.
Zhang
, “
Heavy flags undergo spontaneous oscillations in flowing water
,”
Phys. Rev. Lett.
94
,
094302
(
2005
).
14.
B. S.
Connell
and
D. K.
Yue
, “
Flapping dynamics of a flag in a uniform stream
,”
J. Fluid Mech.
581
,
33
67
(
2007
).
15.
L.
Zhu
and
C. S.
Peskin
, “
Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method
,”
J. Comput. Phys.
179
,
452
468
(
2002
).
16.
L.
Huang
, “
Flutter of cantilevered plates in axial flow
,”
J. Fluids Struct.
9
,
127
147
(
1995
).
17.
T.
Theodorsen
and
W.
Mutchler
, “
General theory of aerodynamic instability and the mechanism of flutter
,”
NACA Report No. 496
,
1935
.
18.
A.
Kornecki
,
E.
Dowell
, and
J.
O’brien
, “
On the aeroelastic instability of two-dimensional panels in uniform incompressible flow
,”
J. Sound Vib.
47
,
163
178
(
1976
).
19.
M.
Argentina
and
L.
Mahadevan
, “
Fluid-flow-induced flutter of a flag
,”
Proc. Natl. Acad. Sci. U. S. A.
102
,
1829
1834
(
2005
).
20.
C.
Eloy
,
C.
Souilliez
, and
L.
Schouveiler
, “
Flutter of a rectangular plate
,”
J. Fluids Struct.
23
,
904
919
(
2007
).
21.
Y.
Watanabe
,
S.
Suzuki
,
M.
Sugihara
, and
Y.
Sueoka
, “
An experimental study of paper flutter
,”
J. Fluids Struct.
16
,
529
542
(
2002
).
22.
T.
Balint
and
A.
Lucey
, “
Instability of a cantilevered flexible plate in viscous channel flow
,”
J Fluids Struct.
20
,
893
912
(
2005
).
23.
G.
Tetlow
and
A.
Lucey
, “
Motions of a cantilevered flexible plate in viscous channel flow driven by a constant pressure drop
,”
Commun. Numer. Methods Eng.
25
,
463
482
(
2009
).
24.
D.
Tang
,
H.
Yamamoto
, and
E. H.
Dowell
, “
Flutter and limit cycle oscillations of two-dimensional panels in three-dimensional axial flow
,”
J. Fluids Struct.
17
,
225
242
(
2003
).
25.
J.
Lighthill
,
Mathematical Biofluiddynamics
(
Society for Industrial & Applied Mathematics
,
USA
,
1975
).
26.
S.
Michelin
and
S. G. L.
Smith
, “
Resonance and propulsion performance of a heaving flexible wing
,”
Phys. Fluids
21
,
071902
(
2009
).
27.
K.
Shoele
and
Q.
Zhu
, “
Leading edge strengthening and the propulsion performance of flexible ray fins
,”
J. Fluid Mech.
693
,
402
432
(
2012
).
28.
K.
Shoele
and
Q.
Zhu
, “
Performance of a wing with nonuniform flexibility in hovering flight
,”
Phys. Fluids
25
,
041901
(
2013
).
29.
L.
Tang
,
M. P.
Païdoussis
, and
J.
Jiang
, “
Cantilevered flexible plates in axial flow: Energy transfer and the concept of flutter-mill
,”
J. Sound Vib.
326
,
263
276
(
2009
).
30.
D. R.
Miller
, “
Critical flow velocities for collapse of reactor parallel-plate fuel assemblies
,”
J. Eng. Gas Turbines Power
82
,
83
91
(
1960
).
31.
G.
Kim
and
D.
Davis
, “
Hydrodynamic instabilities in flat-plate-type fuel assemblies
,”
Nucl. Eng. Des.
158
,
1
17
(
1995
).
32.
C.
Guo
and
M.
Païdoussis
, “
Analysis of hydroelastic instabilities of rectangular parallel-plate assemblies
,”
J. Pressure Vessel Technol.
122
,
502
508
(
2000
).
33.
A.
Tornberg
and
M. J.
Shelley
, “
Simulating the dynamics and interactions of flexible fibers in stokes flows
,”
J. Comput. Phys.
196
,
8
40
(
2004
).
34.
W.
Huang
,
S. J.
Shin
, and
H. J.
Sung
, “
Simulation of flexible filaments in a uniform flow by the immersed boundary method
,”
J. Comput. Phys.
226
,
2206
2228
(
2007
).
35.
B.
Audoly
and
Y.
Pomeau
,
Elasticity and Geometry
(
Oxford University Press
,
2010
).
36.
C.
Guo
and
M.
Paidoussis
, “
Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow
,”
J. Appl. Mech.
67
,
171
176
(
2000
).
37.
C.
Eloy
,
R.
Lagrange
,
C.
Souilliez
, and
L.
Schouveiler
, “
Aeroelastic instability of cantilevered flexible plates in uniform flow
,”
J. Fluid Mech.
611
,
97
106
(
2008
).
38.
S.
Michelin
,
S. G. L.
Smith
, and
B. J.
Glover
, “
Vortex shedding model of a flapping flag
,”
J. Fluid Mech.
617
,
1
10
(
2008
).
39.
D.
Goldstein
,
R.
Handler
, and
L.
Sirovich
, “
Modeling a no-slip flow boundary with an external force field
,”
J. Comput. Phys.
105
,
354
366
(
1993
).
40.
R.
Glowinski
,
T.
Pan
,
T. I.
Hesla
, and
D. D.
Joseph
, “
A distributed lagrange multiplier/fictitious domain method for particulate flows
,”
Int. J. Multiphase Flow
25
,
755
794
(
1999
).
41.
T.
Pan
,
D.
Joseph
,
R.
Bai
,
R.
Glowinski
, and
V.
Sarin
, “
Fluidization of 1204 spheres: Simulation and experiment
,”
J. Fluid Mech.
451
,
169
191
(
2002
).
42.
M.
Uhlmann
, “
An immersed boundary method with direct forcing for the simulation of particulate flows
,”
J. Comput. Phys.
209
,
448
476
(
2005
).
43.
Z.
Wang
,
J.
Fan
, and
K.
Luo
, “
Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles
,”
Int. J. Multiphase Flow
34
,
283
302
(
2008
).
44.
E.
Uddin
,
W.
Huang
, and
H. J.
Sung
, “
Interaction modes of multiple flexible flags in a uniform flow
,”
J. Fluid Mech.
729
,
563
583
(
2013
).
45.
R.
Maniyeri
and
S.
Kang
, “
Numerical study on bacterial flagellar bundling and tumbling in a viscous fluid using an immersed boundary method
,”
Appl. Math. Modell.
38
,
3567
3590
(
2014
).
46.
Y.
Zang
,
R. L.
Street
, and
J. R.
Koseff
, “
A non-staggered grid, fractional step method for time-dependent incompressible Navier-Stokes equations in curvilinear coordinates
,”
J. Comput. Phys.
114
,
18
33
(
1994
).
47.
R.
Mittal
,
H.
Dong
,
M.
Bozkurttas
,
F.
Najjar
,
A.
Vargas
, and
A.
Von Loebbecke
, “
A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries
,”
J. Comput. Phys.
227
,
4825
4852
(
2008
).
48.
C.
Perng
and
R.
Street
, “
A coupled multigrid-domain-splitting technique for simulating incompressible flows in geometrically complex domains
,”
Int. J. Numer. Methods Fluids
13
,
269
286
(
1991
).
49.
B.
Leonard
, “
Simple high-accuracy resolution program for convective modelling of discontinuities
,”
Int. J. Numer. Methods Fluids
8
,
1291
1318
(
1988
).
50.
H.
Choi
and
P.
Moin
, “
Effects of the computational time step on numerical solutions of turbulent flow
,”
J. Comput. Phys.
113
,
1
4
(
1994
).
51.
C. S.
Peskin
, “
The immersed boundary method
,”
Acta Numer.
11
,
479
(
2002
).
52.
D. L.
Gee
and
R.
Webb
, “
Forced convection heat transfer in helically rib-roughened tubes
,”
Int. J. Heat Mass Transfer
23
,
1127
1136
(
1980
).
53.
J.
Tsia
and
J.
Hwang
, “
Measurements of heat transfer and fluid flow in a rectangular duct with alternate attached–detached rib-arrays
,”
Int. J. Heat Mass Transfer
42
,
2071
2083
(
1998
).
54.
V.
Ozceyhan
,
S.
Gunes
,
O.
Buyukalaca
, and
N.
Altuntop
, “
Heat transfer enhancement in a tube using circular cross sectional rings separated from wall
,”
Appl. Energy
85
,
988
1001
(
2008
).
55.
P.
Promvonge
,
T.
Chompookham
,
S.
Kwankaomeng
, and
C.
Thianpong
, “
Enhanced heat transfer in a triangular ribbed channel with longitudinal vortex generators
,”
Energy Convers. Manage.
51
,
1242
1249
(
2010
).
56.
C.
Guo
and
M.
Paidoussis
, “
Stability of rectangular plates with free side-edges in two-dimensional inviscid channel flow
,”
J. Appl. Mech.
67
,
171
176
(
2000
).
57.
S.
Alben
, “
The flapping-flag instability as a nonlinear eigenvalue problem
,”
Phys. Fluids
20
,
104106
(
2008
).
58.
S. C. R.
Dennis
,
J.
Hudson
, and
N.
Smith
, “
Steady laminar forced convection from a circular cylinder at low Reynolds numbers
,”
Phys. Fluids
11
,
933
940
(
1968
).
59.
M. W.
Chang
and
B. A.
Finlayson
, “
Heat transfer in flow past cylinders at Re < 150—Part I. Calculations for constant fluid properties
,”
Numer. Heat Transfer, Part A
12
,
179
195
(
1987
).
60.
H.
Gan
,
J.
Chang
,
J. J.
Feng
, and
H. H.
Hu
, “
Direct numerical simulation of the sedimentation of solid particles with thermal convection
,”
J. Fluid Mech.
481
,
385
411
(
2003
).
61.
Z.
Wang
,
J.
Fan
,
K.
Luo
, and
K.
Cen
, “
Immersed boundary method for the simulation flows with heat transfer
,”
Int. J. Heat Mass Transfer
52
,
4510
4518
(
2009
).
62.
S.
Turek
and
J.
Hron
,
Proposal for Numerical Benchmarking of Fluid-Structure Interaction Between an Elastic Object and Laminar Incompressible Flow
(
Springer
,
2006
).
63.
R.
Bhardwaj
and
R.
Mittal
, “
Benchmarking a coupled immersed-boundary-finite-element solver for large-scale flow-induced deformation
,”
AIAA J.
50
,
1638
1642
(
2012
).
64.
F.
Tian
,
H.
Dai
,
H.
Luo
,
J. F.
Doyle
, and
B.
Rousseau
, “
Fluid–structure interaction involving large deformations: 3D simulations and applications to biological systems
,”
J. Comput. Phys.
258
,
451
469
(
2014
).