There are infinitely many possible types of meanders, and here we consider only a special type of meandering channel, i.e., the sinusoidal meandering channel. These channels are influenced by various agents that cause fluid motion, with the most common being the pressure gradient. In some cases, the pressure gradient may be mechanical in nature. The flows within these channels are predominantly induced by the pressure gradient, and they have numerous practical applications (Webb and Bergles, 1981; Bergles and Webb, 1985; Jensen et al., 1997; Ligrani et al., 2003). The partition of escalation structures can be performed on those that are responsible for necessary changes in the thermal and physical properties of fluids under consideration and the ones marked by perfection in fluid maxing (Fiebig, 1995a; Fiebig, 1995b; Jacobi and Shah, 1995). Mixing with improved mechanics is the point of attention in engineering devices functioning at small Reynolds numbers (Re) with laminar patterns (Fiebig, 1998). Commonly, it is understood that improved mixing can be done via compelling the fluid from the laminar state into a turbulent one or by the formation vortex generators. Vortex generators are functional but have significant pressure in their drag drawback (Patera and Mikic, 1986; Fiebig and Chen, 1999).

For nearly a century, centrifugal instability has been observed in shear layers. Previous studies focused on simple geometries and canonical flows, making it easy to calculate the curvature or meander of flow patterns. One example is the flow of fluid motion between rotating cylinders, where (Floryan, 1991) investigated the stability states of the distribution of rotation in an inviscid mechanism. In (Mohammadi and Floryan, 2013), addressed a highly viscous problem was addressed, and the critical conditions of the secondary flow onset were determined. In (Xu et al., 2016), comparable instability in curved channels was analyzed, and (Rayleigh, 1917) presented a case of boundary layer flows on concave sheets, taking into account centrifugal instability. In (Taylor, 1923), it was found that the special type of instability is prominent for fluid flow over concave and convex sheets, subject to the condition that the velocity components are non-monotonic in nature. In (Ghalambaz et al., 2016; Hayat et al., 2017; Srinivasacharya and Sibanda, 2020; Mahmud and Uddin, 2021), numerical simulation and analysis of different types of nanofluids flowing past various types of surfaces in the context of fluid flow and heat transfer in different channels was performed (Khan and Ahmed, 2015; Alsaedi et al., 2016; Nadeem et al., 2016; Izadi and Pourmehran, 2017; Khan et al., 2017; Sheikholeslami et al., 2017). These studies explore the impact of various factors such as thermal radiation, magnetic fields, bioconvection, entropy optimization, and chemical reactions on flow and heat transfer characteristics. As such, these articles are highly relevant for those interested in studying the behavior of nanofluids. All of these results apply to constant meanders and establish a clear relationship between streamline meanders and the wall meander. Consequently, the demonstration of this special type of stability revealed that it is represented by a single parameter, which is one possible way to describe this meandering with a single parameter.

The main objective of this study is to investigate the flow of a viscous fluid in a meandering channel with wavy walls driven by a pressure gradient. The equations of motion, including the conservation of the mass and momentum equations and the energy equation, are utilized to analyze the fluid flow and temperature distribution within the channel. The boundary conditions are imposed on the channel walls and at the center of the channel to maintain symmetry. Dimensionless variables are used to transform the governing equations into a dimensionless form for easy tracking of units. Additionally, the study explores special types of instabilities and their significance. In the first part, an assumption of fixed pressure gradient is imposed, assuming direct channel flow with parallel plates, and sinusoidal channel flow is driven by an identical pressure gradient. Flow rate variation is used to supplement the flow losses associated with the meandering channel. The generalized model is solved using perturbation method, and the stream function is calculated and examined using other techniques. The temperature distribution behavior is illustrated through various plots depicting different dimensionless parameters, such as Pr, Re, and Ec.

Consider the flow of a viscous fluid in a meandering channel, which consists of two wavy walls separated by a fixed gap. The classical model for laminar flows in channels, tubes, and ducts is provided by Poiseuille. The flow in such a channel is primarily induced by a pressure gradient. There are various types of flows that pass through ducts and channels, and they have significant practical importance. The classification of scaling structures can be conducted based on those that have variations in the effective thermophysical properties of fluids, and the pursuit of optimal fluid mixing can lead to an infinite number of possible meander types. In this study, we will focus on the simplest sinusoidal channel, whose geometry is illustrated in Figure 1, where the superscript H denotes the upper and G denotes the lower walls of the channel. The channel extends along the x-axis from negative infinity to positive infinity, and the flow is driven by a pressure gradient.

The general form of continuity equation for incompressible viscous flow is: div V = 0 (1)

The velocity vector, denoted as V, has three components for three-dimensional flow. Similarly, the general form of the Navier–Stokes equations applies to steady and viscous flows of constant viscosity. ρ V ∙ ∇ V = − ∇ P + μ ∇ 2 V , (2)

The variables ρ, μ, and P represent density, viscosity, and pressure of the fluid, respectively. To define the dimensionless variables, we use asterisks and reference length L and velocity U, as follows: ∇ * = L ∇ , V * = V U , a n d   P * = P ρ U 2 , (3)

After substituting the non-dimensional variables, which are defined in Eq. 3, into Eqs 1, 2, we obtain the dimensionless continuity and Navier–Stokes equations: div V * = 0 , (4) ρ U 2 L V * ∙ ∇ * V * = − ρ U 2 L ∇ * P * + μ U L 2 ∇ * 2 V * , (5)

The above equation is multiplied by L ρ U 2 , and when simplified we get: ⇒ V * ∙ ∇ *   V * = − ∇ * P * + 1 R e ∇ * 2 V * . (6) where R e = ρ U L μ represents the Reynold’s number.

Let us consider that for flow in a straight channel and for two-dimensional flow, the velocity is represented by the velocity V 0 , pressure by P 0 , the stream function of the flow is denoted by ψ 0 , and Q 0 is the flow rate representative. Then, the continuity equation in two dimensional is transformed to: ∂ u 0 ∂ x + ∂ v 0 ∂ y = 0 , (7)

The x - component of the Navier–Stoke equations is: ρ ∂ u 0 ∂ t + u 0 ∂ u 0 ∂ x + v 0 ∂ u 0 ∂ y = − ∂ P 0 ∂ x + μ ∂ 2 u 0 ∂ x 2 + ∂ 2 u 0 ∂ y 2 , (8)

The y - component of the Navier–Stoke equations is: ρ ∂ v 0 ∂ t + u 0 ∂ v 0 ∂ x + v 0 ∂ v 0 ∂ y = − ∂ P 0 ∂ y + μ ∂ 2 v 0 ∂ x 2 + ∂ 2 v 0 ∂ y 2 , (9)

The boundary conditions for the problem are established by utilizing the no-slip condition, along with considering the geometry of the problem. Thus, the boundary conditions can be written as: u 0 = 0   a t   y = ±   1 . (10)

Equations 7–9 represent the continuity and the x   &   y - momentum equations, respectively, for a one-dimensional straight channel. These equations can be greatly simplified by making certain assumptions, namely, that u 0 ≠ 0   a n d   v 0 = 0 . With this assumption, the solution for the governing equations for straight channel flow can be obtained by combining the axial velocity u 0 , normal velocity v 0 , and pressure P 0 , which can be expressed as follows: V 0 = u 0 η , 0 = 1 − η 2 , 0 P 0 = − 2 ξ R e + C 1 , (11)

Note that Eqs 1–11 for the fully developed flow in a straight channel has been reported in F.M. White [31] and Schlichting [32].

The volume flow rate for the channel flow of fixed width is: Q 0 = ∫ − 1 1 u 0 η d η = 4 3 . (12)

The stream function ψ 0 is simply evaluated as: ∂ ψ 0 ∂ η = u 0 η ⇒ ψ 0 = − η 3 3 + η + C 2 . (13)

In a meandering channel, where the fluid flow moves in the positive x-axis direction, the Reynolds number is defined based on the maximum x-velocity and channel half-height.

The velocity field and other related field quantities for the fluid motion in the meandering channel can be expressed as follows: V b = V 0 + V 1 = u 0 η , 0 + u 1 ξ , η , v 1 ξ , η , P b ξ , η = P 0 ξ + P 1 ξ , η , ψ b = ψ 0 η + ψ 1 ξ , η , Q b = Q 0 + Q 1 , (14)

In the context of the meandering channel geometry, the velocity V 1 , pressure P 1 , stream function ψ 1 , and flow rate Q 1 are dimensionless quantities that are adjusted accordingly. Additionally, the velocity vector V b = u b , v b is decomposed. These components are further decomposed as follows: u b = u 0 η , u 1 ξ , η   a n d   v b = 0 , v 1 ξ , η ] , (15)

This is achieved by substituting the new variables for the velocity field in the dimensionless governing equations, i.e., Eqs 4, 6, and eliminating pressure term.

The continuity equation Eq. 4 in the new variable is: ∂ u b ∂ ξ + ∂ v b ∂ η = 0 , (16) and now using the decomposed forms of u b and v b from Eq. 15, we have: ⇒ ∂ u 1 ∂ ξ + ∂ v 1 ∂ η = 0 , (17)

Note that u 0 and v 0 are defined in Eq. 15, and some derivatives vanish: ∂ u 0 ∂ ξ = ∂ v 0 ∂ ξ = 0 ,

The dimensionless x − momentum Eq. 6 for new variables becomes: ∂ ∂ η ( u 0 ∂ u 1 ∂ ξ ) + ∂ ∂ η ( u 1 ∂ u 1 ∂ ξ ) + ∂ ∂ η ( v 1 ∂ u 0 ∂ η ) + ∂ ∂ η ( v 1 ∂ u 1 ∂ η = − ∂ 2 P 1 ∂ ξ ∂ η + R e − 1 ∂ 3 u 1 ∂ ξ 2 ∂ η + ∂ 3 u 1 ∂ η 3 , (18)

The dimensionless y − momentum Eq. 6 for new variables becomes: ∂ ∂ ξ ( u 0 ∂ v 1 ∂ ξ ) + ∂ ∂ ξ ( u 1 ∂ v 1 ∂ ξ ) + ∂ ∂ ξ ( v 1 ∂ v 1 ∂ η = − ∂ 2 P 1 ∂ ξ ∂ η + R e − 1 ∂ 3 v 1 ∂ ξ 3 + ∂ 3 v 1 ∂ ξ ∂ η 2 , (19)

The pressure term is eliminated by expending the terms ∂ ∂ η ( u 0 ∂ u 1 ∂ ξ ) , ∂ ∂ η ( v 1 ∂ u 1 ∂ η ) , ∂ ∂ ξ ( u 0 ∂ v 1 ∂ ξ ) . in Eqs 18, 19, and then by subtracting Eq. 19 from Eq. 18, we get: ∂ u 0 ∂ η   ∂ u 1 ∂ ξ + u 0 ∂ 2 u 1 ∂ ξ ∂ η + ∂ ∂ η ( u 1 ∂ u 1 ∂ ξ + v 1 ∂ u 1 ∂ η ) + ∂ v 1 ∂ η   ∂ u 0 ∂ η + v 1 ∂ 2 u 0 ∂ η 2 − u 0 ∂ 2 v 1 ∂ ξ 2 − ∂ ∂ ξ ( u 1 ∂ v 1 ∂ ξ + v 1 ∂ v 1 ∂ η ) } = R e − 1 ∂ 3 u 1 ∂ ξ 2 ∂ η + ∂ 3 u 1 ∂ η 3 − ∂ 3 v 1 ∂ ξ 3 − ∂ 3 v 1 ∂ ξ ∂ η 2 , (20)

We define a stream function ( ψ b ) such that: u b = ∂ ψ b ∂ η   a n d   v b = − ∂ ψ b ∂ ξ   ⇒ u 0 + u 1 = ∂ ∂ η ψ 0 + ψ 1   a n d   v 1 = ∂ ψ 1 ∂ η . (21)

Substituting in Eq. 20, we get: ⇒ u 0 ∂ ∂ ξ ∇ 2 ψ 1 − ∂ 2 u 0 ∂ η 2 ∂ ψ 1 ∂ ξ − R e − 1 ∇ 4 ψ 1 = ∂ ∂ ξ ( u 1 ∂ v 1 ∂ ξ + v 1 ∂ v 1 ∂ η ) − ∂ ∂ η ( u 1 ∂ u 1 ∂ ξ + v 1 ∂ u 1 ∂ η ) , (22) w h e r e ∇ 2 ψ 1 = ∂ 2 ψ 1 ∂ η 2 + ∂ 2 ψ 1 ∂ ξ 2   & ∇ 4 ψ 1 = ∂ 4 ψ 1 ∂ η 4 + 2   ∂ 4 ψ 1 ∂ ξ 2 ∂ η 2 + ∂ 4 ψ 1 ∂ ξ 4 .

From the above continuity Eq. 22, we have: ∂ u 1 ∂ ξ + ∂ v 1 ∂ η = 0 ⇒ u 1 ∂ u 1 ∂ ξ + ∂ v 1 ∂ η = v 1 ∂ u 1 ∂ ξ + ∂ v 1 ∂ η = 0 . (23)

By substituting in Eq. 22, we have: u 0 ∂ ∂ ξ ∇ 2 ψ 1 − ∂ 2 u 0 ∂ η 2 ∂ ψ 1 ∂ ξ − R e − 1 ∇ 4 ψ 1 = ∂ ∂ ξ ∂ ∂ ξ u 1 v 1 + ∂ ∂ η v 1 v 1 − ∂ ∂ η ∂ ∂ ξ u 1 u 1 + ∂ ∂ η u 1 v 1 . (24)

The no slip boundary conditions: u 0 y U ξ + u 1 ξ , y U ξ = 0 , v 1 ξ , y U ξ = 0 , u 0 y L ξ + u 1 ξ , y L ξ = 0 , v 1 ξ , y L ξ = 0 . (25)

Reconsidering the total stream function: ψ 0 y U ξ + ψ 1 ξ , y U ξ = 3 4 . ψ 0 y L ξ + ψ 1 ξ , y L ξ = 0 . (26)

Floryan (1997) used the Fourier expansion for the simplification of the above equation: ψ 1 ξ , η = ∑ n = − ∞ n = ∞ ϕ n η e i n α ξ , (27) where ϕ n = ϕ − n * represents the conditions, and the star identifies the complex conjugate.

By substituting Eq. 27 into Eq. 24, we get: ⇒ D 4 ϕ 0 + 2 α   R e   I m ∑ k = 1 k = ∞ k D 2 ϕ k * D ϕ k = 0 , (28) where D = ∂ ∂ η , I m represents the imaginary part. Integrating Eq. 28 w.r.t. ‘ η ’ we get: ⇒ D 3 ϕ 0 + 2 α   R e   I m ∑ k = 1 k = ∞ k D ϕ k * D ϕ k = A   R e , w h e r e   A = ∂ p 1 ∂ ξ m e a n (29)

The equations of the model are valid for a constant pressure gradient, where both the straight and sinusoidal channels are driven by the same gradient. The flow rate variation is taken into account and measured due to the extra flow losses that occur in the meandering channel. Additionally, it is assumed that in Eq. 26, A is equal to zero, and the correction for the flow rate, Q 1 , can be easily calculated as per Floryan’s work in 1997.

Let us assume that the wavelength is a large quantity, denoted as α → 0 . We can then obtain a solution for the problem. To regularize the flow domain, transformation θ = η − 2   S ⁡ cos ⁡ σ is applied, where σ = α ξ represents a slow scale. With these new variables, the governing equations in dimensionless form, denoted as Eqs 4, 6, are transformed into the following set of field equations given as: α ∂ u b ∂ σ + 2 S α ⁡ sin ⁡ σ ∂ u b ∂ θ + ∂ v b ∂ θ = 0 , (30) − 1 + 4 S 2 α 2 sin 2 ⁡ σ ∂ 2 u b ∂ θ 2 − α 2 ∂ 2 u b ∂ σ 2 − 4 S α 2   sin ⁡ σ ∂ 2 u b ∂ θ ∂ σ − 2 S α 2 ⁡ cos ⁡ σ ∂ u b ∂ θ + 2 S α R e ⁡ sin ⁡ σ u b ∂ u b ∂ θ + R e   v b ∂ u b ∂ θ + α   R e   u b ∂ u b ∂ σ + α   R e ∂ p b ∂ σ + 2 S α R e ⁡ sin ⁡ σ ∂ p b ∂ θ = 0 , (31) − 1 + 4 S 2 α 2 sin 2 ⁡ σ ∂ 2 v b ∂ θ 2 − α 2 ∂ 2 v b ∂ σ 2 − 4 S α 2   sin ⁡ σ ∂ 2 v b ∂ θ ∂ σ − 2 S α 2 ⁡ cos ⁡ σ ∂ v b ∂ θ + 2 S α R e ⁡ sin ⁡ σ u b ∂ v b ∂ θ + R e   v b ∂ v b ∂ θ + α   R e   u b ∂ v b ∂ σ + R e ∂ p b ∂ θ = 0 , (32)

The no-slip boundary and free-stream conditions with the total flow rate given in Eq. 25 is reduced to the following simplest form: u b = v b = 0   a t   θ = ± 1 , u b = 1   a t   θ = 0 , Q b = ∫ − 1 1 u b σ , θ d η = 4 3 . (33)

In Eqs 30–32 are three unknown quantities, i.e., u b , v b , and P b that we need to determine by using the perturbation method. These unknown quantities are expanded in terms of α by using the following series: u b = u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + O α 4 , v b = v ^ 0 + α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 + O α 4 , P b = α − 1 P ^ − 1 + P ^ 0 + α P ^ 1 + α 2 P ^ 2 + O α 3 . (34)

Substituting the values of u b , v b , and P b from Eq. 34 into Eqs 30–32, we have: α ∂ ∂ σ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + 2 S α ⁡ sin ⁡ σ ∂ ∂ θ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + ∂ ∂ θ v ^ 0 + α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 = 0 , (35) − 1 + 4 S 2 α 2 sin 2 ⁡ σ ∂ 2 ∂ θ 2 u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 − α 2 ∂ 2 ∂ σ 2 u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 − 4 S α 2 ⁡ sin ⁡ σ ∂ 2 ∂ θ ∂ σ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 − 2 S α 2 ⁡ cos ⁡ σ ∂ ∂ θ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + 2 S α R e ⁡ sin ⁡ σ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 ∂ ∂ θ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + R e v ^ 0 + α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 ∂ ∂ θ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + α   R e u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 ∂ ∂ σ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 + α   R e ∂ ∂ σ α − 1 P ^ − 1 + P ^ 0 + α P ^ 1 + α 2 P ^ 2 + 2 S α R e ⁡ sin ⁡ σ ∂ ∂ θ α − 1 P ^ − 1 + P ^ 0 + α P ^ 1 + α 2 P ^ 2 = 0 , (36) − 1 + 4 S 2 α 2 sin 2 ⁡ σ ∂ 2 ∂ θ 2 α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 − α 2 ∂ 2 ∂ σ 2 α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 − 4 S α 2   sin ⁡ σ ∂ 2 ∂ θ ∂ σ α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 − 2 S α 2 ⁡ cos ⁡ σ ∂ ∂ θ α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 + 2 S α R e ⁡ sin ⁡ σ u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 ∂ ∂ θ α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 + R e v ^ 0 + α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 ∂ ∂ θ α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 + α   R e u ^ 0 + α   u ^ 1 + α 2   u ^ 2 + α 3   u ^ 3 ∂ ∂ σ α   v ^ 1 + α 2   v ^ 2 + α 3   v ^ 3 + R e ∂ ∂ θ α − 1 P ^ − 1 + P ^ 0 + α P ^ 1 + α 2 P ^ 2 = 0 . (37)

The analysis of heat transfer in a meandering channel involves studying the heat transfer mechanisms that occur due to meanders in the channel. Heat transfer occurs in three modes: conduction, convection, and radiation. In the present study, a case of heat transfer through conduction and convection is taken into account. Conduction occurs when heat flows from regions of high temperature to regions of low temperature through a solid material. In a meandering channel, heat is conducted through the walls of the channel. The rate of heat transfer through conduction is proportional to the temperature gradient, the thermal conductivity of the material, and the cross-sectional area of the channel. Convection occurs when heat is transferred by the movement of fluids, either liquids or gases. In a meandering channel, heat is transferred through convection due to the flow of fluid through the channel. The rate of heat transfer through convection is proportional to the temperature difference between the fluid and the channel wall, the velocity of the fluid, and the heat transfer coefficient of the fluid. The heat transfer is analyzed by using the energy equation, which relates the rate of heat transfer to the temperature distribution in the channel. The energy equation is taken into account to study the conduction and convection mechanisms of heat transfer and the thermal properties of the meandering channel in the flow of fluids.

The problem is modeled using the Navier–Stokes equations, along with the continuity and energy equations. To obtain a non-linear, simple partial differential equation, a defined set of transformations is applied. The regular perturbation technique is utilized to expand the heat transfer, T b , in terms of " α " in the meandering channel.

The energy equation is: ρ c p D T D t = k ∇ 2 T + Φ , (38) where T is the temperature distribution, ρ is density, c p is specific heat, k is thermal conductivity, and Φ is the dissipation term. The boundaries condition of the problem is at the upper wall, T = T 1 , and at the lower wall , T = T 0 .

Consider the following dimensionless variables in order to transform the energy equation into dimensionless form: T * = T − T 0 T 1 − T 0 , Φ * = L 2 Φ U 2 μ , (39)

The dimensionless numbers that appear in Eq. 38 are the Reynolds and Prandtl numbers, which are defined as: R e = ρ U L μ , Pr ⁡ = μ c p k , R e . Pr ⁡ = ρ U L c p k , 1 R e P r = k ρ U L c p , (40) u * ∂ T * ∂ x * + v * ∂ T * ∂ y * = 1 R e P r ∂ 2 T * ∂ x * 2 + ∂ 2 T * ∂ y * 2 + E c R e Φ * , (41)

Equation 41 is the dimensionless form of the energy equation. On the other hand, the boundaries conditions for the temperature distribution in the dimensionless form are obtained, and the dimensional boundary conditions are given by taking: T * x , y = T − T 0 T 1 − T 0 , (42)

The temperature at the upper wall is T = T 1 , so using definition from Eq. 42: y = 1 , T = T 1 , T * x , y = T − T 0 T 1 − T 0 = T 1 − T 0 T 1 − T 0 , T * x , y = 1 , (43)

The temperature at the lower wall is T = T 0 , so using definition from Eq. 42: y = − 1 , T = T 0 , T * x , y = T − T 0 T 1 − T 0 = T 0 − T 0 T 1 − T 0 , T * x , y = 0 , (44)

The temperature of a flow in a straight channel is denoted by T ^ 0 , and so we have: u 0 ∂ T ^ 0 ∂ x + v 0 ∂ T ^ 0 ∂ y = 1 R e P r ∂ 2 T ^ 0 ∂ x 2 + ∂ 2 T ^ 0 ∂ y 2 + E c R e 2 ∂ u 0 ∂ x 2 + 2 ∂ v 0 ∂ y 2 + ∂ v 0 ∂ x + ∂ u 0 ∂ y 2 , (45)

The boundaries condition for the temperature at the upper wall, T ^ 0 = 1 , and at the lower wall, T ^ 0 = 0 .

Putting Eq. 11 into Eq. 45, after simplification we have: ∂ 2 T ^ 0 ∂ y 2 = − 4 P r E c y 2 , (46)

Integrating Eq. 46 twice w.r.t ” y ” , we get: T ^ 0 y = − P r E c y 4 3 + C 1 y + C 2 , (47) where C 1 and C 2 are the constant of integration, which can be determined by using the boundary conditions imposed at the upper wall of the channel, T ^ 0 1 = 1 , and at lower wall of the channel, T ^ 0 − 1 = 0 . T ^ 0 y = P r E c 3 1 − y 4 + 1 2 1 + y , (48)

As the V b X = V 0 X + V 1 X and T b X = T 0 X + T 1 X , where V b X = u b i ^ + v b j ^ and T b X are the velocity and temperature of the total flow, respectively, and V 1 X = u 1 i ^ + v 1 j ^ and T 1 X are the velocity and temperature induced by the channel geometry, respectively.

It is assumed that α is a small quantity, and for this choice of α we expand the unknown quantities of T ^ b by using the regular perturbations to find T ^ b . The dimensionless form of the energy equation is transformed and expanded using the following transformation and expansion: η = y − 2 s c a s ξ , ξ = α x , (49) T ^ b = T ^ 0 + α   T ^ 1 + α 2   T ^ 2 + O α 3 , (50)

We get: u 0 2 S α sin ξ ∂ T ^ 0 ∂ η + α ∂ T ^ 1 ∂ η + α 2 ∂ T ^ 2 ∂ η + α ∂ T ^ 0 ∂ ξ + α ∂ T ^ 1 ∂ ξ + α 2 ∂ T ^ 2 ∂ ξ + α 2 u 2 2 S α sin ξ ∂ T ^ 0 ∂ η + α ∂ T ^ 1 ∂ η + α 2 ∂ T ^ 2 ∂ η + α ∂ T ^ 0 ∂ ξ + α ∂ T ^ 1 ∂ ξ + α 2 ∂ T ^ 2 ∂ ξ + α v 1 ∂ T ^ 0 ∂ η + α ∂ T ^ 1 ∂ η + α 2 ∂ T ^ 2 ∂ η = 1 R e P r [ 4 S 2 α 2 sin 2 ξ ∂ 2 T ^ 0 ∂ η 2 + α ∂ 2 T ^ 1 ∂ η 2 + α 2 ∂ 2 T ^ 2 ∂ η 2 + α 2 ∂ 2 T ^ 0 ∂ ξ 2 + α ∂ 2 T ^ 1 ∂ ξ 2 1 + α 2 ∂ 2 T ^ 2 ∂ ξ 2 + 4 α 2 S ⁡ sin ξ ∂ 2 T ^ 0 ∂ η ∂ ξ + α ∂ 2 T ^ 1 ∂ η ∂ ξ + α 2 ∂ 2 T ^ 2 ∂ η ∂ ξ + 2 S α 2 ⁡ cos ⁡ ⁡ ξ ∂ T ^ 0 ∂ η + α ∂ T ^ 1 ∂ η + α 2 ∂ T ^ 2 ∂ η + ∂ 2 T ^ 0 ∂ η 2 + α ∂ 2 T ^ 1 ∂ η 2 + α 2 ∂ 2 T ^ 2 ∂ η 2 ] + E c R e [ ∂ u 0 ∂ η 2 + α 2 ( 8 S 2 sin 2 ξ ∂ u 0 ∂ η 2 + 2 ∂ u 0 ∂ η ∂ u 2 ∂ η + 4 S s i n ξ ∂ u 0 ∂ η ∂ v 1 ∂ η + 2 ∂ v 1 ∂ η 2 + 8 S s i n ξ ∂ u 0 ∂ η ∂ u 0 ∂ ξ + 2 ∂ u 0 ∂ ξ + 2 ∂ u 0 ∂ η ∂ v 1 ∂ ξ ) + O α 3 ] , (51)

The boundaries condition on the upper wall is y = h = 1 , and the boundaries condition on the lower wall is y = − h = − 1 for the small α . T ^ 0 1 = 1 , T ^ 1 1 = 0 , T ^ 2 1 = 0 , (52) T ^ 0 − 1 = 0 , T ^ 1 − 1 = 0 , T ^ 2 − 1 = 0 , (53)

With the solution represented as a power expansion in terms of σ and by equating the coefficients, we have: u ^ 0 θ = 1 − θ 2 , u ^ 1 σ , θ = 0 , u ^ 2 σ , θ = 4 3 S θ ⁡ cos ⁡ σ 1 − θ 2 , u ^ 3 σ , θ = 2 315 R e S   θ ⁡ sin ⁡ σ θ 2 − 1 θ 4 + 21 S θ θ 2 − 4 cos ⁡ σ + θ 2 − 34 , v ^ 0 σ , θ = 0 , v ^ 1 = − 2 S ⁡ sin ⁡ σ 1 − θ 2 , v ^ 2 σ , θ = 0 , v ^ 3 σ , θ = 1 3 S θ 2 − 1 θ 2 − 1 + 8 S θ ⁡ cos ⁡ σ sin ⁡ σ , P ^ − 1 σ = − 2 σ R e + C 0 , P ^ 0 = C 1 , P ^ 1 σ , θ = − 4 R e − 1 S 2 S σ − θ ⁡ sin ⁡ σ − S ⁡ sin ⁡ 2 σ , P ^ 2 σ , θ = 2 S ⁡ cos ⁡ σ θ − 2 θ 3 3 + θ 5 5 + 11 15 S 2 ⁡ cos ⁡ 2 σ , T ^ 0 η = P r E c 3 1 − η 4 + 1 2 1 + η , T ^ 1 η = 1 1260 sin ξ P r R e − 1 + S − 1 + η 2 − 105 − 5 + η 2 + 8 E c P r η − 11 − 11 η 2 + 10 η 4 ,

And T ^ 2 = 1 1260 [ 105 2 pr 2 R e 2 − 1 + S cos ξ − η 4 + 1 + 1176 E c P r S − η 5 + η ) cos ξ + 4 E c Pr 3 ⁡ R e 2 − 1 + S η 7 − η cos ξ + 10 9 E c Pr 3 ⁡ R e 2 − 1 + S − η 9 + η cos ξ + 88 E c Pr 3 ⁡ R e 2 − 1 + S 2 η 2 − 1 sin 2 ξ + 420 Pr 2 ⁡ R e 2 − 1 + S 2 η 3 − η sin 2 ξ + 4 3 E c P r 11 Pr 2 ⁡ R e 2 1 + S 2 + 3780 S 2 − η 4 + 1 sin 2 ξ + 56 E c Pr 3 ⁡ R e 2 − 1 + S 2 − η 6 + 1 sin 2 ξ + 20 Pr 2 ⁡ R e 2 − 1 + S 2 η 7 − η sin 2 ξ + 50 E c Pr 3 ⁡ R e 2 − 1 + S 2 η 8 − 1 sin 2 ξ + 112 9 E c Pr 3 ⁡ R e 2 − 1 + S 2 − η 10 + 1 sin 2 ξ + 84 P r R e − 1 + S − η 5 + η sin ξ − S cos ξ + 2 P r R e − 1 + S sin ξ + 7 2 Pr 2 ⁡ R e − 1 + S η 6 − 1 cos ξ R e − 32 E c S sin ξ + 4 3 Pr − η 3 + η cos ξ E c − 11 Pr 2 ⁡ R e 2 − 1 + S − 2100 S +    630 R e − 1 + S S sin ξ +    1 2 − η 2 + 1 cos ξ 105 − 5 Pr 2 ⁡ R e 2 − 1 + S + 24 S − 352 E c Pr 2 ⁡ R e − 1 + S S sin ξ +    20 E c Pr 2 ⁡ R e − 1 + S S − η 8 + 1 sin 2 ξ ] .

Stream function ψ σ , θ

The stream function ψ σ , θ is used to define the velocity components such that: u b = ∂ ψ σ , θ ∂ θ , v b = − 2 α   S ⁡ sin ⁡ σ   ∂ ψ σ , θ ∂ θ + ∂ ψ σ , θ ∂ σ , (54)

By integrating u b with ‘ θ ’ and after finding the constant of integration, we get: ψ σ , θ = 1 1260 ( 420 S α 2 − 1 + θ 2 2 ⁡ cos   σ + r e S α 3 − 1 + θ 2 2 − 67 + 2 θ 2 + θ 4   sin ⁡ σ + 4 1 + θ 2 − 105 − 2 + θ +    r e S 2 α 3 10 − 20 θ + 30 θ 2 − 12 θ 3 − 6 θ 4 + 3 θ 5 sin ⁡ 2 σ ) . (55)