Permeability prediction of fibrous porous media by the lattice Boltzmann method with a fluid-structure boundary reconstruction scheme

Abstract

The D3Q27 lattice Boltzmann method (LBM) combined with the fluid-structure boundary reconstruction (fsBR) scheme and the interpolated bounce back (IPBB) method is extensively evaluated to predict the permeability of nonwoven fibrous porous media. The fsBR-IPBB method transfers digitally defined step-like boundary data, e.g. the three-dimensional structure data obtained by the X-ray computed tomography, to continuous smooth boundary data via level-set functions. It leads to highly accurate calculations despite low lattice resolutions of thin fibers with circular cross-sectional shapes, compared to the conventional half-way bounce back (HWBB) method. The fsBR-IPBB method is first applied to predict the permeability of two different arrays of impermeable circular cylinders and verified by comparing the results with the data in the literature. We then validate the method referring to the numerically and experimentally obtained permeability of six types of nonwoven fabrics prepared by the industrial hydroentanglement process. Finally, the discussion on the applicability and the limitation of the macroscopic correlation models to estimate permeability of porous media is carried out. The results show that although the calculated permeability is in reasonable agreement with the measured one with an error of 8.1–16.3%, analytical or empirical correlation models fail to give the correct trend due to the highly inhomogeneous and anisotropic properties of hydroentangled nonwoven fabrics.

Introduction

Nonwoven fabrics are fibrous porous materials, which are assemblies of fibers bonded together in the form of sheets, webs, or bats. They have been utilized in numerous industrial applications such as materials for absorbency, filtration, car interior, fabric reinforcement, etc. Recently, high-performance nonwoven fabrics have replaced the traditional fabric materials owning to the development of the production technology [1]. One of the remarkable properties of the nonwoven fabrics in such applications is the permeability, which is significantly important to the product performance [2–4]. Accordingly, a numerical tool for accurately analyzing fluid flows through nonwoven fabrics without conducting expensive and time-consuming experiments has been long investigated for evaluating the performance and designing products.

The earlier analytical models were based on calculating the flows through ordered arrays of unidirectional cylinders in the cylindrical unit cell by the Stokes equations. Happel extended the free-surface model developed to predict the flow through the assembly of spherical particles to the flows through ordered arrays of cylinders [5]. Sangani and Acrivos determined the solutions for the flows through the regular and staggered arrays of circular cylinders as the function of the volume fraction of the cylinders [6]. Drummond and Tahir extended the permeability prediction of Sangani and Acrivos by utilizing a matching technique in the cylindrical unit cell approach [7]. More recently, Tomadakis and Robertson investigated the permeability of the flows through random fiber structures as the function of their formation factor, specific surface area and porosity [8]. Most of analytical models above were based on the assumption of homogeneous and isotropic geometry. The inhomogeneous and anisotropic nature of nonwoven fabrics though is a key factor for designing high-performance products, the relationship between the permeability and complex geometric properties is difficult to represent by these models.

Accurate prediction of the permeability can be obtained from the direct numerical simulation of the fluid flow, which generally includes four key processes: obtaining 2D image stacks by the X-ray micro-computed tomography (CT), image processing, digital 3D reconstruction of the structures and computation of the flows [2,3,9]. The lattice Boltzmann method (LBM) has been successfully applied to investigate the flows through highly complex and anisotropic porous media. Compared to continuum-based approaches, the LBM has following intrinsic properties: simple fluid/structure representation of lattice nodes by binarized on/off treatment, easy implementation of no-slip boundary condition and easy parallelization using graphics processing units (GPU) due to its discrete nature and locality of calculation [10,11]. Gao et al. calculated the permeability of the carbon paper gas diffusion layers and glass beads packing by the LBM combined with the X-ray CT [12]. They indicated that the LBM was efficient to calculate the pore-scale flow in porous media and capable of estimating the anisotropic permeability. Eshghinejadfard et al. calculated the permeability of the real porous media probe reconstructed from a set of CT images by the LBM [13]. They reported that the accuracy of the predicted permeability depended strongly on the domain resolution. Referring to the literature, the following two requirements are to be satisfied in order to accurately calculate the flows through nonwoven fabrics: (i) the fiber shape with a circular cross section of the order of several micrometers is sufficiently resolved, and (ii) computational domains are sufficiently representative to smoothen the effect of local inhomogeneities of nonwoven fabrics of the order of several millimeters.

The conventional non-slip wall boundary condition of the LBM, the half-way bounce-back (HWBB) method, is very easy to implement and its numerical results are second-order accuracy [14]. With the HWBB method, though, curved surfaces are represented by regular cubic grids, resulting in high computational costs to satisfy the above two requirements due to the requirement of fine lattice resolutions to represent fibers and to calculate in high accuracy [15]. On the other hand, in the case with mathematically described fluid-structure boundary data, the improved boundary condition, the interpolated bounce-back (IPBB) method, can be applicable for more accurate calculations with relatively low lattice resolutions to present curved boundaries [16,17]. Pan et al. investigated the several fluid-structure boundary conditions by calculating the flows through a body-centered-cubic array of spheres and a random-sized sphere pack media using the LBM [16]. They concluded that the IPBB method yielded more accurate calculation at coarser resolution and with shorter iteration times, compared to the standard bounce-back (SBB) method. Ahrenholz et al. calculated the flows through porous media reconstructed from highly resolved tomographic binarized data by the LBM with the IPBB method [18]. The Marching Cube algorithm was utilized to reconstruct the surface of the porous media as a set of planar triangles, and the distances between the lattice nodes and the planar triangle surfaces were obtained to apply the IPBB method. They concluded that the IPBB method with fine resolution X-ray CT data was more accurate and led to faster convergence, compared to the SBB method. However, it is not possible to obtain mathematically described fluid-structure surface data of actual fibrous porous media, while the digitized 3D structure dataset acquired by the X-ray CT is the most probable one in the modern engineering.

Therefore, the objective of this study is to provide an improved computational method to accurately predict the permeability of nonwoven fabrics with relatively coarse resolution X-ray CT data, and to establish the relationship between the resolution and the prediction error for the practical engineering usage. We apply the LBM with the combination of the IPBB method and the fluid-structure boundary reconstruction (fsBR) scheme proposed by Suga et al. [19]. The fsBR scheme translates digitized step-like boundary shapes to continuous level-set functions, which can be utilized instead of mathematically described fluid-structure boundary data for the IPBB method. In this study, we first verify the fsBR-IPBB method by calculating the permeability of two different arrays of impermeable circular cylinders. The results are compared with the data in the literature. The fsBR-IPBB method is then applied to the fluid flow computation to obtain the permeability of six types of nonwoven fabrics produced by the industrial hydroentanglement process. The numerical permeability is compared with the experimental permeability to validate the method. Finally, we discuss the applicability and limitation of the macroscopic correlation models to estimate the permeability of porous media.

Methodologies

The Lattice Boltzmann method

The LBM is a computational fluid dynamics technology based on the discretized Boltzmann equation, which treats a fluid flow by calculating the time-evolution of the distribution function for the groups of particles that stream and collide on lattice nodes. See [20–22] for detailed instructions of the LBM. In the streaming process, the particles are propagated to adjacent lattice nodes in certain directions. In the collision process, the particles at each lattice node are redistributed toward a local equilibrium state [10]. This evolution of the distribution function is described by the lattice Boltzmann equation discretized in time and space, which is expressed as

$$f_{\alpha }\left(\mathit{x}+{\mathit{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)=f_{\alpha }\left(\mathit{x},t\right)+{\mathrm{\Omega }}_{\alpha }\left(\mathit{x},t\right)$$

(1)

where $f_{\alpha }\left(\mathit{x},t\right)$ is the density distribution function at time $t$ and location $\mathit{x}$ , subscript $\alpha$ denotes the finite velocity component of the distribution function, ${\mathit{e}}_{\alpha }$ is the discrete velocity which is shown in Table 1 with the particle velocity $c=\Delta /{\delta }_t$ , $\Delta$ is the lattice spacing, ${\delta }_t$ is the time step and ${\mathrm{\Omega }}_{\alpha }\left(\mathit{x},t\right)$ is the collision operator. This represents that particles move with their discrete velocity ${\mathit{e}}_{\alpha }$ to the neighboring lattice node $\mathit{x}+{\mathit{e}}_{\alpha }{\delta }_t$ at the next time step $t+{\delta }_t$ , followed by the collision process under the control of the collision operator ${\mathrm{\Omega }}_{\alpha }$ . While variety of collision operator models have been proposed such as the single-relaxation-time (SRT) model and the multiple-relaxation-time (MRT) model [23,24], the simple SRT Bhatnagar-Gross-Krook (BGK) model [25,26] is employed in this study, which is expressed as

$${\mathrm{\Omega }}_{\alpha }\left(\mathit{x},t\right)=-\frac{1}{\tau }\left[f_{\alpha }\left(\mathit{x},t\right)-{f\hspace{0.17em}}_{\alpha }^{eq}\left(\mathit{x},t\right)\right]$$

(2)

where $\tau$ is the relaxation time and ${f\hspace{0.17em}}_{\alpha }^{eq}\left(\mathit{x},t\right)$ is the equilibrium distribution function. The BGK model describes that the distribution function approaches toward an equilibrium condition at a time scale represented by the relaxation time $\tau$ . The equilibrium distribution function determined by the fluid density and momentum is given as

$${f\hspace{0.17em}}_{\alpha }^{eq}\left(\mathit{x},t\right)=w_{\alpha }\rho \left[1+\frac{{\mathit{e}}_{\alpha }·\mathit{u}}{c_s^2}+\frac{{\left({\mathit{e}}_{\alpha }·\mathit{u}\right)}^2}{2c_s^4}-\frac{\mathit{u}·\mathit{u}}{2c_s^2}\right]$$

(3)

where $w_{\alpha }$ is the weight coefficient, $\rho$ is the density and $c_s$ is the sound speed. The weight coefficient $w_{\alpha }$ and the sound speed $c_s$ are associated to the lattice discretizing models [27]. In this study, the three-dimensional standard square lattice model of D3Q27 is employed. Figure 1 shows the schematic illustration of the D3Q27 lattice model. Table 1 lists $c_s$ , ${\mathit{e}}_{\alpha }$ and $w_{\alpha }$ of the D3Q27 discrete velocity model.

Table 1.

Parameters of the D3Q27 discrete velocity model.

$c_s/c$	${\mathit{e}}_{\alpha }/c$	$w_{\alpha }$
$1/\sqrt{3}$	$\left(0\text{,}0\text{,}0\right)$ $\left(\pm 1\text{,}0\text{,}0\right)$ , $\left(0,\pm 1\text{,}0\right),\left(0\text{,}0,\pm 1\right)$ $\hspace{0.17em}\left(\pm 1,\pm 1\text{,}0\right),\left(\pm 1\text{,}0,\pm 1\right),\left(0,\pm 1,\pm 1\right)$ , $$ $\hspace{0.17em}\left(\pm 1,\pm 1,\pm 1\right)$	$8/27$ $\left(\alpha =0\right)$ $2/\mathrm{27\hspace{0.17em}}\left(\alpha =1, \dots , 6\right)$ $1/\mathrm{54\hspace{0.17em}}\left(\alpha =7, \dots , 18\right)$ $1/\mathrm{216\hspace{0.17em}}\left(\alpha =19, \dots , 26\right)$

Figure 1.

The schematic representation of the D3Q27 discrete velocity model.

The fluid density $\rho$ , velocity $\mathit{u}$ and pressure $p$ are obtained from the conservation law and the equation of state as follows

$$\rho =\sum _{\alpha }{f}_a$$

(4)

$$\rho \mathit{u}=\sum _{\alpha }{f}_a{\mathit{e}}_{\alpha }$$

(5)

$$p=\rho c_s^2$$

(6)

The kinematic viscosity is given as

$$\nu =c_s^2\left(\tau -\frac{1}{2}\right){\delta }_t$$

(7)

Bounce-back boundary conditions

The conventional HWBB method treats curved surfaces as stepwise shaped surfaces on regular cubic grids and assumes that any particle that hits a solid wall during the streaming step is reflected back to the initial position of streaming. Thus, it has fundamental difficulty to calculate the fluid flow in fibrous porous media correctly, typically in the case of calculating large domain, due to the requirement of fine lattice resolution for representing fibers with a smooth circular cross-sectional shape [28]. The improved boundary treatment: the IPBB method with mathematically described fluid-structure interfaces can be applicable for more accurate calculations in low lattice resolutions to present curved boundaries. Hence, this study applies the linear IPBB method [16]. The linear IPBB formulae for $f_L\left({\mathit{x}}_A,t_{n+1}\right)={\stackrel{\sim }{f}}_R\left({\mathit{x}}_C,t_n\right)$ is given as

$$f_L\left({\mathit{x}}_A,t_{n+1}\right)=\left\{\begin{array}{c}\left(1 - 2q\right)f_R\left({\mathit{x}}_A,t_{n+1}\right)+2q{\stackrel{\sim }{f}}_R\left({\mathit{x}}_A,t_n\right), q<1/2\\ \left(1-\frac{1}{2q}\right){\stackrel{\sim }{f}}_L\left({\mathit{x}}_A,t_n\right)+\frac{1}{2q}{\stackrel{\sim }{f}}_R\left({\mathit{x}}_A,t_n\right), q\ge 1/2\end{array}\right.$$

(8)

$$q=\frac{\left||{\mathit{x}}_A-{\mathit{x}}_W|\right|}{\left||{\mathit{x}}_A-{\mathit{x}}_B|\right|}$$

(9)

where $f$ and $\stackrel{\sim }{f}$ respectively denote the post- and pre-collision states of the distribution function. Position ${\mathit{x}}_A$ , ${\mathit{x}}_B$ , ${\mathit{x}}_C$ and ${\mathit{x}}_W$ represent the boundary fluid node, the structure node, the off-lattice location involved in interpolations and fluid-structure interface, respectively. The subscripts L and R indicate left- and right-bound directions, respectively. Figure 2(a) and (b) illustrates the IPBB procedures of $q<1/2$ and $q>1/2$ in a one-dimensional setting, respectively. In the case that ${\mathit{x}}_W$ is located at the half-way between ${\mathit{x}}_A$ and ${\mathit{x}}_B$ , i.e., $q=1/2$ , the incoming distribution function is simply equal to the corresponding outgoing distribution with the opposite momentum, hence it becomes identical to the HWBB method. The IPBB method, thus, is a generalized form of the HWBB method. The linear IPBB formulae requires only local information of the distribution functions and mathematically described fluid-structure interface data.

Figure 2.

The schematic representation of the IPBB procedures of (a) $q<1/2$ and (b) $q>1/2$ .

Fluid-structure boundary reconstruction scheme

The fsBR scheme developed by Suga et al. utilizes the level set method to describe surface boundaries [19]. The level set method is a well-known method for surface capturing and has been applied to a variety of evolving interface problems such as computational geometry, fluid mechanics, computer vision, material science, etc. [29]. The level set function $\varphi$ indicates normal distance from an interface between two phases, whose positive and negative values respectively correspond to one phase and the other [30]. For the following, we define the interface as the positions of $\varphi =0$ , the solid region as the region of $\varphi >0$ and the liquid region as the region of $\varphi <0$ . The magnitude of the level set function $\varphi$ at any location represents the distance from this location to the interface, thus $\left|\nabla \varphi \right|=1$ . In order to maintain its nature as a distance function without changing the zero level set, the reinitialization process is employed by solving the following partial differential equation [31]

$$\frac{\partial \varphi }{\partial t_{\mathit{init}}}=\mathrm{sgn}\left({\varphi }_0\right)\left(1-\left|\nabla \varphi \right|\right)$$

(10)

which can be rewritten as

$$\frac{\partial \varphi }{\partial t_{\mathit{init}}}+\mathit{\omega }·\nabla \varphi =\mathrm{sgn}\left({\varphi }_0\right)$$

(11)

where $t_{\mathit{init}}$ is the time step of the reinitialization process and $\mathit{\omega }$ is the unit normal to the interface and given as

$$\mathit{\omega }=\mathrm{sgn}\left({\varphi }_0\right)\frac{\nabla \varphi }{\left|\nabla \varphi \right|}$$

(12)

with initial conditions $\varphi \left(\mathit{x},0\right)={\varphi }_0\left(\mathit{x}\right)$ . The reinitialization process conserves the shape of the interface defined implicitly by the zero level set of ${\varphi }_0\left(\mathit{x}\right)$ . In the numerical procedure, the above equation is treated with a modified sign function

$$\mathrm{sgn}\left({\varphi }_0\right)=\frac{{\varphi }_0}{\sqrt{{\varphi }^2+{ϵ}^2}}$$

(13)

where $ϵ$ is a small value, which represents a thickness of the interface. This modification smoothens the exact sign function. Thus, the exact location of the interface is located between the range of $-\mathrm{ϵ}\le \varphi \le ϵ$ .

With the digitized boundary surface data, the fsBR scheme utilizes the reinitialization process with the modified sign function to reconstruct the original smooth interface. In this study, initial values of the level set functions are determined as arbitrary negative and positive values at the nodes in the fluid and solid regions, respectively. Since the actual surface is located within the distance of $\Delta$ which is the lattice spacing, $ϵ=\Delta$ is applied. In the case that the converged solution of the reinitialization process fails to reconstruct sufficiently recovered boundary surfaces, the additional relaxation step is applied. Due to the nature of the distance function, the distribution of the level set function away from an interface tends to be relaxed. Thus, the fsBR scheme performs the second reinitialization process after shifting the zero level set position to a certain level. Figure 3 illustrates the application of the fsBR scheme to the two-dimensional circular cylinder with the resolution of 10 pixels in diameter. Figure 3(a) shows the original circular cylinder described by mathematic expressions (black line) and digitized lattice nodes (voxels), respectively. Figure 3(b) shows the contour map of the level set function obtained by the fsBR scheme. The white line represents the reconstructed cylinder surfaces. The fsBR scheme well reconstructs the circular cylinder shape.

Figure 3.

The schematic representation of the fsBR scheme applied to the two-dimensional circular cylinder with the resolution of 10 pixels in diameter: (a) the original circular cylinder described by mathematic expressions (black line) and the digitized lattice nodes (voxels) and (b) the contour map of the level set function obtained by the fsBR scheme. The white line represents the reconstructed cylinder surfaces.

Results and discussions

Verification of the simulation method

In this section, the fsBR-IPBB method is first verified in terms of the lattice resolution by calculating the permeability of the regular arrays of circular cylinders using three types of bounce-back methods: the conventional HWBB method with the digitized boundary surfaces, the fsBR-IPBB method with the reconstructed boundary surfaces described by the level set function and the IPBB method with mathematically defined ideal boundary surfaces. A cubic domain consisting of a cylinder in the middle whose side is twice the cylinder diameter is considered. The lattice resolutions are set to be $D/\Delta =6,\mathrm{ }9,\mathrm{ }12,\mathrm{ }30,\mathrm{ }50$ and $100$ with the cylinder diameter $D$ . The flow perpendicular to the cylinder axis is imposed by applying the pressure gradient in the flow direction with the periodic boundary conditions to calculate the permeability [32]. The simulated Reynolds numbers based on the bulk velocity and the cylinder diameter are in the range of 1.17 × 10⁻³–3.54 × 10⁻¹. (The Mach numbers are thus lower than 10⁻⁴ and they are low enough to regard the flow to be incompressible. This is true for all simulations in this study.) Based on Darcy’s law, the permeability is obtained as

$$K=-\frac{\rho \nu u_D}{\nabla P}$$

(14)

where $\rho$ is the density, $\nu$ is the kinematic viscosity, $u_D$ is the Darcy velocity and $\nabla P$ is the pressure gradient in the direction of $u_D$ [33,34]. Figure 4 shows the variation of calculated permeability of the regular arrays of circular cylinders against the lattice resolution. The permeability is normalized by the cylinder radius $R\left(=D/2\right)$ . It is reasonably seen that the permeability obtained by the HWBB method is asymptotic to the ideal one obtained by the IPBB method with mathematically (ideally) given surface data, corresponding to the increase of the lattice resolution. Although the permeability obtained by the fsBR-IPBB method and the IPBB method with mathematically given surface shape data decreases for the coarse lattice resolutions, the agreement between both methods is very good even with the coarse lattice resolutions of $D/\Delta =6–12$ . It is thus confirmed that we can calculate the flow around the cylinder array with the fsBR scheme as accurate as with the mathematically given ideal surface shapes.

Figure 4.

The variation of the calculated permeability of the regular arrays of circular cylinders against the lattice resolution. The permeability is normalized by the cylinder radius $R$ .

Then, the permeability of two test cases, which are the regular and staggered arrays of circular cylinders represented in Figure 5, is investigated to verify the fsBR-IPBB method. A cubic domain consisting of a cylinder in the middle is considered for the regular cylinder arrays. The lattice resolution is chosen to be six in the cylinder diameter $\left(D/\Delta =6\right)$ . Dimensions of the computational domain vary from $Lx/\Delta \times Ly/\Delta \times Lz/\Delta =9\times 9\times 9$ to $24\times 24\times 24$ corresponding to the porosity from 0.65 to 0.95. A rectangular domain consisting of a cylinder in the middle and four quarter cylinders located at corners is considered for the staggered cylinder arrays. The lattice resolution is the same as that of the regular cylinder arrays. Dimensions of the computational domain vary from $Lx/\Delta \times Ly/\Delta \times Lz/\Delta =10\times 18\times 10$ to $40\times 72\times 40$ corresponding to the porosity from 0.69 $\mathrm{ }$ to 0.98.

Figure 5.

The schematic representations: The regular cylinder arrays of (a) the top view and (b) the 3D view. The staggered cylinder arrays of (c) the top view and (d) the 3D view.

The flows perpendicular and parallel to the cylinder axis are imposed by applying the pressure gradient in the flow directions with the periodic boundary conditions to calculate the transverse and longitudinal permeability, $K_{\perp }$ and $K_{\parallel }$ [32]. The Reynolds numbers based on the bulk velocity and the cylinder diameter are in the range of 2.95 × 10⁻⁴–1.76 × 10⁻² for the regular arrays of circular cylinders and 2.18 × 10⁻³–3.93 × 10⁻² for the staggered arrays of circular cylinders. Figure 6 shows the variation of the calculated permeability of the regular and staggered arrays of circular cylinders against the porosity. The permeability indicated is normalized by the cylinder radius. The present results are compared with earlier simulated permeability and the theoretical permeability in the literature [5,7,35–37]. The calculated permeability in two test cases shows good agreement with the theoretical data and the other simulated data in the literature, and therefore we confirm that the accuracy of the fsBR-IPBB method is good enough in the lattice resolution of $D/\Delta \ge 6$ .

Figure 6.

Variation of the transverse and longitudinal permeability of (a) the regular and (b) the staggered arrays of circular cylinder against the porosity. The permeability indicated is normalized by the cylinder radius.

Validation of permeability of nonwoven fabrics

Materials and properties

Nonwoven fabrics are prepared by the hydroentanglement process, while it is an important industrial polymer process to produce highly entangled nonwoven fabrics. In this process, first a carded web of polymer fibers is transported to the hydroentangle unit by a water-permeable mesh conveyor. Then, the fiber web is exposed to fine and high-pressure water jets, resulting in the formation of a mechanically entangled and bonded fiber structure. Therefore, the mechanical properties of fibers and the water jet pressure are key factors to determine the characteristics of nonwoven fabrics [38–41]. In this study, three types of fibers, rayon with the linear density of 2.2 and 3.3 dtex and polyethylene terephthalate (PET) of 2.2 dtex, are hydroentangled using two different water jet pressures (5 and 8 MPa), resulting in six types of nonwoven fabric samples as shown in Table 2. The porosities of nonwoven fabric samples are calculated as the ratio of the volume in the fiber and total region from the reconstructed smooth surface data. The detailed process is described in the following subsection. Suga et al. [19] reported that the error of this porosity calculation was 0.2% from the ideal value in the case of the body-centered-cubic structure. The pore scales are measured by a porometer (Perm Porometer, PMI). The average fiber diameters are measured from scanning electron microscope (SEM) images of the nonwoven fabrics using voltage of 10 kV at magnifications of 200×. Carded webs of 60 g∕m² basis weight are prepared and oriented in the cross-machine direction. These webs are then transported to the hydroentangle unit containing three manifolds of water jets by a conveyor at the speed of 10 m∕min. The water jet pressure of the first manifold is kept constant at 2 MPa to prewet the webs. The water jets whose diameters are 0.1 mm are arranged to form a single row with the spacing of 1 mm. The webs pass through the hydroentangle unit twice in order to be subjected to the water jets on the front and back sides. The webs are dried and cured at a temperature of 120 °C in an oven.

Table 2.

Design factors of nonwoven fabrics.

Sample ID	Polymer	Linear density[dtex]	Porosity, $\mathrm{\phi }$	Pore scale[µm]	Average fiber diameter[µm]	Water jet pressure [MPa]
						Manifold
						#1	#2	#3
1	Rayon	2.2	0.954	105.3	13.7	2	5	5
2	Rayon	2.2	0.966	98.5	13.7	2	8	8
3	Rayon	3.3	0.926	66.0	17.9	2	5	5
4	Rayon	3.3	0.918	77.6	17.9	2	8	8
5	PET	2.2	0.955	91.4	14.6	2	5	5
6	PET	2.2	0.949	92.2	14.6	2	8	8

The permeability measurement experiments are conducted following the standard method for the measurement of hydraulic conductivity of saturated porous materials (ASTM D 5084-03 method A or JIS A 1218). In this method, a flow in the through-plane of the pre-wetted nonwoven fabric is generated by the constant head difference. Then, the permeated water volume within 30 seconds is measured to obtain the flow rate to calculate the permeability by Darcy’s law. The experimental permeability is obtained as the averaged value with the standard deviation from 5–11 measurements for each nonwoven fabric sample. The standard deviations are within ±15% of the average permeability. All experiments are conducted at a room temperature of 24 °C at the standard atmospheric pressure.

Image Processing of 3D nonwoven fabric structure data

Three-dimensional structure data of nonwoven fabrics are obtained by the X-ray micro-CT (nano3DX, Rigaku). The scanning is performed at a resolution of 2.16 µm∕pixel, resulting in the scanning domain of 3.18 mm × 1.59 mm × 2.70 mm, which includes three periodic regions in the cross-machine direction according to the spacing of the water jet alignment. The lattice resolutions are $D/\Delta =6–8$ with the fiber diameter $D$ . The gray-scale images of nonwoven fabrics are exported as 2D image stacks in the TIFF file format for a further analysis. The image post-processing is performed using the open-source software ImageJ and the bundle of plugins distributed in the FIJI [42]. A 3D median filter with a rolling ball radius of 10 pixels is applied to subtract the background noise in the images. The segmentation process is then performed to binarize the images by means of the threshold algorithm of Otsu’s method with the optimal threshold value determined by Otsu’s algorithm [43]. The reconstructed smooth surface data are then obtained by supplying the binarized image data to the fsBR scheme. Figure 7 shows the schematic illustrations of the image processing procedures and the reconstructed 3D nonwoven fabric structure. Figure 7(a) shows the gray-scale 2D image stacks. Figure 7(b) shows the result of the background reduction process. Figure 7(c) shows the result of the image segmentation process. The black and white pixels represent fibers and void space, respectively. Figure 7(d) shows the 3D structure of a nonwoven fabric reconstructed by stacking of binarized images. Figure 7(e) shows the contour map of the level set function at a certain cross section. The zero level set represents the surfaces of filaments. Figure 7(f) and (g) show the examples of fiber surfaces before and after applying the fsBR scheme, respectively. The smooth fiber surfaces are recovered from the digitized voxel data by the fsBR scheme.

Figure 7.

The schematic illustrations of the image processing and the reconstructed 3D nonwoven fabric structure: (a) the gray-scale 2D image stack, (b) the result of the background reduction process, (c) the result of the image segmentation process, (d) the 3D structure of a nonwoven fabric reconstructed by stacking of binarized images, (e) the contour map of the level set function at a certain cross section, (f) the fiber surfaces before applying the fsBR scheme and (g) the fiber surfaces after applying the fsBR scheme.

Permeability prediction

The binarized data of nonwoven fabric samples obtained by the image post-processing are transferred by the fsBR scheme to the level set functions representing the smooth fiber surfaces. Then, the IPBB method is applied to describe the fiber surface boundaries for fluid flow calculations. The computational domain size is $Lx\times Ly\times Lz=3.02\hspace{0.17em}\mathrm{mm}\times 0.91\hspace{0.17em}\mathrm{mm}\times 2.59\hspace{0.17em}\mathrm{mm}\left(1400\times 420\times \mathrm{1200\hspace{0.17em}}\mathrm{lattices}\right)$ . As the computational domain size is greater than the Brinkman screening length, which is the square root of permeability (that is shown later in this section), it is confirmed that computational domains are sufficiently representative to smoothen the effect of local inhomogeneities of nonwoven fabrics [44,45]. The pressure gradient boundary conditions are applied in the through-plane direction and the symmetry boundary condition is applied to the remaining four sides of the computational domain. The Reynolds numbers based on the bulk velocity and the pore scale are 0.29, 0.29, 0.09, 0.12, 0.18 and 0.18 for sample IDs of 1–6, respectively. The present calculations are performed by the in-house MPI based multi-GPU code written in CUDA-Fortran programming.

The relation between the flow patterns and the bending stiffness is considered for three types of polymer fibers used to produce the nonwoven fabrics. The bending stiffness according to the standard mechanics is defined as $B=EI$ where $E$ is the fiber’s Young’s modulus and $I$ is the area moment of inertia of the fiber’s cross section. Assuming that the fibers have circular cross section, the area moment of inertia is defined as $I=\pi D^4/64$ where $D$ is the fiber diameter [46]. Thus, the averaged bending stiffness of rayon (2.2 and 3.3 dtex) and PET (2.2 dtex) can be calculated as $\left(1.98\hspace{0.17em}\pm \hspace{0.17em}0.42\right)\times {10}^{-11}$ , $\left(5.76\hspace{0.17em}\pm \hspace{0.17em}1.22\right)\times {10}^{-11}$ and $\left(2.11\hspace{0.17em}\pm \hspace{0.17em}0.53\right)\times {10}^{-11}$ N·m², respectively, referring to the fiber diameter obtained by SEM images and the wet condition Young’s modulus data from the literature [47–51]. Note that Young’s modulus of the rayon fiber decreases in the wet condition [51]. The bending stiffness represents the degree to which fibers are deformed by the impact of the water jets. Figure 8(a) to (c) shows velocity contour maps in a certain cross section of three representative nonwoven fabric samples, which are sample IDs of 1, 3 and 5, respectively. Nonwoven fabric samples of ID 1 and 5 have low bending stiffness and high porosity, and the corresponding high-velocity areas are larger while due to voids generated by the water jets, the nonwoven fabric sample of ID 3 has high bending stiffness and low porosity, resulting in high flow resistance due to the densely entangled fibers.

Figure 8.

Velocity contour maps in a certain cross section of the nonwoven fabric sample IDs of (a) 1, (b) 3 and (c) 5.

The fabric permeability is calculated by Darcy’s law described in equation (14). Note that wettability of the nonwoven fabrics is not considered in the present calculation. Since the contact angles of rayon and PET are approximately 40 and 80 degrees, respectively, the effect of the wettability of those relatively hydrophilic materials on the laminar drag is negligible [52]. The simulated results are compared with the permeability data measured by the standard method for the measurement of hydraulic conductivity of saturated porous materials (ASTM D 5084-03 method A or JIS A 1218). Figure 9 shows the calculated and experimental through-plane permeability of the nonwoven fabrics. The indicated permeability is normalized by the average fiber diameter $D$ , which is measured by SEM images of the nonwoven fabrics. The experimental permeability is represented as the averaged value with the standard deviation obtained by 5–11 measurements. The trend of the calculated permeability shows good agreement with that obtained by the experiment. Figure 10 shows the difference between the experimental and calculated permeability of the fsBR-IPBB and the HWBB methods. The results clearly indicate that the fsBR-IPBB method is more accurate than the HWBB method with the difference of 8.1–16.3% while the difference for the HWBB method is 13.7–23.4%. For the computational error, the lattice resolution is to be considered. As the lattice resolutions are 6–8 in the fiber diameter $D$ , referring to Figure 4 the calculation error of the fsBR-IPBB method is estimated to 3.3–9.1% compared with the result of $D/\Delta =100$ , which we consider to be almost the real value, while the calculation error of the HWBB method is estimated to 15.8–28.5%. Accordingly, we consider that the difference between the measured and simulated data mainly comes from this lattice resolution issue. The lattice resolution of twelve in the fiber diameter though is estimated to be sufficient to accurately calculate the permeability by using the fsBR-IPBB method from Figure 4. To achieve this, it is required the scanning resolution of X-ray CT to be finer. However, we confirm that with the fsBR-IPBB method, the predictive accuracy of the fabric permeability can be estimated by the simple verification in Figure 4. This is a very good feature of the predictive method for the practical usage.

Figure 9.

Comparison between the experimental and calculated values of the fabric permeability in the through-plane direction of the nonwoven fabric samples. Permeability is normalized by the average fiber diameter.

Figure 10.

Difference from the experimental values of the fabric permeability.

Macroscopic correlation models

The experimental and calculated values of the fabric permeability are compared with the profiles of the correlation models of the permeability [5,7,37,45,53,54]. The models of Happel and Drummond and Tahir comprise of a cylindrical unit cell involving flow through ordered arrays of unidirectional cylinders [5,7]. The model of Gebart was derived from first principles for flow in unidirectional reinforcement consisting of regularly ordered fibers [37]. The model of Woudberg was derived by applying a weighted average equation to the existing analytical pore-scale model [53]. Jackson and James empirically extended the model of Drummond and Tahir to regular arrays [54]. The model of Clague et al. was derived from an empirical fitting to their numerical results of ordered arrays [45]. Figure 11 compares the variation of the permeability as a function of porosity $\phi$ . The permeability is normalized by the averaged fiber radius $R$ . The model of Woudberg over-predicts the permeability of the present nonwoven fabrics, while the models of Drummond and Tahir, Happel, Jackson and James and Clague et al. under-predict it. The model of Gebart though is overall in acceptable agreement with the experimental permeability, none of the models sufficiently represent the tendency between the porosity and the permeability of the present nonwoven fabrics. Those analytical or empirical correlation models to predict permeability of 3D fibrous porous media are based on the assumption of microstructures to be either of unidirectional, layered or isotropic random fashion [55]. As hydroentangled nonwoven fabrics prepared in this study have local inhomogeneity due to the water jet nozzle arrangement, fiber entanglement induced by the high-pressure water jet and bending stiffness of polymer fibers, the above assumption leads to the prediction error. Hence, we confirm that it is difficult to predict the fabric permeability by those models [39,56]. Therefore, the computational method demonstrated in this study is more useful for permeability prediction. For further improvement to decrease the error between the experimental and the calculated permeability, a finer scanning resolution of X-ray CT is required. Referring to Figure 4, the accuracy of the fsBR-IPBB and the ideal shapes-IPBB, which discretize the space in voxels, depend on the lattice resolution at the same level. Therefore, we consider that preparing the lattice resolution of at least twelve in the fiber diameter is required to accurately calculate the permeability, which leads to additional computational costs. Although the LBM with unstructured grids adapted to the fiber shapes may more accurately calculate the flow with a lower resolution, such a method in general requires very sophisticated grid generation and much higher computational costs than those of the LBM with voxel-based grids [57–59].

Figure 11.

Comparison of the experimental, calculated and analytical profiles of the fabric through-plane permeability. Permeability is indicated as non-dimensional value normalized by the average fiber radius.

Concluding summary

In this study, computations of the fluid flows in fibrous porous media are presented using the LBM with the fsBR-IPBB method employing a level-set function. The present calculations are performed by in-house MPI based multi-GPU code written in CUDA-Fortran programming language. The fsBR-IPBB method is first verified by computing the permeability of the well-known geometries of cylinder arrays. The results of two test cases show that the permeability calculated by the fsBR-IPBB with the lattice resolution of six in the cylinder diameter are in good agreement with the data reported in the literature. The fsBR-IPBB method is then applied to the fluid flows through nonwoven fabrics produced by the industrial hydroentanglement process. The fabric permeability of six types of nonwoven fabric samples are obtained by the experimental measurements to validate the fsBR-IPBB method. The comparison shows that the fsBR-IPBB method is accurate enough in the permeability prediction of nonwoven fabrics. The prediction error mainly comes from the scan resolution though the error level can be pre-estimated. For further improvement of the accuracy, it is inevitable to perform the computation with finer-resolution X-ray CT scanning. It obviously requires additional costs. The relationship between the permeability and the porosity of the nonwoven fabrics obtained by the experiments and the LBM is found to be different from those obtained by the analytical or numerical correlation models due to the highly anisotropic geometry of the nonwoven fabrics produced by the hydroentanglement process. The numerical method demonstrated in this study is more reliable to predict the fabric permeability of the nonwoven fabrics and can be applied to variety of porous media such as gas diffusion layers of fuel cells, etc.