# Optimal heat transfer enhancement in plane Couette flow

###### Abstract

Optimal heat transfer enhancement has been explored theoretically in plane Couette flow. The velocity field to be optimised is time-independent and incompressible, and temperature is determined in terms of the velocity as a solution to an advection-diffusion equation. The Prandtl number is set to unity, and consistent boundary conditions are imposed on the velocity and the temperature fields. The excess of a wall heat flux (or equivalently total scalar dissipation) over total energy dissipation is taken as an objective functional, and by using a variational method the Euler–Lagrange equations are derived, which are solved numerically to obtain the optimal states in the sense of maximisation of the functional. The laminar conductive field is an optimal state at low Reynolds number . At higher Reynolds number , however, the optimal state exhibits a streamwise-independent two-dimensional velocity field. The two-dimensional field consists of large-scale circulation rolls that play a role in heat transfer enhancement with respect to the conductive state as in thermal convection. A further increase of the Reynolds number leads to a three-dimensional optimal state at . In the three-dimensional velocity field there appear smaller-scale hierarchical quasi-streamwise vortex tubes near the walls in addition to the large-scale rolls. The streamwise vortices are tilted in the spanwise direction so that they may produce the anticyclonic vorticity antiparallel to the mean-shear vorticity, bringing about significant three-dimensionality. The isotherms wrapped around the tilted anticyclonic vortices undergo the cross-axial shear of the mean flow, so that the spacing of the wrapped isotherms is narrower and so the temperature gradient is steeper than those around a purely streamwise (two-dimensional) vortex tube, intensifying scalar dissipation and so a wall heat flux. Moreover, the tilted anticyclonic vortices induce the flow towards the wall to push low- (or high-) temperature fluids on the hot (or cold) wall, enhancing a wall heat flux. The optimised three-dimensional velocity fields achieve a much higher wall heat flux and much lower energy dissipation than those of plane Couette turbulence.

Key words: Mixing enhancement, variational methods

## 1 Introduction

The improvement of the performance of heat exchangers has always been crucial to efficient energy use. What is a flow field with the ability of the most efficient heat transfer? This may be a naive question in the development of flow control technique aimed at heat transfer enhancement. If we can find an optimal state in which heat transfer will be enhanced while suppressing energy loss, it will serve as the target of flow control.

For buoyancy-driven convection, Malkus (1954) raised a question of an upper limit to heat transfer and an optimal state realising it, and Howard (1963) formulated a variational problem to derive a theoretical upper bound on a heat flux. In the case of velocity fields satisfying the equation of continuity, Howard obtained the upper bound under the assumption that an optimal state has a single horizontal wavenumber, and subsequently Busse (1969) considered the same problem by using a multiple-boundary-layer approach to derive an asymptotic upper bound in a large-Rayleigh-number limit. Moreover, by using the same method, Busse (1970) found an upper bound on momentum transfer (or equivalently energy dissipation and so enstrophy), and conjectured an optimal velocity field consisting of hierarchical streamwise vortices in plane Couette flow with the aid of the so-called attached-eddy hypothesis (Townsend 1976). Later Doering and Constantin developed a new variational approach called ‘the background method’ and obtained a rigorous upper bound on energy dissipation in plane Couette flow (Doering & Constantin 1992, 1994) and a heat flux in Rayleigh–Bénard convection (Doering & Constantin 1996). The development of this method has triggered remarkable advancements in identification of optimal momentum transfer in shear flows (Nicodemus et al. 1997, 1998a, 1998b; Plasting & Kerswell 2003) and optimal heat transfer in thermal convection (Kerswell 2001; Otero et al. 2002; Doering et al. 2006). Recently, Hassanzadeh et al. (2014) have reported optimal velocity fields for two-dimensional thermal convection without taking account of the Navier–Stokes equation as a constraint on fluid motion. They considered two-dimensional divergence-free velocity fields under one of the two types of constraints on velocity fields (fixed kinetic energy or fixed enstrophy). On the other hand, Sondak et al. (2015) have taken another approach to optimal heat transfer in the two-dimensional thermal convection, in which they optimised the aspect ratio of the two-dimensional domain for their steady solutions to the Navier–Stokes equation to achieve maximal heat transfer. To our knowledge, however, the velocity field for optimal heat transfer has not been determined in shear flows as yet.

Shear flow turbulence has an ability of significantly high heat transfer in comparison with laminar flow, but it results in friction loss as a consequence of simultaneous promotion of momentum transfer. This is well known as the similarity between heat and momentum transfer in engineering (Reynolds 1874; Chilton & Colburn 1934). Especially, when the Prandtl number, the ratio between thermal and momentum diffusivity, is close to unity, the heat and momentum transfer exhibit strong similarity (Dipprey & Sabersky 1963). Although one of the ultimate goals in turbulence control is more heat transfer and less friction, their achievement has been extremely difficult due to the similarity. Recently, however, control strategies to break the similarity have been proposed. Momentum transfer by viscous flow obeys the Navier–Stokes equation, whereas heat transfer is governed by the energy equation, i.e. an advection-diffusion equation for temperature. By applying the feedback control using the adjoint equations (Bewley et al. 2001; Kasagi et al. 2012) to heat transfer in turbulent channel flow with blowing and suction on the wall, Hasegawa & Kasagi (2011) and Yamamoto et al. (2013) have numerically achieved the dissimilarity even when the Prandtl number is equal to unity. Experimentally, several practicable passive or active control techniques have been developed using riblet (Sasamori et al. 2014), permeable wall (Suga et al. 2011), micro actuator and sensor (Kasagi et al. 2009), and so on. Hence, the accomplishment of dissimilar heat transfer enhancement by flow control might be expected in the near future.

The aim of this study is to find an optimal velocity field which will provide new insight into such flow control. We focus on optimal heat transfer in plane Couette flow. In this flow, the temperature and velocity field exhibit complete similarity in a laminar and conductive state, when a constant temperature difference is imposed between the two parallel walls. Furthermore, when the Prandtl number is set to unity, a turbulent state shows strong similarity between heat and momentum transfer. By taking the excess of the scalar dissipation (or equivalently, a total wall heat flux) over the energy dissipation as an objective functional, an optimal velocity field is determined by maximising the functional under the constraint of the continuity equation for the velocity and the advection-diffusion equation for the temperature.

This paper is organized as follows. In §2 we describe the flow configuration and the mathematical relation between scalar (or energy) dissipation and a wall heat (or momentum) flux, and then present the formulation of a variational problem. In §3, the details of numerical procedures and the parameter dependence of optimal velocity fields are explained. The spatial structures and the statistics of the optimal states are presented in §4, and the hierarchical structure of the optimal fields is discussed in §5. In §6 we discuss the significant effects of three-dimensional vortical structures on heat transfer enhancement in the optimal state. Section 7 is devoted to summary and discussion. Properties of the external body force to be added for the achievement of the optimal velocity field are described in Appendix A. An analytical interpretation of heat transfer intensification mechanisms through the spanwise inclination of a quasi-streamwise tubular vortex is given in Appendix B.

## 2 Mathematical formulation

### 2.1 Flow configuration

Figure 1 shows the configuration of the velocity and temperature fields. The flow is driven by the two parallel plates moving in the opposite directions at a constant speed . The upper (or lower) wall surface is held at higher (or lower) constant temperature (or ). The coordinates, , and (or , and ), are used for the representation of the streamwise, the wall-normal and the spanwise directions, respectively, and their origin is on the midplane of the two plates separated by a distance . The corresponding components of the velocity are given by , and (or , and ) in the streamwise, the wall-normal and the spanwise directions, respectively.

Let us consider heat transfer in an incompressible time-independent velocity field fulfilling the continuity equation

(2.0) |

We suppose that the temperature field is determined as a solution to an advection-diffusion equation

(2.0) |

where denotes thermal diffusivity. The velocity and temperature fields are supposed to be periodic in the - and -directions with the periods, and . The boundary conditions

(2.0) |

are imposed on the walls.

### 2.2 Scalar dissipation and heat flux

We introduce the -average and the -average given respectively by

(2.0) | |||||

(2.0) |

We take the -average of equation (2.1) to obtain

(2.0) |

the -integration of which gives us the total heat flux

(2.0) |

at any wall-normal position . By further taking the average of equation (2.2) in the -direction we have

(2.0) |

Now we decompose the temperature into a conductive state and a fluctuation about it as

(2.0) |

The boundary conditions of the temperature fluctuation are

(2.0) |

The substitution of decomposition (2.2) in (2.2) yields

(2.0) |

We substitute decomposition (2.2) for equation (2.1) to have

(2.0) |

The product of equation (2.2) with gives us

(2.0) |

By further taking the -average of equation (2.2) and taking into account boundary conditions (2.2), we obtain

(2.0) |

where the right-hand side represents the dissipation of variance , hereafter referred to as ‘scalar dissipation’. As a consequence, we have the expression of the wall heat flux

(2.0) |

### 2.3 Energy dissipation and momentum flux

The steady motion of a viscous fluid is described by the Navier–Stokes equation

(2.0) |

where , , and are pressure, mass density, kinematic viscosity, and the -th component of external body force per unit mass, respectively. The inner product of equation (2.3) with the velocity yields the local energy budget equation

(2.0) |

By taking the -average of equation (2.3) and taking into account boundary conditions (2.1), we obtain the total energy budget equation

(2.0) |

where we have supposed a null mean pressure gradient in the wall-parallel (- and -) directions. The left-hand side of equation (2.3) is mean wall shear stress, while the first and second terms in the right-hand side originate from energy dissipation and an energy input by the external force, respectively.

### 2.4 Objective functional and Euler-Lagrange equation

Let us now consider the optimisation problem of a functional to find an optimal state for the dissimilarity between momentum and heat transfer. As can be seen from equation (2.2) and (2.3), the wall heat flux corresponds to the total scalar dissipation, and the wall momentum flux (i.e. the skin friction) is related with the total energy dissipation. Therefore, we shall search for an optimal state which maximises the functional given by

(2.0) |

where () is the dimensionless weight of the contribution from the heat transfer against the momentum transfer, is specific heat at constant pressure, and and represent Lagrangian multipliers. We suppose that satisfies boundary conditions on the walls,

(2.0) |

By introducing the following dimensionless quantities

we rewrite the objective functional (2.4) in a dimensionless form as

(2.0) |

where we have dropped the primes attached to the dimensionless quantities,

(2.0) |

is new dimensionless weight, and

(2.0) |

are the Reynolds number, the Prandtl number, and the Eckert number, respectively. Usually , and thus we suppose that so that we may take . We decompose and into references and infinitesimal perturbations as

(2.0) |

Substituting (2.0) into (2.4), using integration by parts and recalling boundary conditions (2.1), (2.2) and (2.0), we obtain the first variation of as

(2.0) |

A stationary point, at which the first variation of vanishes, i.e. , for any perturbation and , is determined by the Euler–Lagrange equations

(2.0) | |||||

(2.0) | |||||

(2.0) | |||||

(2.0) |

## 3 Numerical procedures

### 3.1 Finding initial guesses

We find velocity fields which give a maximal value of the objective functional with a combination of the steepest ascent method and the Newton–Krylov method. To globally seek a maximal point of , we start with finding initial guesses in the following procedure.

Since (3.0) implies that

(3.0) |

in order that is incompressible we invoke the continuity for ,

(3.0) |

Taking the double curl, the single curl and the -average of equation (2.0), we have the equations for the wall-normal velocity , the wall-normal vorticity , and the averaged velocities,

(3.0) | |||||

(3.0) | |||||

(3.0) | |||||

(3.0) |

with

(3.0) | |||||

(3.0) |

from which has been eliminated. The corresponding boundary conditions are

(3.0a) | |||||

(3.0b) | |||||

(3.0c) | |||||

(3.0d) | |||||

(3.0e) | |||||

(3.0f) |

In the periodic streamwise and spanwise directions, , vorticity , , and are approximated by truncated Fourier series

(3.0) |

where and are positive integers. Note that the Fourier coefficients of and for are given by and through

(3.0) | |||||

(3.0) |

In the wall-normal direction, , , , , and are expanded, in terms of the polynomials which satisfy boundary conditions (3.1), as

(3.0a) | |||||

(3.0b) | |||||

where represents the -th order Chebyshev polynomial. The nonlinear terms are evaluated using a spectral collocation method. Aliasing errors are removed with the aid of the rule for the Fourier transform and the rule for the Chebyshev transform.

### 3.2 Newton–Krylov iteration

In nonlinear dynamical systems stemming from the Navier–Stokes equation, a combination of the Newton–Raphson method and the Krylov subspace method has been used to efficiently obtain nonlinear equilibrium and periodic solutions (Viswanath 2007, 2009). In this study, employing the initial guesses obtained as mentioned in §3.1, we perform the Newton–Krylov iteration to find nonlinear solutions to the Euler–Lagrange equations (2.0)–(2.0).

Introducing pseudo-time derivative terms into equations (3.0)–(3.0), the -average of equation (2.0), equations (2.0) and (2.0), we have the following nonlinear evolution equations

(3.0) | |||||

(3.0) | |||||

(3.0) | |||||

(3.0) | |||||

(3.0) | |||||

(3.0) | |||||

(3.0) | |||||

(3.0) |

‘steady’ solutions to which satisfy the Euler–Lagrange equations (2.0)–(2.0). By using a spectral Galerkin method, we can convert equations (3.0)–(3.0) into a nonlinear autonomous dynamical system

(3.0) |

where represents a state point in the phase space spanned by the Fourier–Chebyshev coefficients of and . The time advancement is carried out by the Crank–Nicholson scheme for the diffusion terms and the 2nd-order Adams–Bashforth scheme for the others.

Let be an initial guess to find a maximal point of . Integrating equation (3.0), we define the flow map as

(3.0) |

for any ‘time’ , where . Supposing that is sufficiently close to a stationary point and thus is a small correction vector, we linearise equation (3.0) at to obtain

(3.0) |

where

(3.0) |

being an identity matrix with the same dimension as the Jacobian matrix . By solving the linear equation (3.0) the correction vector is obtained to renew as

(3.0) |

We iterate the procedure until a termination criterion

(3.0) |

is satisfied, where denotes an -norm. To solve equation (3.0) we use the GMRES method (Saad & Schultz 1986) which is one of the Krylov subspace methods. In this method, an approximation of a correction vector can be got without calculating matrix , by minimising the norm of a residual vector in a Krylov subspace.

### 3.3 Direct numerical simulation of turbulent plane Couette flow

For comparison with optimal states obtained by the maximisation of the functional , direct numerical simulations are also performed for turbulent momentum and heat transfer in plane Couette flow. In the simulations we numerically integrate the incompressible Navier–Stokes equation and the advection-diffusion equation for a passive scalar, i.e. temperature, supplemented by boundary conditions (2.1). The evolution equations are discretised employing a spectral Galerkin method based on the expansions (3.0) and (3.1). Time marching of the equations is done with a combination of the Crank–Nicholson and 2nd-order Adams–Bashforth scheme. The Prandtl number is set to unity. Turbulence data are accumulated in the range of the Reynolds number, –, in a computational domain of .

### 3.4 Global convergence and dependence on domain size

Label | |||
---|---|---|---|

A | |||

B |

We perform the maximisation of the functional starting with 30 different initial guesses extracted from turbulent flow fields in order to examine the dependence of a computed local maximum on the initial guesses. Figure 2 shows the functional against the number of the iteration mentioned in § 3.1 for the distinct 30 initial guesses at for and . For all the initial guesses, the value of rapidly approaches the vicinity of the two local maxima to be mentioned just below. We have eventually obtained the two different local maxima (the larger local maximum is labeled as A, and the smaller one as B). Table 1 compares the values of together with the total scalar dissipation and the total energy dissipation ,

(3.0) |

between the local maxima A and B. The present objective functional has multimodality, but there appear only two local maxima for many initial guesses.

The optimal velocity fields exhibit three-dimensional structures which consist of large-scale streamwise circulation rolls and smaller-scale quasi-streamwise vortex tubes near the wall (see figures 5 and 6). The size of the large rolls corresponds to half the spanwise period . Table 2 summarises each maximal state obtained in the various domains at for and . It is found that takes the largest value for among , , ; , , , , suggesting that there exist finite values of the streamwise and spanwise periods which lead to the global optimal. However, our goal in this paper is not to identify the optimal domain for all the parameters including . As will be described below, the velocity fields, which are realised as a solution to the Euler–Lagrange equations (2.0)-(2.0), have extremely high performance for heat transfer and exhibit qualitatively similar flow structures being roughly independent of . In the following, therefore, we shall choose the domain size to discuss characteristic flow structures and the mechanism of heat transfer enhancement. The effects of the domain size on the scalar dissipation and the energy dissipation can be seen to be minor.

1000 | 1.0 | 0.1 | (32,128,256) | |||||

1000 | 1.0 | 0.1 | (32,128,256) | |||||

1000 | 1.0 | 0.1 | (64,128,256) | |||||

1000 | 1.0 | 0.1 | (32,128,256) | |||||

1000 | 1.0 | 0.1 | (64,128,256) | |||||

1000 | 1.0 | 0.1 | (64,128,512) |

### 3.5 Effects of weight on optimal state

We next consider the effects on the optimal state of the weight of the scalar dissipation in the objective functional .

Let us first note that in the special case of a streamwise-independent two-dimensional velocity field the weight can be removed from the Euler–Lagrange equations (2.0)–(2.0) as follows. Taking the curl of equation (2.0), we have

(3.0) |

where

(3.0) |

is the Jacobian determinant. The streamwise-independent velocity field can be expressed, in terms of the streamfunction , as

(3.0) |

where the streamwise velocity has been uniquely determined from the -component of equation (2.0) and its independence of the -direction (). The streamfunction is related to the streamwise vorticity via the Poisson equation

(3.0) |

Using , equations (2.0) and (2.0) are respectively rewritten as

(3.0) | |||

(3.0) |

The boundary conditions of are

(3.0) |

Introducing the new variable and the new parameter as

(3.0) |