Heat transfer enhancement by asymmetrically clamped flexible flags in a channel flow

International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Heat transfer enhancement by asymmetrically clamped flexible flags in a channel flow Jae Bok Lee, Sung Goon Park, Hyung Jin Sung ⇑ Department of Mechanical Engineering, KAIST, 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea

a r t i c l e

i n f o

Article history: Received 16 May 2017 Received in revised form 5 September 2017 Accepted 24 September 2017 Available online 28 September 2017 Keywords: Fluid–structure–thermal interaction Penalty immersed boundary method Heat transfer enhancement Flexible flags Thermal mixing

a b s t r a c t Two flexible flags clamped in a heated channel were numerically modeled to investigate the dynamics of the flexible flags and their effects on heat transfer enhancement. The penalty immersed boundary method was adopted to analyze the fluid–structure–thermal interaction between the surrounding fluid and the flexible flags. A system comprising the thermally conductive flags in an asymmetric configuration (FAC) with respect to the channel centerline is described for the first time in the present study. The effect of the resulting vortices on heat transfer enhancement was investigated. The FAC generated a reverse Kármán vortex street that encouraged a greater degree of thermal mixing in the wake compared to the vortical structures generated by the flags in a symmetric configuration (FSC). The ratio of FAC occupying a cross-section to the channel height decreased, resulting in a decrease in the pressure drop compared to FSC. The FAC significantly improved the thermal efficiency compared to the FSC. The effects of the gap distance between FAC (G/L) and the ratio of the channel height to the flag length (H/L) on the thermal enhancement were characterized to identify the parameters that optimized the thermal efficiency. The relationship between the flapping dynamics and the heat transfer properties was examined in detail. The presence of the FAC with the optimal parameters increased convective heat transfer by 207% and the thermal efficiency factor by 135% compared to the baseline (open channel) flow. The thermal efficiency factor obtained in the present study was compared with that obtained in the previous studies. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Efficient heat transfer is an essential function of electrical devices, heat sink systems, heat exchangers, and a variety of industrial applications. Researchers have enhanced heat transfer in pipe flow systems by introducing surface roughness and various types of vortex generators. In previous studies [1–3], rib-roughened surfaces were found to enhance convective heat transfer in pipe flows at a given friction compared to sand-grain roughened surfaces. An alternative method of enhancing heat transfer involves using vortex generators. Zhu et al. [4] investigated the effects of wing-type and winglet-type vortex generators on heat transfer enhancement in turbulent channel flows. They reported that longitudinal vortex generators, combined with roughness elements, could increase heat transfer by 450%. Biswas et al. [5] conducted experimental and numerical investigations to determine the flow structure generated by a winglet-type vortex generator. The generated vortices disrupted the growth of the thermal boundary layer and enhanced convective heat transfer at the channel walls. The transverse ⇑ Corresponding author. E-mail address: [email protected] (H.J. Sung). https://doi.org/10.1016/j.ijheatmasstransfer.2017.09.094 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

vortex generators induced swirling flows that significantly improved convective heat transfer in the internal flows. Rigid vortex-generating tabs, such as delta tabs, rectangular tabs, and trapezoidal tabs, were mounted on the inner walls of a pipe or a duct to enhance fluid mixing and enhance heat transfer [7–11]. The rigid vortex generators, however, blocked the streamwise flows, resulting in the significant pressure drop penalty, high manufacturing cost, and more complex thermal systems [4,6]. The augmentation of heat transfer and the pressure drop must be considered simultaneously in an assessment of the thermal efficiency. Much attention has focused on the use of flexible structures in thermal systems to minimize the pressure drop and simplify the system. Fernandez and Poulter [12] enhanced convective heat transfer by inserting a rectangular sheet metal flag that was allowed to flap in a tubular turbulent flow. The flag significantly increased the heat transfer, and the pressure drop was not significant enough to detract from the heat transfer enhancement. An actuated flexible reed was installed into an air-cooled heat sink channel, leading to a remarkable improvement in the thermal performance [13]. A self-oscillating reed was used to increase the thermal performance in an air-cooled heat sink channel system [14,15]. The fluttering motion of the self-oscillating reed disrupted

1004

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Nomenclature A c1, c2 cp Eloss loss E

flapping amplitude constants in the feedback law heat capacity net energy loss time-averaged value of the net energy loss F Lagrangian momentum force Fh hydrodynamic force Fr restoring force f Eulerian momentum force fr friction factor flapping frequency ff fV vortex shedding frequency H channel height k thermal conductivity coefficient ha overall convective heat transfer coefficient L flag length Nu = haL/k Nusselt number p pressure Pr = cp l/k Prandtl number Qnet net heat flux qi heat sources from immersed bodies Re = q0UL/l Reynolds number si curvilinear coordinates T fluid temperature T0 fluid temperature at inlet flag temperature Ts t time

the momentum and thermal boundary layers at the channel walls and enhanced fluid mixing between the thermal boundary layers and the channel core flow. Numerical simulations have been used to investigate the kinematics of flexible structures as they affect fluid mixing and heat transfer, while optimizing the parameters for a given flow condition [16–24]. Arbitrary Lagrangian Eulerian (ALE) and immersed boundary method (IBM) techniques have been widely used to simulate fluid–flexible structure–thermal interactions. Ali et al. [18,19] adopted an ALE method that required remeshing procedures at each time step to simulate the kinematics of flexible flaps clamped at the channel walls and their effects on fluid mixing and heat transfer. The flexible flaps increased both the mixing between scalars and the thermal efficiency compared to the baseline flow. IBM has been preferentially used to simulate fluid–flexible structure–thermal interaction problems because it mitigates the difficulties associated with re-meshing at deforming or moving bodies immersed in a fluid domain [25–29]. Previous IBM studies [20–23] revealed that a self-oscillating reed clamped longitudinally at the channel centerline enhanced convective heat transfer; however, clamping the edge of the reed at the channel centerline required considerable attention and a need for additional devices that fixed the flexible reed [13,14]. To overcome these shortcomings, Lee et al. [24] proposed using a pair of flexible flags clamped transversally at the channel walls. They clamped the thermally conductive flags in a symmetric configuration (FSC) onto the channel walls to explore the effects of the clamped flags on heat transfer enhancement at various bending rigidities, channel heights, and Reynolds numbers. The presence of the clamped flags increased the net heat flux by 185% and the thermal efficiency by 106% compared to the baseline flow. The 6% increase in the thermal efficiency was insufficient, however, and the proposed thermal system will require additional improvements to the thermal efficiency. The FSC acted as an obstacle to the streamwise flow and generated significant mechanical energy loss by introducing a pressure drop that reduced the thermal efficiency.

U u Ui,ib Xd Xi Xi,ib ytip yw

bulk mean velocity at the inlet fluid velocity velocity at the immersed boundary length of the computational domain flag position immersed boundary position y-position of the flag tip from the center of channel y-position of the flag tip from the wall

Greek symbols d delta function c bending rigidity g thermal efficiency factor l dynamic fluid viscosity r tension force q density ratio q0 fluid density q1 flag density temperature on the massive boundary Ui Ui, ib temperature on the massless boundary Subscripts in out w 0

inlet output plane wall channel without vortex generator/reference

In the present study, the system comprising the thermally conductive flags in an asymmetric configuration (FAC) with respect to the channel centerline was newly proposed in an effort to enhance the mean heat flux and alleviate the mechanical energy loss in the heat sink channel system. We adopted the penalty IBM and simulated the kinematics of flexible flags clamped asymmetrically with respect to the channel centerline. The effects of the vortical structures generated by the FSC and FAC on the heat transfer were analyzed. A parametric study was performed to obtain an optimal parameter set that maximized the thermal efficiency as a function of the channel height (H/L) and the gap distance between flexible flags (G/L). Based on the optimal parameter set proposed in the previous study [24], the bending rigidity (c) and the Reynolds number (Re) were set to 0.04 and 600, respectively. The heat transfer as a function of the flapping dynamics was examined for various G/L and H/L values. The FAC defined over the optimal parameter set provided a higher thermal efficiency than the FSC and provided a superior thermal efficiency compared to other heat transfer enhancement techniques proposed in previous studies.

2. Problem formulation 2.1. Problem description Fig. 1 presents a schematic diagram of the computational domain. The length of the flexible flags is L. The leading edges of the upper and lower flags were clamped vertically at the top and bottom walls, respectively. The gap distance between the leading edges is G. The trailing edges of the flexible flags were free. The clamped position of the upper flag was fixed at 6L from the inlet (x = 0). The clamped position of the lower flag was G/L upstream from that of the upper flag. The initial positions of the clamped leading edge were ð0; H=2LÞ and ðG=L; H=2LÞ for the upper and lower flexible flags, respectively. Lout =L indicates the distance

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

1005

Fig. 1. Schematic diagram of the computational domain.

between the leading edge of the upper flag and the output plane. A velocity profile based on the Poiseuille channel flow with a mean velocity (U) of 1 and a constant temperature of T ¼ T 0 ¼ 0 was applied at the inlet boundary ðx ¼ 6LÞ. No-slip conditions ðu ¼ v ¼ 0Þ and a constant wall temperature ðT ¼ T w ¼ 1Þ were imposed at the channel walls. The velocity and temperature gradients at the outlet boundary (x = 26L) were set to zero with respect to the x direction. The motion of the fluid flow was defined on the Eulerian coordinate. The reference points in the x and y directions were defined as the clamped position of the upper flag and the channel centerline, respectively. The motions of the flexible flags were described on the curvilinear coordinate, s. The incompressible fluid flow was governed by the continuity equation and the Navier–Stokes equations,

r  u ¼ 0;

ð1Þ

@u 1 þ u  ru ¼ rp þ r2 u þ f ; @t Re

ð2Þ

where u is the velocity vector, p is the pressure, and f is the momentum forcing term that enforced the no-slip conditions on the flag. The Reynolds number was defined by Re ¼ q0 UL=l; where l is the dynamic viscosity. The fluid density q0 was used to define the characteristic density, and the flag length L defined the characteristic length. The bulk mean fluid velocity at the inlet U defined the characteristic velocity. The governing equations were nondimensionalized using the following non-dimensional scales: L=U 2

@X @X  ¼ 1: @s @s

@T 1 þ u  rT ¼ r2 T þ q; @t Re  Pr

ð3Þ

where Pr is the Prandtl number Pr ¼ lcp =q0 , cp is the heat capacity, k is the thermal conductivity, q is the heat source of the immersed body. The dimensionless temperature was defined as T ¼ ðT  T 0 Þ=ðT w  T 0 Þ, where T w and T 0 indicate the wall temperature and the fluid temperature at the inlet, respectively. The temperatures of all flags were set equal to the constant wall temperature so that the thermally conductive flags were analogous to ideal fins. The motions of the flexible flags were governed by

!  @X i @2 @2Xi  2 c 2  F i; ri @si @si @si




ð4Þ

where the subscript i denotes the lower (i = 1) and upper flags (i = 2), s is the curvilinear coordinate on the flexible flag, X i ¼ ðX i ðs; tÞ; Y i ðs; tÞÞ is the position of the flag, ri is the tension force along the curvilinear coordinate, c is the bending rigidity, and F i is the Lagrangian momentum force exerted on the flags by the surrounding fluid. The governing equation of the flexible flags was non-dimensionalized according to the following characteristic scales: q1 U 2 L2 for the bending rigidity, q1 U 2 for the tension force,

ð5Þ

The fluid–structure–thermal interaction between a pair of flags and the surrounding fluid was resolved by adopting the penalty IBM. The interaction between the flexible flags and the fluid flow was derived using a feedback forcing law [27,30],

Z

t

F i ¼ c1

0

ð6Þ

0

ð7Þ

ðU i;ib  U i Þdt þ c2 ðUi;ib  U i Þ; 0

Z Q i ¼ c1

t

0

ðCi;ib  Ci Þdt þ c2 ðCi;ib  Ci Þ;

where c1 and c2 are large negative free constants (a = 1.74  107, b = 2.60  102) [23,31], Uib is the fluid velocity on the immersed boundary obtained by interpolation, U is the velocity of the inverted flags expressed by U = dX/dt. The massive boundary of the flexible flags was linked strongly to a massless boundary through a virtual stiff spring with damping. In Eqs. (6) and (7), U i and Ci are the velocity and temperature of the massive boundary, and U i;ib and Ci;ib are the velocity and temperature of the massless boundary,

Z

U i;ib ðs; tÞ ¼

XF

2

for time, q0 U for pressure, and q0 U =L for the Eulerian momentum forcing f . The fluid temperature described on the Eulerian coordinate was governed by the energy equation,

@2 Xi @ ¼ 2 @s @t i

and q1 U 2 =L2 for the Lagrangian momentum force. The inextensibility condition was imposed on the flexible flags,

uðx; tÞdðX i ðs; tÞ  xÞdx;

ð8Þ

Tðx; tÞdðX i ðs; tÞ  xÞdx;

ð9Þ

Z

Ci;ib ðs; tÞ ¼

XF

where XF is the fluid domain. Similarly, the Lagrangian momentum forces and heat sources were calculated according to

Z

f i ðx; tÞ ¼ q

XS

F i ðx; tÞdðx  X i ðs; tÞÞds;

ð10Þ

Z qi ðx; tÞ ¼

XS

Q i ðx; tÞdðx  X i ðs; tÞÞds;

ð11Þ

where XF is the domain of the immersed boundary, q is the density ratio q1 =ðq0 LÞ. The Eulerian momentum force f and the heat source of the immersed body q were determined by summing the values obtained from each flag (f 1 ; f 2 and q1 ; q2 ). The procedure used to solve the fluid-flexible structure-thermal interaction is described elsewhere [23,24]. The interaction Eqs. (8)–(11) were obtained by adopting a four-point smoothed delta function [25]. In the present study, the fluid governing equation was solved using the fractional step method [32,33]. 2.2. Output variables The net heat flux (Q net ), mechanical energy loss (Eloss ), and thermal efficiency factor (g) were defined to examine the heat transfer enhancement. The net heat flux in the thermal energy was expressed as

1006

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

"Z Q net ¼ qcp

Z

H=2

H=2

uTdyjout 

#

H=2

H=2

uTdyjin :

ð12Þ

The mechanical energy loss due to the pressure drop between the inlet and output planes was calculated according to

Z Eloss ¼

Z

H=2

H=2

updyjin 

H=2

H=2

updyjout :

ð13Þ

The thermal efficiency factor was defined as the ratio of the mean heat flux of the system comprising the flexible flag to that of the system without a secondary apparatus at a constant pumping power,

g¼

Q net Q 0;net



E0;loss Eloss

!1=3 ;

ð14Þ

where the subscript 0 indicates the value measured in the absence of secondary devices. The bar denotes a time-averaged value in Eq. (14). A thermal efficiency factor exceeding 1 was expected to increase the heat transfer for the same pumping power. The mean heat flux, mechanical energy loss, and thermal efficiency were calculated over a cross-sectional profile of the output plane. For the purpose of determining an appropriate value of Lout =L, we studied the dependence of the thermal enhancement on x/L by considering the flexible flags in a symmetric configuration (FSC). The mean heat flux as a function of Lout =L is presented in Fig. 2. The output plane was fixed at x = 25L, and the mean heat flux was examined by varying Lout =L. As Lout =L increased, the mean heat flux converged to the optimal value of 0.673 at Lout ¼ 25L  26L and then decreased from 27L. In this study, the minimum length required to enhance heat transfer, Lout , was determined to be 25L, and Lout =L ¼ 25 remained constant. 3. Results and discussion 3.1. Heat transfer enhancement by varying the flag configuration In a previous study [24], the heat transfer enhancement by the presence of flexible flags in a symmetric configuration (FSC) with respect to the channel centerline was explored as a function of the bending rigidity, channel height, and Reynolds number. The presence of the flexible flags, under the optimal parameter set (c = 0.04, H/L = 2.5, Re = 600), increased the net heat flux by 185% and the thermal efficiency factor by 106% compared to the baseline flow. In the present study, we proposed using flexible flags in an asymmetric configuration (FAC) to further enhance the thermal efficiency compared to the FSC system reported previously. The bending rigidity (c = 0.04) and Reynolds number (Re = 600) remained fixed unless otherwise stated. Fig. 3 plots the vorticity and temperature fields obtained in the FSC and FAC systems. The

gap distances between the leading edges of the upper and lower flags (G/L) were 0 and 1.0 for FSC and FAC, respectively. The time-averaged values of the net heat flux, mean drag coefficient, mechanical energy loss via pressure drop, and thermal efficiency are listed in Table 1. The vortices shed from the flapping flags were symmetric with respect to the channel centerline, as shown in Fig. 3(a). Two staggered rows of vortices comprising positive circulations near the top wall and negative circulations near the bottom wall passed through the wake, and the structure of the reverse Kármán vortex street was maintained in the wake, as shown in Fig. 3(b). The angular momentum of the vortical structures refreshed the thermal boundary layer and induced the heated fluid near channel walls into the channel centerline. Fig. 4(a) plots the absolute value of  j(x, y, t), the time-averaged local vorticity jx

jxjx ¼

L Hðtn  t0 Þ

Z

tn

t0

Z

H=2L

H=2L

jxjðx; y; tÞdydt;

ð15Þ

which indicated that the vortex strength was correlated with the degree of heat transfer from the wall to the center of the channel, at each streamwise position. The vortex strength near the trailing edges of the flapping flags increased and then decreased as the generated vortices moved downstream and dissipated into the surrounding fluid, as shown in Fig. 4(a). The mean absolute value of   R 26L 1 the vorticity jxj ¼ 32L jxjx dx for the FAC exceeded that for 6L the FSC, as shown in Table 1, indicating that heat transfer in the FAC system was better than that in the FSC system. Switching the flag configuration to the asymmetric configuration increased the mean heat flux due to the global augmentation of the vortex strength. Flexible flags in the FSC were clamped at a given cross-section. Flexible flags in the FAC were clamped at different cross-sections, respectively. The flexible flags in the FSC provided a narrower path of fluid flow between the flags than those in the FAC. In the FSC, flapping flags clamped at a given cross-sectional position acted as obstacles within the streamwise flow and raised the form drag further than those in the FAC. The drag coefficient in the FSC system was larger than that in the FAC system, as shown in Fig. 4 (b) and (c) due to the increased form drag. The mechanical energy loss in the FAC system was smaller than that in the FSC system, as shown in Table 1. The asymmetric arrangement of flapping flags produced a high vorticity, improved heat transfer, and reduced the mechanical energy loss. The asymmetric arrangement of the thermally conductive flags is expected to play an important role in enhancing the thermal efficiency. A parametric study across flow conditions is required to optimize the thermal efficiency of this system. Fig. 5 plots the mean heat flux, mechanical energy loss, and thermal efficiency factor for various G/L and H/L values. As the channel height decreased, the thermal mixing region occupying the wake region expanded, and the blockage ratio increased, which increased both the mean heat flux and the mechanical energy loss, as shown in Fig. 5(a) and (b). The thermal efficiency trend was similar to the mean heat flux trend shown in Fig. 5(c). For comparison purpose, we considered an alternative thermal enhancement measure based on the pressure drop penalty and the convective heat transfer enhancement,

Nu ¼

ha L ; k

ð16Þ

where Nu is the Nusselt number, ha is the overall convective heat transfer coefficient, and k is the thermal conductivity. The friction factor fr was computed using Fig. 2. Mean heat flux as a function of Lout =L for the symmetric configuration.

fr ¼

2Dh Dp ; 4X d qf U 2

ð17Þ

1007

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Fig. 3. Vorticity and temperature fields for the (a) symmetric and (b) asymmetric configurations (G/L = 1.0).

Table 1 Quantitative comparison of the mean heat flux, mean drag coefficient, mechanical energy loss, and thermal efficiency factor for the (a) symmetric and (b) asymmetric configurations.

Symmetric Asymmetric

j jx

 net Q

 C d

 E loss

g

1.529 1.581

0.670 0.688

2.037 1.917

1.300 1.221

1.06 1.11

Fig. 4. (a) Streamwise variations of the time-averaged local vorticities for the symmetric (blue line) and asymmetric (red line) configurations (G/L = 1.0). Time history of the drag coefficients for the (b) symmetric and (c) asymmetric configurations. Dashed lines denote the mean values of the drag coefficients. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

where Dh ð¼ 2HÞ is the hydraulic diameter of the channel, Xd is the length of the computational domain, and Dp is the pressure drop between the inlet and outlet. The Fanning friction factor was applied to a laminar channel flow without secondary devices,

f r;0 ¼ 24=ReDh :

ð18Þ

The alternative thermal efficiency was calculated by defining the thermal enhancement factor ðg0 Þ,


0

g ¼

Nu Nu0



fr f r;0

!1=3 :

ð19Þ

FAC displayed a thermal enhancement factor ðgmax Þ increase of 129% and an alternative thermal enhancement factor ðg0max Þ increase of 135%. The heat sink channel system devised in the present study provided a higher thermal efficiency factor than was obtained in previous studies, as shown in Table 2.

1008

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Fig. 5. (a) Mean heat flux, (b) mechanical energy loss, and (c) thermal efficiency factor for Re = 600, c ¼ 0:04.

Table 2 Maximum values of the thermal efficiency factor.

Lee et al. [24] Park et al. [22] Shoele and Mittal [20] Ali et al. [19] Present

gmax

g0max

1.06 1.22 1.25 – 1.29

– – – 1.34 1.35

3.2. Effect of the gap distance between flexible flags We explored the kinematics of flexible flags and their effects on the heat transfer enhancement as a function of G/L. The flag movement at any given time depended on the relationship between the hydrodynamic force (Fh) and the restoring force (Fr). Fig. 6 shows a schematic diagram of the flapping flag clamped at the bottom wall. The trailing edges farthest and closest from the bottom wall were denoted by (1) and (2) in Fig. 6, respectively. The hydrodynamic force exerted on the flexible flag acted as a driving force for flag deflection, as shown in Fig. 6(a). The elastic energy stored in the flexible flag during deflection provided a force that restored the flag to its initial position, as shown in Fig. 6(b). As the hydrodynamic force and the restoring force alternately predominated over the counterpart, the trailing edges of the flags periodically moved up and down. The flapping behavior displayed a quasi-periodic oscillation and a continuous fluttering motion. Fig. 7 plots the time-averaged position of the trailing edges and the flapping amplitude as a function of G/L. In Fig. 7(a), the values of yw for the upper and lower flags indicated the time-averaged values of the gap distances with respect to the top and bottom walls, respectively. Fig. 7(a) shows that the value of yw for the lower flag asymptotically converged to that for the upper flag. Because the flapping behavior was quasi-periodic, the flapping amplitude was defined by

A¼

1 Npeak

X

N peak

jyjpeak  ytip j;

ð20Þ

j¼1

where Npeak is the number of peaks in the quasi-steady state, yjpeak are the values of the local peaks, and ytip is the time-averaged posi-

tion of the trailing edge with respect to the channel centerline. In Fig. 7(b), the flow field between the upper and lower flags became unsteady due to the swirling flow induced by the lower flag. Fig. 7 (b) reveals that the difference between the flapping amplitudes of the upper and lower flags was irregular due to the interactions between the flexible flags. This result indicated that the irregular flapping amplitude trend resulted from the irregularly generated flow field between the flexible flags at G/L < 1.6. As G/L increased beyond 1.6, the difference between the values of A/L for the upper and lower flags decreased over the range 1.6 < G/L < 2.6 and then became constant at G/L > 2.6. The interactions between the flexible flags for G/L > 2.6 were not as significant as for G/L < 2.6. Fig. 8 plots the frequency spectra of the net heat flux and trailing edges for Re = 600, c = 0.04, and H/L = 1.75. The insets of Fig. 8 (a) and (b) plot the time histories of the trailing edges of the upper flags in a quasi-steady state as displaying a pattern of recurrence with a quasi-periodic time interval. For G/L = 0.8 and 1.4, the dominant frequencies in the net heat flux ðf Q G=L¼0:8 ¼ 0:032, f Q G=L¼1:4 ¼ 0:023Þ corresponded to the inverses of the quasiperiodic time interval. f Q G=L¼0:8 and f Q G=L¼1:4 corresponded to the minor frequency in the trailing edge of the upper flag ðf f G=L¼0:8 ¼ 0:032, f f G=L¼1:4 ¼ 0:023Þ in Fig. 8(a) and (b), indicating that the heat transfer enhancement trend was correlated with the quasi-periodic flapping motion at G/L  1.4. Fig. 9 plots the time history of the difference between the trailing edges ðY  ¼ jytip;u j  jytip;l jÞ. The dashed lines denote the upper and lower envelopes of the time history of Y  . The difference between the upper and lower envelopes was maximized when the flags oscillated out of phase with the phase angle of p. The difference between the envelopes was zero when the flag oscillations were in phase. Fig. 10 plots the vorticity and temperature fields corresponding to the instants a and b marked in Fig. 9. For G/L = 0.8 and 1.4, the position at which the vortices merged, point b, shifted downstream relative to a, as shown in Fig. 10(a) and (b). At a in Fig. 10 (a) and (b), the upper and lower flags reduced the convection velocity in the upper and lower half of the channel, respectively, by blocking the main flow whereas the convection velocity near the channel centerline increased at x = 0 based on the mass convergence. The counter-clockwise (CCW) vortices were close enough to

Fig. 6. Schematic diagram of the flapping flag clamped at the bottom wall. (a) The hydrodynamic force predominated over the restoring force. (b) The restoring force predominated over the restoring force.

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

1009

Fig. 7. (a) Time-averaged gap distance between the trailing edges and the channel walls, and (b) flapping amplitude at Re = 600, c ¼ 0:04, H/L = 1.75. Red and blue lines indicate the upper and lower flags, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Frequency spectrum of the net heat flux and the trailing edges of the flexible flags at (a) G/L = 0.8, (b) G/L = 1.4, (c) G/L = 1.6, and (d) G/L = 2.6. The insets indicate the time history of the trailing edge of the upper flag.

Fig. 9. Time history of the difference between the trailing edges of the upper and lower flags for (a) G/L = 0.8, (b) G/L = 1.4, (c) G/L = 1.6. Dashed lines in (a) and (b) indicate the envelopes of the time history.

interact with the same-signed vortices and merge due to a decrease in the convection velocity in the upper half of the channel, as shown at point a of Fig. 10(a) and (b). The angular momentum of the merged vortex induced a clockwise (CW) vortex in the lower half of the channel to shift toward the channel centerline, forging a reverse Kármán vortex street in the wake region. The reverse Kármán vortex street disrupted the growth of the thermal boundary layer and enhanced the thermal mixing between the wall region and the channel centerline. At b in Fig. 10(a) and (b), the

flapping flags were in-phase, leading to the increased velocity near the channel centerline based on the mass convergence. The trailing edges of the flexible flags oscillating in phase moved toward the channel walls after leaving the vortices generated near the channel centerline at b in Fig. 10(a) and (b). The augmented convection velocity at the channel centerline shifted the vortices downstream at a faster velocity than the velocity at which vortices were generated from the out-of-phase flapping flags. The position at which the in-phase vortices merged was downstream of the corresponding position for the out-of-phase vortices. In-phase flapping forged the reverse Kármán vortex street further downstream compared to out-of-phase flapping. A downstream shift in the reverse Kármán vortex street generation decreased the thermal mixing region in the wake. Fig. 10(a) and (b) show that the wake region for G/ L = 0.8 included a thermal mixing region that was larger than the region obtained for G/L = 1.4. Fig. 11 shows the mean heat flux, mechanical energy loss via the pressure drop, and the thermal efficiency factor for various G/L values. Fig. 11(a) shows that the mean heat flux decreased as G/L increased from G/L = 0.8 to G/L = 1.4. The ratio of the cross-sectional area occupied by the flags to the channel height decreased with increasing G/L, which decreased the drag. The mechanical energy loss decreased as G/L increased from G/L = 0.8 to G/L = 1.4, as shown in Fig. 11(b). A previous study [24] demonstrated that the vortices generated from the flexible flags merged in the wake. The number of merged vortices corresponded to the ratio between the flapping frequency

1010

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Fig. 10. Vorticity and temperature fields at the instants a and b denoted in Fig. 8 for (a) G/L = 0.8, (b) G/L = 1.4, and (c) G/L = 1.6.

Fig. 11. (a) Mean heat flux, (b) mechanical energy loss, and (c) thermal efficiency factor at Re = 600, c = 0.04, H/L = 1.75.

and the vortex shedding frequency. In this study, the number of merged CCW vortices corresponded to the ratio between the flapping frequency of the upper flag and the vortex shedding frequency. Two CCW vortices merged in the upper half of the channel. As G/L increased from 0.8 to 1.6, the flapping frequency increased, as shown in Fig. 8(a)–(c). The increase in the flapping frequency increased the number of shedding vortices that detached from the trailing edge per unit time. For G/L = 1.6, the vortex shedding frequency was half of the dominant flapping frequency in the trailing edge of the upper flag ðf f G=L¼1:6 ¼ 0:486; f V G=L¼1:6 ¼ 0:243Þ. The vortex shedding frequency corresponded to the dominant frequency in the net heat flux ðf Q G=L¼1:6 ¼ 0:243Þ, as shown in Fig. 8(c), indicating that the heat transfer was correlated with the vortex shedding frequency. Unlike the systems for G/L  1.4, the inverse of the quasi-periodic time interval on the frequency spectrum (Fig. 8(c) negligibly affected the vortex shedding and net heat flux for G/L = 1.6. Vortex merging in the G/L = 1.6 system at b began upstream of a as compared to G/L = 0.8 and 1.4 in Fig. 10(c). As the dominant frequency corresponding to the net heat flux variations had a correlation with the vortex shedding

frequency, in-phase flapping resulted in active thermal mixing near the channel centerline, as observed at a and b in Fig. 10(c). For G/L = 1.6, the reverse Kármán vortex street was generated upstream of the corresponding position for G/L = 1.4, as shown in Fig. 10(b) and (c). The thermal mixing was more active at the channel centerline because thermal mixing began further upstream than in G/L = 1.4. The mean heat flux increased as G/L increased from 1.4 to 1.6 in Fig. 11(a). As G/L increased to 2.6, the fluid momentum in the upper half of the channel between the flexible flags decreased compared to the corresponding value obtained for G/L  1.6. The decreased fluid momentum in the upper half of the channel reduced the hydrodynamic force exerted on the upper flag. The resistance created by the upper flag relative to the hydrodynamic force increased, which increased yw in the upper and lower regions, as shown in Fig. 7(a). The convection velocity of vortex decreased with decreasing fluid momentum in the upper half of the channel. In the region immediately behind the upper flag, the CCW vortices in the upper half of the channel passed through the wake with a slower convection velocity than the CW vortices in the lower half of the channel. The CCW vortex immediately behind the upper flag caught up with

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

the previously generated CCW vortex due to the decrease in the convection velocity. The vortices generated for G/L = 2.6 merged, and the reverse Karman vortex street was formed at a position (x = 2.8) upstream of the corresponding positions in the G/L = 0.8, 1.4, and 1.6 systems due to the reduction in the convection velocity. As the merged vortex was forged upstream and close enough to affect the flapping behavior, a second dominant frequency at the trailing edge of the upper flag was shown for G/L = 2.6. The second dominant frequency ðf f G=L¼2:6 ¼ 0:24Þ that appeared in Fig. 8(d) was half of the first dominant frequency ðf f G=L¼2:6 ¼ 0:48Þ and corresponded to the vortex shedding frequency ðf V G=L¼2:6 ¼ 0:24Þ: The vortex shedding frequency corresponded to the dominant frequency for the net heat flux ðf f G=L¼2:6 ¼ 0:24Þ, indicating that the vortex shedding behavior was closely correlated with the thermal mixing behavior in the wake region. Fig. 12 plots the vorticity and temperature fields for Re = 600, c = 0.04, H/L = 1.75, and G/L = 2.6. 1 (V1) – 5 (V5) denote the vortical structures. As the trailing edge of the upper flag moved toward the top wall, V1 and V2 were detached and then shed into the wake at t = 269.8, as shown in Fig. 12(a). V2 caught up with V1, and V1 and V2 coalesced to form a single vortical structure (V3) at t = 271.6. The CW vorticity (V4) was generated by the angular momentum of V3 and affected the flapping motion and vortex generation near the trailing edge of the upper flag. The angular momentum of V3 entrained the clockwise vortex (V5) near the lower half of the channel into the channel centerline. V5 merged with a same-signed vortex to form a single vortical structure (V6). V3 and V6 formed a single large vortex that was characterized by V30 and V60 during the subsequent vortex merging step. V30 and V60 generated a shear layer that stretched from the walls into the channel centerline. V30 and V60 induced the shear layers from the bottom and top walls, respectively, as observed in Fig. 12(a). The generated shear layer passing through the reverse Kármán vortex street induced the thermal mixing region to diffuse throughout the channel centerline, as shown in Fig. 12(a) and (b). For Re = 600, c = 0.04, H/L = 1.75, and G/L = 2.6, the efficiency exceeded 129% of the efficiency obtained from the channel flow without flags, as shown in Fig. 11(c). The amount of heat transferred

1011

increased to 199% of the heat transferred in the baseline flow, as shown in Fig. 11(a). 3.3. Effect of the channel height The effects of the channel height on the flapping dynamics and heat transfer were investigated. Fig. 13 plots the time-averaged gap distances between the trailing edges and the channel walls and the flapping amplitude as a function of H/L. As H/L increased, the streamwise velocity gradient near the channel walls with respect to the y axis decreased. The fluid momentum exerted on the flexible flags decreased with increasing H/L. The resistance introduced by the flexible flags increased compared to the hydrodynamic force, which increased yw , as shown in Fig. 13(a). The lower hydrodynamic force reduced the driving force for the flapping motion, which decreased the flapping amplitude at higher values of H/L, as shown in Fig. 13(b). Fig. 14 plots the frequency spectrum of the normalized flapping amplitude for various values of H/L. The insets of Fig. 14(b) and (c) indicate the frequency spectra of the net heat flux for H/L = 2.0 and 2.25, respectively. The time history of the trailing edge at H/L = 1.5–2.75 is shown in Fig. 15. The red lines indicate the distance between the upper flag and the top wall, and the blue lines indicate the distance between the lower flag and the bottom wall. The time histories of yw for the lower flag were multiplied by 1 in the plot to facilitate a comparison with the time histories of the upper flag in Fig. 15. As H/L decreased, the hydrodynamic force increased. The augmented hydrodynamic perturbations within certain values of H/L produced flapping instabilities, and the flapping frequency of the upper flag differed from that of the lower flag. For H/L = 1.5–2.0, the flapping behavior was quasi-periodic, and the phase angle for each trailing edge varied in time, as shown in Fig. 15(a)–(c). For H/L = 2.25–2.75, the flapping frequency of the upper flag agreed well with that of the lower flags in Fig. 14(c) and (d), and the phase angle was constant in Fig. 15(d)–(f). Note that the second dominant frequency of the upper flag corresponded to the vortex shedding frequency, as shown in Fig. 8d). In-phase flapping generated vortical structures that promoted active thermal mixing near the channel centerline when the

Fig. 12. Snapshots of (a) the vorticity and (c) temperature fields at Re = 600, c ¼ 0:04, G/L = 2.6.

1012

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Fig. 13. (a) Time-averaged gap distance between the trailing edges and the channel walls, and (b) flapping amplitude at Re = 600, c ¼ 0:04, G/L = 2.6. Red and blue lines indicate the upper and lower flags, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Frequency spectrum of the normalized flapping amplitude for (a) H/L = 1.5, (b) H/L = 2.0, (c) H/L = 2.25, and (d) H/L = 2.5.

Fig. 15. Time history of the trailing edges of the flexible flags at (a) H/L = 1.5, (b) H/L = 1.75, (c) H/L = 2.0, (d) H/L = 2.25, (e) H/L = 2.5, (f) H/L = 2.75.

dominant frequency corresponding to the net heat flux variations was correlated to the vortex shedding frequency, as shown in Figs. 9(c) and 10(c). For H/L = 1.5–2.0, the flapping frequency of the upper flag deviated from that of the lower flag in Fig. 8d) and Fig. 14(a) and (b). For H/L = 2.0, the flapping frequency ðf f H=L¼2:0 ¼ 0:480Þ was twice the vortex frequency ðf V H=L¼2:0 ¼ 0:239Þ, indicating that the merged vortex consisted of two single vortices. The vortex shedding frequency corresponded to the dominant frequency of the net heat fluxðf Q H=L¼2:0 ¼ 0:239Þ, as shown in Fig. 14(b). Because the net heat flux was correlated with the vortex shedding frequency, the flexible flags oscillating in phase increased heat transfer in the wake. The flexible flags oscillated in phase near t = 250–275 in Fig. 15(c). As the in-phase and out-of-phase flapping motions occurred alternately, the time history of the net heat flux displayed a quasi-periodic signal in Fig. 16(a). In Fig. 16(a), the net heat flux tended to increase near t = 275–310, indicating that the vortices generated from the flags oscillating in phase reached the output plane at around 25 s. In

Fig. 16(b), the CCW vortices fully merged near x = 4, and the reverse Kármán vortex street formed immediately behind the vortex merging position. As the merged vortices passed through the wake, they produced shear layers that stretched from the walls to the channel centerline. In Fig. 16(c), the shear layers induced diffusion of the thermal mixing region throughout the channel centerline and enhanced heat transfer in the wake. For H/L = 2.25–2.75, the flapping frequency of the upper flag has the same value of the flapping frequency of the lower flag, as shown in Fig. 14(c) and (d). Unlike H/L = 1.5–2.0, the flags oscillated out of phase continuously with a constant phase angle for H/L = 2.25–2.75. The flapping frequency ðf f H=L¼2:25 ¼ 0:471Þ was twice the vortex frequency ðf V H=L¼2:25 ¼ 0:235Þ. Single vortices merged and coalesced to form a single vortical structure in the wake. The vortex shedding frequency corresponded to the dominant frequency of the net heat flux ðf Q H=L¼2:25 ¼ 0:235Þ in Fig. 14 (c). The flexible flags oscillated with a constant frequency; therefore, the variations in the net heat flux displayed a periodic trend,

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

1013

Fig. 16. (a) Time history of the net heat flux, snapshots of the (b) vorticity and (c) temperature fields at Re = 600, c ¼ 0:04, G/L = 2.6, H/L = 2.0.

Fig. 17. (a) Time history of the net heat flux, snapshots of the (b) vorticity and (c) temperature fields at Re = 600, c ¼ 0:04, G/L = 2.6, H/L = 2.25.

as shown in Fig. 17(a). The vortex strength and magnitude in the H/L = 2.25 system decreased compared to that in the H/L = 2.0 system due to a decrease in the flapping amplitude, as shown in Fig. 13(b). As H/L increased to 2.25, the block ratio decreased,

and the channel core region adjacent to the flexible flags became more open than that for H/L = 2.0. The inlet flow passed through along the channel centerline with decreasing the block ratio and dragged the generated vortices into the wake. The vortices merged

1014

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015

Fig. 18. (a) Mean heat flux, (b) mechanical energy loss via the pressure drop, and (c) thermal efficiency at Re = 600, c ¼ 0:04, G/L = 2.6.

further downstream in the H/L = 2.25 system than those observed in the H/L = 2.0 system, as shown in Figs. 16(a) and 17(a). The merged vortices formed a single vortical structure that swept up the thermal boundary layer and enhanced thermal mixing between the near-wall fluid and the channel core flow. As H/L decreased from 2.25 to 2.0, the distance between the merged vortices and the walls decreased sufficiently that the angular momentum of the merged vortices stretched the shear layer from the wall to the channel centerline. In Fig. 16(b) and 17(b), the vortices in the H/L = 2.25 system merged farther from the upper flag than in the H/L = 2.0 system. The reverse Kármán vortex street for H/L = 2.25 formed further downstream than it did for H/L = 2.0. Fig. 17c) shows that the thermal mixing region along the channel centerline was less intense than for H/L = 2.0. Fig. 18 plots the mean heat flux, mechanical energy loss, and thermal efficiency factor as a function of H/L. Fig. 18(a) shows that the mean heat flux at H/L = 2.0 was 117% of the value obtained at H/L = 2.25. As H/L decreased from 2.25 to 2.0, the flapping instability that increased the hydrodynamic force affected the flapping frequency. For H/L = 1.5–2.0, the flapping frequencies differed, and the mechanical energy loss was higher than it was for H/L = 2.25–2.75. The thermal efficiency factor was 114% at H/L = 2.0 and had the maximum value of 1.29 at H/L = 1.75.

position at which the reverse Kármán vortex street was initially forged and the presence of a shear layer that stretched from the wall to the channel centerline. Under the optimal parameter set, the CCW vortices in the upper half of the channel merged immediately behind the upper flag. The merged CCW vortices comprising the reverse Kármán vortex street induced the formation of a shear layer that began at the bottom wall. The shear layer, which was identified to CW vorticity at the lower half of the channel, extended from the wall to the channel centerline. The shear layer, which was associated with the CCW vorticity in the upper half of the channel, was generated in the same manner. As the shear layers passed through the reverse Kármán vortex street, the thermal mixing region diffused broadly across the channel centerline, which enhanced the thermal efficiency. The asymmetric configuration of the thermally conductive flags provided a superior thermal efficiency compared to other heat transfer enhancement techniques proposed in previous studies. The heat sink channel system proposed in the present study offered a convective heat transfer of 207% and a thermal efficiency of 135% of the values obtained from the symmetrically configured flag system. Conflict of interest None declared.

4. Conclusions

Acknowledgments

We examined the kinematics of flexible flags and the effects of their configurations on the flow behavior and heat transfer in a channel flow by adopting the penalty immersed boundary method. Numerical analysis indicated that a FAC system should enhance thermal efficiency. The thermal efficiency of the FAC system was compared with that of the FSC system described previously [24]. The vortical structures detached from the FSC displayed a symmetric configuration with respect to the channel centerline, whereas the FAC generated vortical structures in an asymmetric configuration. The mean vorticity in the FAC was higher than that in the FSC. An asymmetric flag configuration increased the mean vorticity, which encouraged more active thermal mixing in the wake than was observed in the FSC. The ratio between the cross-sectional area occupied by the flags and the channel height decreased. For the FAC, the increase in the mean vorticity and the decrease in the block ratio increased the thermal efficiency compared to the FSC system. A parametric study was performed to assess the effects of G/L and H/L on the mean heat flux, the mechanical energy loss, and the thermal efficiency. The optimal range of G/L and H/L values was determined by evaluating the thermal efficiency. Two key factors for obtaining a high thermal efficiency were identified: the

This study was supported by the Creative Research Initiatives program (No. 2017-013369) of the National Research Foundation of Korea (MSIP) and the Global Leading Technology program (N10042433) of the Ministry of Trade, Industry and Energy. References [1] R.L. Webb, E.R.G. Eckert, R.J. Goldstein, Heat transfer and friction in tubes with repeated-rib roughness, Int. J. Heat Mass Transfer 14 (1971) 601–617. [2] J.C. Han, L.R. Glicksman, W.M. Rohsenow, An investigation of heat transfer and friction for rib-roughened surfaces, Int. J. Heat Mass Transfer 21 (1978) 1143– 1156. [3] D.L. Gee, R.L. Webb, Forced convection heat transfer in helically rib-roughened tubes, Int. J. Heat Mass Transfer 23 (1980) 1127–1136. [4] J.X. Zhu, M. Fiebig, N.K. Mitra, Numerical investigation of turbulent flows and heat transfer in a rib-roughened channel with longitudinal vortex generators, Int. J. Heat Mass Transfer 38 (1995) 495–501. [5] G. Biswas, K. Torii, D. Fujii, K. Nishino, Numerical and experimental determination flow structure and heat transfer effects of longitudinal vortices in a channel flow, Int. J. Heat Mass Transfer 39 (1996) 3441–3451. [6] M. Fiebig, Embedded vortices in internal flow: heat transfer and pressure loss enhancement, Int. J. Heat Fluid Flow 16 (1995) 376–388. [7] M. Reeder, M. Samimy, The evolution of a jet with vortex-generating tabs: realtime visualization and quantitative measurements, J. Fluid Mech. 311 (1996) 73–118.

J.B. Lee et al. / International Journal of Heat and Mass Transfer 116 (2018) 1003–1015 [8] J.K. Foss, K.B.M.Q. Zaman, Large- and small-scale vortical motions in a shear layer perturbed by tabs, J. Fluid Mech. 382 (1999) 307–329. [9] S. Dong, H. Meng, Flow past a trapezoidal tab, J. Fluid Mech. 510 (2004) 219– 242. [10] C. Habchi, T. Lemenand, D.D. Valle, H. Peerhossaini, Turbulent mixing and residence time distribution in novel multifunctional heat exchangers-reactors, Chem. Eng. Proc. 49 (2010) 1066–1075. [11] H.M. Kaci, C. Habchi, T. Lemenand, D.D. Valle, H. Peerhossaini, Flow structure and heat transfer induced by embedded vorticity, Int. J. Heat Mass Transfer 53 (2010) 3575–3584. [12] J.L. Fernandez, R. Poulter, Heat transfer enhancement by means of flag-type insert in tubes, Int. J. Heat Mass Transfer 30 (1987) 2603–2609. [13] P. Hidalgo, F. Herrault, A. Glezer and M. Allen, Heat transfer enhancement in high-power heat sinks using active reed technology, in: Proceedings of the 16th International Workshop on Thermal Investigations of ICs and Systems, Barcelona, Spain, 2010, pp. 1–6. [14] F. Herrault, P.A. Hidalgo, C.-H. Ji, A. Glezer, M.G. Allen, Cooling performance of micromachined self-oscillating reed actuators in heat transfer channels with integrated diagnostics, in: Proceedings of the IEEE 25th International Conference on MEMS, 2012, pp. 1217–1220. [15] P. Hidalgo, S. Jha, A. Glezer, Enhanced heat transfer in air cooled heat sinks using aeroelastically fluttering reeds, in: Proceedings of the 21st International Workshop on Thermal Investigations of ICs and Systems, Paris, France, 2015, pp. 1–6. [16] T.S. Park, H.J. Sung, A nonlinear low-Reynolds-number k-e model for turbulent separated and reattaching flows – (I) Flow field computations, Int. J. Heat Mass Transfer 38 (1995) 2657–2666. [17] T.S. Park, H.J. Sung, Development of a near wall turbulence model and application to jet impingement heat transfer, Int. J. Heat Fluid Flow 22 (2001) 10–18. [18] S. Ali, C. Habchi, S. Menanteau, T. Lemenand, J.-L. Harion, Heat transfer and mixing enhancement by free elastic flaps oscillation, Int. J. Heat Mass Transfer 85 (2015) 250–264. [19] S. Ali, S. Menanteau, C. Habchi, T. Lemenand, J.-L. Harion, Heat transfer and mixing enhancement by using multiple freely oscillating flexible vortex generators, Appl. Therm. Eng. 105 (2016) 276–289.

1015

[20] K. Shoele, R. Mittal, Computational study of flow-induced vibration of a reed in a channel and effect on convective heat transfer, Phys. Fluids 26 (127103) (2014) 1–24. [21] A.K. Soti, R. Bhardwaj, J. Sheridan, Flow-induced deformation of a flexible thin structure as manifestation of heat transfer enhancement, Int. J. Heat Mass Transfer 84 (2015) 1070–1081. [22] S.G. Park, B. Kim, C.B. Chang, J. Ryu, H.J. Sung, Enhancement of heat transfer by a self-oscillating inverted flag in a Poiseuille channel flow, Int. J. Heat Mass Transfer 96 (2016) 362–370. [23] S.G. Park, C.B. Chang, B. Kim, H.J. Sung, Simulation of fluid-flexible body interaction with heat transfer, Int. J. Heat Mass Transfer 110 (2017) 20–33. [24] J.B. Lee, S.G. Park, B. Kim, J. Ryu, H.J. Sung, Heat transfer enhancement by flexible flags clamped vertically in a Poiseuille channel flow, Int. J. Heat Mass Transfer 107 (2017) 391–402. [25] C.S. Peskin, The Immersed boundary method, Acta Numer. 11 (2002) 479–517. [26] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. [27] W.-X. Huang, S.J. Shin, H.J. Sung, Simulation of flexible filaments in a uniform flow by the immersed boundary method, J. Compu. Phys. 226 (2007) 2206– 2228. [28] W.-X. Huang, H.J. Sung, An immersed boundary method for fluid-flexible structure interaction, Compu. Methods Appl. Mech. Eng. 198 (2009) 2650–2661. [29] D. Kamensky, M.-Chen. Hsu, D. Schillinger, J.A. Evans, A. Aggarwal, Yuri Bazilevs, M.S. Sacks, T.J.R. Hughes, An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves, Compuy. Methods Appl. Mech. Eng. 284 (2015) 1005–1053. [30] D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comput. Phys. 105 (1993) 354–366. [31] S.J. Shin, W.-X. Huang, H.J. Sung, Assessment of regularized delta functions and feedback forcing schemes for an immersed boundary method, Int. J. Numer. Meth. Fluids 58 (2008) 263–286. [32] K. Kim, S.-J. Baek, H.J. Sung, Animplicit velocity decoupling procedure for the incompressible Navier-Stokes equations, Int. J. Numer. Meth. Fluids 38 (2002) 125–138. [33] W.-X. Huang, H.J. Sung, Three-dimensional simulation of a flapping flag in a uniform flow, J. Fluid Mech. 653 (2010) 301–336.