Chemical Engineering Science 207 (2019) 556–580

Contents lists available at ScienceDirect

Chemical Engineering Science

journal homepage: [Link]/locate/ces

Enhancement of heat and mass transfer in a microchannel via passive

oscillation of a flexible vortex generator
Abdolrahman Dadvand a,⇑, Soheil Hosseini a, Saharnaz Aghebatandish a, Boo Cheong Khoo b
a
Faculty of Mechanical Engineering, Urmia University of Technology, Urmia, Iran
b
Department of Mechanical Engineering, National University of Singapore, Singapore

h i g h l i g h t s

Heat and mass transfer increase in a microchannel by an elastic VG is studied.

Total Nu number has an increase of 18.46% with respect to the rigid case.

Darcy friction factor decreases 42.33% with respect to the rigid case.

Thermal performance factor increases 42% with respect to the rigid case.

Mixing index increases 16.86% with respect to the rigid case.

a r t i c l e i n f o a b s t r a c t

Article history: The advantages of both the passive and active heat and mass transfer enhancement techniques are used.
Received 4 March 2019 A beam (rigid or flexible) as vortex generator (VG) is placed downstream of a cylindrical obstacle on the
Received in revised form 23 June 2019 lower wall of a microchannel. The governing equations are solved using ALE approach. The elastic beam
Accepted 25 June 2019
oscillates due to the forces exerted by the periodic flow behind the obstacle on it leading to the formation
Available online 26 June 2019
and separation of periodic vortices from the tip of the beam. These vortices disrupt the thermal boundary
layer and prevent it from re-growing and hence increase the heat transfer rate dramatically. The total Nu
Keywords:
number increases 18.46%, the Darcy friction factor experiences a decrease of 42.33%, the thermal perfor-
Fluid-structure interaction
Arbitrary Lagrangian-Eulerian
mance factor increases 42% and the mixing index increases 16.86% with respect to the rigid beam case.
Passive oscillation Finally, the flexible beam used here would not experience failure in the laminar flow regime.
Flexible vortex generator Ó 2019 Elsevier Ltd. All rights reserved.
Neo-Hookean model

1. Introduction (1) Vortices disturb the boundary layer growth which acts as a
thermal insulator.
Heat and mass transfer phenomena have numerous and impor- (2) Vortices circulate the heated flow in the vicinity of the
tant applications ranging from every day used home appliances to heated wall to penetrate the cool core of the main flow
advanced/high-tech industries such as automotive, aerospace, air and vice versa.
conditioning and refrigeration, cooling of the electronic devices, (3) Vortices increase the local velocity in the vicinity of the
nuclear reactors and so on (Babar and Ali, 2019; Khattak and Ali, heated wall resulting in heat transfer enhancement.
2019; Rehman et al., 2019; Sajid and Ali, 2019; Webb and Kim,
2005). Therefore, heat and mass transfer enhancement has Vortex generation techniques are categorized into two groups,
received considerable attention from many researchers. For this namely, the passive and active techniques (Kakaç et al., 1987). In
purpose, different methods have been developed among which the passive techniques, the vortex is generated via surface or geo-
vortex generation in the fluid flow is deemed the most successful metrical modifications to the flow channel by incorporating inserts
(Webb and Kim, 2005). Vortex generation can lead to heat transfer or additional devices without any external power source. In the
enhancement for the following reasons (Ahmed et al., 2012): active techniques, the vortex is generated by applying an external
power source which can be either mechanical or electromagnetic.
Although the active techniques give rise to higher heat and mass
transfer enhancement than the passive ones, these techniques,
⇑ Corresponding author. nonetheless, have not shown much applicable potential as it is
E-mail address: [Link]@[Link] (A. Dadvand).

[Link]
0009-2509/Ó 2019 Elsevier Ltd. All rights reserved.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 557

Nomenclature

(a,b) position of center of cylindrical obstacle (m) TMAC tangential momentum accommodation coefficient
A area (m2) TV transverse vortex
c concentration (mol m3) TVG transverse vortex generator
C Right Cauchy-Green deformation tensor (m) VG vortex generator
Cf fanning friction factor (1)
Cp specific heat in constant pressure (J kg1 K1) Dimensionless numbers
d diameter of cylindrical obstacle (m) Kn Knudsen number
D mass diffusivity (m2 s1) Le Lewis number
e relative error (1) Ma Mach number
E elasticity modulus (Pa) Nu Nusselt number
f Darcy friction factor (1) Pe Peclet number
F deformation gradient tensor (1) Pr Prandtl number
h beam height (m), grid size (m) Re Reynolds number
H channel height (m) Sc Schmidt number
I unity tensor (1)
j Colburn factor (1) Greek symbols
J ratio of volumes (1) a thermal diffusivity (m2 s1)
k thermal conductivity (W m1 K1)
aV tangential momentum accommodation coefficient (1)
L beam distance from channel inlet (m) b thermal expansion coefficient (K1)
L channel length (m) c ratio of specific heats (1)
Lc characteristic length (m) e thermal expansion (1)
N total number of cells (1)
fT temperature jump coefficient (1)
p pressure (Pa), apparent order of convergence (1) g thermal performance factor (1)
00
q heat flux (W m2) k lame constant (Pa), mean free path (m)
r grid refinement factor (1) l lame constant (Pa), dynamic viscosity (Ps s)
S second Piola-Kirchhoff stress tensor (Pa)
m kinematic viscosity (m2 s1)
t time (s), beam thickness (m) q density (kg m3)
T temperature (K) rV viscos slip coefficient (1)
u solid displacement (m)
rT thermal slip coefficient (1)
U fluid velocity (m s1) s shear stress (Pa)
U mean inflow velocity (m s1)
V volume (m3)
W strain energy density (Pa m) Subscripts
(x,y) Cartesian coordinate (m) 0 initial
cr critical
D diameter
Acronyms E elastic
3D three dimensional
F fluid
ALE arbitrary Lagrangian-Eulerian G grid
CCW counter-clockwise H heat transfer
CFD computational fluids dynamics M mass transfer
CV coefficient of variations
R rigid
CVD counter-rotating vortex pair S solid
CW clockwise th thermal
FFT fast Fourier transformation tot total
FSI fluid-structure interaction
W wall
FVG flexible vortex generator
GCI grid convergence index
IB immersed boundary Superscripts
LBM lattice boltzmann method T transpose
LDA laser doppler anemometry * dimensionless
spatial average
LV longitudinal vortex
LVG longitudinal vortex generator total (spatial and time) average
00
MI mixing index flux
RVG rigid vortex generator

difficult and sometimes impossible to provide external power and hence are more effective on the heat transfer enhancement
input in many cases especially in the micro/nano scales. The vortex than the TVs. Due to their 3D nature the LVs are identified only
generated by these two techniques is mainly divided into trans- in 3D simulations.
verse vortex (TV) and longitudinal vortex (LV) based on the direc- The obstacle such as a flap, a tab, a winglet, a block or a blade
tion of its rotating axis. The axes of TVs and LVs are respectively placed in a channel for vortex generation is called vortex generator
perpendicular and parallel to the main flow direction. The LVs (VG). There are abundant numerical and experimental studies on
are stronger than the TVs and affect a broader region of the flow the heat and mass transfer enhancement using VGs. These studies
558 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

focus mainly on the effects of different parameters of VGs such as in regime. Similar results have been obtained by Fu and Tong
their geometry, their inter-distance, their arrangement and their (2002). Mirzaee et al. (2012) studied numerically the forced convec-
angle of attack on heat and mass transfer enhancement. Fiebig tion heat transfer in a microchannel with an oscillatory elastic VG
(1998) presented a fairly comprehensive review on different mounted inside the channel with possible application in heat
geometries of the VGs. The general consensus is that the rectangu- exchangers. They found that for a fixed Reynolds (Re) number, the
lar and triangular VGs have essentially the same effect on the heat elastic VG results in a smaller pressure drop along the channel and
transfer enhancement. In addition, the LVGs are more effective at the same time a higher heat transfer enhancement as compared
than TVGs in heat transfer enhancement. Moreover, under the to the rigid VG. Inspired by the motion of ciliated micro-organisms
same conditions the overall heat transfer increase due to VGs in occurring in nature, Khatavkar et al. (2007) proposed a design for
laminar flow regime is higher than that in turbulent flow regime. an active micro-mixer consisting of two artificial cilia in the form
Mohebbi et al. (2018) used Lattice Boltzmann Method (LBM) to of micro-actuators covering the channel wall. They modeled the
study the forced convection flow in a channel with extended sur- fluid and solid interactions using the fictitious-domain method
faces (blocks) attached to its both walls on the effect of different and found that the mixing effectiveness is higher when the actuators
blocks’ height and spacing between them. According to their find- are arranged on the same wall. Lambert and Rangel (2010) explored
ings, the Nusselt number increases with the height of the blocks. In the capacity of a flexible flap to increase mixing in a microchannel
contrast, it decreases with the spacing between the blocks. Higher numerically. The flap was actuated using a distributed external
Re number leads to formation of bigger and stronger vortexes force. Their results showed that mixing is enhanced for larger flap
which occupy the whole space between the blocks as well. Ma displacements and for dimensionless frequencies between 1 and 2.
et al. (2010) investigated experimentally the single-phase heat In addition, they concluded that optimal mixing would occur when
transfer enhancement using four pairs of rectangular block as LVGs the flap length is 2/3 times the microchannel height.
in a narrow rectangular channel. They observed that in the laminar Thus far, it was declared that the oscillatory VGs are more effec-
flow regime the VGs lead to an increase of 100.9% in heat transfer tive than the fixed ones and that the elastic VGs are more efficient
and an increase of 11.4% in pressure loss. These values are, respec- than the rigid ones in both the heat transfer and mass transfer
tively, 87.1% and 100.3% in the turbulent flow regime. Using both enhancement. The primary restriction with the active techniques
numerical simulations and laboratory experiments, Habchi et al. is the need for an external energy source and its practical imple-
(2010) studied the heat and mass transfer in a circular pipe by a mentation in small scales making the active techniques to be lim-
row of VGs made up of four diametrically opposed trapezoidal tabs. ited to the theoretical level. Recently, a new concept has been
They showed that the hairpin-like vortices as well as the counter- proposed that takes the advantages of both the active and passive
rotating vortex pair (CVP) generated behind the tabs had remark- techniques simultaneously. The idea is to oscillate the VGs by the
able effects on the heat and mass transfer enhancement. Habchi hydrodynamic forces of the fluid rather than by a costly external
et al. (2010) conducted computational fluid dynamics (CFD) simu- force. This method has no limitations in mini/micro size applica-
lation and laser doppler anemometry (LDA) measurements to tions such as micro heat exchangers and micro reactors. Using
investigate the arrangement effects of trapezoidal tabs on turbu- the immersed boundary method (IBM) Shoele and Mittal (2014)
lent mixing in multifunctional heat exchangers-reactors. They studied the flow-induced vibration of a reed and its effect on con-
studied three different tab arrangements, namely, classical aligned, vective heat transfer in a channel and found that heat transfer
alternating and reverse tab arrays. They found that the last enhances due to the vibration of the reed. Soti et al. (2015) simu-
arrangement could generate higher and more powerful CVP and lated the heat transfer enhancement in a heated channel laminar
hence give superior heat transfer efficiency than the first two, flow due to the flow-induced deformation of an elastic thin plate
but its pressure drop was about 1.4 times higher. The structure attached to lee side of a rigid cylinder using FE-IBM. Ali et al.
of flow field around three-dimensional obstacles was investigated (2015) investigated numerically the role of freely oscillating elastic
by Becker et al. (2002). They found that as the angle of attack is flaps on the mixing process and heat transfer in a two-dimensional
increased the vortices generated behind the obstacle are changed laminar flow. The oscillations of flaps generated an unsteady lam-
from TV to LV and thereby heat transfer increases. Depaiwa et al. inar flow with complex coherent vortices detaching periodically
(2010) investigated experimentally the forced convection heat from the tip of the flaps that, in turn, caused the flexible flaps
transfer and friction loss characteristics for turbulent airflow located downstream of the main flaps to oscillate. It was shown
through a constant heat flux channel solar air heater with rectan- that the mixture quality and the overall heat transfer were
gular winglet VG. Ten pairs of the winglet VGs with various attack enhanced up to 98% and 134% relative to those in the rigid cases.
angles of 60°, 45° and 30° placed in both the upstream direction Ali et al. (2017) studied numerically the heat transfer and mixing
and in the downstream direction. They found that the best thermal performances of a three-dimensional heat exchanger/reactor con-
performance was achieved when VGs are placed in the upstream sisting of a circular pipe in which five arrays of four equally spaced
direction and at the attack angle of 30°. trapezoidal VGs are inserted with an angle of 45° with respect to
All the aforementioned works are associated with the passive the pipe wall. They examined both the elastic VGs and rigid VGs
vortex generation techniques. In the following, we shall refer to and found that the elastic VGs and the rigid VGs could respectively
the works that have used active techniques. In almost all of these improve the overall heat transfer of about 118% and 97% with
works, a periodic mechanical or electromagnetic force is applied to respect to an empty pipe. In addition, the elastic VG configuration
oscillate the VG with controlled frequency and amplitude. Yang shows an extreme improvement of the mixture quality with an
and Chen (2008) investigated numerically the effect of transient increase of 195% with respect to the rigid VG. Li et al. (2017) devel-
flow field structures, and the heat transfer characteristics of heated oped a novel self-agitator to enhance the thermal performance of a
blocks in the channel with a transversely oscillating cylinder using a plate fin heat exchanger via its air-side heat transfer enhancement
Galerkin finite element (FE) formulation with Arbitrary Lagrangian– and evaluated it both numerically and experimentally. The oscilla-
Eulerian (ALE) method. They found that the oscillatory cylinder tion of the self-agitator could enhance the heat transfer in plate fin
enhances the heat transfer as compared to the fixed cylinder. Fur- and save the pumping power of up to 40% with the same rejected
thermore, they found that due to the phenomenon of resonance in heat. Park et al. (2017) developed the IBM to study the forced con-
the channel flow (i.e., when the oscillating frequency of the cylinder vective heat transfer around a flexible circular cylinder placed in a
approaches to the vortex shedding frequency), the heat transfer rate square cavity. They observed two different stable states of
is enhanced more remarkably. This phenomenon is called the lock- stretched-stable and self-sustained flapping depending on the Re
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 559

number and the flexible property of the cylinder. The convective For ease of comparison, the problem is subdivided into four
heat transfer decreased for the stretched-stable state due to the cases: (i) simple channel, (ii) channel with a rigid beam, (iii) chan-
elastic deformation, while enhancement of up to 7% for the self- nel with a rigid beam and a cylindrical obstacle and (iv) channel
sustained flapping state as compared to the heated rigid cylinder. with a hyperelastic beam and a cylindrical obstacle. The results
Lee et al. (2017) investigated the effects of two flexible flags in a of the last two cases in terms of Nusselt (Nu) number, Colburn fac-
symmetric configuration on heat transfer enhancement in a heated tor, thermal performance factor and mixing index are compared.
channel using the penalty IBM. They found that the flexible flags Finally, the failure probability of the hyperelastic beam is
significantly enhanced the thermal efficiency compared to rigid discussed.
flags. In addition, the flexible flags lead to an increase of about
185% in the net heat flux and 106% in the thermal efficiency factor 3. Assumptions
compared to the baseline (open channel) flow provided the opti-
mal parameter set are used. Lee et al. (2018) also investigated The Re number based on the mean velocity at the channel inlet
the dynamics of two flexible flags in an asymmetric configuration U and the obstacle diameter d (Red ¼ qUd=lÞ is set to 100 to ensure
and their effects on heat transfer enhancement in a heated channel generation of periodic flow behind the cylinder. The flow is consid-
using the penalty IBM. They found that the flexible flags lead to an ered incompressible (Ma < 0:3). In addition, the Re number based
increase of about 207% in the convective heat transfer and 135% in
on the channel’s hydraulic diameter ðReH ¼ qU2H=lÞ is set to
the thermal efficiency factor compared to the baseline (open chan-
1000, which is sufficiently smaller than the critical Re number
nel) flow provided the optimal parameter set are used. Eventually
(Recr ¼ 1500 Bejan and Kraus, 2003) to guarantee the flow in the
Yu et al. (2018) investigated the influence of the ratio of the length
channel to remain in the laminar flow regime. Furthermore, it is
of an inverted flag L to the channel width W (C* = L/W) on the gen-
assumed that the diffusion of momentum, energy and mass occur
erated vortex dynamics and heat transfer experimentally using
at the same rate (i.e. Pr ¼ Sc ¼ Le ¼ 1). Therefore, as Pe ¼ Re  Pr,
time-resolved particle image velocimetry technique. Three differ-
the Pe number for heat and mass transfer are obtained to be equal
ent values of C* were examined depending on which three distinct
to 1000 (i.e. PeH ¼ PeM ¼ 1000), which states that the diffusion is
vortex dynamics, i.e., the von Kármán vortex street, the G mode,
negligibly small. The material of beam is assumed to be silicone
and the singular mode, were observed. It was found that increase
rubber that is widely used in MEMS and NEMS technologies.
in C* strongly increased heat transfer from the channel wall.
Table 1 gives the mechanical properties of this material (https://
In the present work, a fixed rigid cylindrical obstacle is placed at
[Link]/[Link]?ArticleID=920).
a distance from the channel inlet to generate periodic flow behind
the cylinder (Re ¼ 100). This periodic flow causes a flexible hyper-
4. Governing equations
elastic beam placed downstream of the cylinder to oscillate. The
results are compared with the corresponding case of the rigid
4.1. Fluid domain
beam. Simulation of the motion and deformation of hyperelastic
beam and its interaction with the surrounding fluid is called
Several methods have been proposed for the solution of FSI
fluid-structure interaction (FSI) and needs coupling of fluid and
problem involving moving boundaries such as Arbitrary Lagran-
solid equations (Morand and Ohayon, 1995). We use the finite ele-
gian Eulerian (ALE) method and immersed boundary (IB) method
ment based commercial software COMSOL Multiphysics version
(Van Loon et al., 2007; Peskin, 2002). In the present work, the
5.4. Time integration is carried out with backward differentiation
ALE method is used. In this method the fluid-solid interface is
(Euler) formulae. In addition, P2 =P1 Taylor–Hood finite elements
tracked by mesh and hence the mesh velocity Ug appears in the
are used to discretize the velocity components and pressure, and
P2 finite elements are chosen for temperature and concentration. advection term of the governing equations (for more information
the reader is referred to (Donea et al., 2004). The continuity,
momentum, heat and mass transfer equations governing the fluid
2. Problem description flow are as follows:

The main problem under consideration involves laminar and

rU¼0 ð1Þ
incompressible flow passing through a microchannel which is
@U

  1
3000 lm long and 500 lm high. The schematic of the channel þ U  Ug  r U ¼ mr2 U  rp ð2Þ
along with corresponding nomenclature is given in Figs. 1 and 2. @t qf
The microchannel includes a hyperelastic beam with height
h ¼ 250 lm and thickness t ¼ 20 lm placed vertically on the lower @T

þ U  Ug  rT ¼ ar2 T ð3Þ
wall of the channel at a distance l ¼ 1000 lm from the inlet. The @t
tip of the beam is filleted with radius 10 lm. A cylindrical obstacle
@C

of diameter d ¼ 100 lm is located at a distance a ¼ 250 lm from þ U  Ug  rC ¼ Dr2 C ð4Þ
@t
the channel inlet on the centerline of the channel.

Fig. 1. Schematic of the channel with nomenclature.

560 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 2. Fluid and solid domains and boundary conditions.

Table 1 where I1 C is the first invariant of C tensor, l and k are the Lame con-
Mechanical properties of silicone rubber ([Link]
stants that are related to the elasticity modulus and Poisson’s ratio
ArticleID=920).
as follows,
Property Symbol value Unit
E
Density qs 2500 kg=m3 l¼ ð9Þ
Elasticity modulus Es Pa
2ð 1 þ m Þ
5  106
Poison ratio ms 0.48 
Thermal conductivity ks 2.5 W=ðm  KÞ mE
k¼ ð10Þ
Specific heat Cp s 1200 J=ðkg  KÞ ð1 þ mÞð1  2mÞ
Thermal expansion bs 275 1=K
The ratio of volumes in deformed state (V) and at rest (V0) con-
figurations J is defined as the determinant of deformation gradient
where U denotes the fluid velocity vector, p is the pressure, qf is the tensor F, i.e.,
fluid density, m is the kinematic viscosity, t is time, T is temperature, V
a is the thermal diffusivity, C represents the mass concentration and J¼ ¼ det½F  ð11Þ
V0
D is the mass diffusivity.
8
< J > 1; expansion
>
4.2. Solid domain
J ¼ 1; isochoric
>
:
The equation of motion of the solid domain (hyperelastic beam) J < 1; contraction
is as follows, The continuity equation (1), momentum equation (2) and mass
2
@ u transfer equation (4) are solved for the fluid domain and equation
qs ¼ r  FS ð5Þ (5) is solved for the solid domain. However, the energy equation (3)
@t2
is solved in both the fluid and solid domains thereby the tempera-
where u denoted the solid displacement, F is the deformation gra- ture distribution in hyperelastic beam and its thermal expansion
dient tensor, S is second Piola-Kirchhoff stress tensor defined as, are obtained. The thermal strain due to the temperature difference
F ¼ I þ ru ð6Þ DT is defined as,
eth ¼ bDT ð12Þ
@W
S¼2 ð7Þ where b denotes the thermal expansion coefficient.
@C
In these two relations I is the unity tensor, W is the strain
4.3. Initial and boundary conditions
energy density function and C represents the right Cauchy-Green
T T
deformation tensor, which is defined as C ¼ F F, where F is the
4.3.1. Initial conditions
transpose of F.
The initial conditions of the problem are:
Different hyperelastic material models have been developed
based on the different strain energy density functions W. The sim- Uðx; y; 0Þ ¼ 0 m=s; Pðx; y; 0Þ ¼ 0 Pa; Tðx; y; 0Þ
plest hyperelastic material model is the Saint Venant–Kirchhoff
¼ 300 K; Cðx; y; 0Þ ¼ 1 mol=m3 and uðx; y; 0Þ ¼ 0 m:
model (Bower, 2009) which is just an extension of the linear elastic
material model to the nonlinear regime. Among other models one
may refer to the Neo-Hookean model (Attard, 2003), Mooney-
4.3.2. Hydrodynamic boundary conditions
Rivlin model (Rivlin and Saunders, 1951), Ogden model
The problem under consideration is microfluidic system as at
(Holzapfel, 2000); Varga model (Holzapfel, 2000); Arruda-Boyce
least one of the dimensions of the system is of micro size (Bruus,
model (Arruda and Boyce, 1993); Storakers model (Storåkers,
2008). In addition, the Knudsen number defined as Kn ¼ k=Lc
1986), Gent model (Gent, 1996) and Gao model (Gao, 1997). In
(where k is the mean free path and Lc is the characteristic length
the present work, the Neo-Hookean model (Attard, 2003) is
in micro size) lies between 0.1 and 0.001, and hence the modified
employed where the strain energy density functions W is defined
boundary conditions (Maxwell’s slip condition for velocity and Von
as,
Smoluchowski temperature jump for temperature (Karniadakis
1
C  1 et al., 2006) are enforced on the channel walls and on the beam
W¼ l I1  3  l ln J þ kðln JÞ2 ð8Þ
2 2 boundary.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 561

The Maxwell’s slip condition is used for the velocity on the solid The zero-heat flux (adiabatic) boundary condition (q00 H ¼ 0) is
walls, applied on the circular obstacle. As the energy equation is solved
for the fluid and solid domains, the beam wall does not need ther-
dU x d
U ¼ U w þ rV k þ rT m ðln T w Þ ð13Þ mal boundary condition.
dy dx
where Uw ¼ 0 m=s is the wall velocity, and rV and rT are the vis- 4.3.4. Mass transfer boundary conditions
cous slip and temperature slip coefficients, respectively, defined At the inlet of the channel two fluids with the same properties
as rV ¼ ð2  aV Þ=aV and rT ¼ 3=4. The coefficient aV denotes the but with different mass concentrations of 0 mol/m3 and 1 mol=m3
tangential momentum accommodation coefficient (TMAC) indicat- enter the channel. At the outlet of the channel the concentration is
ing the portion of the incident molecules that reflect diffusely. It extrapolated. At all other boundaries the mass flux is set to zero
takes a value between 0.85 and 1.0 (Karniadakis et al., 2006; (q00 M ¼ 0).
Sharipov, 2011), which is considered equal to 0.9 in the present
work. 4.3.5. Boundary conditions for solid domain
The flow in the channel entrance is considered fully developed. On the hyperelastic beam-fluid interface the normal stresses are
At the inlet of the channel a parabolic velocity profile is given as in equilibrium, i.e., rf  n ¼ rs  n. In addition, the base of the beam
follows, has always zero displacement (u ¼ 0).
 
3 4y
U inlet ð0; y; t Þ ¼ U 2 ðH  y Þ ð14Þ
2 H 5. Results and discussion

where U ¼ 1 mm=s is the mean velocity at the inlet and H is the 5.1. Validation
channel height.
Finally, the outlet gage pressure is set to zero, i.e., 5.1.1. Validation I
Poutlet ðL; y; tÞ ¼ 0. The numerical benchmark proposed by Turek and Hron (2006)
is used for validation of the present FSI problem. The benchmark is
4.3.3. Heat transfer boundary conditions associated with the laminar incompressible flow around an elastic
The temperature at the channel inlet is set to 300 K, i.e., plate attached to a cylindrical obstacle. The elastic plate oscillates
Tinlet ð0; y; tÞ ¼ 300 K. At the outlet the temperature is not specified due to the hydrodynamic forces exerted by the flow upon it. The
but extrapolated. According to the Von Smoluchowski temperature dimensions (in mm) of the problem are given in Fig. 3. The elastic
jump, the fluid temperature near the top and bottom walls is given plate is slightly shifted from the horizontal centerline in order to
as, make the flow become axisymmetric.
dT The boundary conditions considered for this problem include
T ¼ T w þ fT k ð15Þ the parabolic velocity profile (fully developed flow) at inlet, zero
dy
gage pressure at outlet, no-slip conditions at the walls. The Re
where Tw ¼ 360 K is the wall temperature. In addition, fT is the number is considered to be equal to 100. The physical properties
temperature jump coefficient defined as follows, of the problem, which is identical to Case FSI2 of (Turek and
2  aV 2c Hron, 2006), are given in Table 2.
fT ¼ ð16Þ Fig. 4 depicts the displacement of the tip of the elastic plate in x
aV ð1 þ cÞPr
and y directions. After approximately 11 s, the elastic plate shows
Here c is the ratio of specific heats taken to be equal to 1 and Pr harmonic oscillations with constant amplitude and frequency.
denotes the Prandtl number, which takes a value of 1 (see Fig. 5 shows the transient drag and lift forces exerted by the flow
Section 3). on the cylinder and plate.

Fig. 3. Schematic of the benchmark proposed by Turek and Hron (2006).

Table 2
Physical properties of the problem, which is identical to Case FSI2 of (Turek and Hron, 2006).
   
qs kg=m3 Es ½Pa ms ½ qf kg=m3 lf ½Pa  s U ½m=s Re½

10,000 1:4  10 6 0.4 1000 1 1 100

562 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 4. Displacement of the tip of the elastic plate.

Fig. 5. Drag and lift forces exerted by the flow on the cylinder and the elastic plate behind it.

The average and amplitude of a parameter at the final period of flow wise inside the channel (the blockage ratio is b ¼ 1=4Þ. Walls
the oscillation can be obtained using the following relations (Turek are considered adiabatic except in a portion of the bottom wall
and Hron, 2006); located downstream of the obstacle where the temperature is
maintained at 353 K, while the temperature of the incoming flow
1
Average ¼ ðmax þ minÞ ð17Þ is 293 K. The viscosity and thermal conductivity are considered
2 quadratic functions of temperature. For more details on the prob-
lem one may refer to (Meis et al., 2010). The Results are depicted
1
Amplitude ¼ ðmax  minÞ ð18Þ in Fig. 8 in terms of temperature contours (Fig. 8a), vorticity con-
2 tours (Fig. 8b) and time-averaged local Nu number (Fig. 8c). A good
In addition, the oscillations frequency is obtained by applying agreement is observed between the current simulation and that of
the fast Fourier transformation (FFT) to the displacement signal Meis et al. (2010).
(see Figs. 6 and 7).
The results are given in Table 3 and are compared with their 5.2. Grid independence test
counterparts given in (Turek and Hron, 2006). There is good agree-
ment between the current results and those in (Turek and Hron, In order to determine the grid size for which the results are grid
2006). independent, we have used the method proposed by Celik et al.
(2008). In this method the cell size h is defined as:
5.1.2. Validation II " #1=2
In the second validation, following Meis et al. (2010), a 2D, 1 XN

unsteady, laminar flow (ReH = 600) of water in a non-isothermal h¼ DA i

N i¼1
microchannel is considered. A cylindrical obstacle is placed cross-
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 563

Fig. 6. FFT of displacement signal of the tip of the plate.

Fig. 7. FFT of the signal of the drag and lift forces.

Table 3
Comparison between present results and ref (Turek and Hron, 2006).

X-component of displacement of the beam tip Y-component of displacement of the beam tip
Present results Turek and Hron (2006) Error (%) Present results Turek and Hron (2006) Error (%)
Mean ½mm 14.54 14.58 0.27% 1.3 1.23 5.7%
Amplitude ½mm 12.54 12.44 0.8% 80.69 80.06 0.8%
Frequency ½Hz 3.74 3.8 1.58% 1.87 2 6.5%
Drag force Lift force
Present results Turek and Hron (2006) Error (%) Present results Turek and Hron (2006) Error (%)
Mean ½N=m 203.95 208.83 2.3% 0.96 0.88 9.1%
Amplitude ½N=m 69.48 73.75 5.8% 225.13 234.2 3.87%
Frequency ½Hz 3.72 3.8 2.1% 1.87 2 6.5%

where DAi denotes the area of ith cell and N is the total number of In Section 4.1 it was mentioned that in the ALE framework the
cells. Then three different grid sizes are chosen so that the grid fluid/solid mesh deforms following the movement/deformation of
refinement factor r ¼ hcoarse =hfine > 1=3: For more details of the the solid body (see Fig. 10).
method, the reader is referred to Celik et al. (2008). Table 4 gives The total Nusselt number Nutot and the mixing index MI are
the specifications of three different grid sizes. Fig. 9 depicts the important parameters that are used here for grid independence
unstructured mesh #3. It can be seen that the grid has been refined study. The method proposed by Celik et al. (2008) is employed.
near all the walls of the problem. The results are given in Table 5, where P represents the apparent
564 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 8. Comparison of the Instantaneous temperature contours (a), Instantaneous vorticity contours (b) and time averaged local Nusselt number (c) for laminar flow in a non-
isothermal microchannel with a circular obstacle with blockage ratio b ¼ 1=4 and Re = 600.

x  y
Table 4
x ¼ ;y ¼ ð19Þ
Specifications of the grids used for grid independence test. H H
Mesh #1 Mesh #2 Mesh #3 Fig. 11 depicts the variations of the time averaged Nusselt num-
Total number of cells N 8527 14,609 31,924 ber (Nu) along the channel for different Re numbers. The time aver-
h P i12
1:32  108 1:01  108 0:68  108 aged Nusselt number (Nu) is defined as follows,
Grid sizeh ¼ N1 Ni¼1 ðDAi Þ

Grid refinement factor r ij ¼ hj =hi  1.307 1.485

Z t
1
NuðxÞ ¼ Nuðx; tÞdt ð20Þ
order of convergence, e21 denotes relative error, GCI21 is the grid t 0

convergence index. For a fixed Re number, the thickness of the thermal boundary
layer along the channel increases and hence the heat transfer coef-
5.3. Simple channel ficient and/or the Nu number decrease. This is clearly observed in
Fig. 11. In addition, the Nu number increases with the Re number.
First the variations of Nu number along the channel for different The ratio of the Re number based on the channel height ReH and
Re numbers of 100–1000 are investigated. The Re number is the Re number based on the diameter of the cylindrical obstacle
defined as ReH ¼ qUð2HÞ=l . Red is considered equal to 10, i.e.,ReH =Red ¼ 2H=d ¼ 10. So, since
The special dimensions are normalized by the height of the we need Red ¼ 100 for the flow behind the cylindrical obstacle to
channel H. be oscillatory, ReH is set to be equal to 1000.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 565

Fig. 9. Non-uniform unstructured mesh, which is refined near the walls.

Fig. 10. Tracking the fluid-solid interface by moving mesh in ALE method.

Table 5
Results of grid independence test.
nel with the rigid beam (see Section 5.4)) becomes steady after a
P e21 GCI21 certain time. For the cases involving the cylindrical obstacle (see
Nutot 1.45 0.67% 1.78% Sections 5.5 and 5.6), however, the flow is unsteady for
MI 1.32 1.03% 2.26% Red ¼ 100. Nevertheless, for all the cases studied in the present
work including both the channel with cylindrical obstacle (see Sec-
tions 5.5 and 5.6) and without it (see Sections 5.3 and 5.4), the sim-
It should be noted that the flow in the channel without cylindri- ulations are carried out based on the unsteady flow assumption.
cal obstacle (i.e., the simple channel (see Section 5.3) and the chan- Therefore, all the diagrams associated with the variations of Nu
566 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 11. Variations of Nu number versus the length of the channel.

number along the channel length are depicted in terms of time jsw j
Cf ¼ 1 ð21Þ
averaged Nu number. q U2
2 f

5.4. Effect of rigid beam on heat transfer where sw is the shear stress at the lower wall of the channel. Fig. 13
depicts the variations of C f along the lower wall. It can be seen that
In this section the effect of rigid beam attached to the lower at x ¼ 5:15, the friction factor C f ¼ 0:
wall on the heat transfer is investigated. The channel considered It is very important here to investigate the variations of the Nu
in this section does not include a circular obstacle. Fig. 12 demon- number along the top wall of the channel. In the case of simple
strates the variations of Nu number along both the top and bottom channel due to the symmetry of the channel with respect to the
walls of the channel. A long fixed vortex with negligible effect on centerline of the channel, the Nu number along the bottom and
the Nu number is generated behind the beam. In fact, the fixed vor- top walls is the same. However, when the rigid beam is added to
tices act as insulation against heat transfer and hence reduce the the channel at x ¼ 2, the geometry becomes asymmetric. At the
heat transfer. The major increase in the Nu number occurs around location of the beam the flow passage becomes narrower. As a
the reattachment point where the cold core fluid is injected to the result, the local flow velocity near the top wall increases and hence
hot fluid near the lower hot wall. Thereafter the thermal boundary the Nu number along the top wall of the channel is enhanced.
layer grows again and hence the Nu number decreases. At the reat- The time and spatial averaged Nu number is defined as,
tachment point the Fanning friction factor C f becomes zero. The Z L
1
Fanning friction factor is defined as, Nu ¼ NuðxÞdx ð22Þ
L 0

100

Fig. 12. Effect of rigid beam on the Nu number.

 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 567

Fig. 13. Variations of the Fanning friction factor along the lower wall of the channel.

According to Fig. 12, the total Nu number (i.e., due to the increased local flow velocity near the walls as a result of
Nutot ¼ NuTop þ NuBottom ) associated with the second case the obstacle presence. In this case, the total Nu number is equal to
(Nutot ¼ 30:23) is greater than its counterpart for the simple chan- 32.40 (see Fig. 14).
nel case (Nutot ¼ 26:42).
5.6. Effect of cylindrical obstacle and hyperelastic beam on heat
transfer
5.5. Effect of cylindrical obstacle and rigid beam on heat transfer
In this section, the rigid beam is replaced by a hyperelastic one.
In this section the effect of cylindrical obstacle on the heat The problem becomes a time-dependent one. The dimensionless
transfer is investigated. The Re number based on the cylinder’s time is defined as,
diameter is set to be equal to 100, for which the flow behind the
tU
obstacle becomes periodic (Zdravkovich, 1997). The periodic vor- t ¼ ð23Þ
tices (von-Karman Vortex Street) lead to better mixing and hence H
mass transfer in the microchannel. However, they have less effect The vortices generated behind the cylindrical obstacle, cause
on the heat transfer rate as they are far enough from the walls. At the hyperelastic beam to oscillate. These oscillations are formed
the location of the obstacle ðx ¼ 0:5Þ, a slight increase in the heat due to the periodic forces which are exerted by the flow on the
transfer from both the top and bottom walls is observed, which is beam without any external power source. The oscillation mecha-

Fig. 14. Effect of cylindrical obstacle on the Nu number in a channel with RVG.
568 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

nism is described as follows. When a vortex shedded from the saved in the beam and its springback property, the beam moves
cylindrical obstacle hits the beam, it causes a forward displace- backward. This scenario repeats periodically causing the beam to
ment of the beam. As the vortex leaves the beam, due to the energy oscillate.

Fig. 15. Velocity field at different times for the case of elastic beam.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 569

In the case of rigid beam, the vortex generated behind the beam boundary layer starts re-growing which decreases Nu number.
is stationary, which has not sensible effect on the heat transfer However, when the beam becomes hyperelastic, the vortex gener-
enhancement. In this case after the reattachment point the thermal ated behind the beam is oscillatory. This is due to the fact that in

Fig. 16. Pressure contours at different times for the case of elastic beam.
570 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

the case of hyperelastic beam, whenever the beam moves back- of a vortex, which is separated from the tip during the forward
ward in the opposite direction of the main flow, a vacuum region motion of the beam. At the next period of the oscillation, another
is generated behind the tip of the beam. This leads to the formation vortex is generated and shedded. Therefore, unlike the rigid beam

Fig. 17. Vorticity contours at different times for the case of elastic beam.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 571

where a stationary vortex is generated behind the beam, in the The time-averaged Nu number (Nu) along both the top wall and
hyperelastic beam a train of the vortices form and move toward bottom wall of the channel is depicted in Fig. 20. It can be seen that
the outlet of the channel destructing the thermal boundary layer, the Nu number is high near the channel inlet because the temper-
mixes the cold fluid of the core and the hot fluid close to the heated ature gradient is high there. It then decreases due to the growth of
bottom wall and hence increases the heat transfer dramatically. the thermal boundary layer. Then due to the presence of the beam
The time-dependent flow regime in terms of velocity, pressure
the Nu number decreases further for the bottom wall because of
and vorticity fields is shown in Figs. 15–17, respectively.
the low flow velocity but increases for the top wall because of
From the pressure contours (Fig. 16) it is clearly evident that the
increase in the flow velocity. After the location of the beam, due
pressure just upstream of the hyperelastic beam is greater than
to the sudden expansion, the flow velocity decreases near the top
that just downstream of it. This leads to a negative pressure gradi-
wall of the channel and hence the Nu number decreases. For the
ent that causes the beam to bend in the flow direction. This pres-
bottom wall, however, the formation of the train of vortices behind
sure gradient becomes oscillatory after t ¼ 10 due to the
the beam destruct the thermal boundary layer and also combine
oscillatory (periodic) nature of the flow behind the cylindrical
the cold fluid in the core of the channel with the hot fluid adjacent
obstacle. The oscillatory pressure gradient causes the hyperelastic
beam to oscillate accordingly. to the wall. Therefore, the Nu number increases again. These find-
In Fig. 17 the vortices shown in blue color are clockwise (CW) ings can also be deduced from the temperature contours (Fig. 19)
and those shown in red color are counter-clockwise (CCW). The so that the Nu number is related directly to the temperature gradi-
important point is the formation of CCW vortex near the top wall ent on the corresponding wall. The total average Nu number asso-
of the channel. This vortex is stationary for the rigid beam case, ciated with the top wall NuTop is equal to 20.79, which is about the
but it is oscillatory for the case with hyperelastic beam. Although same as that associated with the rigid beam case. However, its
such a vortex was reported in the previous works using rigid beam,
value for the bottom wall NuBottom equals 17.59 showing remark-
the reason behind it was not revealed. This may be attributed to
able increase as compared with the rigid beam case. This may be
the fact that at the location of the beam the cross section of the
attributed to the oscillatory nature of the vortices generated
flow passage decreases and hence the flow velocity increases there
behind the hyperelastic beam which prevent the boundary layer
(according to the conservation of mass). However, after the beam, a
from being reconstructed after being destroyed by the vortices.
sudden expansion occurs and hence the flow velocity decreases
The total Nu number Nutot is equal to 38.38 showing an increase
near the top wall of the channel. This leads to an increase in the
of 18:46% with respect to its rigid case counterpart.
pressure (due to the Bernoulli effect) and hence to an undesirable
The Colburn factor represents the ratio of thermal energy trans-
pressure gradient near the top wall. This positive pressure gradient
ferred to the mechanical energy consumed, i.e.,
results in flow separation and vortex generation near the top wall
of the channel. This is also deduced from the variations of the Nutot
j¼ 1
ð24Þ
dimensionless time-averaged pressure along the top wall of the ReH Pr 3
channel (see Fig. 18). It is observed that after the point x ¼ 3:4
The value of Colburn factor for the case of hyperelastic beam is
the pressure gradient along the top wall is positive.
0.03838, which is 18.46% greater than its counterpart of the rigid
The instantaneous snapshots of the temperature field are
beam case.
depicted in Fig. 19. The oscillation of the flexible beam leads to
The hyperelastic beam represents lower resistance against flow
the formation of a train of vortices near the bottom wall behind
than the rigid one. The variation of thermal performance factor of
it that destruct the thermal boundary layer and also combine the
the channel from the elastic beam case (e) to the rigid beam case
cold fluid in the core of the channel with the hot fluid adjacent
(r) is defined as (Bejan and Kraus, 2003);
to the bottom wall. These would result in the increase in the tem-
perature gradient and hence the Nu number as compared with the   13
Nue fe
rigid beam case. g¼ ð25Þ
Nur fr

Fig. 18. Variations of the dimensionless time-averaged pressure along the top wall of the channel.
572 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 19. Temperature field at different times for the case of elastic beam.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 573

Fig. 20. Time-averaged Nu number in the case of hyperelastic beam.

Table 6 field in the case of rigid beam at t ¼ 30. It is observed that the two
Comparison between the thermal performance of the channel with the rigid beam fluids have been mixed quite well.
and hyperelastic beam cases.
Fig. 22 depicts the concentration field associated with the
Rigid beam Hyperelastic beam Variations (%) hyperelastic beam case. It is observed that the concentration is
Nutot 32.40 38.38 18.46%" quite smooth, i.e., the two fluids have been mixed quite well. This
j 0.0324 0.03838 18.46%" can be deduced quantitatively from Fig. 23. This figure depicts the
f 1.89 1.09 42.33%# variations of concentration along two horizontal lines of y ¼ 0:25
g   42%"
and y ¼ 0:75 at time t  ¼ 30. The concentration approaches 0.5 at
the outlet of the channel. The velocity vectors related to the con-
centration field (Fig. 22) are demonstrated in Fig. 24.
The criterion by which the mixing quality in any cross section of
where f is the Darcy friction factor, which is defined as,
the microchannel is evaluated is the mixing index (MI) defined as
  (Lambert and Rangel, 2010);
Dh Dp
f ¼ ð26Þ
L 12 qU 2 MI ¼ 1  CV ð27Þ

The value of g experiences 42% increase when the rigid beam is where CV denotes the coefficient of variations, which is defined as
replaced by the hyperelastic one in the channel (see Table 6 for the ratio of standard deviation to the average,
more detail). "P #12
n
It should be noted that the CPU time is about 3 h and 2 min 1 i¼1 ðc i  cmean Þ2
CV ¼ ð28Þ
using a computer with the specifications of Intel Core i5-2430M cmean n
@ 2.40 GHz and RAM 4.00 GB.
Here cmean is the mean concentration in a specified cross section,
ci is the concentration of the ith grid cell and n denotes the total
5.7. Mass transfer and mixing between two fluids number of grid cells on that section. The mixing index MI takes a
value between 0 and 1. The closer the value of MI is to 1, the better
In this section the mass transfer and the mixing performance of the mixing quality. Usually MI > 0.95 is desirable for most of appli-
the microchannel in the case of cylindrical obstacle and hyperelas- cations (Thakur et al., 2003). Fig. 25 shows the concentration pro-
tic beam is compared to the cylindrical obstacle and rigid beam file and MI value at the outlet cross section for the rigid beam case
case. The channel inlet is divided into two identical sections. The and hyperelastic beam case.
fluid concentration in the upper and lower half sections are set It is observed that the MI = 0.963 for the hyperelastic beam case
to be 1and 0 mol=m3 , respectively. Fig. 21 shows the concentration showing an increase of 16.86% with respect to the rigid beam case

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Fig. 21. Concentration field in the case of cylindrical obstacle and rigid beam at time t ¼ 30.
574 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 22. Concentration field at different times for the case of elastic beam.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 575

Fig. 23. Variations of concentration along two horizontal lines of y ¼ 0:25 and y ¼ 0:75 at time t  ¼ 30.

for which MI = 0.824. Therefore, FVG’s are more efficient than is equal to 0.69 MPa. This vale is also much less than the yield
RVG’s in mixing of fluids and more importantly they reduce the stress of 2.4 MPa. Therefore, it can be anticipated that in the lami-
mixing length. This especially would help the engineers design nar flow regime (ReH < 1500) considered in the present work, the
more compact and smaller mixers. employed silicon rubber beam is in the safe design zone and it
would not experience failure.
5.8. Deformation of the hyperelastic beam Fig. 31 depicts the x-component of the beam’s tip displacement
for different Re numbers. It is found that as the Re number
In this section, first the deformation of the beam for ReH ¼ 1000 increases, the beam oscillates with shorter amplitude and higher
is investigated in more detail. Then the effects of Re number on the frequency. This can also be deduced from Fig. 32, which is FFT of
beam behavior is studied by comparing the results for ReH ¼ 1000, Fig. 31. These findings reflect well the flow physics because by
1100 and 1200. increasing the Re number, the frequency of vortices that are shed-
The von Mises stress along with the beam displacement for ding from the circular obstacle placed upstream of the beam
ReH ¼ 1000 is depicted in Fig. 26. It is observed that the highest increases causing the beam to oscillate faster. Table 7 provides
stress is associated with the base of the beam where it is attached more detailed information about Fig. 31.
to the bottom wall of the channel. Based on the distortion energy
theory, the failure occurs when the von Mises stress exceeds the
yield strength of the material. The yield strength of the silicone 6. Conclusions
rubber is equal to 2:4MPa ([Link]
aspx?ArticleID=920). Fig. 27 demonstrates the maximum von Heat and mass transfer enhancement using passive and active
Mises stress along the time indicating that it does not exceed vortex generators (VGs) has many scientific and engineering appli-
0.62 MPa. This would guarantee that the beam would not fail at cations. The present work employs the advantages of both the pas-
this Re number. sive and active VGs. A hyperelastic beam as a VG is placed in a
Fig. 28 depicts the displacement of the beam tip in both the x microchannel which oscillates due to the periodic vortices gener-
and y directions. After approximately t  ¼ 10, the beam oscillates ated by a cylindrical obstacle placed further upstream in the chan-
harmonically with constant frequency and amplitude. According nel. In comparison with the rigid beam, the elastic beam has the
to the relations (17) and (18) the average and amplitude of the har- following advantages.
monic oscillations of the beam in the x and y directions are
134:45  7:6lm and 43:75  4:95lm, respectively.
In the case of rigid beam a fixed vortex is generated behind the
The frequency of the harmonic oscillations of the beam is beam which has very small effect on the heat and mass transfer
obtained by applying the fast Fourier transform (FFT) to the dis- enhancement. After flow reattachment which occurs at a cer-
placement signal of the beam tip. This is found to be 2:6Hz (see tain distance downstream of the beam, the thermal boundary
Fig. 29). layer starts re-growing and heat transfer rate decreases again.
In order to investigate the effect of Re number on the oscillatory In the case of hyperelastic beam, the oscillation of the beam
motion of the hyperelastic beam, three different ReH numbers of leads to the formation and separation of periodic vortices which
1000, 1100 and 1200 are considered. Higher values of Re number generate a train of vortices downstream of the beam that may
are not considered because of the restriction of the ALE method extend towards the channel outlet. These vortices prevent the
in modeling large deformations. Fig. 30 indicates the variations thermal boundary layer from re-growing and hence increase
of the maximum von Mises stress with time for the three Re num- the heat transfer rate dramatically. The total Nu number (the
bers. It is obviously observed that as the Re number increases the sum of Nu number on the top and bottom walls) associated with
fluid flow can exert more stress on the elastic beam. The highest the hyperelastic beam case has an increase of 18.46% with
stress exerted on the beam is associated with ReH ¼ 1200, which respect to the rigid beam case.
576 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 24. Velocity vectors at different times for the case of elastic beam.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 577

The hyperelastic beam represents lower resistance against the

flow and hence a lower pressure loss along the channel occurs.
So, the Darcy friction factor in the case of hyperelastic beam
experiences a decrease of 42.33% as compared to the rigid case.

The thermal performance factor (defined as the ratio of the
change in the heat transfer rate to change in friction factor) in
the case of hyperelastic beam is 42% more than its counterpart
in the rigid beam case.

The mixing index, which is a criterion measuring the mixture
quality and concentration homogeneity increases about
16.86% with respect to the rigid beam case.

According to the distortion energy theory, the hyperelstic beam
fall in the safe design range and it would not fail in the laminar
flow regime. In addition, increase of Re number causes the beam
to oscillate with shorter amplitude but higher frequency.

Fig. 25. Time-averaged concentration profile and MI value at the outlet cross
section of the channel.

Fig. 26. von Mises stress and displacement of the hyperelastic beam.

Fig. 27. Maximum von Mises stress in the hyperelastic beam.

578 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Fig. 28. Displacement of the tip of the hyperelastic beam.

Fig. 29. FFT of the displacement signal of the beam tip.

Fig. 30. Maximum von Mises stress in the hyperelastic beam for different Re numbers.
 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580 579

Fig. 31. Horizontal component of the displacement of the tip of the hyperelastic beam for different Re numbers.

Fig. 32. FFT of the displacement signal of the beam tip for three different Re numbers.

Table 7
Effect of Re number on the characteristics of the hyperelastic beam. Ali, S., Habchi, C., Menanteau, S., Lemenand, T., Harion, J.-L., 2015. Heat transfer and
mixing enhancement by free elastic flaps oscillation. Int. J. Heat Mass Transf. 85,
ReH ¼ 1000 ReH ¼ 1100 ReH ¼ 1200 250–264.
Ali, S., Habchi, C., Menanteau, S., Lemenand, T., Harion, J.-L., 2017. Three-
Mean [lm] 134.45 140.32 146.26
dimensional numerical study of heat transfer and mixing enhancement in a
Amplitude [lm] 7.60 4.92 3.31 circular pipe using self-sustained oscillating flexible vorticity generators. Chem.
Frequency [lm] 2.60 2.86 3.58 Eng. Sci. 162, 152–174.
Arruda, E.M., Boyce, M.C., 1993. A three-dimensional constitutive model for the
large stretch behavior of rubber elastic materials. J. Mech. Phys. Solids 41 (2),
Declaration of Competing Interest
389–412.
Attard, M.M., 2003. Finite strain––isotropic hyperelasticity. Int. J. Solids Struct. 40
None (17), 4353–4378.
Babar, H., Ali, H.M., 2019. Towards hybrid nanofluids: preparation, thermophysical
properties, applications, and challenges. J. Mol. Liq.
Appendix A. Supplementary material Becker, S., Lienhart, H., Durst, F., 2002. Flow around three-dimensional obstacles in
boundary layers. J. Wind Eng. Ind. Aerodyn. 90 (4–5), 265–279.
Bejan, A., Kraus, A.D., 2003. Heat transfer handbook. John Wiley & Sons.
Supplementary data to this article can be found online at Bower, A.F., 2009. Applied Mechanics of Solids. CRC Press.
[Link] Bruus, H., 2008. Theoretical Microfluidics. Oxford University Press, Oxford.
Celik, I.B., Ghia, U., Roache, P.J., 2008. Procedure for estimation and reporting of
uncertainty due to discretization in CFD applications. J. Fluids Eng.-Trans. ASME
References 130 (7).
Depaiwa, N., Chompookham, T., Promvonge, P., 2010. Thermal enhancement in a
solar air heater channel using rectangular winglet vortex generators. In: Energy
Ahmed, H., Mohammed, H., Yusoff, M., 2012. An overview on heat transfer
and Sustainable Development: Issues and Strategies (ESD), 2010 Proceedings of
augmentation using vortex generators and nanofluids: approaches and
the International Conference on, IEEE, pp. 1–7.
applications. Renew. Sustain. Energy Rev. 16 (8), 5951–5993.
580 A. Dadvand et al. / Chemical Engineering Science 207 (2019) 556–580

Donea, J., Huerta, A., Ponthot, J.-P., Rodriguez-Ferran, A., 2004. Encyclopedia of Mirzaee, H., Dadvand, A., Mirzaee, I., Shabani, R., 2012. Heat transfer enhancement
Computational Mechanics, Vol. 1: Fundamentals., Chapter 14: Arbitrary. in microchannels using an elastic vortex generator. J. Enhanced Heat Transfer
Lagrangian-Eulerian Methods. Wiley & Sons. 19 (3).
Fiebig, M., 1998. Vortices, generators and heat transfer. Chem. Eng. Res. Des. 76 (2), Mohebbi, R., Rashidi, M., Izadi, M., Sidik, N.A.C., Xian, H.W., 2018. Forced convection
108–123. of nanofluids in an extended surfaces channel using lattice Boltzmann method.
Fu, W.-S., Tong, B.-H., 2002. Numerical investigation of heat transfer from a heated Int. J. Heat Mass Transf. 117, 1291–1303.
oscillating cylinder in a cross flow. Int. J. Heat Mass Transf. 45 (14), 3033–3043. Morand, H.J.-P., Ohayon, R., 1995. Fluid Structure Interaction-Applied Numerical
Gao, Y., 1997. Large deformation field near a crack tip in rubber-like material. Theor. Methods. Wiley.
Appl. Fract. Mech. 26 (3), 155–162. Park, S.G., Chang, C.B., Kim, B., Sung, H.J., 2017. Simulation of fluid-flexible body
Gent, A., 1996. A new constitutive relation for rubber. Rubber Chem. Technol. 69 (1), interaction with heat transfer. Int. J. Heat Mass Transf. 110, 20–33.
59–61. Peskin, C.S., 2002. The immersed boundary method. Acta Numer. 11, 479–517.
Habchi, C., Lemenand, T., Della Valle, D., Peerhossaini, H., 2010. Turbulent mixing Rehman, T.-U., Ali, H.M., Janjua, M.M., Sajjad, U., Yan, W.-M., 2019. A critical review
and residence time distribution in novel multifunctional heat exchangers– on heat transfer augmentation of phase change materials embedded with
reactors. Chem. Eng. Process.: Process Intensification 49 (10), 1066–1075. porous materials/foams. Int. J. Heat Mass Transf. 135, 649–673.
Habchi, C., Lemenand, T., Valle, D.D., Peerhossaini, H., 2010. Turbulence behavior of Rivlin, R.S., Saunders, D., 1951. Large elastic deformations of isotropic materials VII.
artificially generated vorticity. J. Turbul. (11), N36 Experiments on the deformation of rubber. Phil. Trans. R. Soc. Lond. A 243 (865),
Holzapfel, A.G., 2000. Nonlinear Solid Mechanics II. 251–288.
[Link] in. Sajid, M.U., Ali, H.M., 2019. Recent advances in application of nanofluids in heat
Kakaç, S., Shah, R.K., Aung, W., 1987. Handbook of single-phase convective heat transfer devices: a critical review. Renew. Sustain. Energy Rev. 103, 556–592.
transfer. Sharipov, F., 2011. Data on the velocity slip and temperature jump on a gas-solid
Karniadakis, G., Beskok, A., Aluru, N., 2006. Microflows and Nanoflows: interface. J. Phys. Chem. Ref. Data 40, (2) 023101.
Fundamentals and Simulation. Springer Science & Business Media. Shoele, K., Mittal, R., 2014. Computational study of flow-induced vibration of a reed
Khatavkar, V.V., Anderson, P.D., den Toonder, J.M., Meijer, H.E., 2007. Active in a channel and effect on convective heat transfer. Phys. Fluids 26, (12) 127103.
micromixer based on artificial cilia. Phys. Fluids 19, (8) 083605. Soti, A.K., Bhardwaj, R., Sheridan, J., 2015. Flow-induced deformation of a flexible
Khattak, Z., Ali, H.M., 2019. Air cooled heat sink geometries subjected to forced flow: thin structure as manifestation of heat transfer enhancement. Int. J. Heat Mass
A critical review. Int. J. Heat Mass Transf. 130, 141–161. Transf. 84, 1070–1081.
Lambert, R.A., Rangel, R.H., 2010. The role of elastic flap deformation on fluid mixing Storåkers, B., 1986. On material representation and constitutive branching in finite
in a microchannel. Phys. Fluids 22, (5) 052003. compressible elasticity. J. Mech. Phys. Solids 34 (2), 125–145.
Lee, J.B., Park, S.G., Kim, B., Ryu, J., Sung, H.J., 2017. Heat transfer enhancement by Thakur, R., Vial, C., Nigam, K., Nauman, E., Djelveh, G., 2003. Static mixers in the
flexible flags clamped vertically in a Poiseuille channel flow. Int. J. Heat Mass process industries—a review. Chem. Eng. Res. Des. 81 (7), 787–826.
Transf. 107, 391–402. Turek, S., Hron, J., 2006. Proposal for numerical benchmarking of fluid-structure
Lee, J.B., Park, S.G., Sung, H.J., 2018. Heat transfer enhancement by asymmetrically interaction between an elastic object and laminar incompressible flow. pp. 371–
clamped flexible flags in a channel flow. Int. J. Heat Mass Transf. 116, 1003– 385.
1015. Van Loon, R., Anderson, P., Van de Vosse, F., Sherwin, S., 2007. Comparison of various
Li, Z., Chen, Y., Xu, X., Li, K., Ke, Z., Zhou, K., Chen, H.-H., Huang, G., Chen, C.-l., Chen, fluid–structure interaction methods for deformable bodies. Comput. Struct. 85
C.-H., 2017. Air-side heat transfer enhancement with a novel self-agitator. In: (11–14), 833–843.
ASME 2017 Heat Transfer Summer Conference, American Society of Mechanical Webb, R.L., Kim, N., 2005. Enhanced Heat Transfer. Taylor and Francis, NY.
Engineers, pp. V002T010A008–V002T010A008. Yang, Y.-T., Chen, C.-H., 2008. Numerical simulation of turbulent fluid flow and heat
Ma, J., Huang, Y.P., Huang, J., Wang, Y.L., Wang, Q.W., 2010. Experimental transfer characteristics of heated blocks in the channel with an oscillating
investigations on single-phase heat transfer enhancement with longitudinal cylinder. Int. J. Heat Mass Transf. 51 (7–8), 1603–1612.
vortices in narrow rectangular channel. Nucl. Eng. Des. 240 (1), 92–102. Yu, Y., Liu, Y., Chen, Y., 2018. Vortex dynamics and heat transfer behind self-
Meis, M., Varas, F., Velázquez, A., Vega, J., 2010. Heat transfer enhancement in oscillating inverted flags of various lengths in channel flow. Phys. Fluids 30, (4)
micro-channels caused by vortex promoters. Int. J. Heat Mass Transf. 53 (1–3), 045104.
29–40. Zdravkovich, M., 1997. Flow around circular cylinders. Fundamentals, vol. 1.