Heat Transfer in Non-Newtonian Flows by a Hybrid Immersed Boundary–Lattice Boltzmann and Finite Difference Method

Abstract

A hybrid immersed boundary–lattice Boltzmann and finite difference method for fluid–structure interaction and heat transfer in non-Newtonian flow is presented. The present numerical method includes four parts: fluid solver, heat transfer solver, structural solver, and immersed boundary method for fluid–structure interaction and heat transfer. Specifically, the multi-relaxation time lattice Boltzmann method is adopted for the dynamics of non-Newtonian flow, with a geometry-adaptive technique to enhance the computational efficiency and immersed boundary method to achieve no-slip boundary conditions. The heat transfer equation is spatially discretized by a second-order up-wind scheme for the convection term, a central difference scheme for the diffusion term, and a second-order difference scheme for the temporal term. The structural dynamics is numerically solved using a finite difference method. The major contribution of this work is the integration of spatial adaptivity, thermal finite difference method, and fluid flow immersed boundary-lattice Boltzmann method. Several benchmark problems including the developing flow of non-Newtonian fluid in a channel, non-Newtonian fluid flow and heat transfer around a stationary cylinder and flow around a stationary cylinder with a detached filament are used to validate the present method and developed solver. The good agreements achieved by the present method with the published data show that the present extension is an efficient way for fluid–structure interaction and heat transfer involving non-Newtonian fluid. The heat transfer around an oscillating cylinder in non-Newtonian fluid flow at Reynolds number of 100 is also numerically studied using the present solver, considering the effects of the oscillating frequency and amplitude. The results may be used to expand the currently limited database of fluid–structure interaction and heat transfer benchmark studies.

1. Introduction

Experiments are expensive and challenging with increasing complexity. As an alternative research method, numerical simulations have attracted much attention in the past twenty years [1]. As an alternative method of computational fluid dynamics (CFD), the lattice Boltzmann method (LBM) has received remarkable attention since its origin [2,3]. After the proposal of the Bhatnagar–Gross–Krook (BGK) collision operator, the efficiency and flexibility of the LBM have been enhanced significantly. Consequently, the BGK-based LBM has been successfully applied to various fluid dynamics problems, such as multiphase flow [4,5,6,7], complex flow in porous media [8], and non-Newtonian flow [9,10]. To address the numerical instability in single-relaxation-time LBM due to low viscosity, entropic [11,12], two-relaxation-time (TRT), and multi-relaxation-time (MRT) LBM methods have been developed [4]. The major advantages of the LBM are its simplicity, easy implementation, explicit calculations, and intrinsic parallel nature. However, the uniform mesh used in the standard LBM limits its applications in engineering [13]. Many efforts have been made to remove this drawback, including differential type and interpolation type LBM. Grid refinement is another way to retain the simplicity of the LBM and enhance computational efficiency. The mesh is locally refined where high resolution is needed. The multi-block technique proposed by Yu et al. [14] divides the fluid field into several blocks with different mesh sizes. As the mesh is initialized before the computation in this method, the mesh refinement is not involved with the fluid field. Wu and Shu proposed a solution-adaptive refinement scheme, where the refinement procedure is decided by the fluid field [13]. By using the refinement technique, the computational efficiency can be highly enhanced. Fluid–structure interaction (FSI) problems extensively exist in engineering. Thanks to the representation of the forcing term introduced by Guo et al. [15], the involved external or internal force can be considered approximately, which also enables the successful combination of the LBM and the immersed boundary method (IBM) (know as IB-LBM) [16,17,18]. The IB-LBM proved to be an efficient method for FSI and has been applied to various problems successfully [10,19,20,21,22,23,24,25,26].

Non-Newtonian fluid is extensively involved in many industrial substances, i.e., blood flow, polymer, and melt solutions. The high viscosity levels in conjunction with the complexity in the behaviour of the stress terms of momentum equations for non-Newtonian fluids generally give rise to the limitation of fluid processing under laminar flow conditions [27,28]. Due to the potential of rheological characteristics (shear-thinning, shear-thickening, etc.) in the applications of engineering, remarkable efforts have been made to develop efficient numerical methods to handle non-Newtonian flows and heat transfer over the last few decades [9,27,29]. As the LBM has proven to be an efficient method for CFD [3], it is thought to be a good idea to employ the LBM to handle non-Newtonian flow and heat transfer. In terms of the simulation of thermal fluid systems, most of the previous thermal LBM models adopt the multi-speed approach [30,31,32], in which additional particle speeds are needed to obtain the energy equation at the macroscopic level [33]. The additional particle speeds decrease the computational efficiency; the stability is also limited. In order to keep the simplicity property of the LBM, Peng et al. proposed a simplified thermal LBM, neglecting the compression work done by the pressure and viscous heat dissipation [34]. In this study, an alternative method to solve the heat transfer in non-Newtonian fluid flow is presented, conserving the simplicity and efficiency of the LBM. The fluid dynamics and heat transfer are solved by the LBM and finite difference method, respectively, with multi-relaxation time and a geometry-adaptive mesh technique to achieve good stability and computational efficiency. In addition to the thermal process, FSI is also achieved using an penalty immersed boundary method. In conclusion, the present method can handle problems involving FSI with complex boundaries and heat transfer in non-Newtonian fluid flow.

In this paper, a hybrid immersed boundary–lattice Boltzmann and finite difference method for FSI and heat transfer in non-Newtonian flow is presented, with the geometry-adaptive Cartesian mesh technique to enhance the computational efficiency. In this contribution, we present for the first time an integration of spatial adaptivity, thermal finite difference method, and fluid flow IB-LBM. The organization of the present paper is as follows. The numerical method is presented in Section 2. Validations of the developing flow of power-law fluid in a channel, power-law fluid flow and heat transfer around a stationary cylinder are presented in Section 3. Section 4 presents the flow around a stationary cylinder with a detached filament. The numerical simulation of heat transfer around an oscillating cylinder is presented in Section 5. Final conclusions are given in Section 6.

2. Numerical Method

In this paper, the two-dimensional non-Newtonian flow and heat transfer involving complex boundaries are considered, where the LBM is adopted for the fluid solver, the finite difference method is used to solve the heat transfer equation and structural dynamics, and the complex no-slip boundaries are achieved by an IB method. The details of the current method will be given in this section.

2.1. Fluid Solver

In the MRT-based IB-LBM, the evolution equation of the velocity distribution function $g_i$ along the i-th direction at position $\mathit{x}$ with BGK approximation is expressed as [35]

$$g_i(\mathit{x}+{\mathit{e}}_i\Delta t,t+\Delta t)=g_i(\mathit{x},t)-{\mathsf{\Omega }}_i(\mathit{x},t)+\Delta t{G}_i,i=0,1,\cdots ,8$$

(1)

where $\Delta t$ is the time step. The collision operator ${\mathsf{\Omega }}_i$ and $G_i$, which represent the body force effects on the distribution function, are defined as

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_i& =& -{\left(M^{-1}SM\right)}_{ij}[g_j(\mathit{x},t)-g_j^{eq}(\mathit{x},t)]\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill G_i& =& {[M^{-1}(I-S/2)M]}_{ij}{F}_j\hfill \end{array}$$

(3)

where M is a $9\times 9$ transform matrix for the D2Q9 model (two dimensional with 9 directions) used here, and S is a non-negative diagonal matrix related to the fluid viscosity. The details for the determination of S can be found in [36]. The lattice speed ${\mathit{e}}_i$ is defined as

$${\mathit{e}}_i=\left\{\begin{array}{cc}(0,0),\hfill & i=0\hfill \\ (\cos (\pi (i-1)/2),\sin (\pi (i-1)/2))\frac{\Delta x}{\Delta t},\hfill & i=1,2,3,4\hfill \\ \sqrt{2}(\cos (\pi (i-9/2)/2),\sin (\pi (i-9/2)/2))\frac{\Delta x}{\Delta t},\hfill & i=5,6,7,8\hfill \end{array}\right.$$

(4)

where $\Delta x$ is the lattice spacing. Then, the macro density and momentum are given as follows,

$$\rho =\sum _{i=1}^8{g}_i,\rho \mathit{u}=\sum _{i=1}^8{g}_i{\mathit{e}}_i+\frac{1}{2}\mathit{f}\Delta t.$$

(5)

The local equilibrium distribution function $g_i^{eq}$ and the force term $G_i$ are calculated by

$$g_i^{eq}={\omega }_i\rho [1+\frac{{\mathit{e}}_i·\mathit{u}}{c_s^2}+\frac{\mathit{u}\mathit{u}:({\mathit{e}}_i{\mathit{e}}_i-c_s^2\mathit{I})}{c_s^4}]$$

(6)

$$G_i={\omega }_i[\frac{{\mathit{e}}_i-\mathit{u}}{c_s^2}+\frac{{\mathit{e}}_i·\mathit{u}}{c_s^4}{\mathit{e}}_i]·\mathit{f}$$

(7)

where the weights ${\omega }_i$ are given by ${\omega }_0=4/9$, ${\omega }_i=1/9$ for $i=1,2,3,4$, and ${\omega }_i=1/36$ for $i=5,6,7,8$. The sound speed $c_s=\Delta x/(\sqrt{3}\Delta t)$, and $\mathit{f}$ is the force acting on the fluid. The relaxation time is related to the kinematic viscosity $\nu$ in the Navier–Stokes equations in terms of $\nu =(\tau -0.5)c_s^2\Delta t$. The rheological equation of state for power-law non-Newtonian fluid considered in this paper is defined by [27]

$$\nu =m{\left(\frac{I_2}{2}\right)}^{(n-1)/2}$$

(8)

where m is the power-law consistency index, n is the power-law fluid behaviour index, and $I_2$ is the second invariant of the rate of strain tensor. Additionally, the multi-block technique developed by Yu et al. [14] is combined with the geometry-adaptive method to enhance the computational efficiency.

2.2. Heat Transfer Solver

The governing equation for the heat transfer in the fluid is given as [37]

$$\rho c_p(\frac{\partial T}{\partial t}+\mathit{u}·\nabla T)=\kappa {\nabla }^{2}T+q$$

(9)

where T is the temperature, $c_p$ and $\kappa$ are the specific heat and thermal conductivity coefficient of the fluid, respectively. q is the source term.

The second-order upwind scheme is used to discretize the convection term in Equation (9) explicitly, expressed as

$$\begin{array}{cc}\hfill \frac{\partial T}{\partial x}& =\frac{1}{2\Delta x}(3T_i-4T_{i-1}+T_{i-2}),u>0\end{array}$$

(10)

$$\begin{array}{cc}\hfill \frac{\partial T}{\partial x}& =\frac{1}{2\Delta x}(-3T_i+4T_{i+1}-T_{i+2}),u<0.\end{array}$$

(11)

A second-order internal central difference scheme is adopted for the spatial discretization on the boundaries. The subscripts indicate the position.

The diffusion term in Equation (9) is discretized using a second-order explicit scheme, expressed as

$$\frac{{\partial }^{2}T}{\partial x^2}=\frac{1}{\Delta x^2}(T_{i-1}-2T_i+T_{i+1}),$$

(12)

and the internal difference scheme is also used for the boundary nodes. All the discretizations of T in the y direction is similar to that in the x direction. In addition, a second-order explicit scheme is used for the temporal discretization. Therefore, the heat transfer solver has a second-order accuracy spatially and temporally, which is consistent with the accuracy of the LBM. The heat transfer equation is also solved on the geometry-adaptive Cartesian mesh that used for the fluid dynamics.

2.3. Structural Solver

The geometrically nonlinear motion for the filament considered here is governed by [38,39,40,41,42,43]

$$m_s\frac{{\partial }^2\mathit{X}}{\partial t^2}-\frac{\partial }{\partial s}[K_s\left(\right|\frac{\partial \mathit{X}}{\partial s}|-1)\frac{\partial \mathit{X}}{\partial s}]+K_b\frac{{\partial }^4\mathit{X}}{\partial s^4}={\mathit{F}}_f$$

(13)

where $m_s$ is the linear density of the filament, s is the Lagrangian coordinate along the length, $K_s$ is the stretching coefficient, $T\left(s\right)$ is the tensile stress defined by $T\left(s\right)=K_s\left(\right|\frac{\partial \mathit{X}}{\partial s}|-1)$, $\mathit{X}$ is the position vector of a point on the filament, $K_b$ is the bending rigidity, and ${\mathit{F}}_f$ is the hydrodynamic stress exerted by the fluid.

The filament is discretized by $N_f$ initially equally spaced nodal points, and the position of the m-th node at time level n is denoted by ${\mathit{X}}_m^n$. The tensile force at the m-th node is calculated by a finite-difference scheme expressed as

$$\frac{\partial }{\partial s}[T\left(s\right)\frac{\partial \mathit{X}}{\partial s}]=\frac{1}{\Delta s}(T_{m+1/2}{\psi }_{m+1/2}-T_{m-1/2}{\psi }_{m-1/2})$$

(14)

where $\Delta s$ is the grid spacing, the tension T and tangent vector, $\psi =\partial \mathit{X}/\partial s$, at the segment center, $m+1/2$, are both computed using a second-order central difference scheme. The bending force are computed using a finite difference scheme expressed as

$$\frac{{\partial }^4{\mathit{X}}_l}{\partial s^4}=\frac{1}{\Delta s^4}\sum _{m=2}^{N_f-1}({\mathit{X}}_{m+1}+{\mathit{X}}_{m-1}-2{\mathit{X}}_m)(2{\delta }_{m,l}-{\delta }_{m+1,l}-{\delta }_{m-1,l})$$

(15)

where ${\delta }_{k,l}$ is a Kronecker symbol, defined by ${\delta }_{k,l}=0$ when $k\ne l$ and by ${\delta }_{k,l}=1$ when $k=l$. The details of the schemes can be found in [44].

2.4. The IB Method for Fluid–Structure Interaction and Heat Transfer

The pIB method developed by Kim and Peskin [45] is used to handle the no-slip boundaries between the rigid body and the fluid. The interaction force between the fluid and the structure can be determined by the feedback law [45,46]:

$$\mathit{F}=\alpha {\int }_0^t({\mathit{U}}_{ib}-\mathit{U})dt+\beta ({\mathit{U}}_{ib}-\mathit{U})$$

(16)

where ${\mathit{U}}_{ib}$ is the boundary velocity obtained by interpolation at the IB, $\mathit{U}$ is the structure velocity, and $\alpha$ and $\beta$ are large positive free constants. The force acting on the Lagrange structure from the ambient fluid can be taken as a concentrated force acting on the corresponding nodes; thus, it can be added to the body force $\mathit{f}$ in Equation (7). Compared to the sharp-interface method [47,48,49,50], the pIB method is particularly suitable here, since in this approach all the grid points within the computational domain are treated with a unified equation.

The transformation between the Euler and Lagrange variables can be realized by the Dirac delta function. The interpolation of velocity and the spreading of the Lagrange force to the adjacent grid points are expressed as

$${\mathit{U}}_{ib}(s,t)={\int }_V\mathit{u}(x,t){\delta }_h(\mathit{X}(s,t)-\mathit{x})d\mathit{x}$$

(17)

$$\mathit{f}(\mathit{x},t)=-{\int }_{\mathrm{\Gamma }}\mathit{F}(s,t){\delta }_h(\mathit{X}(s,t)-\mathit{x})ds$$

(18)

where $\mathit{u}$ is the fluid velocity, $\mathit{X}$ is the coordinates of structural nodes, $\mathit{x}$ is the coordinates of fluid, s is the arc coordinate, and V is the fluid domain, and $\mathrm{\Gamma }$ is the structure domain.

The position of the filament is updated explicitly by

$${\mathit{X}}^{n+1}={\mathit{X}}^n+{\mathit{U}}_{ib}\Delta t$$

(19)

where ${\mathit{U}}_{ib}$ is calculated using (17). The internal force density of the filament can be calculated using (14) and (15). The details of the computation of the inertial force defined in (13) can be found in [17]. The ${\mathit{F}}_f$ can then be calculated by Equation (13). The force $\mathit{F}$ acting on the fluid in Equation (18) is the reaction force of ${\mathit{F}}_f$ defined in Equation (13). Equation (18) is used to distribute the force to the ambient fluid.

The smooth function ${\delta }_h$ is used to approximate the Dirac delta function

$${\delta }_h(x,y)=\frac{1}{h^2}\lambda (\frac{x}{h})\lambda (\frac{y}{h}).$$

(20)

In this paper, the four-point delta function introduced by Peskin [51] is used

$$\lambda \left(r\right)=\left\{\begin{array}{cc}\frac{1}{8}(3-2|r|+\sqrt{1+4\left|r\right|-{4\left|r\right|}^2}),\hfill & 0\le \left|r\right|<1\hfill \\ \frac{1}{8}(5-2|r|-\sqrt{-7+12\left|r\right|-{4\left|r\right|}^2}),\hfill & 1\le \left|r\right|<2\hfill \\ 0,\hfill & 2\le \left|r\right|.\hfill \end{array}\right.$$

(21)

The heat transferred from the immersed boundary to the fluid can be written as [52]

$$q={\int }_{\mathrm{\Gamma }}Q(s,t){\delta }_h(\mathit{X}(s,t)-\mathit{x})ds$$

(22)

where Q is the virtual boundary heat flux. When the temperature boundary is used, then Q can be calculated using the penalty immersed-boundary method, expressed as

$$Q={\alpha }_T(T_{ib}-T_0)$$

(23)

where ${\alpha }_T$ is a large factor, $T_0$ is the boundary condition on the immersed boundary, and $T_{ib}$ is the temperature of the virtual boundary, interpolated by

$$T_{ib}(s,t)={\int }_{V}T(x,t){\delta }_h(\mathit{X}(s,t)-\mathit{x})d\mathit{x}.$$

(24)

3. Validations

3.1. The Developing Flow of Non-Newtonian Power-Law Fluid in a Channel

The developing flow through a confined channel has been extensively used as a benchmark validation for non-Newtonian fluid. In this section, the two-dimensional steady laminar developing flow is considered to validate the in-house code. In the numerical simulation, the flow of power-law fluid develops with a uniform inlet velocity ($U_0$) through a rectangular channel of h in height and L in length, as shown in Figure 1. The macro density and velocities on the boundaries are set as $\mathit{u}=(U_0,0)$ at $x=0$, $\mathit{u}=(0,0)$ at $y=0$ and h for the channel wall, and $\rho ={\rho }_0$ at $x=L$. The non-equilibrium extrapolation presented by Guo et al. [53] is adopted to calculate the velocity distribution functions on the boundaries.

The numerical simulation is performed in a computational domain with the size of $15h\times 1h$ and the uniform mesh size is $h/50$. The dimensionless time step $\Delta t^{*}=\Delta t{U}_0/h=0.001$. The Reynolds number (Re) that controls this problem is defined as Re = ${\rho }_0{h}^n{U}_0^{2-n}/m$. The Reynolds number of 10 is used in the current steady simulation, and five power-law indexes n ranging from $0.3$ to $1.5$ (including shear-thinning $n<1$, shear-thickening $n>1$, and Newtonian fluid $n=1$) are considered. The analytical solution of the velocity profiles of the fully developed flow can be expressed as [27]

$$\frac{u_f(y,n)}{U_0}=\frac{2n+1}{n+1}(1-|1-\frac{2y}{h}{|}^{(n+1)/n}).$$

(25)

The fully developed velocity profiles from present numerical simulations and that predicted by Equation (25) are presented in Figure 2. According to the comparisons in Figure 2, excellent agreements (discrepancies $<2\%$) are achieved by the present numerical method. The velocity peaks presented in

Table 1. The velocity peaks at Re = 10.

Sources	$\mathit{n}=\mathbf{0.3}$	$\mathit{n}=\mathbf{0.5}$	$\mathit{n}=\mathbf{0.7}$	$\mathit{n}=\mathbf{1.0}$	$\mathit{n}=\mathbf{1.5}$
Present	1.242	1.333	1.411	1.501	1.591
Analytical	1.231	1.333	1.412	1.500	1.600

also show the good accuracy of the current method.

3.2. Non-Newtonian Power-Law Fluid Flow and Heat Transfer around a Stationary Cylinder

Non-Newtonian power-law fluid flow and heat transfer from a stationary cylinder is considered to validate the accuracy of the present numerical method in handling FSI and heat transfer. The two-dimensional power-law fluid inlets with a uniform flow $U_0$ over a stationary circular cylinder (of diameter D), with the computational domain extending from $(-15D,-5D)$ to $(25D,5D)$, as shown in Figure 3.

In order to improve the computational efficiency and validate the geometry-adaptive technique utilized in present in-house code, four blocks are used in the simulations, as shown in Figure 4, and the most refined mesh size of the fluid domain is $D/80$. The dimensionless time step $\Delta t^{*}=\Delta t{U}_0/D=6.25\times {10}^{-4}$.

The non-dimensional parameters that control this problem include the Reynolds number (Re) and the Prandtl number (Pr), which are defined as

$$\mathrm{Re}=\frac{{\rho }_0{D}^n{U}_0^{2-n}}{m},\mathrm{Pr}=\frac{c_{p}m{(U_0/D)}^{n-1}}{\kappa }.$$

(26)

The drag coefficient, lift coefficient, and Strouhal number are defined as

$$C_D=-\frac{\sum F_x}{0.5\rho U_0^{2}D},C_L=\frac{\sum F_y}{0.5\rho U_0^{2}D},\mathrm{St}=\frac{D}{U_{0}T}$$

(27)

where T is the vortex shedding period, and $F_x$ and $F_y$ are the horizontal and vertical components of $\mathit{F}$ acting on the cylinder calculated by Equation (16).

Firstly, we consider the unsteady flow of non-Newtonian power-law flow over a cylinder at Re $=100$ (at which the flow patterns show von-Kármán periodic vortex shedding behind the cylinder) to validate the geometry-adaptive fluid solver in terms of handling FSI problem. In order to achieve the early von-Kármán vortex shedding, a vertical velocity perturbation is introduced into the initial inlet boundary. Four different indexes $n=0.6,1.0$, $1.4$, and $1.8$ are calculated.

Table 2. A uniform flow over a stationary cylinder at Re = 100. St: Strouhal number; $C_{D,m}$: drag coefficient; $C_L$: lift coefficient.

n	Sources	St	${\mathit{C}}_{\mathit{D},\mathit{m}}$	${\mathit{C}}_{\mathit{L}}$
0.6	Present	0.182	1.258	0.375
	Patnana et al. [54]	0.180	1.180	–
	Tian et al. [27]	0.188	1.179	0.367
1.0	Present	0.164	1.415	0.349
	Patnana et al. [54]	0.164	1.341	0.325
	Tian et al. [27]	0.160	1.430	0.360
	Wang et al. [55]	0.161	1.450	0.310
	Xu and Wang [56]	0.171	1.423	0.340
	Tian et al. [17]	0.166	1.43	–
1.4	Present	0.159	1.546	0.345
	Patnana et al. [54]	0.150	1.497	–
	Tian et al. [27]	0.161	1.523	0.356
1.8	Present result	0.152	1.661	0.327
	Patnana et al. [54]	0.139	1.630	–
	Tian et al. [27]	0.155	1.657	0.356

shows the present results of St, $C_{D,m}$ (mean drag coefficient) and $C_L$ as well as results from the literature. It shows that the present results agree well with these other values. A snapshot of the vortex contour for Re = 100 and $n=1.4$ after periodic shedding are presented in Figure 5.

After the validation of the fluid solver in computing FSI problems involving non-Newtonian flow, we further consider the steady flow and heat transfer around a cylinder. The average Nusselt number ($Nu$) used for quantitative comparison with the available data from the literature are defined as

$$Nu=\frac{1}{2\pi }{\int }_0^{2\pi }Nu\left(\theta \right)\mathrm{d}\theta$$

(28)

where $Nu\left(\theta \right)$ is the local Nusselt number on the surface of the cylinder, it can be evaluated using the temperature field according to

$$Nu\left(\theta \right)=-\frac{\partial T}{\partial {\mathit{n}}_s}$$

(29)

where ${\mathit{n}}_s$ is the unit vector normal to the surface of the cylinder.

The average Nusselt number from the current computation along with the available data from the literature of two Reynolds number (Re = 10 and 40) and three indexes ($n=0.8$, $1.0$, and $1.4$) are presented in

Table 3. Averaged Nusselt number ($Nu$) for forced convection heat transfer from a stationary cylinder to power-law fluids at Pr = 1.0.

Re	Sources	$\mathit{n}=\mathbf{0.8}$	$\mathit{n}=\mathbf{1.0}$	$\mathit{n}=\mathbf{1.4}$
10	Present	2.089	2.038	1.973
	Bharti et al. [57]	2.123	2.060	1.973
	Tian et al. [27]	2.208	2.150	2.075
	Soares et al. [58]	2.116	2.058	1.973
40	Present	3.714	3.588	3.401
	Bharti et al. [57]	3.830	3.653	3.400
	Tian et al. [27]	3.923	3.769	3.554
	Soares et al. [58]	3.736	3.570	3.325

. It is found that the present results are in good agreement with other published data. With the increase of n from $0.8$ to $1.4$, the average Nusselt number decreases slightly at both Reynolds numbers. Comparisons between the Reynolds numbers of 10 and 40 show that the average Nusselt number changes more remarkably with larger Reynolds numbers. It is because of the enhancement of the heat convection effects under the larger Reynolds number. The snapshots of the temperature field from current simulations are presented in Figure 6. As is readily seen from Figure 6, the length of the vortex attached behind the cylinder stretches significantly with n increasing from 0.8 to 1.4, which also qualitatively agrees with the analysis in [54].

4. Heat Transfer around a Stationary Cylinder with a Detached Filament

A stationary cylinder with a detached flexible filament is considered in this section, as shown in Figure 7, where D is the diameter of the fixed cylinder, L is the length of the filament, and l is the distance from the origin of the cylinder to the fixed end of the filament. This model has been extensively used to study fish behaviors in a street wake [17,59,60,61]. Simulations are performed at $l/D=3.0$, $m_s=0$, and Re = 100. The non-dimensional stretching coefficient defined by $K_s^{*}=K_s/\left(\rho U^2\right)$ is 1000 (the filament is approximately inextensible), and the non-dimensional bending rigidity defined by $K_b^{*}=K_b/\left(\rho U^2{D}^3\right)$ is ${10}^{-4}$. The mesh size around the cylinder and filament is D/160. The dimensionless time step $\Delta t^{*}=\Delta t{U}_0/D=6.25\times {10}^{-4}$. Three power-law indices $0.6$, $0.8$, and $1.0$, and two filament lengths $2.5D$ and $4.0D$ are considered.

The averaged drag coefficient ($C_d$), the Strouhal number ($St$), the averaged Nusselt number ($Nu$) of the cylinder, and the vertical flapping amplitude ($A/D$) of the trailing end of the filament are presented in

Table 4. The averaged drag coefficient ($C_d$), the Strouhal number ($St$), the averaged Nusselt number ($Nu$) of the cylinder, and the vertical flapping amplitude $A/D$ of the trailing end of the filament.

Sources	n	$\mathit{L}/\mathit{D}$	${\mathit{C}}_{\mathit{D}}$	$\mathit{St}$	$\mathit{Nu}$	$\mathit{A}/\mathit{D}$
Present	0.6	2.5	1.25	0.175	5.881	1.10
		4.0	1.24	0.175	5.872	1.28
Present	1.0	2.5	1.38	0.156	5.407	1.11
Sui et al. [60]		2.5	1.41	0.156	–	1.14
Present		4.0	1.37	0.152	5.401	1.30
Tian et al. [17]		4.0	1.39	0.153	–	1.34
Present	1.4	2.5	1.51	0.143	5.077	1.10
		4.0	1.48	0.141	5.069	1.30

. The good agreements of the drag coefficient, Strouhal number, and flapping amplitude between the present results and the data from [17,60] show the excellent accuracy of the present solver. The longer filament generates a larger flapping amplitude, but the filament length does not significantly affect other fluid and temperature characteristics, as shown in Table 4. The results from non-Newtonian fluid flow show that the shear-thickening fluid generates larger drag comparing with shear-thinning fluid, but $Nu$ presents an opposite tendency, as analysed in Section 3.2. Additionally, the current results show that the power index does not have significant influence on the flapping amplitude. The instantaneous vortex and temperature contours are presented in Figure 8. The non-Newtonian results presented here can be used to extend the limited database involving FSI.

5. Heat Transfer around an Oscillating Cylinder in Non-Newtonian Fluid Flow

5.1. Physical Problem

In this section, the heat transfer around an oscillating cylinder in non-Newtonian fluid flow is considered, which can be taken as validations for heat transfer involving moving boundaries in the future. In this problem, the center of the cylinder undergoes a translational motion governed by

$$H=\frac{H_0}{2}\sin \left(2\pi ft\right)$$

(30)

where $H_0$ and f are the oscillating amplitude and frequency, respectively. The non-dimensional parameters governing this problem including oscillating frequency and amplitude can be written as

$$f^{*}=fD/\left(\mathrm{St}{U}_0\right),H^{*}=H_0/D.$$

(31)

The Reynolds number and Prantdl number are defined as those in Equation (26). In the current simulation, the influences of the oscillating amplitude, frequency, and power index are considered.

5.2. Effects of the Oscillating Amplitude

In the current simulations, the cylinder oscillates in the frequency of the vortex shedding, so the value of the non-dimensional oscillating frequency ($f^{*}$) is 1.0. The Reynolds number and the Prantdl number defined in Equation (26) are 100 and 1.0, respectively. Three indexes ($n=0.6,1.0$, and 1.4) are numerically computed to include shear-thinning, Newtonian, and shear-thickening fluid. The Strouhal numbers defined in Equation (27) for these cases have already been given in Table 2 above. Three non-dimensional oscillating amplitudes ($H^{*}$) $0.25$, $0.5$, and $1.0$ are considered to identify the effects of the amplitude. The time-averaged Nusselt numbers and their amplitudes ($\Delta Nu$) from the present computations are presented in

Table 5. Time-averaged Nusselt number ($Nu$) and its amplitude ($\Delta Nu$) for forced convection heat transfer from an oscillating cylinder to power-law fluids at Re = 100, Pr = 1.0, and $f^{*}=1.0$.

${\mathit{H}}^{*}$	$\mathit{n}=\mathbf{0.6}$	$\mathit{n}=\mathbf{1.0}$	$\mathit{n}=\mathbf{1.4}$
1.0	6.677, 0.431	5.993, 0.295	5.563, 0.262
0.5	6.368, 0.115	5.786, 0.086	5.394, 0.067
0.25	6.296, –	5.683, –	5.312, –

, where the amplitudes for $H_0=0.25$ are negligible (less than 1.0%) and thus not presented. The results show that the time-averaged Nusselt number is much higher in the shear-thinning ($n<1.0$) fluid than in the shear-thickening ($n>1.0$) fluid, which agrees with that in the previous steady cases. It also shows that the increasing oscillating amplitude enhances the heat convection effects, and the oscillation of the averaged Nusselt number tends to be more remarkable with increasing $H_0$. Figure 9 presents time histories of the averaged Nusselt numbers for $H=1.0$ in two periods. The averaged Nusselt number has a frequency of $2f$. The snapshots of the temperature fields for $n=0.6$ and $H_0=1.0$ in a period are further presented in Figure 10.

5.3. Effects of the Oscillating Frequency

Here, in order to study the effects of the oscillating frequency on the heat transfer around an oscillating cylinder in a uniform flow, the oscillating amplitude is fixed at $H^{*}=1.0$ with Re = 100. Two more oscillating frequencies $f^{*}=0.8$ and $1.2$ are considered for comparison. The time-averaged Nusselt numbers and amplitudes are presented in

Table 6. Time-averaged Nusselt number ($Nu$) and its amplitude ($\Delta Nu$) for forced convection heat transfer from an oscillating cylinder to power-law fluids at Re = 100, Pr = 1.0, and $H^{*}=1.0$.

${\mathit{f}}^{*}$	$\mathit{n}=\mathbf{0.6}$	$\mathit{n}=\mathbf{1.0}$	$\mathit{n}=\mathbf{1.4}$
0.8	6.326, 0.283	5.776, 0.224	5.408, 0.163
1.0	6.677, 0.431	5.993, 0.295	5.563, 0.262
1.2	6.901, 0.591	6.152, 0.421	5.567, 0.351

. The results show that the increasing oscillating frequency of the cylinder enhances the heat convection, resulting in larger averaged Nusselt numbers and amplitudes for both shear-thinning and shear-thickening fluids. However, the enhancement in the shear-thinning fluid is much more significant than that in the shear-thickening fluid.

The instantaneous vorticity and temperature contours for $n=0.6$ are presented in Figure 11. Similar Kármán-street type wake is observed in all cases. It is found that the increasing oscillating frequency of the cylinder tends to increase the vertical distance between the two adjacent vortices. Specifically, the vortices shed at $f^{*}=0.8$ are almost in a line (as marked in Figure 11a) as oscillating frequency increases, and the vortices tend to separate vertically (as marked in Figure 11c). The temperature contours show that the heat transfer in the fluid exactly follows these vortices.

6. Conclusions

In this paper, a hybrid immersed boundary–lattice Boltzmann and finite difference method for fluid–structure interaction and heat transfer in non-Newtonian flow is presented, with the geometry-adaptive Cartesian mesh technique to enhance the computational efficiency. The present numerical method includes four parts: the fluid solver, heat transfer solver, structural solver, and immersed boundary method for fluid–structure interaction and heat transfer. Specifically, the multi-relaxation-time-based LBM is adopted for the computation of non-Newtonian flow, with the geometry-adaptive Cartesian mesh technique to enhance the computational efficiency. The heat transfer equation is spatially discretized by a second-order up-wind scheme for the convection term, a central difference scheme for the diffusion term, and a second-order difference scheme for the temporal term. The structural dynamics is numerically solved using a finite difference method. The no-slip boundary is achieved using an immersed boundary method.

The proposed method and developed solver is validated with standard benchmark problems involving the non-Newtonian fluids flow, i.e., the developing flow of non-Newtonian fluid in a channel, heat transfer around a stationary cylinder and flow around a stationary cylinder with a detached filament. The good agreements achieved by the present method with the published data show that the present extension is an efficient way for fluid–structure interaction and heat transfer involving non-Newtonian fluid.

The heat transfer around an oscillating cylinder in non-Newtonian fluid flow is also numerically studied using the present method, considering the effects of the oscillating frequency and amplitude. The results show that the increasing frequency and amplitude tend to enhance the heat convection, resulting in larger time-averaged Nusselt numbers. The enhancement that occurs in the shear-thinning fluid is much more significant than that in the shear-thickening fluid. The present results may be used to expand the currently limited database of FSI and heat transfer benchmark studies.