2.1. Triply periodic minimal surface porous media

Among various TPMS porous media (e.g. Lord & Mackay Reference Lord and Mackay2003; Al-Ketan & Abu Al-Rub Reference Al-Ketan and Abu Al-Rub2019), we choose three common geometries: (i) gyroid – represents a complex internal architecture typical of TPMS porous media; (ii) ‘primitive’ – has a distinctive converging–diverging geometric shape; (iii) TPMS-based body-centred cubic (BCC) – owing to its open cell structure. Note that primitive is somewhere between gyroid and BCC in terms of complexity and open cell type.

2.1.1. The TPMS design

There are various ways to design TPMS topologies, and we adopt nodal surface approximation (e.g. Zhao et al. Reference Zhao, Liu, Fu, Zhang, Zhang and Zhou2018; Al-Ketan & Abu Al-Rub Reference Al-Ketan and Abu Al-Rub2019) because of the available mathematical expression, which enables us to produce customised topologies using MATLAB. In this method, either a solid or a sheet type topology is designed by choosing a level-set constant $c$ , which allows us to control the porosity $\phi$ by controlling the sheet thickness with different values of $\pm c$ . The following equations, respectively, for gyroid, primitive and BCC, which use the implicit method to create zero-valued surfaces solving nodal equations, are the most commonly used mathematical expressions to create TPMS porous media:

(2.1a) \begin{align} \cos (k_xx)\sin (k_yy)+\cos (k_yy)\sin (k_zz)+\cos (k_zz)\sin (k_xx)-c=0, &{\,\rm \quad[gyroid]} \end{align}

(2.1b) \begin{align} \cos (k_xx)+\cos (k_yy)+\cos (k_zz)-c=0, &{\,\rm \quad[primitive]} \end{align}

(2.1c) \begin{align} \begin{split} \cos (2k_xx)+\cos (2k_yy)+\cos (2k_zz)-2[\cos (k_xx)\cos (k_yy)&\\ +\cos (k_yy)\cos (k_zz)+\cos (k_zz)\cos (k_xx)]-c=0. &{\,\rm \quad[BCC]} \end{split} \end{align}

Here $x$ , $y$ and $z$ are spatial coordinates in the 3-D Cartesian system, with $k_x$ , $k_y$ and $k_z$ defined by $k_i =2\pi {n_i}/{L_i}$ ( $i=x,y,z$ ) where $n_i$ is the number of unit cells and $L_i$ is the total lattice length in the $i$ direction. Hence, the length of each unit cell is $l_i = L_i/n_i$ , where we keep the cell size the same in all directions, $l=l_i$ (e.g. Peng & Tran Reference Peng and Tran2020). In practice, (2.1a ), (2.1b ) and (2.1c ) are solved using MATLAB, and its built-in functions of isosurface and isocaps are used to generate desired surfaces of TPMS porous media satisfying the above relationship. Figure 1(a), 1(b) and 1(c) show the unit cell geometry of gyroid, primitive and BCC, respectively, whereas the corresponding figure 1(d), 1(e) and 1( f) show a $4\times 4$ cell structure.

2.1.3. Manufacturing of TPMS porous media for experiments

Three different 3-D printing manufacturing processes are used to obtain TPMS geometries of the different $\phi$ and $l$ as described in table 1. We divide the experiments into two sets: Set 1 and Set 2. Set 1 has $l=12$ mm, and is manufactured using a specialised 3-D printer, Stratasys Object 260 Connex3. Using this 3-D printer, the minimum unit cell size that can be obtained without being damaged from harvesting and removal of supporting material is approximately 12 mm. With this optimised length, a unit cell of each TPMS porous media, 12 mm $\times$ 12 mm $\times$ 12 mm is designed, and based on this unit cell size, the values of $c$ are optimised for each porous media to create targeted porosity namely, $\phi$ = 0.85, 0.70 and 0.55.

Table 1.

Details of the TPMS porous media used in this study. Note that the design of G4 and Metal-G4 are the same; however, the former is manufactured using SLA and the latter by SLM as shown in figure 14(b) and 14(c). The permeabilities $K_{\!Q1}$ and $K_{\!Q2}$ (cf. (1.1)) are discussed in § 3.2. The top Set 1 has $l=12$ mm fixed, whereas the bottom Set 2 has $\phi \approx 0.85$ and varying $l$ . Note that the first and last row are the same.

Set 2 cases in table 1 have a nominally fixed $\phi \approx 0.85$ and varying $l$ , and are denoted G1 to G4. Stereolithography (SLA), a fine-resolution technology, is used to obtain $l=$ $4$ , $6$ and $9$ mm with $\phi \approx 0.85$ . Computational drawings of the three TPMS geometries are shown in figure 2, and photographs of each sample for both sets (cf. table 1) are presented in figure 3. The geometry ‘Metal-G4’ (with $l=9$ mm) is an aluminium TPMS that is 3-D printed with selective laser melting (SLM), and has the same bulk parameters as G4 that is printed in plastic with SLA. The geometry of Metal-G4 has a slightly rougher surface finish than G4 owing to the different 3-D printing process, which becomes apparent when pressure drop is measured. The roughness is estimated to be approximately 50 $\mu \textrm {m}$ . Case G5 in Set 2 is with $l=12$ mm, and is the same as the top row in Set 1.

Figure 2.

All dimensions are in millimetres. The TPMS porous media in the cylinder shape with a diameter of 20.54 mm and a length of 120 mm. Panels (a), (b) and (c) are gyroid, primitive and BCC, respectively, in isometric, side and front views. This corresponds to Set 1 (see table 1) where the unit cell size $l=12$ mm, and the TPMS are placed inside a pipe of inner diameter 20.6 mm, i.e. little less than two unit cells within the pipe cross-section. Note that for Set 2 (not shown here), the smallest $l=4$ mm, i.e. over five periodic cells within the pipe cross-section.

Figure 3.

Photographs of one of each 3-D-printed TPMS porous media (of a diameter of 20.54 mm and a length of 120 mm) with different porosities. See table 1 for details of G1, G2, etc.

Figure 4.

Schematic diagram of the experimental set-up for the pressure drop measurement. Numbers represent corresponding parts. Here $\unicode{x2460}$ Constant head tank, $\unicode{x2461}$ Flow conditioner, $\unicode{x2462}$ Entrance length, $\unicode{x2463}$ and $\unicode{x2464}$ TPMS porous media are inserted in the entire measurement section. Pressures are measured at three locations as indicated $P_1$ , $P_2$ , and $P_3$ in the inset. $\unicode{x2465}$ Needle valve at the outlet, $\unicode{x2466}$ Load cell.

To characterise fluid flow through the TPMS porous media, the porous media is inserted in an acrylic pipe of inner diameter 20.6 mm with thickness 1.98 mm. (We will describe the full set-up in § 2.2.) To fill the length of 1.56 m acrylic pipe for Set 1 (and a shorter pipe for Set 2, as described later), a cylinder-shaped porous media is required. This full 1.56 m length of the porous media is manufactured by dividing into 13 smaller pieces of length 120 mm due to the limited size of a built plate in Stratasys Object 260 Connex3 (cf. figures 2 and 3). During the design process, we particularly pay attention to maintaining the complete unit cell as multiple 120 mm pieces are inserted into the acrylic pipe. These 13 pieces are aligned in the pipe to make the 1.56 m of a single porous medium. To cover the first set, nine separate acrylic pipes are used for three different TPMS porous media with three different porosities. For the second set of experiments G1 to G5 with $l$ = $4$ , $6$ and $9$ mm (cf. table 1), which is conducted after finishing first set of experiments, a shorter pipe of 36 cm is used. The results in different length pipes make hardly any difference (discussed in § 4.3), and the shorter pipe is allowed for a more comprehensive parametric study of the gyroid geometry.

2.2. Experimental set-up

A schematic diagram of the experimental set-up for the pressure drop measurements is shown in figure 4. The set-up consists of three main sections: the inlet, the measurement section and the outlet. The inlet has a constant pressure head gravity driven system allowing water to be fed to a nozzle which includes honeycomb as flow straighteners. The fluid travels 1.3 m of empty acrylic pipe with an inner diameter of 20.6 mm before entering the pipe with the porous media inserted (for a pipe length of 1.56 m in Set 1 and 36 cm for Set 2). The entrance length of 1.3 m is such that an approximately fully developed laminar flow enters the porous media section, except for the highest flow rates that we encounter. We note that it is more important to achieve a fully developed flow in the porous section beyond the inlet pipe, which will be confirmed later in the paper.

For each porous media section of the pipe, the first pressure tap is located 2 cm from the start of the porous media section, the second and the third taps are 75 and 150 cm from the first tap. To accurately measure pressure changes at specific locations, five holes – each with a diameter of 4 mm – are made around the upper-half of the circumference to connect pressure taps, allowing for the measurement of the averaged pressure at that location. When the porous media are inserted into the acrylic pipe, extra care is taken to avoid the inserted porous media obstructing the holes. To ensure that the holes are not blocked, the silicone tubes connected to the taps are inspected and confirmed to be completely filled with water.

A manometer fabricated in-house (cf. figure 4) is carefully designed to minimise the effect of surface tension (with a relatively large manometer diameter of 35 mm), and height measurements are carried out with submillimetre resolution using a digital camera with a ruler for calibration. At the outlet, a calibrated load-cell is used to record the change of water mass with time. To detect any abnormality of the load-cell during experiments, a beaker with known weight and volume is used to collect the outlet water with a timer and later compared with the load-cell output.

For Set 2, we have considered a shorter 36 cm pipe. Tests with gyroid $l=12$ mm and $\phi =0.85$ show that the results are almost identical for the longer and 36 cm pipe. For each of the five pipes used in the Set 2 measurements, the first pressure tap is located 3 cm from the start of the gyroid, the second tap is 15 cm from the first and the last one 15 cm from the second tap. Thus, the last tap is 3 cm from the end of the porous medium, as shown in Appendix B (figure 14), where we also describe further details of the experimental set-up for Set 2. This includes special connectors for pressure tapping required for the metal gyroid – case Metal-G4 in table 1.

Figure 5.

Pressure drop within the measurement section in figure 4 without porous media. Three different symbols are three repeat experiments. (a) Pressure drops over $L=1.5\rm \,m$ at various superficial velocities ( $u_s$ ). The dashed line is the linear fitting. (b) Normalised pressure drop versus $\textit{Re}_{\textit{pipe}}$ on log–log axes and its comparison with analytical $64/\textit{Re}_{ { pipe}}$ .

2.3. Testing of pressure measurement system

Before conducting experiments on pressure drop in our porous media, we first measure the pressure drop in the acrylic pipe without the porous media. This measurement can then be compared with the analytical solution for laminar pressure drop in smooth pipe flows. By controlling the needle valve at the outlet, shown in figure 4, a range of flow rates is obtained, and the pressure drop at each flow rate is measured. The measurement is repeated three times and all the results are plotted together. Figure 5(a) shows pressure drop per unit length versus superficial velocity (which is the same as bulk velocity in this empty pipe). As the flow rates increase, the pressure drop also increases as shown in figure 5(a), and this increasing trend satisfies the expected linear relationship. Figure 5(b) shows the same data, but plotted on a log–log axes with $ C_f = (-\Delta\! P/\Delta x)\,{d_{\textit{pipe}}}/(( {1}/{2})\rho u_s^2)$ against $\textit{Re}_{\textit{pipe}} = u_s d_{\textit{pipe}}/\nu$ , where $d_{\textit{pipe}}$ is the empty pipe diameter. The data follows the laminar solution of $64/\textit{Re}_{\textit{pipe}}$ reasonably well, providing confidence in our measurement system. Note that the pressure difference across $L=1.5$ m is less than $5$ Pa m⁻¹, whereas the pressure gradient in our porous media is at least an order of magnitude higher, which further supports the adequacy of the present system.

Figure 6.

Pressure drops of TPMS porous media with three different porosities – Set 1 in table 1. (a) All pressure drop results. The fitting to extract $K_{\!Q1}$ and $K_{\!Q2}$ from (1.1) are shown separately: (b) gyroid; (c) primitive; (d) BCC. Red dashed lines () are the least square fitted (1.1) to data.

3.1. Results of pressure drop measurement

Figure 6(a) presents the experimental results of the net pressure drop ( $-\Delta\! P/\Delta x$ ) for the three TPMS geometries of Set 1 (gyroid, primitive and BCC, cf. table 1), which have a fixed unit cell size $l=12$ mm, and three different porosities of $\phi =0.85$ , 0.70 and 0.55. The colour of symbols in the figure fades as the $\phi$ decreases. The measurement is repeated five times for each porous media, and the results are plotted together to show the reliability of the results and the repeatability of the measurements. Note that in the repeat experiments we do not try to keep $u_s$ constant. The pressure drop within each TPMS porous medium increases as $\phi$ decreases. This is because a reduction in porosity implies a reduction in the void space for fluid to flow, which increases the viscous and pressure drag for a given $u_s$ and, in turn, increases flow resistance. As such, the largest pressure drops are observed at $\phi =0.55$ . In general, pressure drops in BCC are lower than those of gyroid and primitive. Compared with primitive at $\phi =$ 0.85 and 0.70 in figure 6(a), gyroid shows a larger pressure drop. However, at $\phi =0.55$ , both gyroid and primitive (and ) have almost identical pressure drops at lower $u_s$ , and then primitive has a slightly larger pressure drop as the superficial velocity increases beyond $u_s\approx 0.01$ m s⁻¹.

Figure 7.

(a) Pressure drop in gyroid with $\phi \approx 0.85$ and different $l$ – Set 2 in table 1. (b) Fitted second-order polynomial equations (1.1) to data in red dashed lines ().

Furthermore, comparing the results of each porous medium at, say $\phi =0.85$ (in figure 6 a), a gradual linear increase in pressure drop is observed until $u_s\approx 0.01$ m s⁻¹ (or less) for gyroid and primitive, and $u_s\approx 0.04$ m s⁻¹ for BCC. As flow rate increases further, pressure drop increases faster than linear, reflecting the expected change from a viscous to a pressure drag dominated regime as suggested by (1.1). A similar trend is found in other open-cell type porous media, such as metal foams (e.g. Dukhan Reference Dukhan2006; Oun & Kennedy Reference Oun and Kennedy2014) and pack beds of spherical particles (e.g. Ergun Reference Ergun1952; Lovreglio et al. Reference Lovreglio, Das, Buist, Peters, Pel and Kuipers2018).

Another way to interpret figure 6(a) is to hold a constant $-\Delta\! P/\Delta x$ , say 400 Pa m⁻¹, and move horizontally as $u_s$ increases. We notice that for the same pressure drop, gyroid and primitive with $\phi =0.55$ (and ) have the smallest bulk velocities, and then gyroid and primitive at $\phi =0.70$ (and ) and subsequently at $\phi =0.85$ (and ), which represents increasing bulk velocity with increasing porosity. The BCC structures allow the highest flow rates among the three geometries, with $\phi =0.85$ () resulting in the highest in the BCC category.

The pressure drop results for Set 2, where $\phi \approx 0.85$ and the unit cell size $l$ vary, are presented in figure 7(a). Not surprisingly, increasing $l=$ 4, 6, 9, 12 mm corresponding to G1, G2, G3 (or G4) and G5 results in a reduction of $-\Delta\! P/\Delta x$ . Note that the abscissa of figure 7(a) is larger compared with that in figure 6(a) to accommodate the increased pressure resulting from smaller $l$ . Interestingly, comparing the pressure drop results of porous media manufactured using SLA technique with plastic (G4) with SLM in metal (Metal-G4) with the same porosity and unit cell size, the pressure drop in Metal-G4 (SLM) is much larger than G4 (SLA). When $u_s \lessapprox 0.02$ m s⁻¹, the trend of pressure drop is similar. However, as velocity increases, the pressure drop in the SLM sample increases dramatically. It is known that metal 3-D printing with SLM cannot produce surfaces that are as smooth as in SLA using plastic, and this increased surface roughness is likely the cause of the elevated pressure drop even though the bulk characteristics of both Metal-G4 and G4 are the same. In figure 7(a), in a manner similar to figure 6(a), the linear trend at smaller $u_s$ gives way to a quadratic-like dependence at higher $u_s$ . This linear and quadratic variation in net pressure drop is quantified next following (1.1) in terms of $K_{\!Q1}$ and $K_{\!Q2}$ .

Figure 8.

Pressure drops of TPMS porous media with three different porosities. The black solid line represents the Ergun equation (1.4). For symbols, see table 1. (a) With $d \mapsto d_{\!p1}$ , (b) $d \mapsto d_{\!p2}$ in (1.4).

3.2. Permeability of TPMS porous media

Two permeability coefficients in (1.1) are usually obtained by fitting the experimental results with a second-order polynomial equation in $u_s$ . In the fitted equation, the coefficients of linear and quadratic terms correspond, respectively, to $\mu /K_{\!Q1}$ and $\rho /K_{\!Q2}$ . Since dynamic viscosity and density of the working liquid are known, $K_{\!Q1}$ and $K_{\!Q2}$ can be calculated. The least-square fitted equations to the pressure drop data of gyroid are shown in red dashed lines () in figure 6(b) for Set 1 and in figure 7(b) for Set 2, and all the data had an $R^2\gt 0.99$ . The two determined coefficients $K_{\!Q1}$ and $K_{\!Q2}$ for gyroid are tabulated in table 1, with the common trend that the values of $K_{\!Q2}$ decrease as porosity decreases. In other words, ‘form drag’ coefficient $C_Q\equiv 1/K_{\!Q2}$ in table 1 increases significantly as the internal structure of the gyroid acts as a blockage or resistance against the fluid flow. The results from Set 2, which consists only of gyroid with varying $l$ and similar $\phi$ , are presented in figure 7(b) and table 1. The table shows that with increasing $l$ (i.e. larger void regions) Darcy permeability $K_{\!Q1}$ increases, whereas $C_Q$ reduces. The case with Metal 3-D printing (Metal-G4) exhibits a higher $C_Q$ (i.e. ‘drag’) compared with the case G4 that has the same bulk properties but better surface finish. This again reinforces the comments made earlier about the increased drag owing to enhanced surface roughness in metal 3-D printed porous media.

Figure 9.

Friction factor $f$ (see (4.1)) distributions based on the hydraulic diameter $d_{\!{H}}$ in (a), and with $d_{\!{H}}$ replaced with $d_{{H\hbox{-}\textit{equ}}}$ in (b). The black solid line represents the Ergun equation (4.2) written in $d_{\!{H}}$ , whereas the blue dashed line is the linear part (i.e. only first term on the right-hand side) of (4.2).

Interestingly, and rather surprisingly, an attempt to fit (1.1) to the data of primitive and BCC in the least-square sense results in an erratic and sometimes a negative coefficient for  $u_s$ (i.e. a negative $K_{\!Q1}$ ) that is non-physical. An attempt at forcing a least-square fit with a non-negative coefficient for $u_s$ results in a value of zero for $K_{\!Q1}$ , again non-physical. This mathematical anomaly points to a physical difference between the form of pressure drop data for gyroid and the other two geometries (i.e. primitive and BCC). A careful observation of some of the primitive and BCC data in figure 6(a) or 6(c) and 6(d) shows that the pressure drop functions ‘bend outwards’ for higher velocities, making a linear plus a quadratic fit less viable. This becomes clear later (cf. figures 8 or 9) where the friction factor versus $Re/\phi$ trend for primitive and BCC differs considerably from gyroid. As will be discussed later, primitive and BCC have a behaviour that has similarities to a pipe flow (i.e. connected voids), which has a substantial ‘transitions’ region apart from the laminar to turbulent region. An additional issue encountered when fitting (1.1) to all cases is related to the number of data points used in fitting. As shown in Appendix A, for gyroid ( $\phi =0.85$ ), changing the number of data points does not change $K_{\!Q1}$ and $K_{\!Q2}$ beyond 10 %–20 %, whereas for primitive (and BCC), one could obtain a positive $K_{\!Q1}$ for a reduced and a negative $K_{\!Q1}$ for an increased number of data points. Appendix A further attempts to highlight these variations using one example from each gyroid and primitive by successively increasing the number of data points in fitting (1.1), starting from $u_s=0$ . In fact, the role of these large connected void regions presented in primitive and BBC remains to be fully investigated and it seems that this is a potential reason for the negative values of $K_{\!Q2}$ .

Although empirically determined $K_{\!Q1}$ and $K_{\!Q2}$ using (1.1) are widely used in the study of metal foams (and other porous media) to provide valuable information for modelling a specific geometry or configuration, it has several limitations. (i) For cases where a fitting (1.1), such as for gyroid in our experiments or multitude of other porous media like metal form or packed bed of spheres, dimensional $K_{\!Q1}$ and $K_{\!Q2}$ are less satisfactory, which also limits comparison between different geometries. (ii) There are other geometries, such as primitive and BCC in our case, where (1.1) cannot fit the data, owing to a different physical mechanism at play. (iii) Furthermore, a naive fitting of (1.1) produces different results for varying data count, especially for the ‘non-canonical’ porous media like primitive and BCC, and to a lesser extent in gyroid.

Some of these issues will be resolved in the next section, and in doing so, we will also demonstrate a procedure that will allow a unified comparison of all porous media. The non-dimensional pressure drop ( $f$ ) versus Reynolds number is the starting point. The difficulty, however, is the definition of a length scale required for non-dimensionalisation. We will begin with the ‘pore diameters’ defined in § 2.1.2 as the length scale, before moving to the hydraulic diameter ( $d_{\!{H}}$ ). As we shall see, $f$ versus $Re$ based on $d_{\!{H}}$ by itself is also not completely satisfactory, which will lead us to define an ‘equivalent hydraulic diameter’ ( $d_{{H\hbox{-}\textit{equ}}}$ ) to characterise and compare various porous media. We will commence with the data collected for Set 1 (for fixed $l=12$ mm, and varying $\phi$ and geometry), and towards the end show results from Set 2 (fixed $\phi \approx 0.85$ and varying $l$ ).

4.1. Hydraulic diameter and resulting $f$ versus $Re$ distribution

The concept of hydraulic diameter usually works well for similar geometries in diverse fluid flow situations. As discussed in § 1 hydraulic diameter $d_{\!{H}} \equiv 4V_\textit{fluid}/({\rm wetted\,area})$ , an estimation of $d_{\!{H}}$ requires void volume available for flow ( $V_\textit{fluid}$ ) that is related to porosity, the wetted surface area that includes the total surface area of TPMS porous media, and the surface area of the acrylic pipe that does not contact TPMS porous media.

Appendix C provides details of extracting $d_{\!{H}}$ . Briefly, Autodesk Netfabb computer software is used to estimate $V_\textit{fluid}$ and the wetted surface area. The values of $V_\textit{fluid}$ are experimentally estimated by immersing samples in water and calculating the increased volume. The difference between computer-estimated and experimental values is $\lessapprox 1.5\,\%$ . The final computed $d_{\!{H}}$ values for all samples are reported in table 2. Consistent with our expectation that smaller $\phi$ should result in a higher pressure drop and hence equivalent to a packed bed of smaller spheres, Set 1 data shows that reduction in $\phi$ results in a reduced  $d_{\!{H}}$ . Likewise, in Set 2, with a fixed $\phi \approx 0.85$ , an increase in unit cell size $l$ shows an increased $d_{\!{H}}$ .

Table 2.

Equivalent hydraulic diameter $d_{{H\hbox{-}\textit{equ}}}$ and equivalent sphere diameter $d_{{equ}}$ ; $d_{{equ}} = (3/2)d_{{H\hbox{-}\textit{equ}}}(1-\phi )/\phi$ ; the Darcy permeability $K_1$ from (4.5); and $K_2$ estimated in § 5.

To plot $f$ versus Reynolds number, the Ergun equation (1.4) is modified in-terms of the hydraulic diameter $d_{\!{H}}$ by using, $d=(3/2) d_{\!{H}}\,(1-\phi )/\phi$ , and with

(4.1) \begin{align} &f \equiv \frac {-\Delta\! P/\Delta x}{\rho u_s^2/d}\frac {\phi ^3}{(1-\phi )} = \frac {-\Delta\! P/\Delta x}{\rho u_s^2/d_{\textrm { H}}}\phi ^2 \left (\frac {3}{2} \right )\!, \end{align}

(4.2) \begin{align} & f = 150\left (\frac {2}{3}\right )\frac {\phi }{\textit{Re}_{d_{\!{H}}}} + 1.75, \end{align}

where, $\textit{Re}_{d_{\!{H}}} = u_s\,d_{\!{H}}/\nu$ . Note that the derivation of the Ergun equation begins with $d_{\!{H}}$ (rather than $d$ ); so, in that sense (4.2) is more broadly applicable than (1.4). Recall that the constants 150 and 1.75 in both (4.2) and (1.4) are for random packed spheres.

Using $d_{\!{H}}$ from table 2 and $f$ evaluated as in (4.1), figure 9(a) shows Set 1 data on a log–log $f$ versus $\textit{Re}_{d_{\!{H}}}$ plot using appropriate symbols for each case. The Ergun equation (4.2) is presented with a black solid line. The net pressure drop, composed of viscous and pressure drag, dominates, respectively, at lower and higher $\textit{Re}_{d_{\!{H}}}$ . Data from TPMS porous media start at different ordinate positions, but most porous media show a linear relationship that is suggested by the Ergun equation for low $\textit{Re}_{d_{\!{H}}}/\phi$ values; $\textit{Re}_{d_{\!{H}}}/\phi \lessapprox 10$ for gyroids; $\lessapprox 50$ for $\phi =0.55$ primitive and $\lessapprox 100$ for other primitives; $\lessapprox 200$ for BCC. This can be qualitatively explained by noticing that gyroid has the most ‘contorted’ path of all the three geometries (cf. figure 2), whereas primitive and BCC have an increasingly ‘easier’ flow path. An approximate analogy can be drawn between these cases and rough-walled pipe flows, where a reduction in the wall-roughness (i.e. a reduced flow contortion) increases the linear dependence region. Enhanced flow contortions result from flow separations around roughness or strongly curved pore regions (with a subsequent adverse pressure gradient), and hence, a dominance of the pressure drag over the viscous drag.

It is clear that the majority of the data in figure 9(a) do not follow the Ergun equation, and this is expected because the constants in the Ergun equation are tuned only for pressure drop through packed spheres. Nevertheless, $f$ versus $\textit{Re}_{d_{\!{H}}}$ is not uncommon to be plotted for a variety of porous media, such as metal forms (Richardson, Peng & Remue Reference Richardson, Peng and Remue2000; Dietrich et al. Reference Dietrich, Schabel, Kind and Martin2009), or packed rods, cylinder, plates, etc. (Richardson Reference Richardson2002), or more recently even for a TPMS ‘diamond’ geometry (Hawken et al. Reference Hawken, Reid, Clarke, Watson, Fee and Holland2023), and the aim then is to characterise the constants $K_{\!Q1}$ and $K_{\!Q2}$ , or their equivalent numerical values, such as 150 and 1.75 in (4.2) as a function of different parameters such $\phi$ , $l$ and type of geometry. As alluded in the Introduction § 1, there are two essential deficiencies: first, different functional fitted curves for each specific case make comparison of different cases difficult. Furthermore, in many cases, an estimation of $d_{\!{H}}$ is not even possible, and an arbitrary ‘pore diameter’ is chosen to plot $f$ versus $Re$ , making comparison of different cases even more difficult. Second, a closer examination of the data in figure 9(a) shows that at higher $\textit{Re}_{d_{\!{H}}}/\phi$ , especially for primitive and BCC, $f(d_{\!{H}})$ increases before becoming a constant. This ‘dip’ in the $f$ versus $\textit{Re}_{d_{\!{H}}}$ implies that the two-constant quadratic model (1.1) or the Ergun equation (4.2) cannot completely capture the flow physics associated with variation of $f$ as the flow ‘transitions’ from the laminar regime to a fully turbulent regime.

Towards a resolution to the first deficiency, where no single existing ‘diameter’ can be uniquely related to a case with specific $\phi$ , $l$ and geometry, we will define an ‘equivalent hydraulic diameter’ $d_{{H\hbox{-}\textit{equ}}}$ , that is possible for all porous media (and not just TPMS which is the focus here) and is unique to each porous medium.

4.2. Equivalent hydraulic diameter, and non-dimensional pressure drop

The essential idea is to find the equivalent hydraulic diameter that moves the laminar drag law region of each geometry to that given by the laminar region in the Ergun equation (for packed spheres). Thus, for each geometry, we are trying to find the equivalent hydraulic diameter that gives the same laminar viscous drag as packed spheres. Said differently, we plan to identify the ‘linear’ or laminar region for each geometry in figure 9(a) and move it towards the Ergun equation (black line in figure 9 a) by varying the hydraulic diameter as a free parameter. The diameter that best matches (4.2) in the linear regime is $d_{{H\hbox{-}\textit{equ}}}$ . As an aside, we note that this is similar to defining an ‘equivalent sand-grain roughness’ in rough wall-flows in the fully turbulent regime, whereas here we focus on the laminar regime. Below, we will present two methods to estimate $d_{{H\hbox{-}\textit{equ}}}$ that should provide the same $d_{{H\hbox{-}\textit{equ}}}$ within numerical and experiential data errors. With such a $d_{{H\hbox{-}\textit{equ}}}$ , that is listed in table 2 for different cases, $f$ versus $\textit{Re}_{d_{{H\hbox{-}\textit{equ}}}}$ can be obtained and is plotted in figure 9(b) for Set 1. As required, the ‘laminar’ region of the data follows the linear part of the Ergun equation: $f = 150 ({2}/{3} ){\phi }/{\textit{Re}_{d_{{H\hbox{-}\textit{equ}}}}}$ shown in the blue dashed line. Before discussing figure 9(b) or features of $d_{{H\hbox{-}\textit{equ}}}$ , we discuss the two methods for estimating $d_{{H\hbox{-}\textit{equ}}}$ .

Method 1 (this is a two-step process). In Step 1, the ‘laminar’ region or the range of $\textit{Re}_{d_{\!{H}}}/\phi$ where $f_H$ varies linearly with $\textit{Re}_{d_{\!{H}}}/\phi$ (in figure 9 a) is determined. Step 1 starts by fitting a linear line with slope –1 (on log–log axes) to the experimental data points in an increment of one data point at each time, and stops where the error is minimum. In Step 2, we vary the value of $d_{\!{H}}$ such that the linear region in the data matches the linear part of the Ergun equation ( $f = 150 ({2}/{3} ){\phi }/{\textit{Re}_{d_{{H\hbox{-}\textit{equ}}}}}$ ). The details of this two-step process is illustrated in Appendix D.1 using figure 15 by taking the gyroid at $\phi =0.85$ as an example. Note that although we start the process from a known $d_{\!{H}}$ , one could start from any reasonable but arbitrary diameter value, and the procedure will converge to the same $d_{{H\hbox{-}\textit{equ}}}$ . This is further evidenced by the second method for estimating $d_{{H\hbox{-}\textit{equ}}}$ described below, which does not rely on the knowledge of $d_{\!{H}}$ .

Method 2 (this is a one-step process). Similarly to Method 1, we again determine the linear drag region, but here in the $(-\Delta\! P/\Delta x)$ versus $u_s$ distribution, as those presented in figures 6 and 7. We start fitting a straight line passing through the origin, with slope ( $m$ ) as a free parameter, starting at $u_s=0$ , one data point at a time, until the error becomes a minimum. Using the dynamic viscosity $\mu$ , the fitted line becomes

(4.3) \begin{equation} \frac {-\Delta\! P}{\Delta x}=\frac {1}{K_1}\mu u_s, \end{equation}

where, the permeability strictly within the linear region: $K_1 \equiv \mu /m$ . Now, the Ergun equation (4.2) can be rearranged to retain only $(-\Delta\! P/\Delta x)$ on the left-hand side,

(4.4) \begin{equation} \frac {-\Delta\! P}{\Delta x}=\left [\left (\frac {2}{3}\right )^2 \frac {150}{d_{\!{H}}^2 \phi } \right ]\mu u_s + \left [ \frac {1.75}{d_{\!{H}} \phi ^2 (2/3)} \right ]\rho u_s^2. \end{equation}

Comparison of the linear part of (4.4) and (4.3) allows us to estimate $d_{{ H\hbox{-}\textit{equ}}}$ , such that

(4.5) \begin{equation} d_{{H\hbox{-}\textit{equ}}} = \left (\frac {2}{3}\right ) \left [ \frac {150\,K_1}{\phi } \right ]^{1/2} \!. \end{equation}

Further details of Method 2 by considering gyroid with $\phi =0.85$ data are presented in Appendix D.2 using figure 16 (where, a different symbol, $\tilde {d}_{{H\hbox{-}\textit{equ}}}$ is used to represent values obtained from (4.5) to distinguish it from those estimated by Method 1).

Although theoretically one would expect the same $d_{{H\hbox{-}\textit{equ}}}$ by both Methods 1 and 2, practically, we obtain the same values within $\pm 0.05$ mm. As such, we present only $d_{{H\hbox{-}\textit{equ}}}$ values obtained from Method 1 in table 2. The slight difference is owing to the different fitting procedures carried out on an experimental data set with inherent measurement errors. Once $d_{{H\hbox{-}\textit{equ}}}$ is found, we can also define an equivalent sphere diameter $d_{{equ}} \equiv (3/2)d_{{H\hbox{-}\textit{equ}}}(1-\phi )/\phi$ , that is also presented in table 2. So, physically an arbitrary porous media with $d_{{equ}}$ will produce the same drag as packed spheres of diameter $d_{{equ}}$ in the laminar region. We note that our direct estimation of $K_1$ (leading to $d_{{H\hbox{-}\textit{equ}}}$ by Method 2) or $d_{{H\hbox{-}\textit{equ}}}$ (leading to $K_1$ as in Method 1) is empirically estimated from the available data within linear region around at the smallest Reynolds number. In the limit of $Re=0$ (that can never be achieved in an experiment, unlike in a numerical simulation of Stokes flow), the $K_1$ will correspond to the true viscous laminar permeability; this draws attention to the distinction between the viscous–inertial laminar and purely viscous laminar regimes (e.g. Lasseux et al. Reference Lasseux, Abbasian, Ali and Ahmadi2011; Ahmed & Bottaro Reference Ahmed and Bottaro2023). At this point, it is also worth comparing the present experimentally obtained data with the numerical simulations of Ahmed & Bottaro (Reference Ahmed and Bottaro2023) for gyroid and primitive. Towards this, in table 2 we present $K_1$ normalised with the unit cell length, i.e. $K_1/l^2$ , which for gyroid and primitive range between $0.1\times 10^{-3}$ and $2\times 10^{-3}$ that is reasonably close to the range of values presented by Ahmed & Bottaro (Reference Ahmed and Bottaro2023) (in their figure 4). Furthermore, we note the differing trend of gyroid and primitive as $\phi$ changes; for $\phi =0.85$ , the $K_1/l^2$ value is higher for primitive compared with gyroid, whereas gyroid exhibits a higher value of $K_1/l^2$ compared with primitive for $\phi =0.55$ , and these trends are again consistent with the simulations of Ahmed & Bottaro (Reference Ahmed and Bottaro2023). We, however, note that a closer comparison shows that our values of $K_1$ are smaller than $Re=0$ , Stokes limit, presented by the numerical simulations of Ahmed & Bottaro (Reference Ahmed and Bottaro2023). Leaving aside the inherent measurement uncertainties and differences that can arise by a comparison between uncoordinated experiments and numerical simulations, this lower $K_1$ may point to the finite inertial effect and indicates that we are likely not at the Stokes limit, especially for the larger $l$ . The Reynolds number used by them is $ \rho M l^3/\mu ^2 \approx 1880 M$ for $l=12$ mm and $ 70 M$ for $l=4$ mm with our experimental $\rho$ and  $\mu$ values, where $M=-\Delta\! P/\Delta x$ in Pascals per metre. For this Reynolds number to be ‘small’ (say, less than $\approx$ 5000 from their figure 4), we should have $-\Delta\! P/\Delta x$ below a few Pascals per metre and 71 Pa m⁻¹ for $l=12$ and 4 mm, respectively. Within our linear fitting regime, the majority of cases lie above this limit, and we therefore anticipate a finite inertial contribution to the estimates of $K_1$ . Another observation, now from Set 2 (G1 to G5) where $\phi \approx 0.85$ and  $l$ is changing, is that the values of $K_1/l^2$ are relatively similar (at least lesser variations compared with Set 1). This shows that although the larger length scale of the flow set-up, i.e. the pipe diameter and the total length of the porous media or the number of unit cells affect the estimated $K_1$ , the effect is likely minor.

Once $d_{{H\hbox{-}\textit{equ}}}$ is estimated, it is now possible to compare different configurations as the Reynolds number increases. Figure 9(b) shows the pressure drop results in terms of $d_{{H\hbox{-}\textit{equ}}}$ . The functional forms of $f$ for gyroid at the three different porosities are reasonably close to the Ergun equation, unlike those of primitive and BCC. Figure 10(a) shows the same data as in figure 9(b) but only for gyroid removing the clutter. The closeness of gyroid data to the Ergun equation suggests the ‘equivalence’ of flow dynamics within the gyroid geometry and randomly packed spheres because of the intricate geometric structure of gyroid. It is evident that $\phi =0.85$ resembles packed spheres closest, and decreasing $\phi$ results in an increased drag compared with random spheres. For example, gyroid with $\phi =0.55$ will result in a higher turbulent drag than packed spheres of an equivalent diameter, whereas $\phi =0.85$ shows a slightly lower drag (compared with spheres) in the fully turbulent regime.

Figure 10.

Friction factors with the equivalent hydraulic diameters for (a) gyroid; (b) primitive; (c) BCC. Black solid and blue dashed lines are the same as in figure 9. The dashed lines are models based on $K_1$ and $K_2$ , whereas the dotted lines are based on $K_{\!Q1}$ and $K_{\!Q2}$ obtained by a least square fitting to (1.1), which are discussed in § 5.

Figure 10(b) and 10(c) show corresponding data for primitive and BCC. Beyond the linear region, primitive and BCC show different relationships compared with the Ergun equation, but have a similar pattern for the different porosities as the flow rates increase. The converging–diverging shape of the primitive geometry (see figures 1 b and 1 e and 2 b) and the open structure with the node and strut in the middle of the BCC cases (see figures 1 c, 1 f and 2 c) deviates from random spheres, and highlights the large-scale inhomogeneities in these porous media. In primitive (cf. figure 10 b), $\phi =0.85$ seems to show a slightly reduced drag in the fully turbulent region (corresponding to highest $\textit{Re}_{d_{{H\hbox{-}\textit{equ}}}}/\phi$ ), whereas smaller porosities appear to reach a slightly higher asymptotic drag compared with packed spheres. In BCC, on the other hand, it shows a reduced turbulent drag at all $\phi$ , and this low-drag behaviour is hardly affected by changing $\phi$ . It is important to highlight that the comparison of drag between different parameters of the porous media in fully turbulent regime is only possible after the definition of $d_{{H\hbox{-}\textit{equ}}}$ . Specific to the TPMS geometries considered here, we observe a lower turbulent drag for $\phi =0.85$ in all three geometries, and for all $\phi$ values for BCC. It is quite likely that TPMS with almost zero mean curvature seems favourable to drag reduction when compared with packed spheres that have a fixed and finite curvature. One can speculate that it is likely that lower pressure drag, which is usually related to flow separation, is somewhat curtailed in these lower curvature geometries.

Also, note that the insets in figure 10(b) and 10(c) highlight the start of the transitional region where inertial effects begin to manifest on a larger scale in the primitive and BCC cases, and there is a ‘dip’ behaviour in $f$ distribution as $Re$ increases. That is, after the laminar region (where $f \sim (\textit{Re}_{d_{{H\hbox{-}\textit{equ}}}}/\phi )^{-1}$ ), $f$ increases before levelling towards an approximately constant fully turbulent region. This increase in $f$ with $Re$ cannot be expressed as the usual linear–quadratic drag law (cf. (1.1) or the Ergun equation (1.4)) with a sum of laminar and fully turbulent flow drag contained in the existing formulations of porous media.

This complex behaviour of $f$ is, however, reminiscent of flow through rough pipes presented on a ‘Moody’ diagram whereas the flow transitions from laminar to turbulent, an increase in $f$ with $Re$ is quite common. In pipes, as the flow transitions from a laminar regime caused by the roughness-induced instabilities, the increased momentum transport owing to almost exponentially growing instabilities (e.g. Philip & Cohen (Reference Philip and Cohen2010), with sufficient space to grow in an empty pipe, i.e. possibility of long-range correlations) results in a rapid increase in $f$ . As $Re$ is further increased, the nonlinear transfer of energy (e.g. Cohen, Philip & Ben-Dov Reference Cohen, Philip and Ben-Dov2009) arrests this growth of $f$ and brings to a constant value that is a fixed fraction of the total dynamic pressure in the pipe. We suspect a similar phenomenon happening in primitive and BCC geometries with their straight pipe-like structure (see figure 2 bii and 2 cii), and roughness that is inherent at the edges. The period structure with large empty streamwise expansive regions in primitive and BCC allows the possibility of larger scale ‘transitional’ structures (with possible long-range correlations) to exist within these geometries. This is in contrast to a ‘random porous media’ (such as a packed bed of spheres) where such structures are not allowed owing to the non-periodic nature of the geometry. Note that although gyroid has a periodic structure, the periodicity is hampered in the streamwise direction of the flow (see figure 2 a ii), resulting in gyroid to behave more closely to a random medium. In fact, the effect of large connected void regions forming direct pathways along the streamwise direction in primitive and BCC is reflected in the much larger $K_1$ in table 2, unlike gyroid which has a smaller $K_1$ . Overall, this suggests that a more complex model, unlike the quadratic Ergun equation, is required to model flows in the primitive and BCC geometries that possess period geometries with extended streamwise regions.

The geometric similarity of the gyroid (compared with other geometries considered here) and randomly packed beds is clearly observed in figure 10(a). We now focus on our attention to Set 2 (cf. table 1) for gyroid where $\phi \approx 0.85$ and the unit cell size $l$ is varied.

4.3. Effect of the number of unit cells

Figure 11.

Here $f$ versus Reynolds number based on (a) $d_{\!{H}}$ , and (b) $d_{{H\hbox{-}\textit{equ}}}$ for data for Set 2 (cf. table 1) where $\phi \approx 0.85$ for $l=$ 4 mm G1 (), 6 mm G2 (), 9 mm G3 and G4 (and ), 12 mm G5 () and metal sample $l=$ 9 mm Metal-G4 ().The black solid line represents the Ergun equation (4.2) written in $d_{\!{H}}$ , whereas the blue dashed line is the linear part of (4.2). For comparison purposes, the pressure drop measurement with the gyroid at $l=$ 12 mm and $\phi =0.85$ (case G5) in table 1 is conducted in the Set 2 measurement section (36 cm), and several points are included in (a) with yellow symbols. Insets show zoomed-in view. The dashed green lines in (b) are models based on $K_1$ and $K_2$ as discussed in § 5.

Figure 11(a) and its inset show non-dimensionalised pressure drop in terms of the hydraulic diameter $d_{\!{H}}$ for Set 2, where $l=$ 4, 6, 9 and 12 mm. For the same $\textit{Re}_{{d_{\!H}}}/\phi$ , a reduction in cell-size results in a corresponding increased dimensionless drag $f$ (see (4.1) where $d_{\!{H}}$ is replaced with $d_{{H\hbox{-}\textit{equ}}}$ ) and not just the dimensional drag. The variation in dimensionless drag (similar to figure 9 a) suggests the unsuitability of $d_{\!{H}}$ in furnishing a self-similar profile. Like most practical applications, we recall that our porous media is contained within a casing, i.e. a cylindrical pipe. Therefore, the pipe diameter will also feature as another parameter that can affect the results. Although $d_{\!{H}}$ accounts for the fluid coming in contact with the outer pipe, the collapse of the data for different  $l$ is not satisfactory. We suspect that as the pipe diameter becomes increasingly large compared with $l$ , the effect of outer geometry will become less significant. Nevertheless, as presented in figure 11(b), the use of $d_{{H\hbox{-}\textit{equ}}}$ collapses all the different cases within the laminar region, and shows the expected variations as the flow becomes turbulent. Within the turbulent region (highlighted in the inset of figure 11 b) the data separate into two groups: $l=4$ and 6 mm cases group together, whereas the $l=9$ and 12 mm cases fall on top of each other. Interestingly, the metal $l=9$ mm cases, which we know has a rougher surface finish than the other plastic surfaces, merges with the smaller $l$ value cases. This indicates that surface roughness within the turbulent region acts dynamically as smaller cell sizes. Consistently, $d_{{H\hbox{-}\textit{equ}}}$ for Metal-G4 and G4 are 2.9 and 3.0 mm, respectively (cf. table 2), although both have the same bulk geometric characteristics. Note that the different  $\phi$ data presented in figure 10 separate in the turbulent region, because as $\phi$ decreases $d_{{H\hbox{-}\textit{equ}}}$ decreases (see table 2), indicating larger ‘roughness’ effect with decreasing $\phi$ . Similarly, decreasing $l$ in figure 11(b) evidently results in a reduced value of $d_{{H\hbox{-}\textit{equ}}}$ leading to a ‘rougher’ turbulent flow. Note that there is no expectation that the data would collapse within the turbulent region, as the interaction of multiple turbulent length scales with the porous media is less likely to be represented by $d_{{H\hbox{-}\textit{equ}}}$ defined for the laminar region. However, for approximate engineering calculations, where if we ignore the variation of data in figure 11(b) within the turbulent region, the Ergun equation (1.4) with $d_{{H\hbox{-}\textit{equ}}}$ seems sufficient to estimate the complete drag distribution for the gyroid geometry.