## Abstract

A hybrid numerical–analytical solution is developed for laminar flow development in a parallel plate duct partially filled with porous media. The integral transform method is employed in combination with a single domain reformulation strategy for representing the heterogeneous media within the channel. A novel eigenfunction expansion basis is proposed, including abrupt spatial variations of physical properties due to the domain transitions. The introduction of the new basis allows for a solution with similar convergence rates as in previous applications with simpler formulations, as demonstrated through a careful convergence analysis of the expansions. The inherent automatic error control characteristic of the integral transforms approach then provides benchmark results for the developing velocity profile. Moreover, a physical analysis further verifies the consistency of both the proposed expansion and the mixed symbolic–numerical code developed. A detailed verification with a finite-element commercial code is also performed.

## 1. Introduction

The large number of practical applications involving the interaction between porous media and fluid flow, e.g. groundwater contamination [1,2], enhanced oil recovery [3], redox flow batteries and fuel cells [4,5], gives rise to a need for reliable modelling and numerical simulation schemes for both scientific and engineering usage.

Over the years, several efforts were devoted to develop better modelling strategies for the flow in adjacent fluid and porous regions. Such works have been particularly focused on the non-trivial modelling of the interface between free flow and flow inside the porous medium. Some of the earlier works approached this physical situation with a two-domain model featuring Darcy equation in the porous region and Navier–Stokes equations for the free fluid [6–8]. The coupling between the two equations was accomplished with a semi-empirical slip boundary condition imposed at the interface between the fluid and porous domains.

Adopting the Brinkman viscous correction to Darcy's Law [9] allows for the employment of continuity boundary conditions at the interface between a fluid and a porous medium. Neale & Nader [10] have shown that the Brinkman model coupled with the Navier–Stokes equations through continuity boundary conditions rendered equivalent results to the employment of Darcy model. This equivalency was guaranteed only if the adopted slip coefficient is proportional to the square root of the effective viscosity within the porous medium.

In the context of linear stability of natural convection in superposed fluid and porous layers, a critical comparison between the two main approaches has already been reported in the literature [11]. Additionally, an interesting alternative proposing a single domain formulation of the physical problem was included. Instead of modelling the two regions individually, a single set of equations with different properties defined within the appropriate bounds of the fluid and the porous media was used. The single domain formulation showed a considerable deviation of the stability curves from the two-equation approaches commonly used. However, in a later work, it has been determined that such deviation was due to the localized effects of the abrupt variation of physical properties at the interface between the fluid and the porous medium [12].

The capabilities of discrete numerical methods, such as finite-volume and finite-element methods, to solve a large range of physical problems of high complexity have made these techniques ubiquitous in science and technology calculations. Nevertheless, benchmark results are still needed for the verification of the computational codes developed, in the light of the inherent numerical errors associated with the discrete methodologies. Purely analytical methods can provide benchmark results to any prescribed accuracy, but they are usually restricted to simple linear problems in regular geometries. In this context, with the goal of extending the applicability of analytical-type techniques, the so-called hybrid numerical–analytical methods have been proposed to solve more involved problems along the last few decades. Among these methods, the generalized integral transform technique (GITT) has proved to be a fairly general methodology, able to solve diffusion and convection–diffusion problems with variable properties, moving boundaries, nonlinear source terms, etc. The main advantages of this methodology are the automatic error control and the mild increase in computational cost for an increasing number of independent variables. A complete description of the method and a thorough revision of the applications can be found in a few reference works [13–16]. Applications of GITT to convection–diffusion problems involving porous media are abundant [11,12,17]. Regarding the solution of fluid flows, the GITT has been applied to both the primitive variables [18–21] and streamfunction-only [22–31] formulations. The latter was usually preferred due to the automatic satisfaction of the continuity equation and better convergence characteristics, though directly applicable only to two-dimensional problems. In the case of three-dimensional formulations, the scalar and vector potential formulations were also considered [32], which, however, lack the advantage of collapsing the momentum equations into a single fourth-order ODE system, as achievable in the streamfunction formulation for the two-dimensional case.

In recent years, a new solution strategy within the GITT framework, known as the single domain formulation, has been under development to enable a more straightforward handling of convection–diffusion problems in heterogeneous media and complex geometries [33–38]. The single domain formulation proposes the treatment of different media and irregular domains through abrupt variations of physical properties and source terms in the formulation. As mentioned before, in a certain sense, this idea has already been applied to evaluate the linear stability in cavities partially filled with a porous medium [11,12]. Nevertheless, the single domain formulation has evolved into a more systematic methodology, allowing for analytical solution in various cases and increasing the overall robustness of the integral transforms approach.

The present work aims at advancing the combination of the single domain formulation with the integral transform method to solve the Navier–Stokes equations. In order to achieve such goal, a new eigenvalue problem containing information about the abrupt variation of the physical properties is proposed, whose solution is then used as a basis for the eigenfunction expansion of the velocity field vector. A new interpretation to the eigenfunction expansion is provided, through which the velocity vector field can be determined considering the influence of an infinite number of vortices disturbing a base flow. This proposition automatically recovers the streamfunction formulation as a special case for the two-dimensional case. The methodology is applied to a parallel plate channel partially filled with a porous material, under hydrodynamically developing laminar flow conditions. A detailed verification with a finite-element analysis (FEA) is performed and benchmark results for the physical situation under study are offered. Finally, the new perspectives that are offered through the present proposal for the application of integral transforms are discussed, especially in the application of GITT to three-dimensional flows in a more effective way.

## 2. Formulation and solution methodology

### (a) Fluid flow model

The physical situation proposed to be studied in this work is illustrated in figure 1. The picture shows a parallel plate channel partially filled with a porous medium in a symmetrical manner. The flow direction and its uniformity at the entrance are also indicated. Dimensions of interest and a conveniently positioned Cartesian coordinate system are included to facilitate the understanding of the mathematical model.

The model for the fluid flow inside the porous medium is chosen to be the Darcy–Brinkman [9], while within the free fluid the Navier–Stokes equations for steady-state incompressible flow are employed. The coupling at the interface between the two regions is proposed assuming continuity of the velocity vector, the normal stress component and the tangential stress component [11]. The models for the two regions could be solved separately and coupled by the common boundary between the different domains. However, this coupling can be quite cumbersome when complicated geometries or many different media are involved. With the main advantage of accomplishing this coupling automatically and more conveniently for computational purposes, the single domain formulation is here employed. Instead of designating one model for each different region, a single generalized model with spatially variable physical properties is proposed. The model deemed appropriate to the physics under consideration is shown as
** u** is the dimensionless velocity vector,

*ρ*is the dimensionless density,

*p*is the dimensionless pressure field,

*μ*is the dimensionless dynamic viscosity,

*Re*is the Reynolds number based on the hydraulic diameter and

*Da*is the Darcy number (table 1).

The dimensionless quantities can be obtained from their dimensional counterparts in the following way:
*a*–*g*
where *u*_{0} is the uniform entry longitudinal velocity component, *h* is half the height of the channel, *ρ** is the density, *ρ*_{0} is the fluid density, *p** is the pressure field,

Following the single domain formulation strategy, the properties used in the model of equation (2.2) present spatial variation. Let *V*_{f} and *V*_{p} be the regions occupied by the free fluid and the porous medium, respectively. In consonance with the choice of working with the intrinsic pressure, instead of the superficial pressure, during the volume averaging process, an effective viscosity inversely proportional to the porosity was applied to the porous region [39]. Such pressure can be more closely related to the pressure in the free fluid flow, allowing for a more straightforward single domain formulation. Additionally, neglecting the inertial effects within the porous medium, and disregarding any dissipative term in the free fluid flow, the following definitions are adopted:
*ε* is the porosity of the porous medium.

To enable the solution of the aforementioned model, a set of boundary conditions is needed. The use of the Brinkman correction in the Darcy equation allows the imposition of a no-slip boundary condition at the wall contour illustrated in figure 1. At the entrance, limited to the region indicated with the arrows in figure 1, a uniform velocity profile was imposed. The outlet boundary condition is prescribed assuming that the channel is long enough for the flow to be fully developed. Hence, a dimensionless fully developed velocity profile *F*(*y*), to be detailed in a later section, is considered at the outlet. The adopted single domain approach eliminates the need to explicitly specify boundary conditions at the interface between the free fluid and porous layers. The appropriate continuity boundary conditions naturally arise from the variations of the physical properties defined in equations (2.4)–(2.6). Mathematically, the boundary and interfacial conditions can be summarized as follows:
*a,b*
*c,d*
*e*
*f*
where *x*_{o} is the dimensionless channel length, *f* is the velocity profile at the entry of the channel and *n** _{i}* is a vector normal to the interface between free fluid and porous medium.

### (b) Filter solution

Integral transform methodologies can strongly benefit from the employment of simple analytical solutions to simplified versions of the original model as filters [13–16]. The aim is to extract information from equation and boundary source terms, or even from linearized versions of nonlinear differential operators, so as to reduce their influence on slowing down convergence of the proposed eigenfunction expansions. For the particular physical situation being studied in this work, a fully developed velocity profile naturally appears as the longitudinal coordinate *x* grows, as an equivalent to a steady-state solution in an actually transient problem. The fully developed velocity profile not only identically satisfies equation (2.1) and the r.h.s. of equation (2.2), but also perfectly represents the outflow boundary condition imposed, as described in the previous section.

The fully developed velocity profile was obtained from equation (2.2) in a Cartesian coordinate system consistent with the one illustrated in figure 1, with the simplification of a longitudinally invariant flow. The resulting equation consists of an ordinary differential equation with the constant longitudinal pressure gradient as the source term. The multi-layer characteristic of the physical configuration here being analysed is introduced through a transversally variable viscosity. In order to obtain an analytical solution to the fully developed velocity profile conveniently, different equations with constant viscosities are associated with each domain, coupled by their common boundaries with continuity conditions. Hence, the set of coupled ordinary differential equations to be solved for the filter is shown as follows:
*a,b*
*c,d*
*e,f*
where *w*_{p} is the ratio between the thickness of the porous medium (*h*_{p} in figure 1) and half of the channel height (*h* in figure 1).

The filter velocity vector is then given as follows:

where

The filter velocity vector is then employed to compose the unknown velocity field in the following form:

Finally, the filtered continuity and Navier–Stokes equations that govern the unknown filtered velocity field,

### (c) Eigenvalue problem

The underlying assumption used in the solution of equations (2.12) and (2.13) is that the velocity vector field can be determined by taking into account the influence of an infinite number of vortices disturbing a base flow. This physical interpretation can be translated into the following relation:

The new perspective embodied by the above equation is, in fact, equivalent to the streamfunction-only formulation in two-dimensional cases when a base vector *z*-axis is used [22–31], and shares the same main advantage of identically satisfying the continuity equation. Such equivalency is further explained in electronic supplementary material, S1. However, the interpretation without recurring to the streamfunction definition allows for a more straightforward generalization to three-dimensional problems, provided that a proper base for the expansion (

Equation (2.14) was developed having in mind the integral transformation of the Navier–Stokes equations in all but one of the spatial coordinates. In the present case, the longitudinal coordinate (*x*) is chosen to be handled numerically, due to the convenience of the transformation in the transversal direction (*y*), with its homogeneous boundary conditions. One can reason that the choices for the two functions displayed in the right-hand side of equation (2.14) cannot be arbitrary. From the linearity of the curl operator, the resulting infinite series within the differential operator must be able to span the whole space of possible solutions for the velocity vector. One possibility is to impose that

Besides demanding orthogonality, another desirable feature of the basis is the presence of the largest possible amount of information on the coefficients and operators of the partial differential equations to be solved, which is also a common characteristic in previous GITT solutions [13,15]. In the single domain framework, the solution methodology strongly benefits from the inclusion of the abrupt variation of the physical properties that represent the media transitions. Therefore, some modifications from the classical biharmonic eigenvalue problem [22–31] are here proposed. The eigenvalue problem selected for the present situation is then written as follows:
*a*
*b,c*
*d,e*
*f,g*

The basis for the expansion can then be defined as follows:

With a simple application of an integral operator with the proper kernel in equation (2.15*a*) and the boundary conditions of equation (2.15*b–e*), the following orthogonality property can be derived:

Direct analytical solution to equation (2.15*a–e*) is unlikely to be obtainable in terms of known functions. Therefore, the GITT itself was used to solve this eigenvalue problem, even though other methodologies can be found elsewhere [40], making use of a simpler auxiliary problem, in a similar procedure already applied in previous work on the single domain formulation with integral transforms [33–38]. In order to employ the GITT procedure in solving equation (2.15*a*), the classical biharmonic eigenvalue problem was selected in the following form:
*a*
*b,c*
*d,e*
*f,g*

Considering the transversal symmetry of the physical situation depicted in figure 1, equations (2.18*a–e*) admit the following solution:

The eigenfunction of equation (2.19) bears the orthogonality property as follows:

Such orthogonality property of the auxiliary eigenfunction allows for the following transform-inverse pair to be formed:
*a,b*

The integral transformation process proceeds operating on equation (2.15*a*) with *a*
*b–d*
*e*
*d*

Truncating the algebraic eigenvalue problem of equation (2.23*a*) to a finite order *M* allows the employment of reliable routines available to obtain its eigenvalues and corresponding eigenvectors. The inverse formula of equation (2.22*b*) is then used to obtain the required eigenfunctions, *Mathematica* v. 10.4 platform was used for both the numerical task of solving the eigensystem of equation (2.23*a*) and the symbolic computation of the integrations in equation (2.23*e,f*).

### (d) Transformed problem

The final step in the integral transform methodology consists in operating both the governing equations and the remaining boundary conditions with the appropriate integration kernels. The proposal for the velocity vector presented in equation (2.14) is able to identically satisfy the filtered continuity equation of equation (2.12). Hence, the integral transform needs to be employed only in the handling of the linear momentum conservation equation, i.e. equation (2.13).

Let d*V* be an infinitesimal volume with contour ∂*V* consisting of the whole cross section of the channel depicted in figure 1 and with an infinitesimal longitudinal length d*x*. Applying ** b** and

**are class C**

*c*^{1}vector fields, we then have

Note the disappearance of the pressure gradient term in equation (2.25). The use of the convenient kernel for the integral transform automatically eliminates the pressure term by imposing a curl operation to the pressure gradient. The absence of any pressure-related term avoids the need for proposing an eigenfunction expansion to the pressure field, further simplifying the solution of the flow problem. Further developing the integrations in equation (2.25), as detailed in the electronic supplementary material, S2, leads to the so-called transformed problem as follows:
*a*
with the integral coefficients
*b*
*c*
*d*
*e*
*f*
*g*
*h*
*i*

The boundary conditions at the longitudinal ends of the domain under study must be transformed in a somewhat different way. Let *f*(*y*) be the function representing the entry flow depicted in figure 1. Assuming the linearity of the curl operator holds for equation (2.14), the following must be true:
*a*
and
*b*
where *x*_{o} is the ratio between the length of the channel *L** and *h*. Further derivations detailed in the electronic supplementary material, S3 lead to the following boundary conditions:
*a*
*b*
*c*
*d*

Operating equation (2.28*a–d*) with *a*
*b*
*c,d*

### (e) Solution procedure

Owing to the impossibility of solving infinite sets of equations such as equation (2.23*a–f*) and equation (2.26*a–i*), truncation of the systems must be carried out. The number of terms in both systems is directly governed by how many terms are included in the summations of equations (2.14) and (2.22*b*). Here, the truncation orders are *M* and *N* for the summations in equations (2.22*b*) and (2.14), respectively. These truncation orders are the only parameters to be monitored as far as the convergence is concerned.

Before proceeding with the solution of the transformed problem, the integrations presented in equation (2.26*b–i*) were analytically determined. Substituting the inverse formula of equation (2.22*b*) into the integration formulae and rearranging the results leads to integrations involving the simpler auxiliary eigenfunctions from equation (2.19). These operations were carried out analytically using symbolic computation routines built in the commercial software *Mathematica* v. 10.4. The same process was followed for the integration of equation (2.29*a*). Once numerical values for the coefficients of equation (2.26*a*) and for the boundary condition of equation (2.29*a*) were evaluated, the system was solved using a dedicated boundary value problem solver based on the robust PASVA3 algorithm [41]. This routine was programmed in the *Mathematica* environment, allowing for full integration of the developed code into a single file. However, the PASVA3 algorithm requires a pre-treatment of the system, because the routine has been devised to handle first-order differential equations [41]. A new definition of the dependent variables, described in detail in previous works [13,22,23], was also employed here. In the process, the *N *× *N* system of equation (2.26*a*) becomes a (4*N*) × (4*N*) system of first-order ordinary differential equations. A tolerance of 10^{−4} was fed to the PASVA3 algorithm in every case run to serve as stopping criteria for the routine.

The dimensionless length of the channel was such that the following relation holds for every case studied:

Equation 2.30 was shown to satisfactorily comply with the approximation of invariable velocity profile along the longitudinal direction at the outlet, as will be demonstrated in the results of the next section.

### (f) Test cases

The integral transform solution methodology was tested in two different ways. The first one is aimed at investigating the robustness of the method under different Reynolds numbers. In a second moment, the influence of the Darcy number is tested for a fixed Reynolds number. The values employed for both Darcy and Reynolds numbers follow from typical values found in the application to small-scale (*h* ∼ 100 µm) redox flow batteries [5]. Moreover, the spanning of the Darcy number over 3 orders of magnitude allows for a thorough analysis of the robustness of the proposed solution methodology. All the other relevant parameters remained constant throughout the test cases. Table 2 summarizes the parameters used and also defines the different combinations of Reynolds and Darcy numbers.

## 3. Results and discussion

### (a) Eigenvalue problem

Before proceeding with the solution of the flow equations, the basis for the expansion of equation (2.14) must be analysed. Therefore, a careful convergence analysis was performed *a priori* for the test cases under consideration, by increasing the truncation order *M* until certain convergence criteria were satisfied. Analysing the eigenvalue problem of equation (2.15*a–g*), the only parameter that varies along the different test cases is the Darcy number. Hence, convergence of the eigenvalue problem solution is first demonstrated for all three Darcy numbers used in the calculations.

Firstly, the graphical convergence with increasing *M* was investigated for the most critical case (*Da* = 0.002). Figure 2 illustrates how the 18th eigenfunction behaves, sorted in ascending order of the eigenvalues, as the truncation order is raised with increments of three terms per run. Clearly, almost no difference can be noted for truncation orders larger than 21.

To further evaluate the convergence of the eigenvalue problem solution, table 3 presents the numerical values of the 18th eigenvalue for all three different Darcy numbers. The choice for the 18th eigenvalue is based on the assumption that it is the most critical one among the ones effectively used to converge the velocity vector, as will be discussed further in detail. Convergence to five significant digits was achieved with truncation orders not larger than 42. Lower values for *M* are sufficient to attain similar convergence as the Darcy number is increased. The maximum truncation order calculated for each Darcy number is used to obtain the eigenfunctions and eigenvalues needed for the transformed problem of the flow equations.

### (b) Flow development for different Reynolds numbers

Once the eigenfunctions convergence has been verified, we may now proceed to the solution of the transformed problem for the flow equations, and to the convergence analysis of the eigenfunction expansion proposed in equation (2.14). With a fixed value of *M* = 39 (table 3) and for the Darcy number of 0.02, the convergence of the velocity vector solution is then analysed in terms of the truncation order *N*. The behaviour of the longitudinal velocity component at the channel centreline, for increasing truncation order *N*, is now examined for all three Reynolds numbers. The results are presented for selected longitudinal positions in tables 4–6. Full convergence to at least the third significant digit was attained in all cases for *N* < 18, which is comparable to the convergence rates obtained for the free fluid flow problem between parallel plates with a similar methodology [23]. The longitudinal velocity component converges faster at positions further away from the entry. This fact underscores the importance of the filter solution proposed because the velocity vector tends asymptotically towards the fully developed profile as *x* grows. To verify the criteria for the truncation of the channel length using equation (2.30), results considering a channel twice as long as obtained from this expression are also included in tables 4–6. These results were obtained for the largest truncation order *N* used for the same case. The agreement is perfect within the precision shown, indicating invariance with the channel length beyond the point set by equation (2.30). For co-verification purposes, results from FEA obtained with the commercial software *COMSOL Multiphysics* [42] are also presented. The agreement is almost perfect, with a relative deviation never larger than approximately 0.2% for all selected positions.

Figure 3 depicts the evolution of the centreline longitudinal velocity component with the following modified longitudinal coordinate, in a semi-log scale:

The main reason for the definition of equation (3.1) was the increase of the development length for flows with higher Reynolds numbers. With the introduction of *x*^{+}, all relevant phenomena can be placed into a single interval, facilitating the analysis.

Results for both the integral transforms and the FEA solutions are shown in figure 3. The agreement between the two methodologies for all the range of the longitudinal position and Reynolds numbers, provides additional verification of the present implementation, while also demonstrating the potential of GITT solutions in the establishment of reference results.

The longitudinal velocity component at the centreline undergoes an undershoot for the higher Reynolds number studied (*Re* = 300). This higher energy flow allows for a more pronounced penetration of fluid into the porous medium at the early stages of flow development along the channel, consequently, decreasing the longitudinal velocity component at the centreline. The reverse occurs for the lower Reynolds studied (*Re* = 50), while a somehow neutral behaviour is observed for the case with Reynolds number equal to 150. As *x*^{+} grows, the three curves tend to collapse into one, in accordance with the attainment of the fully developed region. The value of *x*^{+} at which all curves fuse together was identified as being approximately 0.013. This value was considered in the development of equation (2.30) from an *a posteriori* analysis, aimed at minimizing computational cost.

### (c) Flow development for different Darcy numbers

The convergence of the eigenfunction expansion in equation (2.14) was also investigated for different values of the Darcy number. Results for the longitudinal velocity component at the centreline for a Reynolds number equal to 50 are presented in tables 7 and 8, with *Da* = 0.002 and 0.2, respectively, for five different longitudinal positions. Convergence of at least three significant digits was observed, similarly to the behaviour shown in tables 4–6. The same improved convergence as the longitudinal coordinate grows was observed, with the filter solution analytically capturing a significant parcel of the flow field. Results for a channel with twice the length are also presented in tables 7 and 8. The largest truncation order *N* used for each case was set for the test of the channel length. The agreement is perfect for all digits shown, corroborating the adequacy of truncation length expression of equation (2.30). Comparison with the *COMSOL* results, also present in tables 7 and 8, leads to the same conclusion about the adequacy of the integral transform approach in providing reliable benchmark results for the verification of more general codes. Once more, the relative deviations between the FEA and the hybrid approach here proposed fall below approximately 0.2%.

The longitudinal velocity component profiles along the transversal direction, for Reynolds number equal to 50, are shown in figure 4*a–d*. Four different longitudinal positions are studied for three different Darcy numbers. The adherence with the FEA results from *COMSOL Multiphysics* [42] is perfect to the graph scale, further contributing to the co-verification effort. As the longitudinal coordinate grows, the effects of the different Darcy numbers become more evident. As expected, for the larger Darcy number, more fluid tends to penetrate and flow through the porous medium. Clearly, the continuity of the velocity vector is respected via the single domain formulation. The continuity of both the normal and tangential stresses is harder to identify solely based on the graphs of figure 4. With increasing Darcy number, a clear tendency towards the parabolic velocity profile typical of Poiseuille flow between parallel plates is observed in figure 4*d*. Such behaviour is expected, because when the Darcy number tends to infinity (the free fluid flow condition), the corresponding parabolic fully developed profile is recovered.

Results for the transversal velocity component profile along the transversal coordinate are presented in figure 5*a–d*, for the same four different longitudinal positions and the three different values of Darcy number (*Da* = 0.002, 0.02 and 0.2). Again, the *COMSOL Multiphysics* [42] results are in perfect agreement with the integral transform results here obtained. The tendency to further penetrate the porous medium for higher Darcy numbers is again observed when analysing the transversal velocity component profiles. For *Da* = 0.002, the transversal velocity lies majorly in the negative region, consistent with a flow towards the centreline of a developing boundary layer. In other words, with decreasing Darcy numbers, the porous medium becomes increasingly like an impermeable wall. The inverse behaviour is shown by the other extreme among the simulated Darcy numbers. In the results for *Da* = 0.2, the transversal velocity was always positive (towards the wall), indicating little resistance to the flow at the interface between free fluid and porous medium. For the intermediate Darcy number, a somewhat mixed regime was observed, with the transversal velocity values majorly positive, but showing negative values at the early stages of flow development.

## 4. Conclusion

The present work advances the integral transform method, in combination with the single domain formulation strategy, to successfully solve steady-state fluid flow in heterogeneous media. Benchmark results for the physical situation involving developing flow in a parallel plate channel partially filled with a porous medium were obtained, and then employed for co-verification with a well-tested FEA algorithm. Physical reasoning with the velocity profiles variation in both longitudinal and transversal directions offered further evidence of the robustness and reliability of the proposed methodology to solve fluid flow problems in channels partially filled with porous media.

The use of abrupt variations of the physical properties to represent the different material layers within the domain allowed for the avoidance of complicated models involving different partial differential equations coupled through their common domain interfaces. The inclusion of information about the abrupt properties variation into a novel eigenvalue problem proposal rendered fast convergence rates comparable to those previously observed for free fluid flow between parallel plates with a similar integral transform methodology.

In contrast with the more commonly applied streamfunction-only formulation for integral transform solutions of the Navier–Stokes equations, the new interpretation of the velocity eigenfunction expansion offered here allows for a more straightforward generalization of the integral transforms methodology to three-dimensional fluid flow problems governed by the Navier–Stokes equations.

## Data accessibility

This work does not have any experimental data.

## Authors' contributions

All authors conceived the mathematical model, interpreted the computational results and wrote the paper. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

The authors are grateful for the financial support offered by the Brazilian Government agencies CNPq (project no. 401237/2014-1), CAPES and FAPERJ.

## Acknowledgements

K.M.L. is grateful to Prof. Dimos Poulikakos and the Laboratory of Thermodynamics in Emerging Technologies at ETH Zürich, where part of the work was conducted.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3955894.

- Received September 8, 2017.
- Accepted December 1, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.