Effects of topographical wall patterns on flow through porous media

The schematic illustrates a periodic grooved microchannel used for Darcy–Brinkman flow analysis. The flow direction is along the z axis (longitudinal flow). L, a, and h are defined as pattern wavelength, amplitude, and mean channel height, respectively. The top wall is considered flat for both approaches. (a) In the lubrication theory, the lower boundary is defined by an arbitrary shape, represented as $y = a h (x)$ ⁠. (b) In the spectral method, the lower boundary follows a sinusoidal profile, expressed as $y = a cos (\frac{2 π x}{L})$ ⁠. In both approaches, the analysis assumes one periodic cell.

III. LUBRICATION APPROXIMATION

The lubrication theory is an approximation of the Navier–Stokes and continuity equations at low Reynolds numbers for narrow geometries with slow changes in curvature.^57–59 We assumed that the flow occurs along the z axis, while velocities in the x and y directions are negligible. Therefore, the function w depends on x and y. We introduce the non-dimensional variables as follows:

x = L X

⁠,

y = h Y

⁠, and

w = \frac{- p_{z} h L W}{μ_{eff}}

⁠. Accordingly, the governing equations in terms of the non-dimensional variables can be expressed as

δ^{2} \frac{\partial^{2} W_{l}}{\partial X^{2}} + \frac{\partial^{2} W_{l}}{\partial Y^{2}} - κ_{l}^{2} W_{l} = - 1 .

(3)

The aspect ratio of the channel is denoted as

δ = h / L

⁠. The non-dimensional parameter

κ_{l}

=

h \sqrt{\frac{μ}{μ_{eff} K}}

appears in Eq. and that characterizes the porous medium within lubrication framework. The Stokes and Darcian flow limits are given by

κ_{l} \to 0

and

κ_{l} ≫ 1

⁠, respectively. Boundary conditions at both walls are

W_{l} (X, Y = 1) = 0,

(4a)

W_{l} (X, Y = ω H (X)) = 0,

(4b)

where

ω = \frac{a}{h}

is a ratio between the pattern amplitude and the mean channel height of the pattern. Here,

H (X)

is a normalized shape function. Section III A applies the perturbation method to calculate the flow quantities.

A. Perturbation approximation and solution of the hydrodynamic problem

A perturbation method can be employed to obtain an approximate solution for problems involving a wavy channel. The objective of this section is to improve the solution by adding extra terms in the perturbation solution. Correspondingly, we assume that the ratio

δ = h / L = 1 / λ

is small or

δ ≪ 1

⁠, and this allows for the use of

δ

as a perturbation parameter. Therefore, the following expression is used for velocity,

W_{l} = W_{l 0} + δ^{2} W_{l 2} + δ^{4} W_{l 4} + O (δ^{6}) .

(5)

Before proceeding further, the non-dimensional governing equation and boundary conditions need to be determined at different orders of

δ

⁠. Therefore, substituting the above-mentioned expansion in Eq. and collecting terms in powers of

δ

⁠, the coefficients of

δ^{2}

and

δ^{4}

are set to zero. The governing equations and boundary conditions of the problem at leading and second orders will be discussed in Secs. III A 1 and III A 2.

1. The leading order problem: $O (δ^{0})$

The first step involves a standard discussion found in textbook.⁵⁸ Since the problem statement only incorporates

δ^{2}

⁠, which is assumed to be small, the analysis at the leading order

O (δ^{0})

reduces to the well-known classical lubrication theory (CLT) problem. In this context, the primary flow characteristics are governed by the principles of standard lubrication theory. Consequently, the governing equation becomes

\frac{\partial^{2} W_{l 0}}{\partial Y^{2}} - κ_{l}^{2} W_{l 0} = - 1 .

(6)

Boundary conditions are

W_{l 0} (X, Y = 1) = 0

and

W_{l 0} (X, Y = ω H (X)) = 0

⁠. Utilizing these two boundary conditions and solving Eq. for an arbitrary surface,

\begin{matrix} W_{l 0 a r b} = \frac{1}{κ_{l}^{2}} - \frac{1}{κ_{l}^{2}} (\frac{sinh (κ_{l} H (X) + κ_{l} (- Y) + κ_{l}) + cosh (κ_{l} H (X) + κ_{l} (- Y) + κ_{l})}{sinh (κ_{l} H (X)) + cosh (κ_{l} H (X)) + sinh (κ_{l}) + cosh (κ_{l})}) \\ + \frac{1}{κ_{l}^{2}} (\frac{- sinh (κ_{l} Y) - cosh (κ_{l} Y)}{sinh (κ_{l} H (X)) + cosh (κ_{l} H (X)) + sinh (κ_{l}) + cosh (κ_{l})}) . \end{matrix}

(7)

Equation gives the velocity profile for the classical lubrication theory (CLT) for any complex geometries. Additionally, we also focused on a periodically grooved sinusoidal surface, and the shape function is assumed

H (X) = cos (2 π X)

⁠, which yields

\begin{matrix} W_{l 0} = \frac{1}{κ_{l}^{2}} - \frac{sinh (κ_{l} Y) + cosh (κ_{l} Y)}{κ_{l}^{2} (sinh (ω κ_{l} cos (2 π X)) + cosh (ω κ_{l} cos (2 π X)) + sinh (κ_{l}) + cosh (κ_{l}))} \\ - \frac{sinh (ω κ_{l} cos (2 π X) + κ_{l} (- Y) + κ_{l}) + cosh (ω κ_{l} cos (2 π X) + κ_{l} (- Y) + κ_{l})}{κ_{l}^{2} (sinh (ω κ_{l} cos (2 π X)) + cosh (ω κ_{l} cos (2 π X)) + sinh (κ_{l}) + cosh (κ_{l}))} . \end{matrix}

(8)

Equations and are independent of the perturbation parameter (⁠

δ

⁠). It depends on

ω

⁠,

κ_{l}

⁠, X, and Y.

2. The second order problem: $O (δ^{2})$

In the past, many researchers have employed the lubrication approach to examine the effects of micropatterning, typically truncating the analysis at the leading-order solution [

O (δ^{0})

]. The objective of this section is to refine the solution of CLT by incorporating additional terms into the perturbation solution, with the next order being

O (δ^{2})

⁠. This solution is considered an extended lubrication theory (ELT). At the next order,

O (δ^{2}

⁠), the governing equation and boundary condition become

\frac{\partial^{2} W_{l 0}}{\partial X^{2}} + \frac{\partial^{2} W_{l 2}}{\partial Y^{2}} - κ_{l}^{2} W_{l 2} = 0,

(9)

W_{l 2} (X, Y = 1) = 0,

(10a)

W_{l 2} (X, Y = ω H (X)) = 0 .

(10b)

Here, we determine

W_{l 0 X X}

using Eq. , as outlined in Appendix B.

O (δ^{2})

-solution for complex surfaces is provided by

W_{l 2 a r b} = - \frac{(T_{11}) (T_{22}) (H^{″} (X) + κ_{l} H^{'} {(X)}^{2} tanh (\frac{1}{2} (κ_{l} - κ_{l} H (X))))}{4 {cosh}^{2} (\frac{1}{2} (κ_{l} - κ_{l} H (X))) κ_{l}^{2} (- sinh (2 κ_{l} H (X)) - cosh (2 κ_{l} H (X)) + χ_{1})} .

(11)

If a shape function is

H (X) = cos (2 π X)

⁠,

W_{l 2} = \frac{(κ_{l} ω (cos (4 π X) - 1)) (T_{1} + 2 cos (2 π X))}{κ_{l}^{2} T_{8} + T_{9}} .

(12)

The expressions of

χ_{1}

⁠,

T_{11}

⁠,

T_{22}

⁠,

T_{8}

⁠, and

T_{9}

are given in Appendix B. The total flow-field is obtained by

W_{l} = W_{l 0} + δ^{2} W_{l 2} + O (δ^{4})

⁠, where

W_{l 0}

and

W_{l 2}

can be obtained from Eqs. and , respectively.

B. Calculation for scaled flow rate

Equation is utilized to evaluate the flow rate for a single periodic cell of the channel. Based on the current scaling parameters, X varies from

- 1 / 2

to

1 / 2

⁠, while Y ranges from

ω H (X)

to 1. By solving Eq. along the X and Y axes, we obtain the dimensionless flow rate for one periodic cell of the channel, denoted as

Q_{l}

⁠,

Q_{l} = \int_{X = - 1 / 2}^{1 / 2} \int_{Y = ω H (X)}^{1} W_{l} d Y d X .

(13)

By substituting Eq. in Eq. and separating the velocities according to different orders of

δ

⁠, the flow rate can be calculated up to

O (δ^{2})

for

H (X)

=

cos (2 π X)

⁠, as follows:

Q_{l} = Q_{l 0} + δ^{2} Q_{l 2} + δ^{4} Q_{l 4} + O (δ^{6}) .

(14)

The

O (δ^{0})

solution is

Q_{l 0} = \frac{- 4 sinh (\frac{κ_{l}}{2}) + 3 κ_{l} cosh (\frac{κ_{l}}{2}) + κ_{l} cosh (\frac{3 κ_{l}}{2}) - 4 sinh (\frac{κ_{l}}{2}) cosh (κ_{l})}{4 κ_{l}^{3} {cosh}^{3} (\frac{κ_{l}}{2})} + \frac{ω^{2} tanh (\frac{κ_{l}}{2})}{4 κ_{l} {cosh}^{2} (\frac{κ_{l}}{2})} .

(15)

If the bottom wall is smooth or straight (with no topographical patterns,

ω = 0

⁠), then the dimensional flow rate (⁠

Q_{l s}

⁠) is

Q_{l s} = \frac{(- 4 sinh (\frac{κ_{l}}{2}) + 3 κ_{l} cosh (\frac{κ_{l}}{2}) + κ_{l} cosh (\frac{3 κ_{l}}{2}) - 4 sinh (\frac{κ_{l}}{2}) cosh (κ_{l}))}{4 κ_{l}^{3} {cosh}^{3} (\frac{κ_{l}}{2})} .

(16)

The

O (δ^{2})

solution is

\begin{matrix} Q_{l 2} = - \frac{π^{2} (- 4 κ_{l} sinh (2 κ_{l}) + sinh (4 κ_{l}) - 4 κ_{l} cosh (2 κ_{l}) + cosh (4 κ_{l}) - 1)}{κ_{l}^{3} {(1 + sinh (κ_{l}) + cosh (κ_{l}))}^{4}} ω^{2} + ω^{4} \frac{π^{2} (sinh (κ_{l}) + cosh (κ_{l}))}{2 κ_{l} {(sinh (κ_{l}) + cosh (κ_{l}) + 1)}^{6}} (R_{1}) . \end{matrix}

(17)

The expression of

R_{1}

is defined in Appendix B. In order to study confined flow, a key parameter is scaled flow rate or effective hydraulic permeability, which is denoted by

Q_{l} / Q_{l s}

⁠. Here,

Q_{l}

denotes the dimensionless flow rate for a channel with a wavy bottom wall, whereas

Q_{l s}

denotes the dimensionless flow rate for a straight channel of identical width and average thickness (h). The scaled flow rates for both classical lubrication theory (CLT) and extended lubrication theory (ELT) are given by

{(\frac{Q_{l}}{Q_{l s}})}_{CLT} = \frac{3 κ_{l} - 4 sinh (κ_{l}) - 2 sinh (2 κ_{l}) + 4 κ_{l} cosh (κ_{l}) + κ_{l} cosh (2 κ_{l}) + κ_{l}^{2} ω^{2} sinh (κ_{l})}{8 {cosh}^{4} (\frac{κ_{l}}{2}) (κ_{l} - 2 tanh (\frac{κ_{l}}{2}))},

(18)

\begin{matrix} {(\frac{Q_{l}}{Q_{l s}})}_{ELT} = \frac{3 κ_{l} - 4 sinh (κ_{l}) - 2 sinh (2 κ_{l}) + 4 κ_{l} cosh (κ_{l}) + κ_{l} cosh (2 κ_{l}) + κ_{l}^{2} ω^{2} sinh (κ_{l})}{8 {cosh}^{4} (\frac{κ_{l}}{2}) (κ_{l} - 2 tanh (\frac{κ_{l}}{2}))} + δ^{2} \frac{2 π^{2} κ_{l} ω^{2} - π^{2} ω^{2} sinh (2 κ_{l})}{8 {cosh}^{4} (\frac{κ_{l}}{2}) (κ_{l} - 2 tanh (\frac{κ_{l}}{2}))} + O (δ^{4}) . \end{matrix}

(19)

The final integration for scaled flow rate is accomplished using MATHEMATICA software. To assess the impact of a wavy bottom wall on a channel's hydraulic conductance relative to a straight channel, we use the dimensionless flow rate ratio

Q_{l} / Q_{l s}

in the present paper.

IV. SPECTRAL APPROACH

The objective of this section is to develop a model for flow parameters that are accurate for significantly large dependent variables (permeability parameters, wavelengths, and amplitudes). Therefore, the spectral method is applied to investigate the impact of micropatterning on porous media flow as illustrated in panel (b) of Fig. 1. According to the current approach, we introduce the non-dimensional variables as follows:

x = \frac{L}{2 π} X

⁠,

y = \frac{L}{2 π} Y

⁠, and

w = \frac{- p_{z} h L W}{2 π μ_{eff}}

⁠. Note that these dimensionless parameters differ from those used in the lubrication approach (Sec. III). By applying this non-dimensionalization scheme to the governing equation and boundary conditions, i.e., Eqs. , , and , the revised dimensionless governing equation and boundary conditions are as follows:

\frac{\partial^{2} W}{\partial X^{2}} + \frac{\partial^{2} W}{\partial Y^{2}} - κ_{f}^{2} W = - \frac{1}{H},

(20)

W (X, Y = H) = 0,

(21a)

W (X, Y = ε cos (X)) = 0 .

(21b)

Here,

H = \frac{2 π h}{L} = \frac{2 π}{λ}

⁠, where

λ = L / h

represents the relationship between pattern wavelength and the mean channel height. The parameter

ε = \frac{2 π a}{L}

defines the ratio of pattern amplitude to wavelength. Meanwhile,

κ_{f}

=

\frac{L}{2 π} \sqrt{\frac{μ}{μ_{e} K}}

characterizes the porous media in the governing Eq. .

A. Spectral method

The velocity field is assumed in the form of a Fourier series. Due to the channel's periodic nature, the analysis is conducted within a single periodic cell to capture the repeating flow dynamics throughout the channel,

W (X, Y) = F_{0} (Y) + \sum_{n = 1}^{\infty} C_{n} F_{n} (Y) cos (n X) .

(22)

Substituting Eq. in Eq. and collecting the coefficients of

cos (n X)

(⁠

n \geq 0

⁠), we equate them to zero to derive ordinary differential equations for

F_{0} (Y)

and

F_{n} (Y)

⁠. No-slip boundary condition at

Y = H

given by Eq. necessitates that

F_{0} (H) = 0

and

F_{n} (H) = 0

⁠. The resulting equations are

F_{0} (Y) = \frac{1}{H κ_{f}^{2}} - \frac{e^{- H κ_{f}}}{H κ_{f}^{2}} e^{κ_{f} Y} + C_{0} (e^{- κ_{f} Y} - e^{- 2 H κ_{f}} e^{κ_{f} Y}), F_{n} (Y) = e^{- η Y} - e^{- 2 η H} e^{η Y},

(23)

where

η = \sqrt{n^{2} + κ_{f}^{2}}

⁠. Upon substituting Eq. in Eq. , the velocity profile can be expressed as follows:

W = \frac{1}{H κ_{f}^{2}} - \frac{e^{- H κ_{f}}}{H κ_{f}^{2}} e^{κ_{f} Y} + C_{0} (e^{- κ_{f} Y} - e^{- 2 H κ_{f}} e^{κ_{f} Y}) + \sum_{n = 1}^{\infty} C_{n} (e^{- η Y} - e^{- 2 η H} e^{η Y}) cos (n X) .

(24)

Equation is required to adhere to the no-slip boundary condition on the undulating surface [Eq. ], leading to the following form after trigonometric simplification:

\begin{matrix} \frac{1}{H κ_{f}^{2}} - \frac{e^{- H κ_{f}}}{H κ_{f}^{2}} e^{ε κ_{f} cos (X)} + C_{0} (e^{- κ_{f} ε cos (X)} - e^{- 2 H κ_{f}} e^{κ_{f} ε cos (X)}) \\ + \sum_{n = 1}^{\infty} C_{n} (e^{- η ε cos (X)} - e^{- 2 η H} e^{η ε cos (X)}) cos (n X) = 0. \end{matrix}

(25)

Equation is integrated between

- π

and

π

⁠. Using the property of a modified Bessel function that

I_{n} (z) = \frac{1}{π} \int_{0}^{π} e^{z cos (θ)} cos (n θ) d θ

⁠, the resultant equation after integration can be rewritten as

\begin{matrix} C_{0} I_{0} (ε κ_{f}) (1 - e^{- 2 H κ_{f}}) + \sum_{n = 1}^{\infty} C_{n} I_{n} (η ε) ({(- 1)}^{n} - e^{- 2 η H}) \\ = - \frac{1}{H κ_{f}^{2}} + \frac{e^{- H κ_{f}}}{κ_{f}^{2}} I_{0} (ε κ_{f}) . \end{matrix}

(26)

In the above-mentioned equation,

cos (m X)

is multiplied on both sides. Here, m is an element of the set of non-negative integers that will provide a different set of equations with reference to m (⁠

m > 0

⁠). The resultant expressions are integrated from

- π

and

π

[Eq. ]. Equation is considered for

m = 0

⁠, and for

m > 0

⁠, equation is outline as follows:

\begin{matrix} C_{0} I_{m} (ε κ_{f}) ({(- 1)}^{m} - e^{- 2 H κ_{f}}) + \sum_{n = 1}^{\infty} \frac{C_{n}}{2} (I_{m + n} (η ε) ({(- 1)}^{m + n} - e^{- 2 η H}) + I_{m - n} (η ε) \\ ({(- 1)}^{m - n} - e^{- 2 η H})) = \frac{e^{- H κ_{f}}}{H κ_{f}^{2}} I_{m} (ε κ_{f}) . \end{matrix}

(27)

By truncating each summation occurring in Eqs. and to n terms, we obtained a numerical solution for the equation system governing the

n + 1

constants

{C_{02}, C_{1 n}, C_{2 n}}

⁠. We achieved this by truncating the number of equations to

n + 1

and performing matrix inversion.

B. Calculation for scaled flow rate

Equation is used to evaluate the dimensionless flow rate (Q) for one periodic cell of the channel. As per the current scaling parameter, X varies from

- π

to

π

and Y varies from

ε cos (X)

to H,

Q = \int_{X = - π}^{π} \int_{Y = ε cos (X)}^{H} W d YdX,

(28)

which evaluates to

\begin{matrix} Q = (2 π) (\frac{1}{κ_{f}^{2}} - \frac{1}{H κ_{f}^{3}} (1 - e^{- H κ_{f}} I_{0} (ε κ_{f})) - \frac{C_{0}}{κ_{f}} (2 e^{- H κ_{f}} - I_{0} (ε κ_{f}) (1 + e^{- 2 H κ_{f}}))) \\ + (2 π) \sum_{n = 1}^{\infty} C_{n} (\frac{I_{n} (ε η)}{η}) ({(- 1)}^{n} + e^{(- 2 η H)}) . \end{matrix}

(29)

When the pattern amplitude is zero (⁠

ε

= 0), the dimensionless flow rate in a straight channel is given by

Q^{*} = \frac{2 π H κ_{f} - 4 π tanh (\frac{H κ_{f}}{2})}{H κ_{f}^{3}} .

(30)

The ratio of the dimensionless flow rate in the patterned channel to that of the straight channel is denoted as

\frac{Q}{Q^{*}}

and is referred to as the scaled flow rate,

\begin{matrix} \frac{Q}{Q^{*}} = \frac{(\frac{1}{κ_{f}^{2}} - \frac{1}{H κ_{f}^{3}} (1 - e^{- H κ_{f}} I_{0} (ε κ_{f})) - \frac{C_{0}}{κ_{f}} (2 e^{- H κ_{f}} - I_{0} (ε κ_{f}) (1 + e^{- 2 H κ_{f}})))}{\frac{H κ_{f} - 2 tanh (\frac{H κ_{f}}{2})}{H κ_{f}^{3}}} + (\frac{1}{\frac{H κ_{f} - 2 tanh (\frac{H κ_{f}}{2})}{H κ_{f}^{3}}}) \sum_{n = 1}^{\infty} C_{n} (\frac{I_{n} (ε η)}{η}) ({(- 1)}^{n} + e^{(- 2 η H)}) . \end{matrix}

(31)

Equation shows that the scaled flow rate depends on the dimensionless pattern amplitude (⁠

ε

⁠), wavelength (⁠

λ

⁠), and porous parameter (⁠

κ_{f}

⁠). Notably, the expression for the scaled permeability of a channel does not impose any limitations on the values of

ε

⁠,

λ

⁠, and

κ_{f}

⁠. However, to accurately determine the coefficients

C_{n}

⁠, Eq. must be solved numerically. In certain cases, closed-form asymptotic approximations for the coefficients, as given by Eqs. and , can be derived. These cases will be discussed in Secs. V A and V B.

V. NOVEL CONTRIBUTION OF THE SPECTRAL METHOD

The present study explored several novel features of the spectral method, including second-order solutions [ $O (ε^{2})$ ] for the DB flow⁸ and Stokes flow limit.⁵⁵ Equations (26), (27), and (31) are used to calculate an asymptotic solution for the scaled flow rate.

A. Spectral equations are compared with the literature model: DB model

In the past, Ng and Wang⁸ employed the boundary perturbation method to study the DB flow through a sinusoidal grooved wall and developed a mathematical model for flow rate, which is an

O (ε^{2})

-accurate. Here, we assumed the top wall is flat. An alternate form of representing the pattern amplitude is

α

⁠, which is utilized in the literature model.

α

and

ε

are defined as

α

=

a / h

=

ε / H

⁠. A key novel of the spectral method is that we can obtain the literature model using the two terms (⁠

n = 2

⁠) in the Fourier series. Equation is essential for calculating the dimensionless flow rate and also aids in determining the flow rate through a straight channel (⁠

ε = 0

⁠). When we truncate Eqs. and to a

2 \times 2

form, the resulting set of equations provides an

O (α^{2})

accurate solution for the scaled flow rate. Similarly, truncating the series to include additional terms yields an accuracy of

O (α^{n})

⁠. For a

2 \times 2

truncation, the resulting constants are as follows:

C_{0} = - \frac{χ_{2}}{H κ_{f}^{2} (χ_{2} + 1)} + α^{2} \frac{H η_{2} (\frac{κ_{f}}{sinh (\frac{H κ_{f}}{2})} - \frac{2 η_{1}}{tanh (H η_{1})} \frac{1}{cosh (\frac{H κ_{f}}{2})})}{16 κ_{f}} + O (α^{4}) .

(32)

The expressions of

χ_{2}

and

η_{2}

are defined in Appendix B. Here,

η_{1} = \sqrt{κ_{f}^{2} + 1}

and

ε

=

α \times H

⁠,

C_{1} = - \frac{tanh (\frac{H κ_{f}}{2}) (coth (H \sqrt{κ_{f}^{2} + 1}) + 1)}{2 κ_{f}} α + O (α^{3}) .

(33)

Equations and signify the order of perturbation parameter (⁠

α

⁠). For example, the first term in Eq. is an order of

O (α^{0})

⁠, whereas the second term is an order of

O (α^{2})

⁠. The next term will contribute up to

O (α^{4})

⁠. For

C_{1}

⁠, the first term contribution in a series is on

O (ε)

⁠, and the next term will be

O (α^{3})

⁠. The final expression for the scaled flow rate is

\frac{Q}{Q^{*}} = 1 + \frac{(κ_{f}^{2} tanh (\frac{H κ_{f}}{2}) - κ_{f} \sqrt{κ_{f}^{2} + 1} {tanh}^{2} (\frac{H κ_{f}}{2}) coth (H \sqrt{κ_{f}^{2} + 1})) H^{2}}{2 H κ_{f} - 4 tanh (\frac{H κ_{f}}{2})} α^{2} + O (α^{4}) .

(34)

B. The spectral equations as compared with the previous domain perturbation models: Stokes flow

The second key novel contribution of the spectral method is to obtain the scaled flow rate for the Stokes limit (⁠

κ_{f} \to 0

⁠). Stroock et al.⁵⁵ employed the boundary perturbation method to develop a mathematical model for flow rate at Stokes limit, which was accurate at a small amplitude-to-pitch ratio. The extension of this work was carried out by Dewangan and Datta⁴⁷ using the spectral method, which is accurate at a significant amplitude-to-pitch ratio. To obtain the literature model, two terms are considered in Eq. at

κ_{f} \to 0

⁠. The final expression of the scaled flow rate is

\frac{Q}{Q^{*}} = 1 + \frac{3 α^{2}}{2} (2 - H coth H) + O (α^{4}) .

(35)

Equation is similar to Eq. of the literature model⁵⁵ as per the current notation.

VI. RELATION BETWEEN CHARACTERIZING POROUS MEDIUM BY LUBRICATION THEORY AND SPECTRAL THEORY

The present study employs two methodologies to investigate the effect of micropatterning on DB flow. The first approach utilizes lubrication theory (discussed in Sec. III), while the second employs the spectral method (discussed in Sec. IV). Both theories involve distinct scaling parameters. The height of the channel keenly influences the dimensionless parameters in the lubrication theory, whereas, in the spectral method, the wavelength of the patterned surface plays a critical role.

κ_{l} = h \sqrt{\frac{μ}{μ_{e} K}}

is a parameter characterizing the porous medium for the lubrication method, whereas, as per the spectral model,

κ_{f}

is equal to

\frac{L}{2 π} \sqrt{\frac{μ}{μ_{e} K}}

⁠. In this section, we demonstrated a relation for the permeability of porous media between both mathematical approaches, which is defined as

κ = \frac{κ_{l}}{κ_{f}} = \frac{h \sqrt{\frac{μ}{μ_{e} K}}}{\frac{L}{2 π} \sqrt{\frac{μ}{μ_{e} K}}} = \frac{2 π h}{L} = \frac{2 π}{λ} = 2 π δ .

(36)

In Sec. II,

δ

=

\frac{1}{λ}

is defined as a perturbation parameter (⁠

δ

⁠) as per the lubrication theory. The key findings of the present mathematical modeling are demonstrated graphically in Sec. VII.

VII. RESULTS AND DISCUSSION

The present study examines fluid flow through a channel with one grooved and one flat, impermeable wall, with porous material present in the channel. The fluid movement in the porous area is governed by the DB equation, considering only viscous forces. Lubrication equations are thoroughly validated by comparison with numerical simulations (NS) of the DB equation, as discussed in Sec. VII A. A finite-element-based numerical study is also conducted using COMSOL Multiphysics^® software, as discussed in Appendix A. The present NS is validated with the existing literature models^8,55 and the current theoretical models (lubrication approach and grid-free spectral approach). In this paper, numerical values are represented symbolically. This section effectively illustrates the influence of dependent variables on scaled flow rate through graphical representation.

The variation of the scaled flow rate ratio is shown against the non-dimensional pattern amplitude for different values of permeability parameters using CLT and ELT, as discussed in Sec. VII A. Additionally, the prediction from the domain perturbation method⁸ is cross-validated with the present NS, as discussed in Sec. VII B. The aforementioned sections reveal the limitations of existing models in the literature. Finally, Sec. VII C compares predictions from the spectral model (SM) with numerical simulations (NS) and discusses a non-linear behavior of scaled flow rate at different values of permeability parameters. The behavior of flow parameters (e.g., scaled flow rate, velocity, and wall shear stress) is examined with pattern amplitude, wavelength, and permeability of the porous medium.

A. Evaluating lubrication theory against numerical simulation

The impact of dimensionless pattern amplitude (⁠ $ω$ ⁠) on the scaled flow rate for various $κ_{l}$ values and perturbation parameters (⁠ $δ$ ⁠) is depicted in Fig. 2. The solid and dashed curves illustrate the predictions from ELT and CLT, respectively. Equations (18) and (19) are utilized for CLT and ELT, respectively. Panel (a) of Fig. 2 demonstrates the variation of scaled flow rate with pattern amplitude (⁠ $ω$ ⁠) for $κ_{l} = 0.5, 2.0$ at $δ = 0.1$ ⁠. When the channel height (h) is 0.1 times channel length (L), the scaled flow rate ratio (⁠ $Q_{l} / Q_{l s}$ ⁠) increases with pattern amplitude across a range of $κ_{l}$ values. Both the numerical study and theoretical approaches (CLT and ELT) indicate an increasing trend in flow rate with $κ_{l}$ ⁠, as depicted in panel (a) of Fig. 2. This comparison highlights the higher accuracy of ELT over CLT, and a close agreement is found between ELT and numerical values at both the Stokes flow limit (⁠ $κ_{l} \to 0$ ⁠) and the Darcy flow limit (⁠ $κ_{l} ≫ 1$ ⁠). CLT is accurate at smaller pattern amplitudes (⁠ $ω \sim 0.2$ ⁠). Conversely, when the channel height is 0.5 times a channel length (⁠ $δ = 0.5$ ⁠), the scaled flow rate decreases with increasing pattern amplitude, as illustrated in panel (b) of Fig. 2. This decreasing behavior is predicted by both ELT and NS. In contrast, CLT consistently predicts an increasing flow rate with pattern amplitude at both limits, which shows the unphysical behavior of the scaled flow rate. The ratio between h and L is second order in the governing equation [Eq. (3)], and the CLT provides the solution for $O (δ^{0})$ ⁠, which is an incomplete investigation to capture the viscous stress. Consequently, the higher order correction in the CLT, which corresponds to a more accurate evaluation of viscous stresses, improves the applicability of lubrication theory from a significantly higher wavelength. Therefore, the higher-order term plays a significant role in the perturbation analysis to understand fluid behavior. Considering the higher-order term in the CLT accurately captures the viscous stress. Panel (b) indicates that the current theoretical models align closely with numerical values at smaller pattern amplitudes (⁠ $ω \sim 0.1 - 0.2$ ⁠) from the Stokes to Darcian flow limits.

FIG. 2.

Variation of the scaled flow rate with pattern amplitude across different permeability mediums (⁠ $κ_{l}$ ⁠) is presented for $δ = 0.1$ and $δ = 0.5$ ⁠, as illustrated in panels (a) and (b), respectively. Equations (18) and (19) are used to show the prediction from the classical lubrication theory (CLT) and extended lubrication theory (ELT). Numerical simulations (NS) are represented symbolically.

The error percentage of the current theoretical model is assessed against numerical simulations to gauge accuracy. For example, at a perturbation parameter $δ = 0.1$ and a permeability parameter $κ_{l} = 0.5$ ⁠, ELT exhibits an error percentage of 2.644 $%$ ⁠, whereas CLT shows a significantly higher error percentage of 9.54 $%$ ⁠, both calculated at a pattern amplitude $ω = 0.8$ ⁠. When the perturbation parameter is increased to $δ = 0.5$ ⁠, ELT maintains an error percentage below 10 $%$ at $ω = 0.3$ and $κ_{l} = 2$ ⁠. However, the error percentage for ELT surpasses the 10 $%$ threshold at $ω = 0.32$ ⁠. It is observed that the relative difference between NS and ELT models increases as the pattern amplitude $ω$ grows, indicating that the model's accuracy diminishes with larger pattern amplitudes. Section VII B presents a discussion of the results derived from the domain perturbation approach [ $O (α^{2})$ ], spectral method (SM), and numerical simulations (NS).

B. Comparative analysis of domain perturbation method and numerical simulation

Ng and Wang⁸ investigated Darcy–Brinkman flow through a sinusoidal grooved microchannel using the domain perturbation method (DPM). They developed a second-order solution for the scaled flow rate, providing a foundation for understanding the behavior of flow with reference to the pattern surface. This model is accurate at significantly small pattern amplitudes. The present study is an extension of the existing literature model and will develop a model for large pattern amplitudes. In this section, we compare different mathematical models, including DPM and SM, against finite-element-based numerical simulations (NS). Equation (31) represents the prediction from SM, whereas Eq. (34) is a prediction from DPM, which is accurate up to $O (α^{2})$ ⁠. By comparing SM and DPM with NS, we aim to evaluate the accuracy and consistency of these models against numerical results, providing deeper insights into their predictive capabilities.

Figure 3 illustrates the variation of the scaled flow rate (⁠ $Q / Q^{*}$ ⁠) as a function of pattern amplitude (⁠ $α$ ⁠). It compares predictions from the theoretical model with numerical results across various wavelengths (⁠ $λ$ ⁠) and porous media permeabilities (⁠ $κ_{f}$ ⁠). For a smaller wavelength (⁠ $λ = 1$ ⁠), an increase in $α$ leads to the reduced scaled flow rate in the Stokes flow limit $(κ_{f} \to 0)$ ⁠, as shown in panel (a) of Fig. 3. This behavior is captured by all approaches. This compression suggests that the $O (α^{2})$ -solution is accurate at smaller values of pattern amplitude $(α \sim 0.2)$ for $κ_{f} = 0.1, 1.0$ ⁠, and it fails for large values of $α$ ⁠. In contrast, the prediction of scaled flow rate from SM is found to have a good agreement with numerical values at significantly large $α$ and $κ_{f} = 0.1, 1.0$ ⁠. Furthermore, we analyzed for a large wavelength (⁠ $λ = 5$ ⁠), and the permeability of porous media varies from 0 to 1; $\frac{Q}{Q^{*}}$ increases with dimensionless pattern amplitude as shown in panel (c) of Fig. 3. At intermediate wavelength (⁠ $λ = 3$ ⁠) and $κ_{f} = 0.1, 1$ ⁠, the flow rate ratio decreases monotonically with increasing $α$ ⁠. This trend is captured by only $O (ε^{2})$ -solution, as shown in dotted lines in panel (b) of Fig. 3. By contrast, a non-monotonic behavior of $\frac{Q}{Q^{*}}$ is observed using SM [Eq. (31)] and numerical study at both values of $κ_{f} = 0.1, 1.0$ ⁠. This observation suggests that the behavior of the scaled flow rate is non-linear, and the minimum flow rate is noticed at an intermediate wavelength for the Stokes and Darcian flow limits. For significantly large pattern amplitudes (⁠ $α \sim 1$ ⁠), it is noticed that SM and NS exhibit a good agreement at both limiting flow, as indicated by the compression. The case of the large characterizing porous medium $(κ_{f} ≫ 1)$ or Darcian flow limit is discussed in detail in Sec. VII C. We also discussed the increasing and decreasing behavior of the scaled flow rate in the later section.

FIG. 3.

Variation of scaled flow rate ( QQ*) with pattern amplitude ( α) for different permeability mediums ( κf) at different wavelengths ( λ). Panels (a), (b), and (c) represent variations for λ=1, λ=3, and λ=5, respectively. Equations (31) and (34) illustrate the predictions derived from the spectral model (SM) and the asymptotic model [ O(α2)], respectively. Numerical simulations (NS) are denoted symbolically.

Variation of scaled flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) with pattern amplitude (⁠ $α$ ⁠) for different permeability mediums (⁠ $κ_{f}$ ⁠) at different wavelengths (⁠ $λ$ ⁠). Panels (a), (b), and (c) represent variations for $λ = 1$ ⁠, $λ = 3$ ⁠, and $λ = 5$ ⁠, respectively. Equations (31) and (34) illustrate the predictions derived from the spectral model (SM) and the asymptotic model [ $O (α^{2})$ ], respectively. Numerical simulations (NS) are denoted symbolically.

C. Effect of permeability of porous media on scaled flow rate

This section examines the variation in the scaled flow rate with respect to pattern amplitude for large values of permeability parameter (⁠ $κ_{f} ≫ 1$ ⁠), employing both SM and NS. The analysis is also conducted across various wavelengths (⁠ $λ$ ⁠). Equation (31) is employed to plot the scaled flow rate as shown in panels (a) and (b) of Figs. 4 and 5. This compression clearly reveals that $Q / Q^{*}$ initially decreases with $α$ ⁠, followed by an increase as $α$ continues to rise at medium (⁠ $λ = 3$ ⁠) to large (⁠ $λ = 15$ ⁠) wavelengths. A similar trend is also captured by the numerical study. Figures 4 and 5 reveal that the parameter $κ_{f}$ plays a crucial role in determining the transition behavior of the scaled flow rate. At the Stokes limit (⁠ $κ_{f} \to 0$ ⁠) and $κ_{f} \sim 1$ ⁠, the non-monotonic behavior of the scaled flow rate is valid for intermediate wavelength (⁠ $λ = 3$ ⁠). By contrast, a similar trend is also noticed at Darcian flow at any values of wavelength (⁠ $λ \geq 3$ ⁠). As demonstrated in the literature model⁸ and the present lubrication theory model, the dimensional flow rate may exhibit either a monotonic increase or decrease with respect to the wavelength and permeability parameter, as discussed in Secs. VII A and VII B. However, in the spectral model's prediction, the dimensionless flow rate reveals non-linear behavior at both the Stokes and Darcy flow limits. The factors contributing to the observed fluctuations in the dimensional flow rate are examined at the end of this section. Figures 4 and 5 elucidated a good agreement between the numerical values and the prediction from SM. Panels (a) and (b) of Figs. 4 and 5 show that the SM prediction for the scaled flow rate always captured the minimum at large values of $κ_{f}$ ⁠. Each wavelength with respect to $κ_{f}$ has a different pattern amplitude, where flow rate behavior captured the minimum. The transition behavior of scaled permeability occurs at a critical pattern amplitude (⁠ $α_{crt}$ ⁠). Table I shows the values of $α_{crt}$ ⁠; at this point, a non-linear behavior of scaled flow rate is observed. When $α > α_{crt}$ ⁠, the flow rate through a grooved channel is lower than that of standard Poiseuille flow. Conversely, when $α < α_{crt}$ ⁠, the flow rate through a grooved channel exceeds the standard Poiseuille flow rate. When $λ$ is constant, for example, $λ$ = 3, the critical pattern amplitude is ∼0.61 for $κ_{f}$ = 0.5, which is lower than that observed for $κ_{f}$ = 3.0. This finding indicates that $α_{crt}$ increases with $κ_{f}$ ⁠. A similar behavior is also captured for other values of $λ$ ⁠.

FIG. 4.

Variation of scaled flow rate ( QQ*) with pattern amplitude for different permeability ( α) mediums ( κf) at different wavelengths ( λ). Panels (a) and (b) are variations for λ=3 and λ=5, respectively. The solid lines represent the prediction from the spectral model (SM), and Eq. (31) is used to obtain the variation of dimensional flow rate. Numerical simulations (NS) are represented symbolically.

Variation of scaled flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) with pattern amplitude for different permeability (⁠ $α$ ⁠) mediums (⁠ $κ_{f}$ ⁠) at different wavelengths (⁠ $λ$ ⁠). Panels (a) and (b) are variations for $λ = 3$ and $λ = 5$ ⁠, respectively. The solid lines represent the prediction from the spectral model (SM), and Eq. (31) is used to obtain the variation of dimensional flow rate. Numerical simulations (NS) are represented symbolically.

FIG. 5.

Variation of dimensional flow rate ( QQ*) with pattern amplitude for different permeability ( α) mediums ( κf) at different wavelengths ( λ). Panels (a) and (b) are variations for λ=10 and λ=15, respectively. The solid lines represent the prediction from the spectral model (SM), and Eq. (31) is used to obtain the variation of the scaled flow rate. Numerical simulations (NS) are represented symbolically.

Variation of dimensional flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) with pattern amplitude for different permeability (⁠ $α$ ⁠) mediums (⁠ $κ_{f}$ ⁠) at different wavelengths (⁠ $λ$ ⁠). Panels (a) and (b) are variations for $λ = 10$ and $λ = 15$ ⁠, respectively. The solid lines represent the prediction from the spectral model (SM), and Eq. (31) is used to obtain the variation of the scaled flow rate. Numerical simulations (NS) are represented symbolically.

TABLE I.

The transition behavior of dimensional flow rate (⁠ $Q / Q^{*}$ ⁠) is calculated with critical pattern amplitude (⁠ $α_{crt}$ ⁠) at various wavelengths (⁠ $λ$ ⁠) and permeability mediums (⁠ $κ_{f}$ ⁠). The boldface value denotes approximate values of critical pattern amplitude, where the dimensionless flow rate is calculated and transmission behaviour is captured.

$λ$	$κ_{f}$	$α_{crt}$	$Q / Q^{*}$
3	0.5	$\sim$ 0.61	0.9494
3	3	$\sim$ 0.80	0.9298
10	15	$\sim$ 0.58	0.9974
10	25	$\sim$ 0.79	0.9967

$λ$	$κ_{f}$	$α_{crt}$	$Q / Q^{*}$
3	0.5	$\sim$ 0.61	0.9494
3	3	$\sim$ 0.80	0.9298
10	15	$\sim$ 0.58	0.9974
10	25	$\sim$ 0.79	0.9967

We can draw the following conclusions based on the trends observed in Figs. 4 and 5, as well as the aforementioned discussion regarding the impact of $α$ ⁠, $λ$ ⁠, and $κ_{f}$ ⁠. As $α$ increases, there is an increase in the additional area on the lower wall [as shown in panel (b) of Fig. 1] over which viscous forces can enforce the no-slip condition. This increase in area is higher than that of a channel with a flat lower wall at Y = 0. When the wavelength of the pattern is short, the most significant gradients of the flow occur near the lower wall and are inversely proportional to the wavelength. This amplifies the role of the viscous forces acting on the lower surface for the DB flow. Due to this effect, the scaled flow rate decreases with $α$ at short wavelengths. For short waves, the grooves behave as roughness elements that enhance drag. However, in the long-wave limit, the x-gradients in the flow, which are influenced by the wavelength, become less significant than the y-gradients that are determined by the local channel thickness. As a result, when $α$ is small initially and at moderate to large wavelengths with significantly large $κ_{f}$ ⁠, grooves behave as roughness elements that enhance friction. This leads to a negative slope in the flow rate vs amplitude curve, as seen in Figs. 4 and 5. On the other hand, when $α$ is increased, a transition is observed in the flow rate vs amplitude curve. Initially, a minimum is reached, as seen in the figures. However, on further increasing $α$ ⁠, the slope of the curve changes toward a long wave-like positive slope. This is because the y-gradients become more significant than x-gradients due to the formation of progressively narrower constrictions with an increase in $α$ ⁠.

1. Wall shear stress (WSS) distribution near to the grooved wall

This study aims to elucidate the effects of pattern amplitude (⁠ $ε$ ⁠) and permeability parameter (⁠ $κ_{f}$ ⁠) on wall shear stress (WSS) behavior along a grooved wall. We present a graphical representation of WSS distribution in Figs. 6–8. The investigation is performed for small (⁠ $λ = 1$ ⁠) and long (⁠ $λ = 10$ ⁠) wavelengths of the channel. Panels (a) and (b) of Fig. 6 illustrate the variation of WSS across different values of $ε$ at a fixed value of $κ_{f} = 15$ ⁠. If the channel length (L) equals the mean channel height (h) or $λ = 1$ ⁠, WSS reaches its maximum value at the center of the channel (⁠ $X = 0$ ⁠) for all values of $ε$ ⁠. The maximum value of WSS increases with increasing pattern amplitude, as illustrated in panel (a) of Fig. 6. The viscous forces are more prominent near the crest of the grooved surface for a shorter wavelength. For the case where $L ≫ h$ (⁠ $λ = 10$ ⁠) with $κ_{f} = 15$ ⁠, the maximum wall WSS is observed at $X = 0$ for $ε = 0.2$ and $ε = 0.3$ ⁠, as shown in panel (b) of Fig. 6. However, for a slightly larger pattern amplitude (⁠ $ε = 0.4$ ⁠), the maximum WSS decreases, and its location shifts away from $X = 0$ ⁠. This indicates that the viscous effects become less significant as the channel length increases relative to its height.

FIG. 6.

Variation of wall shear wall stress is calculated on the bottom grooved wall for $κ_{f}$ = 15 and various values of $ε$ ⁠. Panel (a) corresponds to $λ = 1$ ⁠, while panel (b) corresponds to $λ = 10$ ⁠. Numerical simulations (NS) are represented symbolically.

FIG. 7.

Variation of wall shear wall stress is calculated on the bottom grooved wall for $ε$ = 0.4 and various values of $κ_{f}$ ⁠. Panel (a) corresponds to $λ = 1$ ⁠, while panel (b) corresponds to $λ = 10$ ⁠. Numerical simulations (NS) are represented symbolically, whereas the solid lines are prediction from the spectral method (SM).

FIG. 8.

Variation of wall shear wall stress is calculated on the bottom grooved wall for $ε$ = 0.4 and various values of $κ_{f}$ ⁠. Panel (a) corresponds to $λ = 3$ ⁠, while panel (b) corresponds to $λ = 5$ ⁠. Numerical simulations (NS) are represented symbolically, whereas the solid lines are prediction from the spectral method (SM).

Figure 7 graphically demonstrates the variation of WSS for different values of $κ_{f}$ (0.1, 1.0, 5.0, 10.0, 15.0, and 20.0) at $ε = 0.4$ for smaller and large wavelength. If $λ = 1$ and $ε = 0.4$ ⁠, the maximum WSS is observed in the Stokes limit (⁠ $κ_{f} ≪ 1$ ⁠), while the minimum WSS is identified in the Darcian limit (⁠ $κ_{f} ≫ 1$ ⁠), as illustrated in panel (a) of Fig. 7. Conversely, for $λ = 10$ ⁠, the maximum WSS is found in the Darcian limit, and the minimum WSS is observed in the Stokes limit, as shown in panel (b) of Fig. 7. The viscous shear effect is highly pronounced in the Stokes flow regime and diminishes with increasing permeability of the porous material. Consequently, the maximum wall shear stress (WSS) is observed in the Stokes flow limit compared to the Darcian flow limit, particularly at small wavelengths. This observation indicates that the maximum and minimum WSS values occur at the narrowest section of the channel, corresponding to $X = 0$ ⁠. Furthermore, we observed that an increase in the permeability parameter significantly affects intermediate wavelengths, such as 3 and 5. Consequently, the influence of WSS on the grooved wall surface becomes less prominent. The distribution of WSS is similar for $λ = 3$ and $λ = 5$ as shown in panels (a) and (b) of Fig. 8. When the channel has a more permeable and shorter wavelength, the viscous effect is less significant. This means there is greater resistance on the grooved wall. As a result, WSS decreases with a more permeable medium. On the other hand, viscous forces are more significant for long wavelengths and larger permeable surfaces. These findings show that amplitude and wavelength strongly depend on the geometrical properties of a porous material. The prediction from the spectral method is perfectly aligned with the numerical values for each case.

2. Variation of scaled flow rate with reference to permeability of porous media (⁠ $κ_{f}$ ⁠)

This section examines the influence of $κ_{f}$ on the variation of scaled flow rate. To elucidate its significance, we utilize both spectral method and numerical simulation approaches. Figure 9 demonstrates the transition from the Stokes flow to the Darcian flow. Equation (31) presents the prediction of the scaled flow rate using the spectral approach, while Eq. (34) provides the $O (α^{2})$ solution for the scaled flow rate. The comparison between the analytical solution and numerical study for $λ$ = 3 and $λ$ = 5 is depicted in panels (a) and (b), respectively. At the Stokes and Darcy flow limits, the results from the current spectral method reveal a good agreement with numerical simulations across all pattern amplitude values. At $λ = 3$ ⁠, the scaled flow rate first decreases and then increases as $κ_{f}$ increases, as shown in panel (a). When $κ_{f}$ varies from 0 to 10 with $λ = 3$ ⁠, the behavior of the scaled flow rate exhibits a decreasing trend, followed by an increasing trend, which can be attributed to the effects of surface roughness and confinement, respectively. Both the numerical simulation and the $O (α^{2})$ solution show a similar trend; however, a clear difference exists between them. In contrast, for $λ = 5$ ⁠, the scaled flow rate steadily decreases as $κ_{f}$ increases. At moderately large wavelengths, viscous forces become more pronounced with increasing values of $κ_{f}$ ⁠. The permeability of the porous medium reduces the scaled flow rate. In the Stokes-to-Darcian flow regime (⁠ $κ_{f} \geq 6$ ⁠), the flow rate through the grooved channel (Q) matches that of a straight channel (⁠ $Q^{*}$ ⁠). This behavior is observed from all approaches, as shown in panel (b). As seen in Fig. 9, larger channel wavelengths are associated with higher flow rates in both Stokes and Darcy flow limits. In contrast, moderate wavelengths show a combination of increasing and decreasing flow rates in these two regimes. A non-linear phenomenon is observed at $κ_{f} =$ 1.885. For $κ_{f} < 1.885$ ⁠, the observed decreasing trend is attributed to the predominant influence of surface roughness. In contrast, for $κ_{f} > 1.885$ ⁠, the increasing trend indicates the confinement effects are becoming more significant. When analyzing DB flow in a grooved channel, it is essential to consider the channel's wavelength. This factor affects flow behavior and its interaction with the patterned surface, influencing parameters such as pressure drop, flow distribution, and overall channel performance.

FIG. 9.

The effect of permeability medium for the dimensionless flow rate ( QQ*) is shown by the spectral method, numerical simulation, and asymptotic results of Ng and Wang.8 Panel (a) corresponds to λ=3, while panel (b) corresponds to λ=5. Numerical simulations (NS) are represented symbolically, whereas the solid lines stand for the prediction from the spectral method (SM).

The effect of permeability medium for the dimensionless flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) is shown by the spectral method, numerical simulation, and asymptotic results of Ng and Wang.⁸ Panel (a) corresponds to $λ = 3$ ⁠, while panel (b) corresponds to $λ = 5$ ⁠. Numerical simulations (NS) are represented symbolically, whereas the solid lines stand for the prediction from the spectral method (SM).

3. Variation of scaled flow rate with reference to the height of the channel

The variation of the scaled flow rate is plotted with scaled minimum spacing between grooved and straight plates for a small amplitude-to-pitch ratio (⁠ $ε = 1$ ⁠) for various permeability medium as shown in Fig. 10. The logarithmic scale is used for the independent variable to capture the low-channel height regime and the minimum in each curve precisely. Similarly, a large amplitude-to-pitch (⁠ $ε = 3$ ⁠) ratio is also investigated here as depicted in Fig. 11. When the flat wall (top wall) is brought closer from a large distance to a patterned surface, two competing phenomena are observed with various amplitude-to-pitch ratios. For example, when the permeability is small (⁠ $κ_{f} = 0.05$ ⁠), the dimensional flow rate is larger for a grooved channel as compared to a straight channel. A similar trend is also observed with different values of permeability as shown in Fig. 10. The flow rate ratio decreases with increasing values of permeability if it is scaled with $h - a$ rather than h. The flow rate ratio is always higher than the standard Poiseuille flow through a straight channel. For a large amplitude-to-pitch, the flow rate ratio is always less than that of the standard Poiseuille flow through a straight channel. In this situation, the non-linear behavior of dimensionless flow rate is identified with various permeability values, as illustrated in Fig. 11. This transition behavior of flow rate can be summarized dimensionally. When the dimensional amplitude-to-pitch ratio (⁠ $ε$ ⁠) is 3 and L = 10 μm, then amplitude (a) is equal to 4.75 μm with different permeability (⁠ $κ_{f}$ ⁠) values such as 0.05, 0.5, 1, and 2. Using these quantities, the scaled flow rate reaches a minimum of 45.67 $%$ ⁠, 56.65 $%$ ⁠, 69.29 $%$ ⁠, and 81.88 $%$ of the flow rate of the plane channel flow with the same h and different $κ_{f}$ ⁠, when h is 30.6 μm. Smaller and larger channel sizes than h = 30.6 μm will result in a larger flow rate ratio to the plane channel flow with the same h. The aforementioned discussion shows that the flow rate strongly depends on the permeability of porous media and the dimension of the pattern amplitude.

FIG. 10.

At the Stokes to Darcian flow limit, variation of dimensionless flow rate ( QQ*) with a scaled minimum spacing between patterned and unpatterned plates for the amplitude-to-pitch ratio ( ε=1). A logarithmic scale is used for the independent variable to visually resolve the low-channel height regime and the minimum in each curve. Equation (31) is used to obtain the variation of the scaled flow rate.

At the Stokes to Darcian flow limit, variation of dimensionless flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) with a scaled minimum spacing between patterned and unpatterned plates for the amplitude-to-pitch ratio (⁠ $ε = 1$ ⁠). A logarithmic scale is used for the independent variable to visually resolve the low-channel height regime and the minimum in each curve. Equation (31) is used to obtain the variation of the scaled flow rate.

FIG. 11.

At the Stokes to Darcian flow limit, variation of dimensionless flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) with a scaled minimum spacing between patterned and unpatterned plates for the amplitude-to-pitch ratio (⁠ $ε = 3$ ⁠). Equation (31) is used to obtain the variation of the scaled flow rate.

4. Variation of scaled flow rate with wavelengths

The plot in Fig. 12 depicts how the scaled flow rate changes with wavelengths (⁠ $λ$ ⁠) for a fixed $α$ value of 0.9. This situation is highly constrained, indicating that the thinnest section of the channel is 19 times bigger than the narrowest portion. An investigation is conducted into the Stokes flow limit, as illustrated in panel (a), and the Darcian flow limit, as depicted in panel (b). Equation (31) is used for the prediction by SM, whereas Eq. (34) is considered for an $O (α^{2})$ -solution. Both theories are compared with NS. A good agreement is found between SM and NS for each value of wavelengths at both limits. For $κ_{f} = 0.5$ ⁠, the flow rate is lower than that of the conventional Poiseuille flow at smaller wavelengths, while it exceeds the Poiseuille flow at larger wavelengths as shown in panel (a). A smaller wavelength is emphasized by the roughness effect, while a longer wavelength is associated with the confinement effect. For $κ_{f} ≫ 1$ ⁠, the flow rate is always higher than that of the conventional Poiseuille flow at any values of wavelength as shown in panel (b). The $O (α^{2})$ model is not accurate at smaller to large wavelengths, whereas the prediction from Eq. (31) is accurate at any value of wavelength for Stokes limit and Darcian limit, as shown in both panels of Fig. 12.

FIG. 12.

The effect of wavelength for the dimensionless flow rate ( QQ*) is shown by the spectral method, numerical simulation, and asymptotic results of Ng and Wang.8 Panel (a) corresponds to κf=0.5, while panel (b) corresponds to κf=2. Numerical simulations (NS) are represented symbolically, whereas the solid lines stand for the prediction from the spectral method (SM).

The effect of wavelength for the dimensionless flow rate (⁠ $\frac{Q}{Q^{*}}$ ⁠) is shown by the spectral method, numerical simulation, and asymptotic results of Ng and Wang.⁸ Panel (a) corresponds to $κ_{f} = 0.5$ ⁠, while panel (b) corresponds to $κ_{f} = 2$ ⁠. Numerical simulations (NS) are represented symbolically, whereas the solid lines stand for the prediction from the spectral method (SM).

5. Variation of longitudinal velocity

To gain a better understanding of the role played by $κ_{f}$ in the current problem, it is essential to analyze the axial velocity profile. Equation (22) is not only helpful to identify the flow rate through a channel, but it is also utilized to understand other flow behavior such as flow profile (velocity). The depth-wise variation of velocity at different transverse locations within the flow cross section is plotted in Figs. 13 and 14. For Darcian flow, the analysis is performed for both small (⁠ $λ = 1$ ⁠) and large (⁠ $λ = 10$ ⁠) wavelengths. The results focus on the thickest part of the channel (⁠ $X = - π$ ⁠), as illustrated in panels (a) and (b) of Fig. 13. However, Fig. 14 presents the flow profile at the thinnest section of the channel (⁠ $X = 0$ ⁠) for various values of $κ_{f}$ ⁠. At various values of $κ_{f}$ ⁠, the comparison in these figures clearly shows that the theoretical distribution matches the numerical distribution in the Y direction when $ε =$ 0.5. At $λ = 1$ ⁠, the velocity profile deviates from a parabolic shape, exhibiting a flattened characteristic in the center of the channel across all values of $κ_{f}$ ⁠, as depicted in panel (a) of Fig. 13. It is observed that the axial velocity decreases with increasing $κ_{f}$ ⁠. Higher values of $κ_{f}$ reduce the effective permeability in the flow direction, thereby increasing the overall resistance to flow. As a result, the fluid velocity diminishes with increasing $κ_{f}$ ⁠. Furthermore, an increase $κ_{f}$ from 2 to 10 results in a reduction of velocity at the thickest section of the channel. With a wavelength of $λ$ = 10 and $κ_{f}$ = 6, the depth-wise velocity distribution displays a parabolic profile. The maximum velocity occurs at the center of the channel, as illustrated in panel (b) of Fig. 13. When $κ_{f}$ increases from 6 to 15, the velocity profile deviates from the parabolic shape. For larger values of $κ_{f}$ ⁠, the axial velocity profile becomes increasingly uniform, with the influence of the wall confined to a thin boundary layer near the surface. As a result, the wall resistance has a minimal contribution to the overall resistance, which is predominantly governed by the porous medium. This results in a decrease in the overall velocity. This comparison signifies the dominance of Darcian flow over the Stokes flow at different wavelengths. To conclude, if the length of the channel (L) is 10 times the mean height (h), it is essential to employ reduced values of $κ_{f}$ in order to achieve a parabolic velocity profile.

FIG. 13.

At various permeability mediums, the variation of velocity profile along the transverse direction is shown at the thickest section of that channel (⁠ $X = - π$ ⁠) for moderately large pattern amplitude (⁠ $ε = 0.5$ ⁠). Equation (22) is used to obtain the variation of the scaled flow rate. Panel (a) corresponds to $λ = 1$ ⁠, while panel (b) corresponds to $λ = 10$ ⁠. Panels (c) and (d) illustrate the variation of axial velocity at $X = - π$ and $X = 0$ for $λ = 10$ ⁠, $κ_{f}$ = 2, 4. The solid lines represent the predictions obtained using the spectral method (SM), while the symbols correspond to the results from the numerical simulations (NS).

FIG. 14.

At various permeability mediums (⁠ $κ_{f}$ ⁠), the variation of velocity profile along the transverse direction is shown at the thinnest section of that channel (⁠ $X = 0$ ⁠) for moderately large pattern amplitude (⁠ $ε = 0.5$ ⁠). Equation (22) is used to obtain the variation of the scaled flow rate. Panel (a) corresponds to $λ = 1$ ⁠, while panel (b) corresponds to $λ = 10$ ⁠. Numerical simulations (NS) are represented symbolically.

Panels (c) and (d) illustrate the variation in the axial velocity for two values of $κ_{f}$ (2 and 4) at two different channel sections (⁠ $X = - π$ and $X = 0$ ⁠) for $λ = 10$ and $ε = 0.5$ ⁠. The results indicate that as $κ_{f}$ increases, the velocity magnitude decreases, particularly at the thickest section of the channel (⁠ $X = - π$ ⁠). In regions with high material porosity, the flow velocity is reduced, and the velocity profile transitions from parabolic to flatter, as shown in panel (c) of Fig. 13. In contrast, at the thinnest section of the channel, the velocity profile remains close to parabolic, and the maximum velocity W values are nearly identical for $κ_{f}$ = 2 and 4, as shown in panel (d). Overall, $κ_{f}$ has a more pronounced effect at shorter wavelengths compared to longer wavelengths. A strong agreement is observed between the spectral and numerical approaches.

At the thinnest section of the channel (⁠ $X = 0$ ⁠), the velocity profile is not parabolic with reference to any values of $κ_{f}$ at $λ$ = 1, as illustrated in panel (a) of Fig. 14. For a large wavelength, the velocity profile is always captured in a parabolic nature or similar to Poiseuille flow as illustrated in panel (b) of Fig. 14. The influence of $κ_{f}$ becomes more pronounced when the length of the channel (L) is equal to its mean height (h). Conversely, when the channel length is ten times the mean height, the effect of $κ_{f}$ diminishes significantly. Overall, the magnitude of the velocity (W) decreases with increasing $κ_{f}$ at any wavelength of the channel. The aforementioned discussion reveals that the flow profiles strongly depend on $λ$ and $κ_{f}$ ⁠. A good agreement is found between the prediction from the spectral model and the numerical simulation. Furthermore, the no-slip boundary condition at both walls is achieved using spectral method and numerical simulations.

The analysis focuses on the behavior of the axial velocity (W) at the Stokes flow limit (⁠ $κ_{f} ≫ 1$ ⁠) for both the thickest (⁠ $X = - π$ ⁠) and thinnest (⁠ $X = 0$ ⁠) sections of the channel, considering both small (⁠ $λ = 1$ ⁠) and large wavelengths (⁠ $λ = 10$ ⁠), as illustrated in Fig. 15. The investigation is conducted for various values of $κ_{f}$ ⁠, including 0.1, 0.3, 0.5, and 1. For $λ = 1$ and each value of $κ_{f}$ ⁠, the maximum velocity is observed at the center of the channel at both sections, as shown in panels (a) and (c). At $κ_{f} ≫ 1$ ⁠, the velocity profile closely resembles standard Poiseuille flow (parabolic profile). As $κ_{f}$ increases from 0.1 to 0.5 and then from 0.5 to 1, the fluid velocity in the porous medium decreases, with the velocity profile becoming progressively flatter at the center of the channel, approaching the Darcian flow limit. Panels (a) and (c) of Fig. 15 indicate that flow resistance increases with permeability at both channel sections, and the influence of the porous material becomes more pronounced slightly away from the wall. For $λ = 10$ and $X = - π$ ⁠, similar trends are observed, as shown in panel (b) of Fig. 15. However, at large wavelengths, the effect of $κ_{f}$ is less pronounced compared to the small wavelength case. The minimum axial velocity is observed at $κ_{f}$ = 1. At $X = 0$ ⁠, the effect of the permeability of the porous media is less significant compared to $X = - π$ ⁠, as shown in panel (d). The maximum velocity for the transition from Stokes to Darcian flow is nearly identical for both sections. For all cases, excellent agreement is found between the spectral method and numerical simulations, as shown in Fig. 15. The above-mentioned analysis demonstrates that the presence of a porous medium effectively suppresses the growth of boundary layers on rigid, impermeable plates, preventing them from extending to the center of the channel. This behavior contrasts sharply with that of plane Poiseuille flow, where boundary layers develop uniformly from the plates and merge at the channel center, as indicated by $\frac{d W (y)}{d y} = 0$ ⁠. Consequently, the introduction of a porous medium distinctly modifies the flow dynamics compared to conventional cases.

FIG. 15.

Variation of axial velocity in the Stokes flow regime is examined at different locations within the channel, specifically at the thinnest (⁠ $X = 0$ ⁠) and thickest (⁠ $X = - π$ ⁠) sections for $λ = 1$ and $λ = 10$ ⁠. The solid lines represent predictions from the special method (SM), while the symbols denote results from the numerical simulations (NS).

Figure 16 presents both the contour and surface plots showing the variation in the magnitude of the velocity (W) with different values of $κ_{f}$ for $ε = 0.5$ and $λ = 10$ ⁠. The figure indicates that as $κ_{f}$ increases from the Stokes to Darcian flow limit, the fluid velocity decreases, particularly at the thickest section of the channel (⁠ $X = - π$ ⁠), due to the influence of permeability in the porous medium. For instance, at the Stokes flow limit $κ_{f}$ = 0, the maximum velocity is ∼0.224, while for $κ_{f}$ = 5, it decreases to around 0.055. In contrast, at the thinnest part of the channel (⁠ $X = 0$ ⁠), the maximum velocity is nearly the same (0.003) for both the Stokes and Darcy limits. This observation suggests that the permeability effect is more pronounced in the thicker sections of the channel than in the thinner sections, particularly for long wavelengths. Additionally, the maximum velocity is observed at the Stokes flow limit, as compared to the Darcy flow limit. In high-permeability conditions, flow resistance in wave channels becomes more pronounced, counteracting the applied pressure gradient. As a result, the velocity profile tends to flatten.

FIG. 16.

For a large wavelength (⁠ $λ = 10$ ⁠) and pattern amplitude (⁠ $ε = 0.5$ ⁠), panels (a)–(h) illustrate the variation in velocity magnitude (W) for different $κ_{f}$ under the Stokes and Darcian flow limits.

VIII. CONCLUSION

In this paper, we present a comprehensive study of hydrodynamic flow through a corrugated channel containing sparse porous material. The primary objective is to investigate Darcy–Brinkman flow in a grooved microchannel, focusing on cases with large pattern amplitudes and high permeability in the porous medium, particularly for complex geometries. The analysis is conducted using both spectral and lubrication approaches: the lubrication theory is applied to complex topographies, while spectral methods are utilized for sinusoidal channels. While previous models in the literature, based on the domain perturbation theory, offer accurate predictions for small pattern amplitudes and low permeability, their applicability diminishes as these parameters increase. Furthermore, studies on Darcy–Brinkman flow have predominantly focused on transverse flow, with limited attention given to longitudinal flow in porous media. To address these challenges, we propose a mathematical model that ensures an accurate scaled flow rate even in scenarios with high permeability and large pattern amplitudes. This study employs a systematic perturbation analysis to derive higher-order terms, extending the traditional lubrication approximation. The analytical results are thoroughly compared with numerical simulations, demonstrating the accuracy and reliability of the proposed model.

The key conclusions are as follows: (i) In both the Stokes and Darcy limits, the higher-order extended lubrication theory (ELT) shows excellent agreement with numerical simulations at $δ = 0.1$ ⁠, demonstrating a significant improvement over the classical lubrication theory (CLT). However, in the case of $δ = 0.5$ ⁠, ELT fails with increasing pattern amplitude. (ii) The spectral model overcomes these limitations of the lubrication approximation at both limiting flows. Based on the Fourier theory, the current proposed model is more accurate than the domain perturbation theory, regardless of dependent parameters. We developed a model that is accurate when the crest of a grooved wall is very close to the upper wall. (iii) According to the literature model, the dimensional flow rate may increase or decrease with dependent parameters, which provides an incomplete behavior for the Darcy–Brinkmann flow through pattern topographies. The present spectral model predicts a non-monotonic behavior of the scaled flow rate at the Darcian limit and provides a broader view of the effect of micropatterning on Darcy–Brinkmann flow. The non-linear behavior of the scaled flow rate is a significant finding of the present work, whereas the literature models are unable to capture this trend. The Darcy viscous drag force is very prominent and depends on the geometry's orientation. (iv) At smaller wavelengths, the maximum shear stress is observed at the crest of the pattern surface for Stokes flow (⁠ $κ_{f} ≪ 1$ ⁠) and Darcian flow limit (⁠ $κ_{f} ≫ 1$ ⁠). In contrast, for a large wavelength, the wall shear stress behavior significantly alters with the permeability medium. The amplitude-to-pitch ratio significantly influences the range of maximum and minimum shear stress values. (v) At the Stokes flow limit, the velocity profile has a parabolic shape, similar to standard Poiseuille flow. However, in the case of a high permeability medium (⁠ $κ_{f} ≫ 1$ ⁠) or the Darcian limit, the profile becomes flatter.

Equations (22) and (31) are key equations to predict the behavior of flow parameters for Darcy–Brinkman flow through a grooved microchannel. We have demonstrated the influence of dependent variables on the scaled flow rate with micropatterning. The findings of the present study are highly significant and have broad applications in multiple fields, such as chemical reactors, filtering equipment, heat exchangers, and flow through rock fractures.

ACKNOWLEDGMENTS

The author would like to thank Dr. Tim Persoons, Department of Mechanical, Manufacturing and Biomedical Engineering, Trinity College Dublin, Dublin, Ireland, for providing access to the software.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Mainendra Kumar Dewangan: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead).

DATA AVAILABILITY

The data that support the findings of this study are available within the article.

APPENDIX A: NUMERICAL SIMULATION

Finite-element analysis using COMSOL Multiphysics is employed to evaluate the present model and assess the limitations of existing literature models. The partial differential equation (PDE) module in COMSOL is utilized to solve the governing equations and boundary conditions. The PARDISO (parallel direct sparse solver) direct solver, a stationary solver with relative tolerance $τ_{r t} = 10^{- 4}$ was used to obtain the numerical values. For the analysis, a single periodic cell of the microchannel is considered, with periodic boundary conditions applied along the y axis within the computational domain. Numerical simulation is performed for the lubrication approach and spectral approach separately. For the lubrication theory, Eq. (3) with Eqs. (4a) and (4b) is solved using the PDEs module in the software, whereas, for the spectral method, Eq. (20) with Eqs. (21a) and (21b) is considered for the solution. No slip boundary conditions at walls (X-edges) and periodic boundary conditions at Y-edges are chosen. The domain is discretized using a triangular mesh, with mesh refinement ensuring accuracy up to 3–4 decimal places to evaluate mesh dependency. A grid sensitivity analysis is performed to examine how variations in the mesh size affect the computed flow quantities. This includes an independent grid study for each wavelength and pattern amplitude. The size of the computational cells depends on the maximum and minimum element sizes, element size growth rate, and maximum allowed iterations. We found that further refining the mesh did not affect the accuracy of the numerical study. The maximum and minimum element sizes are chosen based on the grid independence test as depicted in Table II. In the numerical simulations (NS), we have assessed the limitations of the computational framework. Our findings indicate that the NS are sufficiently robust to handle extreme scenarios, such as (i) cases where the crest of the grooved surface approaches the upper flat surface (⁠ $α \sim 1$ ⁠) and (ii) situations where the permeability of the porous medium is exceptionally high (⁠ $κ_{f} ≫ 1$ ⁠). The spectral model is shown to provide faster predictions compared to the numerical simulations, highlighting its computational efficiency.

TABLE II.

Variation of scaled flow rate at $λ = 10$ ⁠, $κ_{f} = 25$ ⁠, and $ε = 0.5$ for different mesh sizes. The boldface value denotes that the grid sensitivity is accurate up to 4 decimal points.

Element size	Number of cell count	$Q / Q^{*}$
0.2	791	0.991 65
0.1	1 084	0.995 92
0.08	1 637	0.996 47
0.05	4 115	0.996 87
0.03	11 314	0.996 94
0.02	25 255	0.996 95
0.01	101 011	0.996 96

Element size	Number of cell count	$Q / Q^{*}$
0.2	791	0.991 65
0.1	1 084	0.995 92
0.08	1 637	0.996 47
0.05	4 115	0.996 87
0.03	11 314	0.996 94
0.02	25 255	0.996 95
0.01	101 011	0.996 96

APPENDIX B: SUPPORTING EQUATIONS

\frac{\partial^{2} W_{0 a r b}}{\partial X^{2}} = - \frac{sinh (κ_{l} - κ_{l} Y) (H^{″} (X) + κ_{l} H^{'} {(X)}^{2} tanh (\frac{1}{2} (κ_{l} - κ_{l} H (X))))}{2 κ_{l} {cosh}^{2} (\frac{1}{2} (κ_{l} - κ_{l} H (X)))},

(B1)

η_{2} = \frac{sinh (H κ_{f}) + cosh (H κ_{f})}{cosh (\frac{H κ_{f}}{2})},

(B2)

T_{11} = (sinh (κ_{l} (H (X) + 1)) + cosh (κ_{l} (H (X) + 1))),

(B3)

χ_{1} = sinh (2 κ_{l}) + cosh (2 κ_{l}),

(B4)

χ_{2} = sinh (H κ_{f}) + cosh (H κ_{f}),

(B5)

\begin{matrix} T_{22} = ((Y - 2) sinh (κ_{l} (Y - H (X))) - Y sinh (κ_{l} (H (X) + Y - 2)) \\ + H (X) (sinh (κ (Y - H (X))) + sinh (κ_{l} (H (X) + Y - 2)))) . \end{matrix}

(B6)

It is possible to list the coefficients found in Eq. as follows:

\begin{matrix} T_{1} = sinh (κ_{l}) + cosh (κ_{l}) - sinh (κ_{l} ε cos (2 π X)) - cosh (κ_{l} ε cos (2 π X)), \\ T_{2} = sinh (κ_{l}) + cosh (κ_{l}) + sinh (κ_{l} ε cos (2 π X)) + cosh (κ_{l} ε cos (2 π X)), \\ T_{3} = cosh (κ_{l} Y - κ_{l} ε cos (2 π X)) - sinh (κ_{l} Y - κ_{l} ε cos (2 π X)), \\ T_{4} = ε cosh (4 κ_{l}) cos (2 π X) - Y cosh (4 κ_{l}) - Y cosh (2 κ_{l} + 2 κ_{l} Y) + 2 cosh (2 κ_{l} + 2 κ_{l} Y), \end{matrix}

(B7)

\begin{matrix} T_{5} = - ε cos (2 π X) cosh (2 κ_{l} + 2 κ_{l} Y) + Y cosh (2 κ_{l} + 2 κ_{l} ε cos (2 π X)) \\ - 2 cosh (2 κ_{l} + 2 κ_{l} ε cos (2 π X)), \end{matrix}

(B8)

\begin{matrix} T_{6} = Y cosh (2 κ_{l} ε cos (2 π X) + 2 κ_{l} Y) - ε cos (2 π X) cosh (2 κ_{l} ε cos (2 π X) + 2 κ_{l} Y) \\ + ε cos (2 π X) cosh (2 κ_{l} + 2 κ_{l} ε cos (2 π X)), \end{matrix}

(B9)

\begin{matrix} T_{7} = - ε cos (2 π X) sinh (2 κ_{l} + 2 κ_{l} Y) + Y sinh (2 κ_{l} + 2 κ_{l} ε cos (2 π X)) \\ + Y sinh (2 κ_{l} ε cos (2 π X) + 2 κ_{l} Y) - ε cos (2 π X) sinh (2 κ_{l} ε cos (2 π X) + 2 κ_{l} Y) \\ + ε sinh (4 κ_{l}) cos (2 π X) + ε cos (2 π X) sinh (2 κ_{l} + 2 κ_{l} ε cos (2 π X)) \\ - 2 sinh (2 κ_{l} + 2 κ_{l} ε cos (2 π X)) + Y sinh (4 κ_{l}) - Y sinh (2 κ_{l} + 2 κ_{l} Y) + 2 sinh (2 κ_{l} + 2 κ_{l} Y), \end{matrix}

(B10)

T_{8} = (- sinh (κ_{l}) - cosh (κ_{l}) + sinh (κ_{l} ε cos (2 π X)) + cosh (κ_{l} ε cos (2 π X))),

(B11)

T_{9} = sinh (κ_{l}) + cosh (κ_{l}) + sinh (κ_{l} ε cos (2 π X)) + cosh (κ_{l} ε cos (2 π X)),

(B12)

\begin{matrix} R_{1} = - 1 + 8 cosh (κ_{l}) + 4 κ_{l} cosh (κ_{l}) - 12 κ_{l} cosh (2 κ_{l}) - 8 cosh (3 κ_{l}) + 4 κ_{l} cosh (3 κ_{l}) \\ + 4 κ_{l} sinh (κ_{l}) + 8 sinh (κ_{l}) - 12 κ_{l} sinh (2 κ_{l}) + 4 κ_{l} sinh (3 κ_{l}) - 8 sinh (3 κ_{l}) \\ + sinh (4 κ_{l}) + cosh (4 κ_{l}) . \end{matrix}

(B13)

REFERENCES

1.

H.

Darcy

,

The Public Fountains of the City of Dijon

(

Victor Dalmont

,

Paris, France

,

1856

).

https://doi.org/10.1007/BF02120313

2.

H. C.

Brinkman

, “

A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles

,”

Appl. Sci. Res.

1

,

27

–

34

(

1949

).

https://doi.org/10.1021/ie50725a003

3.

J.

Herzig

,

D.

Leclerc

, and

P. L.

Goff

, “

Flow of suspensions through porous media—Application to deep filtration

,”

Ind. Eng. Chem.

62

,

8

–

35

(

1970

).

https://doi.org/10.1115/1.2911418

4.

G.

Hwang

and

C.

Chao

, “

Heat transfer measurement and analysis for sintered porous channels

,”

J. Heat Mass Transfer

116

,

456

–

464

(

1994

).

https://doi.org/10.1016/j.applthermaleng.2016.05.096

5.

P.

Smakulski

and

S.

Pietrowicz

, “

A review of the capabilities of high heat flux removal by porous materials, microchannels and spray cooling techniques

,”

Appl. Therm. Eng.

104

,

636

–

646

(

2016

).

https://doi.org/10.1016/S0096-3003(01)00314-9

6.

E.

Elshehawey

,

E. M.

Elbarbary

, and

N. S.

Elgazery

, “

Effect of inclined magnetic field on magneto fluid flow through a porous medium between two inclined wavy porous plates (numerical study)

,”

Appl. Math. Comput.

135

,

85

–

103

(

2003

).

https://doi.org/10.1016/j.ijheatmasstransfer.2008.08.040

7.

K.

Khanafer

,

B.

Al-Azmi

,

A.

Marafie

, and

I.

Pop

, “

Non-Darcian effects on natural convection heat transfer in a wavy porous enclosure

,”

Int. J. Heat Mass Transfer

52

,

1887

–

1896

(

2009

).

https://doi.org/10.1007/s11242-010-9580-1

8.

C.-O.

Ng

and

C.

Wang

, “

Darcy–Brinkman flow through a corrugated channel

,”

Transp. Porous Media

85

,

605

–

618

(

2010

).

https://doi.org/10.1017/S0022112085003172

9.

J. J.

Higdon

, “

Stokes flow in arbitrary two-dimensional domains: Shear flow over ridges and cavities

,”

J. Fluid Mech.

159

,

195

–

226

(

1985

).

https://doi.org/10.1063/5.0223856

10.

F. Q.

Al-Daamee

,

N. H.

Hamza

,

M. I.

Khan

, and

M.

Al-Dossari

, “

Thermal performance analysis of sinusoidal corrugated channels: A comparative study of thermo-hydraulic and entropy generation

,”

Phys. Fluids

36

,

093622

(

2024

).

https://doi.org/10.1016/j.ijheatfluidflow.2024.109580

11.

F. Q.

Al-Daamee

and

N. H.

Hamza

, “

Hydro-thermal characteristics of flow of different corrugated channels: Experimental and numerical approaches

,”

Int. J. Heat Fluid Flow

110

,

109580

(

2024

).

https://doi.org/10.1016/j.icheatmasstransfer.2024.107814

12.

M.

Fathi

,

M. M.

Heyhat

,

M. Z.

Targhi

, and

A.

Emadi

, “

Semi-porous-fin microchannel heat sinks for enhanced micro-electronics cooling

,”

Int. Commun. Heat Mass Transfer

157

,

107814

(

2024

).

https://doi.org/10.1021/acs.langmuir.3c02150

13.

S. K.

Mehta

and

P. K.

Mondal

, “

AC electrothermal effect promotes enhanced solute mixing in a wavy microchannel

,”

Langmuir

39

,

16797

–

16806

(

2023

).

https://doi.org/10.1063/5.0228640

14.

S.

Yuan

,

X.

Liu

,

X.

Liu

, and

J.

Deng

, “

Effect of induced charge forming active vortex structures in serpentine microchannels on fluid mixing driven by pressure flow

,”

Phys. Fluids

36

,

102014

(

2024

).

https://doi.org/10.1063/5.0194724

15.

A.

Bhattacharya

,

S.

Sarkar

,

A.

Halder

,

N.

Biswas

, and

N. K.

Manna

, “

Mixing performance of t-shaped wavy-walled micromixers with embedded obstacles

,”

Phys. Fluids

36

,

033609

(

2024

).

16.

S.

Sisavath

,

A.

Al-Yaarubi

,

C. C.

Pain

, and

R. W.

Zimmerman

, “

A simple model for deviations from the cubic law for a fracture undergoing dilation or closure

,” in

Thermo-Hydro-Mechanical Coupling in Fractured Rock

(

Birkhäuser

,

Basel

,

2003

), pp.

1009

–

1022

.

https://doi.org/10.1017/jfm.2022.306

17.

M. K.

Dewangan

,

U.

Ghosh

,

T.

Le Borgne

, and

Y.

Méheust

, “

Coupled electrohydrodynamic transport in rough fractures: A generalized lubrication theory

,”

J. Fluid Mech.

942

,

A11

(

2022

).

https://doi.org/10.1016/S0017-9310(01)00335-0

18.

C.-C.

Wang

and

C.-K.

Chen

, “

Forced convection in a wavy-wall channel

,”

Int. J. Heat Mass Transfer

45

,

2587

–

2595

(

2002

).

https://doi.org/10.1088/0022-3727/38/10/005

19.

G.

Croce

and

P.

D'Agaro

, “

Numerical simulation of roughness effect on microchannel heat transfer and pressure drop in laminar flow

,”

J. Phys. D: Appl. Phys.

38

,

1518

(

2005

).

https://doi.org/10.1039/B302569A

20.

I.

Glasgow

and

N.

Aubry

, “

Enhancement of microfluidic mixing using time pulsing

,”

Lab Chip

3

,

114

–

120

(

2003

).

https://doi.org/10.1021/acs.iecr.3c04092

21.

S.

Das

,

V. B.

Vanarse

, and

D.

Bandyopadhyay

, “

Electroosmotic micromixing in physicochemically patterned microchannels

,”

Ind. Eng. Chem. Res.

63

,

5312

–

5325

(

2024

).

https://doi.org/10.1063/5.0218566

22.

M.

Khatibi

,

S. N. A.

Sumit Kumar Mehta

, and

P. K.

Mondal

, “

Surface charge-dependent slip length modulates electroosmotic mixing in a wavy micromixer

,”

Phys. Fluids

36

,

073105

(

2024

).

https://doi.org/10.1115/1.2188507

23.

A.

Haji-Sheikh

, “

Fully developed heat transfer to fluid flow in rectangular passages filled with porous materials

,”

J. Heat Mass Transfer

128

,

550

–

556

(

2006

).

https://doi.org/10.1115/1.2188510

24.

K.

Hooman

and

A. A.

Merrikh

, “

Analytical solution of forced convection in a duct of rectangular cross section saturated by a porous medium

,”

J. Heat Mass Transfer

128

,

596

–

600

(

2006

).

https://doi.org/10.1098/rspa.2017.0193

25.

T.

Karmakar

and

G. R.

Sekhar

, “

A note on flow reversal in a wavy channel filled with anisotropic porous material

,”

Proc. R. Soc. A

473

,

20170193

(

2017

).

https://doi.org/10.1007/s11242-020-01473-2

26.

N. F.

Okechi

and

S.

Asghar

, “

Darcy–Brinkman flow in a corrugated curved channel

,”

Transp. Porous Media

135

,

271

–

286

(

2020

).

https://doi.org/10.1063/5.0158996

27.

G. M.

Moatimid

,

A.

Sayed

, and

M. H.

Zekry

, “

Effects of uniform and periodic magnetic fields at the nonlinear stability of three magnetic fluids in porous media

,”

Phys. Fluids

35

,

074109

(

2023

).

https://doi.org/10.1063/1.1373680

28.

J.-T.

Jeong

, “

Slip boundary condition on an idealized porous wall

,”

Phys. Fluids

13

,

1884

–

1890

(

2001

).

https://doi.org/10.1007/s11242-016-0789-5

29.

M.

Faltas

,

E.

Saad

, and

S.

El-Sapa

, “

Slip–Brinkman flow through corrugated microannulus with stationary random roughness

,”

Transp. Porous Media

116

,

533

–

566

(

2017

).

https://doi.org/10.1007/s11242-009-9496-9

30.

C.

Wang

, “

Darcy–Brinkman flow over a grooved surface

,”

Transp. Porous Media

84

,

219

–

227

(

2010

).

https://doi.org/10.1063/5.0231222

31.

N.

Gupta

and

K.

Vajravelu

, “

Maximal transport of non-Newtonian fluid in an anisotropic rotating porous channel with an inclined magnetic field

,”

Phys. Fluids

36

,

093119

(

2024

).

https://doi.org/10.1063/5.0200187

32.

D.

Bhargavi

,

R.

Aich

, and

N.

Gupta

, “

Thermal enhancement of couple stress fluid flow through anisotropic porous media

,”

Phys. Fluids

36

,

043111

(

2024

).

https://doi.org/10.1016/S0017-9310(02)00032-7

33.

G.

Degan

,

S.

Zohoun

, and

P.

Vasseur

, “

Forced convection in horizontal porous channels with hydrodynamic anisotropy

,”

Int. J. Heat Mass Transfer

45

,

3181

–

3188

(

2002

).

34.

M.

Kaviany

,

Principles of Heat Transfer in Porous Media

, Mechanical Engineering Series Vol.

28

(

Springer

,

New York, NY

,

1991

).

https://doi.org/10.1115/1.3248179

35.

M.

Parang

and

M.

Keyhani

, “

Boundary effects in laminar mixed convection flow through an annular porous medium

,”

J. Heat Transfer

109

,

1039

(

1987

).

https://doi.org/10.1007/s11242-013-0124-3

36.

L.

Yu

and

C.

Wang

, “

Darcy-Brinkman flow through a bumpy channel

,”

Transp. Porous Media

97

,

281

–

294

(

2013

).

https://doi.org/10.1007/s11242-017-0865-5

37.

M.

Faltas

and

E.

Saad

, “

Three-dimensional Darcy–Brinkman flow in sinusoidal bumpy tubes

,”

Transp. Porous Media

118

,

435

–

448

(

2017

).

https://doi.org/10.1016/j.euromechflu.2008.08.003

38.

C.

Wang

, “

Darcy–Brinkman flow past a two-dimensional screen

,”

Eur. J. Mech.-B/Fluids

28

,

321

–

327

(

2009

).

https://doi.org/10.1007/s11242-021-01677-0

39.

N. F.

Okechi

, “

Stokes–Brinkman flow in a rough curved channel

,”

Transp. Porous Media

139

,

513

–

526

(

2021

).

https://doi.org/10.1007/s11242-016-0818-4

40.

E.

Marušić-Paloka

and

I.

Pažanin

, “

On the Darcy–Brinkman flow through a channel with slightly perturbed boundary

,”

Transp. Porous Media

117

,

27

–

44

(

2017

).

https://doi.org/10.1007/s11242-013-0170-x

41.

D. D.

Gray

,

E.

Ogretim

, and

G. S.

Bromhal

, “

Darcy flow in a wavy channel filled with a porous medium

,”

Transp. Porous Media

98

,

743

–

753

(

2013

).

https://doi.org/10.1098/rspa.2017.0234

42.

B.

Tavakol

,

G.

Froehlicher

,

D. P.

Holmes

, and

H. A.

Stone

, “

Extended lubrication theory: Improved estimates of flow in channels with variable geometry

,”

Proc. Math. Phys. Eng. Sci.

473

,

20170234

(

2017

).

43.

M. K.

Dewangan

and

S.

Datta

, “

Flow through microchannels with textured walls: A theory for moderately slow variations

,” in

International Conference on Nanochannels, Microchannels, and Minichannels

(

American Society of Mechanical Engineers

,

2018

), Vol.

51197

, p.

V001T06A001

.

https://doi.org/10.1017/S0022112076000906

44.

L.

Hocking

, “

A moving fluid interface on a rough surface

,”

J. Fluid Mech.

76

,

801

–

817

(

1976

).

https://doi.org/10.1103/PhysRevLett.64.2269

45.

D.

Einzel

,

P.

Panzer

, and

M.

Liu

, “

Boundary condition for fluid flow: Curved or rough surfaces

,”

Phys. Rev. Lett.

64

,

2269

(

1990

).

https://doi.org/10.1063/1.1560925

46.

C.

Wang

, “

Flow over a surface with parallel grooves

,”

Phys. Fluids

15

,

1114

–

1121

(

2003

).

https://doi.org/10.1016/j.euromechflu.2018.02.004

47.

M. K.

Dewangan

and

S.

Datta

, “

Flow through microchannels with topographically patterned wall: A spectral theory for arbitrary groove depths

,”

Eur. J. Mech.-B/Fluids

70

,

73

–

84

(

2018

).

https://doi.org/10.1017/jfm.2020.134

48.

M. K.

Dewangan

and

S.

Datta

, “

Effective permeability tensor of confined flows with wall grooves of arbitrary shape

,”

J. Fluid Mech.

891

,

A12

(

2020

).

https://doi.org/10.1007/s00162-023-00679-6

49.

M. K.

Dewangan

, “

Investigation of stokes flow in a grooved channel using the spectral method

,”

Theor. Comput. Fluid Dyn.

38

,

39

–

59

(

2024

).

50.

I.

Pop

and

D. B.

Ingham

,

Transport Phenomena in Porous Media II

(

Elsevier

,

2002

).

https://doi.org/10.1115/1.2909526

51.

J.

Khodadadi

, “

Oscillatory fluid flow through a porous medium channel bounded by two impermeable parallel plates

,”

J. Fluids Eng.

113

,

509

–

511

(

1991

).

https://doi.org/10.1017/S0022112069000322

52.

C. K.

Tam

, “

The drag on a cloud of spherical particles in low Reynolds number flow

,”

J. Fluid Mech.

38

,

537

–

546

(

1969

).

https://doi.org/10.1017/S002211207200120X

53.

T. S.

Lundgren

, “

Slow flow through stationary random beds and suspensions of spheres

,”

J. Fluid Mech.

51

,

273

–

299

(

1972

).

https://doi.org/10.1017/S0022112074002503

54.

I.

Howells

, “

Drag due to the motion of a Newtonian fluid through a sparse random array of small fixed rigid objects

,”

J. Fluid Mech.

64

,

449

–

476

(

1974

).

https://doi.org/10.1021/ac0257389

55.

A. D.

Stroock

,

S. K.

Dertinger

,

G. M.

Whitesides

, and

A.

Ajdari

, “

Patterning flows using grooved surfaces

,”

Anal. Chem.

74

,

5306

–

5312

(

2002

).

https://doi.org/10.1088/1361-6463/ad2007

56.

M. K.

Dewangan

and

T.

Persoons

, “

Electromagnetohydrodynamic flow through a periodically grooved channel

,”

J. Phys. D: Appl. Phys.

57

,

165002

(

2024

).

57.

W. E.

Langlois

and

M. O.

Deville

,

Slow Viscous Flow

(

Springer

,

1964

), Vol.

173436

.

58.

L. G.

Leal

,

Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes

(

Cambridge University Press

,

2007

), Vol.

7

.

59.

H.

Ockendon

and

J. R.

Ockendon

,

Viscous Flow

(

Cambridge University Press

,

1995

), Vol.

13

.