A spectral/hp element method for thermal convection

We report on a high‐fidelity, spectral/hp element algorithm developed for the direct numerical simulation of thermal convection problems. We consider the incompressible Navier–Stokes (NS) and advection–diffusion equations coupled through a thermal body‐forcing term. The flow is driven by a prescribed flowrate forcing with explicit treatment of the nonlinear advection terms. The explicit treatment of the body‐force term also decouples both the NS and the advection–diffusion equations. The problem is then temporally discretized using an implicit–explicit scheme in conjunction with a velocity‐correction splitting scheme to decouple the velocity and pressure fields in the momentum equation. Although not unique, this type of discretization has not been widely applied to thermal convection problems and we therefore provide a comprehensive overview of the algorithm and implementation which is available through the open‐source package Nektar++. After verifying the algorithm on a number of illustrative problems we then apply the code to investigate flow in a channel with uniform or streamwise sinusoidal lower wall, in addition to a patterned sinusoidal heating. We verify the solver against previously published two‐dimensional results. Finally, for the first time we consider a three‐dimensional problem with a streamwise sinusoidal lower wall and sinusoidal heating which, for the chosen parameter, leads to the unusual dynamics of an initially unsteady two‐dimensional instability leading to a steady three‐dimensional nonlinear saturated state.

higher-order accuracy and preferential convergence (usually exponential) properties of spectral methods by applying a high-order polynomial of order P (hence p-type) in every elemental domain of a coarse finite element type mesh. Spectral/hp element methods can be considered as a high-order extension of classical (traditionally low-order) finite element methods. 2 To the best of our knowledge there have been only a handful of studies which deal with thermal convection using spectral element methods. Saha et al. 3 and Gaddar et al. 4 developed higher-order spectral element models suitable for forced convection without considering thermal body-forcing in the incompressible Navier-Stokes (NS) equations. BniLam and Al-Khoury 5 developed a two-node spectral element model in Fourier space to study a shallow geothermal system consisting of an axisymmetric domain with a borehole heat-exchanger and a soil layer. Most recently, using the spectral/hp element method, Kumar and Pothérat 6 studied a stably stratified convective flow in a semicylindrical cavity with a heated upper boundary bounded by a porous-cold semicircular boundary. They considered through-flow fed in the upper boundary and used a quasi-three-dimensional approach (infinitely extended in the third dimension). They performed direct numerical simulation (DNS) and linear stability analysis of the flow considered.
Certain flow cases such as channel flow, pipe flow, flow with periodic-patterns of geometry (e.g., multibladed turbomachinery or periodic surface roughness) and thermal conditions (spatially periodic wall temperature/heatflux), are often simulated by considering an inflow that is periodic with the outflow. To drive such flows either a constant pressure gradient (essentially a body forcing) or a constant flow rate condition needs to be enforced along the direction of the flow. Although it is easier to implement a body forcing in most numerical simulations, in practice it is difficult to determine the correct magnitude of the body force, especially in flows with complex geometry or turbulent flows. Therefore, it is desirable to implement a constant mass flow rate condition [7][8][9] and in this article we will outline the procedure for this leveraging previous work for nonthermal problems. This is theoretically equivalent to the implementation suggested previously for isothermal cases 7 but we have modified the implementation to consistently manage the addition of the thermal field.
Consideration of the published literature suggests that there is an unmet need to develop a spectral/hp element solver for the DNS of thermal convective flows which incorporates thermal body-forcing, periodic boundary condition, and supports a constant mass flow forcing constraint. While other open-source spectral element codes, such as SEMTEX 10 and NEK5000, 11 support the simulation of thermal convection problems, to the best of our knowledge none support the specific capabilities outlined here. This acts as the primary motivation for this work. Addition of the body forcing term in the incompressible NS equations leads to a coupling of the equations with the energy equation (EE). This contributes an extra challenge in solving the system. In addition to this, the thermal body-forcing term generates a spurious pressure field when one attempts to maintain a constant flow-rate in a periodic boundary setting. Since science research is gradually progressing towards an open-access philosophy, 12,13 we make the source-code of the solver freely available through the open-source package Nektar++ * .
The organization of this article is as follows. Section 2 outlines the governing equations and the temporal and spatial discretization used, with a brief overview of the resulting matrix systems, how boundary conditions are implemented, and how the constant flow-rate condition is achieved. Since the system is unsteady, in Section 3, we discuss the procedure to determine a steady-state solution which is typically necessary for stability analysis, which is to be considered in subsequent work. In Section 4, we validate the solver against previously published results, illustrate some demonstrative flow examples dealing with thermal convection and discuss the novel DNS and saturation state of flow in a heated wavy channel. Finally, we provide a brief summary in Section 5.

Governing equations
We consider the three-dimensional unsteady thermally convective flow of an incompressible, Newtonian fluid in a three-dimensional geometric domain Ω. We use a Cartesian coordinate system with positions denoted by x(x, y, z) and time denoted by t. We assume the fluid has density , kinematic viscosity , thermal diffusion coefficient , and thermal expansion coefficient . Assuming variations of follow the Boussinesq approximation, the flow and temperature fields (with constant material properties) are therefore governed by, Momentum equation: Advection-diffusion equation: Continuity equation: where u(u, v, w) is the velocity, p is the pressure and is the temperature. A body force f b (0, f b , 0) represents the effect of heating imposed on the boundary. The problem is nondimensionalized by the Reynolds number Re = UL , where the characteristic length scale L is given by, for example, half the mean channel height and the characteristic velocity scale U is given by, for example, the maximum inlet velocity. The pressure field is nondimensionalized by the dynamic pressure scale defined as U 2 . The Prandtl number is given by Pr = . For the body force term, we use an explicit form, f b =

Ra
Pr Re 2 , where Ra = g L 3 T s is the Rayleigh number based on the temperature scale T s . One may note that the choice of the temperature scale T s is problem-specific. For example, in the case of a channel flow with uniformly heated walls the temperature scale could be defined as the temperature difference between the hot and cold walls, whereas for the case of a channel flow with spatial periodic heating, the amplitude of the heating could be considered as the temperature scale. The nondimensional temperature is defined as, = (T − T r )∕T s , where T is a dimensional temperature and T r is a reference temperature. For the example cases discussed in this study, the upper wall temperature is assumed to be the reference temperature T r .
To determine a closed-form solution the system of equations (1) is subject to an initial and boundary conditions on the velocity and temperature as Ω  denotes the part of the boundary where Dirichlet boundary conditions are imposed, whilst on Ω  Neumann boundary conditions are imposed.

Temporal discretization
The system (1) is discretized in time using a high-order stiffly stable splitting scheme 14,15 in order to decouple of the pressure and velocity fields. We treat the nonlinear advection terms explicitly which leaves the system coupled only through the body-force term f b . Explicit treatment of this body-force term in the NS equations further decouples the advection-diffusion equation. The temporal scheme used to solve the decoupled system is a projection scheme based on backward-differencing in time. This scheme evaluates the nonlinear advection terms explicitly, solves a Poisson system for the pressure field, and finally solves a series of Helmholtz type problems to enforce viscous forcing and boundary conditions (2). This method is one of a class of velocity-correction projection schemes. 16 A more formal description of the scheme is given below. For the simplicity of subsequent discussions, we write the momentum and the advection-diffusion (energy) equations (1a) and (1b) as a single equation Accordingly, the conditions (2) lead to The time-derivative is discretized using a backward difference formula, while the explicitly treated advection term, F a (V), and body force term, F b (V), are extrapolated from the previous time-steps using an Adams-Bashforth scheme. This leads to the following time-discrete formula for system (3), 1 Δt with where J e and J i are the integration orders of the explicit and implicit terms, respectively, and , , and 0 > 0 are the coefficients of the time-integration scheme given in References 14,17. Adapting their notation, the splitting scheme then proceeds by first advancing the advection and body force termŝ followed byV Here, we note that since the temperature component of is identically zero,̂≡̂and we only need to solve for the pressure p, which couples the velocity components in the momentum system (5). The splitting scheme exploits the incompressibility assumption of the fluid, that is, ∇ ⋅û = 0 in Equation (7), and enables us to obtain the pressure field from the intermediate velocity by solving the Poisson problem whereû is the momentum component ofV. To reduce errors due to the splitting scheme, we employ a consistent (high-order) pressure Neumann boundary condition 14 as where n is the domain unit outward normal, J p is the order of integration used in this boundary condition.
Finally, we calculate the new solution V k + 1 as, leading to the set of four scalar Helmholtz equations: and augmented with appropriate velocity and temperature boundary conditions. One can therefore observe that the solution procedure involves successive evaluation solutions of five elliptic problems: a Poisson problem and four (in 3D) Helmholtz problems.

Spatial discretization
In this section, we give an overview of the spectral/hp element method implemented in Nektar++. The computational domain Ω is discretized, and system (3) solved, by first casting equations (8) and (10) into their weak form using a continuous Galerkin projection, which enforces C 0 -continuity across element boundaries. 17 The five scalar elliptic problems are of the form where ≥ 0. We take the inner product of Equation (11) with respect to a test function v(x), The test functions are defined such that v( Ω  ) = 0. Applying the divergence theorem to the first term of Equation (12) we obtain, which introduces the boundary integral ⟨v, ∇u ⋅ n⟩ Ω = ∫ Ω v ∇u ⋅ n dx, and n is the outward normal to the boundary Ω. Similar to the h-type finite element method, the solution domain Ω is subdivided into N el nonoverlapping physical elements Ω e such that Ω = ⋃ N el e=1 Ω e . For two-dimensional domains, Nektar++ implements quadrilateral and triangular elements whereas for three-dimensional domains tetrahedral, hexahedral, prismatic, and pyramid shaped elements are implemented. Each physical element, Ω e , is mapped to the corresponding standard reference element Ω st , on which the numerical integration and differential operators are constructed, under the mapping x = e ( ) ∶ Ω st → Ω e . As an example, for quadrilateral elements, the standard element Ω st is defined as the square Standard regions for other element shapes are defined similarly, with simplex elements using a Duffy transformation 18 to handle the collapsed coordinate system.
On each standard element the discrete solution u is approximated as, for example, whereû e pq are elemental coefficients, p and q represent the one-dimensional polynomial basis functions of order P (hence p-type) in each of the coordinate directions, with P = 1 corresponding to the well-known linear finite element expansion. The expansion basis consists of either modal (hierarchical) or nodal (point-based) basis functions which are well-suited for element boundary-interior modes decomposition and form a global C 0 continuous expansion. We select the family of Legendre polynomials as basis functions, modified such that interior modes do not have support on the element boundary. Integrations are performed using Gauss-Lobatto-Legendre quadrature points in the standard element Ω st . To avoid polynomial aliasing error 17,19 in the nonlinear (advection) terms, it is noted that the minimum number of quadrature points needed is 3(P + 1)/2. One can appreciate that this polynomial aliasing error is associated with the higher-order spectral element method, in addition to the traditional Fourier aliasing error associated with the use of a complex Fourier series to denote a flow quantity in a homogeneous third dimension (see Section 4.1).
For a Galerkin projection, the test functions are chosen to be the same as the trial functions . The continuous inner-product on a physical element can then be recast in discrete form as where  is the Jacobian of and w i and w j are the Gaussian quadrature integration weights. Derivatives can similarly be expanded in terms of derivatives of basis functions. Expanding u and v in Equation (14), and expressing as a matrix equation, we arrive at the elemental contribution of the system (13), Here, the matrix L e is the discrete Laplacian, M e is the mass matrix, the vector e accounts for the discrete boundary integrals, and the vectorf e is the projection of the forcing functions onto the spectral element basis. For higher orders we exploit the sum-factorization technique 20 to improve computational efficiency. 2,21,22 The C 0 -continuity across element boundaries is enforced through an assembly mapping , with a sparse matrix equivalent A, which sums local contributions to form a global solution vectorû g under the relationû e = Aû g . The local elemental system (15) can then be recast as a global matrix system, where H is the weak Helmholtz matrix, and the system (16) is subject to imposition of boundary conditions. The global solutionû g is computed using a multilevel static condensation approach, 17 with the inner level system being solved using a preconditioned conjugate gradient method. A spectral vanishing viscosity technique, 23 which dampens high-frequency information when the flow is underresolved, is applied to stabilize the convergence property of the solution process.

Boundary conditions
In this section, we provide a very brief overview on the imposition of Dirichlet and Neumann boundary conditions used in Equations (4a)-(4c), since a subtlety discussion is available in Reference 17. For Dirichlet boundary condition, the surface integral vector e = 0 in the local elemental matrix system (15), since by definition v| Ω D = 0. We impose Dirichlet boundary conditions by lifting the known modes out of the system. We decompose the solution into an unknown homogeneous V  (x) contribution and a known nonhomogeneous contribution V  (x) which satisfies the given Dirichlet boundary condition (4b), such that, where, V  ( Ω) = 0. The known degrees of freedom on the Dirichlet boundaries can then be moved to the right-hand side and Equation (13) can be written as, Neumann boundary conditions are enforced through the boundary integral vector e by projecting the function V  onto the quadrature points and using the corresponding surface Jacobian for the contributing elements whose edges touch the Neumann boundary region.

Enforcing a constant flowrate
In addition to the thermal forcing, a constant flowrate forcing is employed as an external means to drive the flow. The volumetric flux, through a given cross-sectional surface C of area (C) is kept constant to a prescribed value Q p . Following the efficient Green's function method, 7 at the first-time step, the linear Stokes equation with unit forcing f s (per unit length along the flow direction) together with homogeneous boundary conditions is solved to obtain a velocity field u s as, It is noted that the solution u s of the linear Stokes system (18) trails the splitting scheme discussed above, and also constitutes a Helmholtz system. A correction factor s at each time-step is calculated as, Then the corrected flowũ k = u k + s u s has the desired flowrate Q p . The linear Stokes system is computed once at the start of the simulation. The complete solution algorithm is depicted in Figure 1. One may note that while solving the linear Stokes system (18) in a thermal convection setting, the presence of nonzero temperature boundary condition(s) and the thermal forcing term f b in system (1a) cause the pressure Neumann boundary condition (Equation (9)) to be nonzero. This is due to the fact that temperature energizes the momentum equations explicit terms (essentially an addition to the nonlinear term) and this sets a nonzero pressure Neumann boundary condition. Since the Stokes solver requires homogeneous pressure Neumann boundary conditions, this pressure boundary condition needs to be explicitly set to zero at the first time step, otherwise resulting spurious pressure fluctuations would lead to divergence of the solution.
One may note that values of the correction factor s is problem specific. For the examples (see Section 4) considered in this study, we have observed that values of s (at steady state) varies between O(10 −2 ) to O(10 −1 ). For a problem with uniformly heated parallel plates channel, s assumes a value of O(10 −2 ) whereas for a two-dimensional corrugated-patterned channel with spatial periodic heating, s assumes a value of O(10 −1 ).

CAPTURING STEADY-STATE SOLUTION
For the purpose of capturing the steady-state solution of the flow system (3), we extend the selective frequency damping (SFD) method developed by Jordi et al. 24 to incorporate the solution of the scalar quantity , briefly describing the method below. The system (3) can be written as, where  contains the NS and the EE operators, that is, the right-hand side of the system (3). In order to determine a target steady-state, a linear stabilization technique with feedback control that contains a low-pass time filter, with filter width Δ and coefficient , is employed as Steady-state is reached when V = V, hence we iterate (19) in time such that ||V − V|| ∞ ≤ , where is the specified tolerance level. This linear SFD solver acts as a wrapper around the original unsteady solver discussed in Section 2.

THERMAL FLOW SIMULATIONS
In this section, we consider both a canonical incompressible flow problem with thermal convection from uniform heating and a more complex case of patterned heating and geometry. Unless otherwise stated, a third-order implicit-explicit (IMEX) time integration scheme proposed by Karniadakis et al. 14 is adopted in the simulations. The inflow and outflow boundaries are considered to be periodic with a constant flowrate forcing being applied along the streamwise (x) direction.

Flow with uniform heating
We consider flow of a fluid in a two-dimensional parallel-plate channel with a fully developed plane Poiseuille profile where the lower plate of the channel is heated uniformly with L = 1 and an isothermal upper plate with U = 0. TA B L E 1 Norm of the L 2 error as a function of polynomial order P for the flow in a two-dimensional channel with lower plate exposed to uniform heating The channel length and height are both considered to be of two units in length. We consider the length scale L as the half-channel height, and the temperature scale T s as the temperature difference between the upper and the lower plates.
No-slip and no-penetration boundary conditions for velocity are imposed on the plates. The initial condition is prescribed as the fully developed velocity profile with linear temperature profile given by: with flowrate forcing (per unit cross-sectional area), Q p = ∫ C u 0 dy = 4∕3. The simulation is performed with the modal expansion basis and 100 quadrilateral elements (mesh with 10 × 10 elements) with flow conditions Re = 100, Ra = 1000, Pr = 0.71, and a time-step size Δt = 0.001. After 10 time units of simulation, no flow bifurcation is observed, and respective velocity and temperature profiles remain unaltered which is consistent with Gage and Reid. 25 Table 1 shows the L 2 error norm of u-velocity as a function of polynomial order P and indicates that spectral convergence is achieved by changing the order of the polynomial P. Identical results are obtained when heating is imposed to the upper plate only.
Next, we consider a quasi-three-dimensional channel with uniform heating at the bottom wall. For this geometry the third dimension z is homogeneous, and the flow quantities are expressed by a Fourier basis 26 (21) where N m is the number of Fourier modes used in the computation. A quasi-three-dimensional simulation is performed for the flow conditions mentioned for the above two-dimensional case with the fully developed initial condition (20) and the initial spanwise velocity w 0 = 0. Considering the constant flowrate forcing in the streamwise (x) direction, we capture a two-dimensional solution that is qualitatively similar to the solution of the problem discussed in the above and the spanwise velocity w remains at zero irrespective to the time units prescribed in the simulation. This conforms with the corresponding analytical solution which suggests that the presence of the third dimension would not amend the two-dimensionality of the solution.

Flow with patterned heating
To observe the performance of the algorithm with respect to the spatial variation of the wall temperature distribution, we first consider a two-dimensional fully developed channel flow with the lower wall exposed to a spatial heating pattern defined as L = 0.5 cos( x) with is the heating wavenumber (such a heating provides alternate heated and cooled zones along the wall) and an isothermal upper wall ( U = 0) as shown in Figure 2(A). The length of the channel is considered to be one-wavelength of the heating. For such a heating pattern, the temperature scale T s is assumed to be the amplitude of the heating, hence the Rayleigh number Ra dictates the amplitude of heating. 27 The computational domain consists of a 20 × 20 mesh with quadrilateral elements distributed such that the mesh is finer near the both walls, as shown in Figure 2(B). We performed a simulation using the modal expansion basis for a flow condition of Re = 5, Ra = 1000, Pr = 0.71, heating wavenumber = 1 with heating wavelength Λ = 2 ∕ . A steady-state solution is determined using the method described in Section 3 with tolerance  Note: max denotes maximum and rms denotes root-mean-square (L 2 norm). Abbreviation: C-S solver, Chebyshev-Spectral solver. a Comparison of the solutions of the current solver using P = 7 with respect to the C-S solver. The current spectral/hp element solution is compared with a Chebyshev-Spectral solver, 28 herewith named as C-S solver, and the flow quantities agree to at least four leading digits (see Table 2). The u-velocity and temperature fields as well as flow topology depicted by streamlines are displayed in Figure 3 and are in qualitative agreement with Reference 28. Similar flow topology was also observed experimentally. 29

Flow with patterned heating and pattered geometry
Surface topography can have profound influence on thermal convection and the dynamics of the flow. The combined effect of the heating pattern (nonuniform heating) and surface corrugation, as shown in Figure 4, is considered to examine the efficacy of the current algorithm, exploiting the advantages of using the spectral/hp element method. In this problem, we consider the geometric pattern of the lower wall as y L = −1 + y b cos( x) whereas the upper wall is straight with y U = −1, and we apply a spatial heating pattern at the lower wall as L = 0.5 cos( x + Φ) and keep the upper wall as isothermal (i.e., U = 0). Here, denotes heating and corrugation wavenumber, y b is the corrugation amplitude, and Φ denotes the phase shift between heating and corrugation patterns. The computational mesh again consists of a 20 × 20 grid with quadrilateral elements and a finer mesh near both walls. The channel length is considered to be one-wavelength of the heating and corrugation (i.e., units) and the mean-height of the channel is to be two units. We performed a simulation with the Gauss-Lobatto-Legendre Lagrangian nodal expansion base since for such a complex flow system nodal expansion base is well suited. [30][31][32] We considered a flow condition of Re = 1, Ra = 1000, Pr = 0.71, amplitude of wall-corrugation y b = 0.05, heating and corrugation wavenumber = 2, and the phase shift between heating and corrugation Φ = ∕2. The steady state condition was enforced using the method described in Section 3 with tolerance level = 10 −6 , and required 36.9 time units of simulation with a time-step size Δt = 0.001. Once again polynomial order of P = 6 and P = 8 were applied to check the computational convergence to steady state. The spectral/hp element solution is compared with an immersed boundary condition (IBC) based spectral solver, 33,34 (here named as IBC-S solver) and the flow quantities agree to at least three significant figures (see Table 3).
The velocity and temperature fields, together with the flow topology depicted by streamlines, displayed in Figure 5 are also in qualitative agreement with those reported by Hossain and Floryan. 33 After ascertaining the veracity of the computational algorithm, we now proceed to investigate a quasi-three-dimensional corrugated channel with patterned heating applied at the bottom wall as shown in Figure 6(A) which has not been previously studied. We investigate a flow case with a higher wall corrugation amplitude of y b = 0.3, for which a converged solution could not be achieved for the IBC-S solver 33 even for a two-dimensional case. The third dimension z is homogeneous and we define a z-wavenumber (spanwise) to denote the periodic z-length as Λ z = 2 ∕ . To resolve this length-scale we prescribe N z equi-spaced planes (with Δz = Λ z ∕N z , as shown in Figure 6(B)), to represent N z /2 complex Fourier modes as N m = 0, 1, 2, … , N z /2 − 1. Figure 7 illustrates the quasi-three-dimensional thermal flow case with heating and corrugation wavenumber = 3, Re = 10, Ra = 3500, Pr = 0.71, y b = 0.3, z-wavenumber = 1.56. Among the various parameters present in this problem, we consider three exemplary phase shifts between heating and corrugation: Φ = 0 so that the heating pattern is in-phase with the wall corrugation to provide the hotter wall temperature zones are located near the crest of the corrugation and the colder wall temperature zones are located near the trough of the corrugation; Φ = ∕2 provides the colder wall temperature zones located at the downward-slopes (crest to trough) of the corrugation whereas the hotter wall temperature zones are located at the upward-slopes (trough to crest) of the corrugation, and Φ = represents heating and corrugation patterns are out-of-phase with the hotter wall temperature zones located at the trough of the corrugation and the colder wall temperature zones are located near the crest of the corrugation. For all the simulations we applied number of z-planes N z = 16, polynomial order P = 7, time-step size Δt = 0.001, and specify an initial condition with the spanwise velocity component w = 0, however with a vertical velocity component v composed of a sinusoidal function of z, and add a similar function in the horizontal velocity component u, that is, where B is the amplitude of the sinusoidal function, Λ z is the spanwise length of the channel and the velocity field follows the incompressibility constraint. One may note that it is necessary in the quasi-three-dimensional approach to have some variation in the z-direction otherwise the semianalytical handling of the Fourier discretization would mean we would not trigger the three-dimensionality of the problem. To identify the flow saturation states, we observe the time evolution of the kinetic energy of each of the Fourier modes defined as E m (t) = 1 2 ∫ Ω ||û m || 2 dx whereû m is the complex Fourier coefficient of the velocity field as defined in Equation (21). The evolution of kinetic energy of the leading Fourier components (N m = 0, 1, 2, 3) is shown in Figure 7. We observe that for the in-phase case of Φ = 0 (see Figure 7(A)), the energy in mean mode (which contains the spanwise average energy of the flow) initially evolves as an unsteady two-dimensional response, as indicated by the oscillatory black line and denoted as mode 0. However, there is then a transfer of energy to higher modes and we first observe mode 1 (as indicated by the green line in the plot) grow exponentially and the initial two-dimensional transients in the mean mode are attenuated. This rather surprisingly leads to a steady three-dimensional saturated state within 50 time units. When Φ = ∕2, we observe that the initial unsteadiness of the mean mode energy is barely detectable, and the growth of mode 1 (also the higher modes) is quite slower. Though the system ultimately settles to a steady three-dimensional saturated state (after 200 time units), this three-dimensionality is weaker compared with the in-phase case, for example, at t = 300 units, the mean mode has an energy of 0.357 units whereas mode 1 has an energy of 0.0098 units, which means the energy of mode 1 is approximately 36.5 times lower than the energy of mean mode. On the other hand, when Φ = , that is, when heating pattern is out-of-phase with the corrugation pattern, mode 1 (and higher) gradually decreases and the system ultimately remains in the two-dimensional state.
The dominant flow structures of respective saturated states are ascertained by the vortex identification method Q-criterion, and are shown in Figure 8 together with the isosurfaces of w-velocity contours at three representative x-locations of the channel: near the inlet, at the middle and near the exit of the channel. The corresponding isosurfaces of the temperature contours along x-and z-plane are shown in Figure 9. For Φ = 0, the flow structure consists of a pair of longitudinal streamwise vortices, where the cross-section of the vortices are smaller in the narrower section of the channel at the crest of the corrugation and gradually expand towards the trough of the corrugation. It is evident from the w-velocity isosurfaces that the nonlinear interactions of the u and v velocities lead the flow to become three-dimensional even though the initial flow was two-dimensional. The temperature contours depict a clear nonlinear variation of respective magnitudes along the spanwise direction. For the phase shifts Φ = ∕2 and Φ = , pairs of transverse spanwise vortices with approximately constant cross-section are observed. For Φ = (out of phase case), zones of colder fluids (shown as negative in Figure 9) dominate while only zones close to the corrugation trough contain hotter (positive ) fluids. Therefore, it can be inferred that the phase shift Φ between the heating and the corrugation patterns can be used as an effective control parameter to produce desired flow-patterns.

SUMMARY
We presented a spectral/hp element algorithm for the DNS of thermal convection problems of incompressible and Newtonian fluids and which can also act as a basis for global linear stability analysis. We demonstrate how the SFD technique can be used as a wrapper function with this algorithm to capture the steady convective flow. The governing equations describing the thermal flow system are the incompressible NS and the EEs with Boussinesq approximation. In addition, we have discussed how a flowrate forcing constraint can be applied. When the nonlinear advection terms, as well as the thermal body-force terms, are treated explicitly we obtain a decoupled system. We then adopt a temporal discretization using a third-order IMEX scheme. The velocity and pressure fields in the momentum equations are decoupled through a higher-order splitting scheme. Illustrative analytical examples of thermal flows demonstrate that the algorithm is able to achieve a rapid convergence. The algorithm is also suitable to simulate nonuniform wall temperature variations as well as geometric nonuniformity. Building on previous research of two-dimensional flow involving patterned heating and geometry, we first perform a validation of the spectral/hp element solver against published results of both a Chebyshev-spectral solver and an immersed boundary solver and excellent agreement is demonstrated. The new spectral/hp element solver is also suitable for three-dimensional problems and using a quasi-three-dimensional discretization of a spectral/hp element-Fourier expansion. For the first time we study the flow in corrugated channel with sinusoidal patterned heating on the lower surface and, in the in-phase case, observe a saturated three-dimensional steady state is present when Re = 10, Ra = 3500. These three-dimensional states are not observed in out-of-phase cases. These flow states deserve further investigation to fully understand the thermal and fluid balance leading to the saturated state as well as the transition process which is likely to be elucidated through global stability analysis.
Finally, we make the algorithm available to the wider community through the open-source package Nektar++ † .

ACKNOWLEDGMENTS
Mohammad Z. Hossain and Spencer J. Sherwin would like to acknowledge the support under the Newton International Fellowship scheme funded by The Royal Society of UK. Grant reference number NIF/R1/180427. Mohammad Z. Hossain would like to acknowledge the generous help received from A. Cassinelli at the early stage of this project.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. Data derived from public domain resources through https://www.nektar.info.