1 Introduction

Membrane distillation (MD) is a desalination technique that uses a membrane to thermally separate potable water from sea or brackish water. It is considered a promising technology as it operates at lower temperatures ( \(< 100\, ^{\circ }\hbox {C}\) ) than conventional thermal desalination processes and can be driven by solar or waste heat Alkhudhiri et al. (2012). Compared to reverse osmosis, the hydrostatic pressure difference across the membrane is reduced in MD and additionally, the process is more resistant to fouling. However, its energy efficiency and economical performance need to be enhanced for large-scale application Saffarini et al. (2012).

In MD, desalination is performed by generating a temperature gradient between the hot feed and cold fresh water (coolant) which are separated by a membrane which retains water but is permeable to water vapors. (s. Fig. 1). The temperature gradient leads to a gradient in partial vapor pressure which drives evaporation of vapor from the hot feed, diffusion through the membrane. In case of direct contact MD, the vapor is condensed directly into the coolant. Salt is not soluble in water vapor and is therefore retained in the feed Khayet and Matsuura (2011). To enhance fresh water production and make MD more economically viable, multidimensional modeling on module scale is needed including at the same time a detailed and accurate representation of the vapor transport through the membrane Shirazi et al. (2016).

Fig. 1
figure 1

Schematic representation of a direct contact membrane distillation module

Full size image

State-of-the-art membranes used in MD are made of polymers like polytetrafluoroethylene (PTFE) or polyvinylidene fluoride (PVDF) which were originally developed for microfiltration processes Khayet (2011). The membrane needs to be hydrophobic but permeable to vapor as it acts as a separating wall for the liquids but at the same time it should allow for a high transport rate of vapor.

The mass transport processes in MD are commonly described by the dusty gas model (DGM) which takes into account Knudsen diffusion, molecular diffusion, viscous flow and surface diffusion of which the last one is usually neglected Mason and Malinauskas (1983). The most general form of the DGM as it is used in MD is stated by Lawson and Lloyd (1997):

$$\begin{aligned}&\frac{N_i^D}{D_{ie}^K} + \sum _{j=1 \ne i}^{n} \frac{p_jN_i^D-p_iN_j^D}{D_{ije}^m} = - \frac{1}{RT} \nabla p_i \end{aligned}$$
(1)
$$\begin{aligned}&N_i^V = - \frac{p_i B_0}{RT\mu } \nabla P \end{aligned}$$
(2)
$$\begin{aligned}&N_i = N_i^V + N_i^D \end{aligned}$$
(3)
$$\begin{aligned}&D_{ije}^m = K_1 P D_{ij}, \qquad D_{ie}^K = K_0 \sqrt{\frac{8RT}{\pi M_i}} \end{aligned}$$
(4)

where \(N_i^D\) is the diffusive molar flux, \(N_i^V\) the viscous molar flux, \(N_i\) the total flux, \(D_{ij}\) the diffusion coefficient of component i in component j, P the total pressure, \(p_i\) the partial pressure of component i, T the temperature, \(\mu\) the fluid viscosity and \(M_i\) the molar mass.

The three constants \(B_0\), \(K_0\), and \(K_1\) depend on the membrane geometry and are recommended to be measured experimentally. In the lack of exact topological pore description, the geometries are conceptually simplified in equivalent structures composed of uniform cylindrical pores. In this case, the aforementioned constants can be approximated and parametrized using the volume averaged membrane pore radius r, the tortuosity factor \(\tau\) and the porosity \(\epsilon\) Lawson and Lloyd (1997):

$$\begin{aligned} B_0 = \frac{\epsilon \, r^2}{8 \, \tau }, \qquad K_0 = \frac{2 \, \epsilon \, r}{3 \, \tau }, \qquad K_1 = \frac{\epsilon }{\tau } \end{aligned}$$
(5)

A high transport flux is thereby achieved by high porosity, large nominal pore size and low tortuosity. However, with increasing nominal pore size, the pressure at which liquid can penetrate the membrane is decreasing thereby limiting the applicable pore diameter Chapter 8 md membrane characterization (2011).

The tortuosity is a measure of the resistance to diffusion because of the existence of a porous structure and relates the length of the actual tortuous paths through the membrane to the membrane thickness Gommes et al. (2009), Tartakovsky and Dentz (2019). Different definitions of tortuosity exist Clennell (1997). Geometric tortuosity is determined in porous media by identifying two parallel planes and comparing the length of the connecting path \(L_{eff}\) with the distance in between the planes L. For the hydraulic tortuosity, the stream length which a particle is traveling through the porous medium is applied as \(L_{eff}\). Throughout this paper, the tortuosity factor is calculated as \(\tau = \left( mean\left( L_{eff}/L\right) \right) ^2\). In MD, this parameter is usually set to 2 or can be used as an adjusting parameter to fit theoretical models to experimentally measured data (or vice versa) Srisurichan et al. (2006), Khayet (2011), Iversen et al. (1997), Khayet et al. (2002), Varsakelis and Papalexandris (2020), Hitsov et al. (2015).

The values of the constants \(B_0\) and \(K_0\) can be measured by gas permeation experiments Lawson et al. (1995), Lawson and Lloyd (1996), Lawson and Lloyd (1997), Khayet et al. (2004), Hitsov et al. (2015). Many MD models apply the simplification for uniform cylindrical pores Alkhudhiri et al. (2012), Khayet (2011), Srisurichan et al. (2006).

A step toward more realistic membrane representation has been done by performing Monte Carlo simulations with a three-dimensional network of interconnected cylindrical pores, which is used to describe the porous membrane Hitsov et al. (2015), Imdakm and Matsuura (2004), Khayet et al. (2010). In these studies, the pore size distribution is set according to a statistical distribution function. The thermodynamic conditions inside the membrane and the vapor flux have been calculated taking into account viscous flow and Knudsen diffusion. Knudsen diffusion in cylindrical pores of varying diameter has been also studied using Monte Carlo simulations Shi et al. (2012).

For other applications of porous media, permeability and diffusivity values were obtained combining three-dimensional (3D) structural mapping with lattice Boltzmann (LB) computation. In a study by Chen et al. (2015), the porous structure of shales for shale gas extraction is reconstructed using the Markov chain Monte Carlo method based on scanning electron microscopy images. Membrane properties are then computed via lattice Boltzmann simulation. Recently, X-ray tomographic microscopy has been used to analyze gas diffusion layers in polymer electrolyte fuel cells, both for the characterization of the transport properties of materials Rosen et al. (2012), Prasianakis et al. (2012), Chen et al. (2013), Hou et al. (2020), as well as for investigations relevant to water evaporation under operating conditions Safi et al. (2017).

Relevant 3D structural mapping techniques for MD membranes include focused ion beam - scanning electron microscopy (FIBSEM) and X-ray tomography. While FIBSEM is usually easier available and offers a higher spatial resolution, X-ray tomography allows structural mapping and material characterization in a non-destructive way without artifacts from milling. Also, tomography allows for in-situ measurements and is more time-effective. The relevant membranes for MD consist of PTFE which is non-conductive, therefore, repeated metal coating would be needed for meaningful SEM imaging. A detailed comparison of 3D imaging techniques is provided by Fam et al. studying catalysts Fam et al. (2018).

In this work, ptychographic X-ray computed tomography Dierolf et al. (2010) is employed to map the 3D structure of three commercial state-of-the-art PTFE membranes in MD. The tortuosity and viscous flow coefficient of the samples are computed using the lattice Boltzmann model by Prasianakis et al. (2009), Prasianakis et al. (2012). The fully resolved fluid flow simulations are complementary to experimental techniques that are used to characterize the membranes. The obtained membrane properties are compared to the viscous flow coefficient from the DGM. Knudsen effects are considered both in the calculated membrane properties and the DGM with the aim of understanding the impact of the assumptions made by the DGM.

2 Experiments

2.1 Sample Preparation

Three relevant state-of-the-art membranes are identified preceding the experiments (Tab. 1). From each of the membranes, cylindrical volumes are prepared that include the full thickness of the membrane as the height (samples 1-3). The diameter of the surface area is chosen to achieve a significant aspect ratio for the numerical simulation while reducing the overall measuring volume to decrease computational and experimental resources. Of the first membrane, a second sample (sample 4) was prepared whose height is only 20 \(\upmu\) m. Instead, it covers a larger surface area.

Table 1 Structural specifications of the imaged membrane samples: Pore size according to manufacturer; height and diameter describe the dimensions of the prepared cylindrical samples such that the height is measured from the PP grid in direction of the membrane thickness
Full size table

The membranes consist of PTFE which is stretched and attached onto a coarser polypropylene (PP) grid. The PTFE itself is not rigid but held in place by the PP grid. Therefore, embedding of the PTFE with epoxy resin is necessary prior to the manufacturing and extraction of cylindrical samples. EPON is used as epoxy resin because of its low viscosity and expansion coefficient to minimize any potential structural changes. The embedding is done in an iterative process to assure the complete filling of all membrane pores. In a first step, membrane pieces of about \(10 \cdot 10\) mm\(^{2}\) are submerged in a mixture of 50% ethanol and 50% epoxy applying vacuum cycling. In a second and third step, this procedure is repeated twice with a pure epoxy bath and vacuum cycling. Afterward, the specimens are set to harden for 24h at \(60~^{\circ }\hbox {C}\).

For the extraction of the cylindrical membrane samples, the epoxy filled membrane pieces are pre-cut with a razor blade to ensure a clean surface and metal-coated on the top-surface. Close to the clean surface, cylindrical samples are milled using plasma focused ion beam cutting. Micro-manipulators are used to transfer the samples onto the pins for ptychography. Meanwhile, gas injection is applied for carbon removal.

2.2 Ptychographic X-ray Computed Tomography

The 3D structures of the samples are obtained using ptychographic X-ray computed tomography. Experiments were conducted at the cSAXS beamline of the Swiss Light Source, PSI, Switzerland. Ptychographic X-ray computed tomography is a “scanning” coherent diffraction imaging technique that provides quantitative 3D density maps Dierolf et al. (2010). For the tomographic reconstructions, iterative phase retrieval algorithms were applied and tomograms of phase and amplitude contrast are obtained by combining multiple projection angles Guizar-Sicairos and Fienup (2008), Thibault et al. (2008). Advantages of this technique lie in the high spatial resolution and the obtained quantitative tomograms Dierolf et al. (2010), Leu et al. (2019).

The experimental setup described by Holler et al. (2014), Holler and Raabe (2015) is extended by a liquid nitrogen cooled and depressurized specimen environment. Thereby, the samples remain in vacuum and at \(-180\,^{\circ }\hbox {C}\) during the entire measurement. The structure of the PTFE at room temperature is preserved as the specimen are embedded in an epoxy matrix.

For the experiments described in this paper, a photon energy of 6.2 keV is used. Ptychographic scans consisted at minimum of 415 diffraction patterns per projection with an exposure time of 0.1 seconds per diffraction pattern. Sampling positions were set using a Fermat spiral scanning grid with an average step size of 2.5 microns resulting in a total acquisition time between 12 and 16.5h per sample depending on the sample volume. The beam size at the sample plane was 4 microns in diameter. The total number of projections of each sample is shown in Table 2.

Tomographic reconstruction was performed as described by Guizar-Sicairos et al. (2011). The resulting voxel size is 38.99 nm. The spatial resolution of the individual tomograms as determined by Fourier shell correlation is found in Table 2.

Larger sample diameters, due to the increased beam exposure, i.e., increased number of required projections to reach a set spatial resolution, enhance the extent and likelihood of beam damage. This favors as such a loss in the achievable spatial resolution which can be seen by comparing sample 1 and sample 4 in Table 2 for example. As such it was favorable to keep the sample diameter as small as possible while increasing the sample height to image a sample representative volume.

2.3 Segmentation

The electron density difference between PTFE and epoxy is sufficient to provide decent contrast in the electron density tomograms (Fig. 2). The PTFE structure is obtained from the tomograms by pixel intensity thresholding aiming at the inclusion of all PTFE fibers for the subsequent flow analysis. A square cuboid is then extracted for the computational domain (Figs. 3, 4). Then, the porosity of the membrane samples can be calculated by relating the sum of all void pixels to the total number of pixels. The segmented tomograms are shown in Figs. 5, 6, 7, and 8 and are openly accessible via https://doi.org/10.5281/zenodo.4467742.

Table 2 Imaging specifications of the membrane samples: Total number of projections; estimation of the resolution according to Fourier shell correlation Guizar-Sicairos et al. (2011), Guizar-Sicairos et al. (2015); domain dimensions (\(x\cdot y\cdot z\)) for the subsequent numerical analysis first in \(\upmu\) m and below in voxels, x normal to membrane plane
Full size table
Fig. 2
figure 2

Histogram of sample 3, each bin contains 256 pixel intensity values

Full size image
Fig. 3
figure 3

A yz orthoslice from sample 2; the membrane material is visible in light gray and the epoxy, i.e., pore space, in dark gray

Full size image
Fig. 4
figure 4

The same slice after thresholding; as computational domain a rectangular area within the membrane material is selected throughout all slices. The length of the rectangle for sample 2 is 15.60 \(\upmu\) m

Full size image
Fig. 5
figure 5

Representation of the segmented tomograms and domains for the computation of the membrane transport properties, Sample 1 (\(59.85 \cdot 22.42\cdot 22.42\) \(\upmu\) m)

Full size image
Fig. 6
figure 6

Representation of the segmented tomograms and domains for the computation of the membrane transport properties, Sample 2 (\(57.70 \cdot 15.60 \cdot 15.60\) \(\upmu\) m)

Full size image
Fig. 7
figure 7

Representation of the segmented tomograms and domains for the computation of the membrane transport properties, Sample 3 (\(62.42 \cdot 21.44 \cdot 21.44\) \(\upmu\) m)

Full size image
Fig. 8
figure 8

Representation of the segmented tomograms and domains for the computation of the membrane transport properties, Sample 4 (\(25.34 \cdot 37.04 \cdot 37.04\) \(\upmu\) m)

Full size image

3 Numerical Model

The lattice Boltzmann method is used to estimate the permeability of the structures. This method is a mesoscopic computational fluid dynamics method stemming from the Boltzmann equation Succi (2001). The velocity field in the membrane samples is computed using the LB-BGK guided equilibrium model, described in detail by Prasianakis et al. (2009), Prasianakis et al. (2012). LB solves the Boltzmann equation for groups of particles on a regular Cartesian grid (lattice). A D3Q27 lattice is applied here, meaning a 3D lattice where in each lattice point 27 discrete velocities are defined that point to different neighboring points. Thereby, the computation consists of two steps: firstly, the populations move to the next lattice points with the respective discrete velocity. Afterward, their interaction with the neighboring populations is computed in the collision step. As the computation is only dependent on the next neighboring computational grid points, its parallelization is straightforward and it allows for efficient computing of large lattices.

The lattice is determined by the segmented tomograms in a way that each pixel corresponds to a lattice grid point. Lattice points where membrane material is present are labeled as solid structure for the computation. At the solid-fluid interface, no-slip flow is assumed which is achieved by half-way bounce back boundary conditions given small Reynolds numbers (\(Re<<1\)) He et al. (1997). On the domain boundaries, periodic boundary conditions are applied.

For each velocity and lattice point, a velocity population is defined which can be viewed as the share of molecules which move in the direction of the specific discrete velocity at this lattice point. These velocity populations \(f_i\) are calculated by the lattice BGK equation (Eq. 6) where \(c_i\) refers to the discrete velocities for \(i = \{1,2...,27\}\) and \(\delta t =1\) in lattice units.

$$\begin{aligned} f_i \left( \mathbf {x} + \mathbf {c}_i \delta t, t+\delta t \right) - f_i\left( \mathbf {x},t\right) = + \frac{2 \delta t}{\delta t+2\tau _{LB}}\left( f_i^{eq} \left( \mathbf {x},t\right) - f_i \left( \mathbf {x},t\right) \right) \end{aligned}$$
(6)

For the calculation of the equilibrium populations \(f_i^{eq}\) the guided equilibrium model is used as it offers optimal accuracy and minimum deviations from the target Navier-Stokes equations Prasianakis et al. (2009). The relaxation parameter \(\tau _{LB}\) is related to the dynamic viscosity by \(\mu = T_{0}\rho \tau _{LB}\), where \(T_{0}=1/3\) is the lattice temperature.

The macroscopic density \(\rho\) and velocity \(\mathbf {u}\) can then be calculated from the velocity populations as shown below:

$$\begin{aligned} \rho = \sum _{i=0}^{26} f_{i} \qquad \rho \mathbf {u} = \sum _{i=0}^{26} f_{i}\mathbf {c}_{i} \end{aligned}$$
(7)

A constant forcing term \(F_i\) accelerates the initially stagnant fluid in the desired direction similar to the effect of an applied pressure gradient, as described in detail in Prasianakis et al. (2012). Driving the flow with a forcing field allows to apply periodic boundary conditions to both ends of the domain, and does not require the numerical implementation of inlet and outlet boundary conditions which could introduce small inaccuracies depending on the boundary condition implementation.

The calculation using fine resolution geometries, which are described by \(\sim\) 0.5 billion computational cells (see Table 2) require a considerable amount of computer memory. For sample 4, 270GB of RAM were needed and the computation was therefore, executed on dedicated large memory nodes on the ULHPC’s high performance cluster.

3.1 Permeability Calculation

In the continuum regime, Darcy’s law can be used to compute the permeability matrix for small Reynolds numbers (\(Re<<1\)) Rosen et al. (2012):

$$\begin{aligned} \kappa _{ab} = \frac{\mu \cdot u_{b,mean}}{\partial p/ \partial a} \end{aligned}$$
(8)

where \(\kappa _{ab}\) is the intrinsic permeability for pressure gradient \(\partial p/ \partial a\) in direction a and mean velocity \(u_{b,mean}\) in direction b. Eq. 8 is applied for the calculation of the membrane samples permeability in the present work using the results from the LB computation. For membrane applications, \(\kappa _{xx}\) is of interest, representing the permeability for acceleration and resulting velocity in x direction. The \(\kappa _{xx}\) represents exactly the parameter \(B_0\) in the viscous flow formulation of the DGM as seen in Eq. 2. The detailed micropore fluid dynamics simulations using the real microscopic structure obtained via tomography, allow to calculate accurately the intrinsic permeability of the structure. In this work, the computationally calculated permeability values obtained for the membrane samples (using Eq. 8) are presented and compared to the ones obtained using the approximation of \(B_0\) for uniform cylinders calculated with Eq. 4.

The permeability obtained for a fluid flowing through the structure in the Darcy flow regime, neglecting any microflow slippage effects is the so-called intrinsic permeability. For gas flows in porous media, however, a higher permeability is found due to slip-velocity effects on the boundary with the solid structure Klinkenberg (1941), Fathi et al. (2012), Li and Sultan (2017). This apparent permeability depends on the flow regime. The standard lattice Boltzmann implementation with the use of the diffusive boundary condition Ansumali and Karlin (2002), Zhang (2011), can successfully describe microflow effects up to a certain Knudsen number Kn \(\sim\) 0.1 Prasianakis and Ansumali (2011), Toschi and Succi (2005). For the present study, it is of interest to evaluate the permeability at Knudsen numbers up to Kn \(\sim\) 1, hence the intrinsic permeability is corrected using correlations from the literature. A correction factor depending on the Knudsen number is proposed by Beskok and Karniadakis (1999) to correlate intrinsic and apparent permeability \(\kappa _{a,ab}\).

$$\begin{aligned} \frac{\kappa _{a,ab}}{\kappa _{ab}} = \left( 1 + \alpha \; Kn \right) \left( 1 + \frac{4 \; Kn}{1-b \; Kn}\right) \end{aligned}$$
(9)

where \(b = -1\) and \(Kn = \frac{\lambda }{L}\) is the Knudsen number relating the mean free path \(\uplambda\) to the characteristic problem length L. The mean free path can be calculated by \(\lambda = \frac{k_B T}{\sqrt{2}\pi d^2 p}\) with \(k_B\) as the Boltzmann constant, T as the temperature and d as the collision diameter.

The parameter \(\alpha\) varies from zero for low Knudsen numbers to an asymptotic value \(\alpha _0\) for large Knudsen numbers (Eq. 10):

$$\begin{aligned} \alpha = \alpha _0 \frac{2}{\pi }\tan ^{-1}\left( \alpha _1 Kn^\beta \right) \qquad \alpha _0 = \left[ \frac{64}{3\pi \left( 1-4/b\right) }\right] \end{aligned}$$
(10)

with \(\alpha _1 = 4.0\) and \(\beta = 0.4\) Beskok and Karniadakis (1999).

Knudsen numbers for MD application lie usually between 0.3 and 0.8. For this range, the apparent permeability of the four membrane samples is computed from the intrinsic permeability according to Eq. 9.

3.2 Tortuosity Calculation

The generation of the connecting paths is facilitated using the LB simulation results stemming from the calculation of the intrinsic permeability as discussed in the previous section. During the simulation, a flow field is generated due to the presence of a forcing field, which acts parallel to the desired measurement direction. The planes perpendicular to the direction of the forcing field, at the sides of the domain are used as the locations where the starting and the end points of the pathways lie. Their distance is defined as L representing the length and thickness of the membrane. As pathways, we use the streamlines of the velocity flowfield which define the \(L_{eff}\) Prasianakis et al. (2012). Following this methodology, the determined tortuosity factors represent the hydraulic tortuosity as defined by Clennell (1997). For these domains, the number of the generated pathways was approximately 20,000.

4 Results

The resulting membrane properties for the four membrane samples are presented in Tab. 3.

Table 3 Resulting porosity, tortuosity and intrinsic permeability values; \(\upkappa _{\mathrm{xy}}\) relates the pressure gradient in x direction to the velocity in y direction
Full size table

The tortuosity factor of the examined membranes is calculated to be slightly larger compared to the reference value of 2 which is commonly assumed for MD studies Khayet (2011). In Fig. 9, the values obtained in this work are compared to common correlations used to calculate the tortuosity from the porosity. For the analyzed membrane samples, the use of the reference value of 2 is enforced compared to the correlations which underpredict the tortuosity. As permeability is inverse-proportional to the tortuosity, permeability and diffusivity are thereby overpredicted.

Fig. 9
figure 9

Comparison of tortuosity factors obtained here and by different correlations: \(\uptau\) obtained in this work, \(\tau = \frac{\left( 2-\epsilon \right) ^2}{\epsilon }\) Iversen et al. (1997), \(\tau = \frac{1}{\epsilon }\) Iversen et al. (1997), \(\uptau\) value commonly used in MD studies Khayet (2011); correlations computed with the porosity reported by the manufacturer

Full size image

The porosity matches the value specified by the manufacturer for sample 2, 3 and 4 (manufacturer values are 85% for sample 1, 3 and 4 and 80% for sample 2). If one takes a closer look at the porosity distribution in axial direction (i.e., the number of void pixels in each yz slice) one finds a heterogeneous distribution (s. Fig. 10). For all samples, a high fluctuation is visible and the porosity close to the PP grid is increased compared to the free surface (\(x = 0\)). This behavior is more profound for samples 1 and 2 while sample 3 and 4 are more homogeneous. A possible explanation is that the investigated yz plane area does not span a representative area of the membrane. Section 4.1 provides a more detailed investigation of the representative elementary volume. The reason for the lower porosity at the free surface might also be found in the composition of the membrane. While the PTFE is stretched and attached on the PP grid, it can bend back on the free surface. The investigation of more membrane samples is needed to fully explain the porosity distribution found in this study.

Fig. 10
figure 10

Porosity values calculated for each yz slice plotted against the axial position: Sample, Arithmetic mean value, Manufacturer value; the x axis starts at the membrane top side and ends at the PP grid side; for comparison the porosity specified by the manufacturer is included in red

Full size image

The values for the intrinsic permeability \(\kappa _{xx}\) lie in the same range as the values reported by Lawson et al. (1995) for other MD membranes. However, compared to the DGM in Fig. 11, two observations can be made: Firstly, the here determined intrinsic permeability values show the same trends for sample 1–3 with sample 3 having the highest permeability of all samples. Sample 4 is manufactured from the same membrane material as sample 1 but differs in the height and aspect ratio. The larger permeability of sample 4 gives reason to believe that also the permeability is increased closer to the PP grid as is the porosity. Secondly, the intrinsic permeability values are consistently smaller than predicted by the DGM. For sample 3, the intrinsic permeability of the mapped volume is even one order of magnitude smaller than predicted by the DGM. Sample 4 shows the closest match. The discrepancy between DGM and computed data needs therefore to be better understood. The LB computation determines the intrinsic permeability very accurately for the investigated geometry by considering all microstructural features and the pore-space topology. The viscous flow contribution of the DGM \(B_0\) is calculated by using the structural properties provided by the manufacturer and stated in Table 1, the standard methodology in MD analysis. Recalling Eq. 5, \(B_0\) is calculated from pore size, porosity and tortuosity factor of which porosity and tortuosity values of the investigated membrane samples are in average agreeing well with the respective manufacturer values as shown in Figs. 9 and 10. The remaining parameter, the pore radius, has a quadratic influence on the viscous flow contribution of the DGM. The pore sizes at which the DGM would match the intrinsic permeability of the membrane samples are shown in Table 4. It can be seen that the permeability of the membrane structure corresponds to a porous medium of roughly a third to half the nominal effective pore size according to the DGM.

Table 4 Fitted pore size values at which the viscous flow contribution of the DGM \(B_0\) equals the calculated intrinsic permeability of the investigated membrane samples in Fig. 11; nominal pore size provided by the manufacturer
Full size table

The differences between \(\upkappa _{\mathrm{xy}}\) and \(\upkappa _{\mathrm{xz}}\) in Table 3 show a preferential ordering of the membrane fibers in the yz plane likely due to manufacturing reasons.

Fig. 11
figure 11

The intrinsic permeability \(\kappa _{xx}\) obtained here is compared to the permeability calculated by the DGM: Intrinsic permeability \(\kappa _{xx}\), \(B_0\) of DGM

Full size image

Correcting the intrinsic velocity for velocity slip effects at the solid boundaries, the apparent permeability as a function of Knudsen number is obtained (Eq. 9). The Knudsen number depends on the physical properties of the diffusing molecules and changes with pore size, temperature and pressure variations. For a porous medium with a given pore size and diffusing species, the variation of Knudsen number is due to changes in temperature. Thereby, for a porous medium with pore size \(d= 0.2 \, \mu m\), the Knudsen number varies between 0.7 and 0.8 for temperature changes between \(40 - 80~^{\circ }\hbox {C}\) at ambient pressure and diffusion of water vapor in saturated air. In Fig. 12, the apparent permeability is presented for Knudsen numbers from 0.3 to 1. For reference purposes, the permeability values \(B_0\) from the DGM are also indicated which scale with geometrical properties but are independent of pressure and temperature. For a given porous medium, \(B_0\) is independent of Knudsen number.

The DGM takes Knudsen effects into account by considering Knudsen diffusion as shown in the introduction, however, not scaling explicitly with Knudsen number. A permeability comparison including Knudsen effects can be performed as a function of temperature when considering the transport of a single species at a given pressure. The partial pressure gradient then equals the total pressure gradient. Under this assumption, the permeability of the DGM including Knudsen diffusion can be calculated as \(B_0 + D_{ie}^K \cdot \mu /p_i\) following Eq. 1, 2 and 3. For water vapor at ambient pressure, the resulting permeability values are shown in Fig. 13 for the common temperature range of MD. It should be noted that in this temperature range boiling conditions are not reached at ambient pressure and this scenario of pure water vapor is therefore thermodynamically not feasible but shown for comparison purposes. A closer fit can be seen for sample 4 while sample 3 with the larger pore size still shows the largest deviation. Due to the fact that the nominal pore diameter of sample 3 is more than twice the pore diameter of sample 1, 2 and 4, the Knudsen number at a given temperature is significantly lower for sample 3. The correction factor for the apparent permeability is therefore also lower for sample 3 to the point that the apparent permeability of sample 4 becomes larger than of sample 3 for all considered temperatures.

Fig. 12
figure 12

The apparent permeability \(\kappa _{a,xx}\) as a function of Knudsen number for Knudsen numbers from 0.3 to 1 and \(B_0\) of the DGM; for a porous medium with pore size \(d= 0.2 \, \mu m\) the Knudsen number varies between 0.7 – 0.8 for temperature changes between \(40- 80~^{\circ }\hbox {C}\) at ambient pressure as indicated by the bar above the graph

Full size image
Fig. 13
figure 13

The apparent permeability \(\kappa _{a,xx}\) and the DGM (\(B_0\) and Knudsen diffusion contribution) as a function of temperature for water vapor at ambient pressure

Full size image

4.1 Homogeneity of the Samples

To evaluate the representative elementary volume (REV) of the membranes, smaller volume sections are analyzed keeping the original aspect ratio and resolution. For each sample, three subsamples with half the original side lengths are calculated spanning an eighth of the total volume (see Fig. 14). Additionally, six subsamples with a quarter of the original side lengths are calculated representing an 64\(^{\mathrm{th}}\) of the original volume. The location of these subsamples is chosen randomly within the complete sample volume. In Fig. 15, the resulting porosity and permeability of these subsamples are presented for each sample. The average values for both subsample seizes are indicated together with the results for the full sample. By comparing the green and red circular marks, it can be seen that the full volume is well approximated by the average of at least three subsamples of an eighth volume. The average of the smaller subsamples (1/64\(^{\mathrm{th}}\)) is only representative for sample 3 and 4. A simulation using (1/64\(^{\mathrm{th}}\)) of the tomogram is at least two orders of magnitude faster compared to the simulation using the full detailed tomogram description. This is due to the fact that for the lattice Boltzmann method, the computational effort scales linearly with the number of grid points. Moreover, steady state can be reached in fewer time steps in smaller geometries. The results for sample 1 and sample 2 indicate not only an inhomogeneous distribution of porosity as already shown above but also permeability values variations of an order of magnitude.

Fig. 14
figure 14

Graphical representation of a subsample with half the original side lengths; from the full volume, subsamples are extracted with constant aspect ratio and randomly chosen location within the full sample

Full size image
Fig. 15
figure 15

Intrinsic permeability and porosity for different volume sizes: Full volume, 1/8 volume, 1/64 volume, 1/8 volume average value, 1/64 volume average value

Full size image

5 Discussion of the Results

The porosity values listed in section 4 agree well with the porosity specified by the manufacturer for three of the samples. The porosity of one sample is below the specified manufacturer value. Considering the porosity across slices perpendicular to the main flow direction, variations are observed. Closer to the PP grid, the porosity is increased compared to the mean value while it is reduced on the free surface.

The tortuosity factor lies for all membrane samples slightly above 2 and justifies the use of the reference value of 2 which is commonly assumed for MD studies Khayet (2011).

The intrinsic permeability is in the same order of magnitude but lower than predicted by the DGM for three samples. For the membrane sample with the larger pore diameter, the DGM predicts the permeability to be one order of magnitude higher which is not reflected in the computations. Similar to the porosity, also the permeability is increased closer to the PP grid compared to the overall membrane height. The permeability obtained in the calculations match the DGM when considering a membrane with reduced effective pore size (a third to half the nominal pore size specified by the manufacturer). Next to tortuosity and porosity which have been analyzed, the effective pore size of the calculated membranes could be compared to the nominal pore diameter provided by the manufacturer as this parameter has a significant impact on the permeability values of the DGM.

The calculated membrane properties show variations along the main flow direction with increasing porosity and permeability at the attachment to the PP grid while decreasing toward the free surface. Sample 4 consists of the same material as sample 1 but does span less height and a larger surface area. It shows the least variations in porosity and the closest fit to the permeability of the DGM. Part of this inhomogeneity along the membrane thickness might be explained by the composition of the membrane as the PP grid preserves the porous shape of the PTFE while back-bending is possible on the free surface. Similar to porosity and permeability, the pore size might be reduced closer to the free surface compared to the nominal pore size and the investigation of the effective pore size distribution could add toward resolving the discrepancy of the permeability values calculated here and predicted by the DGM. If structural analysis proved also an inhomogeneous distribution of pore size, the application of nominal pore size in the DGM should be considered with care in a favor of reduced pore size as transport through smaller pores could have more impact on the permeability. In case of a homogeneous pore size distribution, the question arises whether porosity and tortuosity value alone are sufficient for abstracting these membranes from ideal cylindrical form in the permeability calculation. We finally note that different pore geometries, which are characterized by similar porosity and pore size may differ at micro-structural level and transport properties. In this case, the DGM would not differentiate between the different geometries, while the direct micro-structural fluid dynamics approach does.

For optimal results, the membrane volumes to be mapped with ptychography should have a low aspect ratio and should not exceed the total volume measured in this work. To achieve significant results multiple samples of each membrane material should be mapped. As it is shown here, three volumes of an eighth volume can be averaged to obtain the same results as analyzing the full volume mapped in this work. This procedure also reduces computational requirements. Additionally, special care must be taken to extract the membrane samples from sufficiently large membrane cuts to reduce possible bending of PTFE on the free surface.

Ultra-high resolution fluid dynamics simulations, using the LB methodology, allows to incorporate phase change phenomena to predict evaporation rates at different flow conditions Safi et al. (2017), Yu et al. (2017), Zhang et al. (2015), Qin et al. (2019). Heat transport and humidity transport in the gas phase can be incorporated as passive scalar within the LB framework. At the same time it might the solid material and the heat transport through the solid matrix might play a significant role under circumstances. For that, heat conduction through the solid phase can also be described within this framework. Finally, upon the addition of geochemical reactions and deposition of solids, e.g., due to precipitation Prasianakis et al. (2017), Jiang and Tsuji (2014), the list of necessary elements for predicting membrane degradation and fouling, shall be complete. We shall address these topics in our future research.

6 Conclusion

In this work, three state-of-the-art MD membranes are mapped using X-ray ptychography and their 3D structures are obtained. Using these structures, the flow field is calculated driven by a forcing field using lattice Boltzmann computation. In this way, porosity, tortuosity and permeability values are determined. The porosity values agree well with manufacturer values for three samples. One sample shows some deviations. An inhomogeneous distribution of porosity along the membrane thickness is observed for all samples. The tortuosity factors are slightly higher than the literature value of 2 but justify its usage. The LB computed intrinsic permeability is compared to the values calculated using the DGM. The calculated permeability is smaller than predicted by the DGM, however, for three samples in the same order of magnitude. The results support the application of the DGM for fluid-dynamical simulations of large-scale MD modules in general. As for the porosity, the results indicate permeability variations along the membrane thickness. The permeability obtained in the calculations match the DGM when considering a membrane with reduced effective pore size (a third to half the nominal pore size specified by the manufacturer). Structural analysis and investigation of pore size distribution and mean pore size of the here obtained 3D membrane structures can be performed for a better understanding of the pore network and related membrane properties but is out of the scope of this paper. If structural analysis proved also an inhomogeneous distribution of pore size, the application of nominal pore size in the DGM should be considered with care in favor of a reduced pore size as transport through smaller pores could have more impact on the permeability. In case of a homogeneous pore size distribution, the question arises whether porosity and tortuosity value alone are sufficient for abstracting these membranes from ideal cylindrical form in the permeability calculation.

The results prove the capability of the described method using X-ray ptychography and LB to investigate in depth membrane properties used for membrane distillation processes. They provide understanding of the fluid-structure interaction processes at the microscopic level. Such techniques are the basis for future multiphysics MD simulations, which can include the microstructural heterogeneities to full detail while also including phase change phenomena and heat transfer in the membrane.