Computational Fluid Dynamics Modeling of the Filtration of 2D Materials Using Hollow Fiber Membranes

1. Introduction

Membrane technology has attained promising ubiquity throughout various research areas and industries, e.g., water treatment [1,2,3] , gas separation [4,5] , food processing [6,7,8,9] , pharmaceutical industry [10,11] , etc., due to its substantial advantages over other filtration techniques. The merits of membranes over other conventional filtration technologies are but are not limited to cost and energy efficiency, selectivity, environmental friendliness, and so forth [12]. In addition, the membranes could be tailored for the filtration of specific components by adjustment of their pore size. Despite their significant benefits over other filtration methods, the expansive utilization of membranes is limited by several obstacles. A major impediment in the efficient operation of membranes is the fouling phenomenon caused by the deposition and the accumulation of colloids, salts, macromolecules, etc., on the membrane’s inner wall [13,14,15,16]. The fouling leads to a reduction of the permeate flux, which reduces the efficiency of the membrane performance and significantly increases the maintenance and operating costs, and shortens the membrane lifespan.

The energy efficiency and the environmental friendliness of membrane filtration have enabled them to penetrate into the growing areas of green and sustainable manufacturing, such as the treatment of wastewater effluent in the food industry [17] , recovery of the bio-sourced precursor (succinic acid) of bio-molecules [18] , recovery of the lignin which is the among the most essential bio-fuel compounds [19] , and so forth. Recently, Downing et al. [20] extended the applications of microfiltration membranes to the isolation of 2D materials (graphene) through cross-flow filtration (CFF). Their study revealed that the separation of the graphene flakes from their polydisperse solution, containing graphene flakes and incomplete-exfoliated graphites, significantly improved the energy and cost-efficiency of graphene production compared to other conventional methods. Moreover, as the isolated graphene flakes were collected from the permeate stream, this CFF demonstration avoids the waste of the permeate flow, which is typically discarded in CFF applications. Although the integration of microfiltration technique for the production of 2D materials might be a promising pathway for the advancement of sustainable manufacturing, due to the nascency of this topic, there is a lack of sufficient understanding of the physics of fouling of layered structure materials, e.g., graphite, optimization of the experimental setup, appropriate membrane materials, etc. This lack of knowledge necessitates the development of reliable models to predict the fouling, and its impacts on the CFF performance for the further advancement of microfiltration of 2D materials.

In general, due to the detrimental effects of fouling on the membrane characteristics, modeling the membrane filtration processes to optimize the membrane operation and minimize the adverse fouling effects has become an integral part of membrane studies. Various models have been proposed to simulate the fouling on the membrane, e.g., film theory, shear-induced model, resistance-in-series, Deryaguin-Landau-Verwey-Overbeek (DLVO) theory, etc [21,22]. Among these, the resistance-in-series has received the most notable attention due to its capacity to accommodate different fouling physics, where the transmembrane pressure (TMP) relates to the combination of the resistances posed by the particles deposition and accumulations on the membrane surface through the Darcy relation. For example, different combinations of concentration polarization of the foulants, gel layer formation, pore blockage, cake formation, etc. [23,24,25,26,27] , have been attributed to the fouling physics in the membrane and have been integrated into the computational fluid dynamics (CFD) modeling of microfiltration. However, to the best of our knowledge, no CFD model attempted to simulate the fouling of layered structure materials such as graphite and relate the concentration of 2D solutes in the permeate flow to the fouling behavior.

The fouling of graphite demonstrated different effects on the performance of the membrane compared to fouling of other macromolecules in conventional microfiltration applications, e.g., spherical colloids, proteins, polysaccharides, microorganisms, etc. In the conventional cases, fouling often leads to an increase in the transmembrane pressure (TMP) to maintain the desired permeate flux [28,29,30,31,32]. In contrast, the fouling of graphite was demonstrated to lead to a slight decrement of TMP in some stages of fouling (please refer to Figure S2 of the reference [20]). Moreover, the solute concentration in the permeate stream which is often neglected, is crucial for examining the filtration performance. To fill this gap and set the stage for optimizing filtration performance, we laid the initial step to simulate the filtration of exfoliated graphenes from the graphite-graphene dispersion. This manuscript is structured as follows: 1 - In the methods section, we describe our assumptions and mathematical expressions. 2 - The Results and Discussion section presents the reproduction of experimental results and the impact of key simulation parameters. 3 - The conclusion section outlines the fundamental findings of the paper.

2. Methods

In this study, we employed a 2-dimensional domain to model a hollow fiber membrane. The domain, which represents the membrane’s lumen side, consists of four edges corresponding to the inlet, outlet, symmetric axis, and inner wall of the membrane (Figure 1). The membrane characteristics and the experimental setup parameters are obtained from Table S1 and Table S2 of Downing’s report [20]. We assumed a uniform distribution of the inlet flow to each individual fiber inside the cartridge, resulting in an inlet average velocity ($U_{inlet}$) of 0.65 m/s, calculated using the equation $U_{inlet}$ = $4Q_{inlet}$/(N$\pi$$D_f$), where $Q_{inlet}$ is the inlet flow rate, N is the number of fibers, and $D_f$ is the diameter of the fiber. Furthermore, we assumed that each fiber contributes equally to the permeate flux, leading to a permeate velocity ($U_{perm}$) of 5.5×10${}^{-6}$ m/s quantified by the $Q_{perm}$/($N\pi$$D_f$$L_f$) relation, where $Q_{perm}$ is the experimentally set permeate flow rate, and $L_f$ is the length of the fiber. The numerical simulations were performed using COMSOL Multiphysics, version 5.5, to solve the governing equations in the axis-symmetric cylindrical (r, z) coordinate system.

2.1. Governing equations

The velocity and pressure profiles throughout the simulation time are obtained by the numerically solving of continuity equation (Equation (1)), Navier-stokes in the radial direction (Equation (2)), and z-direction (Equation (3)).

$$\frac{1}{r}\frac{\partial \left(u_{r}r\right)}{\partial r}+\frac{\partial u_z}{\partial z}=0$$

(1)

$$\rho \left(\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+u_z\frac{\partial u_r}{\partial z}\right)=-\frac{\partial p}{\partial r}+\mu \left(\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_r}{\partial r}\right)+\frac{{\partial }^2{u}_r}{\partial z^2}\right)$$

(2)

$$\rho \left(\frac{\partial u_z}{\partial t}+u_r\frac{\partial u_z}{\partial r}+u_z\frac{\partial u_z}{\partial z}\right)=-\frac{\partial p}{\partial z}+\mu \left(\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_z}{\partial r}\right)+\frac{{\partial }^2{u}_z}{\partial z^2}\right)$$

(3)

The equations of mass conservation for the permeating (exfoliated graphene) and fouling (non-exfoliated graphite) components were coupled with the Navier-Stokes equation to determine the concentration profile of each material. The subscript i is used to denote either the permeating or fouling material.

$$\frac{\partial C_i}{\partial t}+u_r\frac{\partial C_i}{\partial r}+u_z\frac{\partial C_i}{\partial z}=D_i\left(\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial C_i}{\partial r}\right)+\frac{{\partial }^2{C}_i}{\partial z^2}\right)$$

(4)

The continuity and Navier-Stokes equations were subjected to the following boundary conditions. At the inlet, the flow was assumed to be fully developed.

$$\mathrm{At}z=0:u_z=u_{z,max}\left(1-\left(\frac{r}{R_f}\right)\right),u_r=0$$

(5)

At the outlet, the atmospheric pressure assumption was applied.

$$\mathrm{At}z=L_f:P=P_{atm}$$

(6)

The symmetric condition was set at the centerline of the fiber.

$$\mathrm{At}r=0:\frac{\partial u_z}{\partial r}=0,u_r=0$$

(7)

And, at the surface of the membrane, the no-slip condition was set, and the radial velocity was fixed to the permeate velocity ($U_{perm}$)

$$\mathrm{At}r=R_f:u_z=0,u_r=U_{perm}$$

(8)

In addition, for both the permeating (graphene) and fouling (non-exfoliated graphite) components, the feed concentration was set to the inlet, a non-diffusive flux condition was assumed at the outlet, and the symmetric condition was applied to the centerline of the fiber.

$$\mathrm{At}z=0:C_i=C_{i,feed}$$

(9)

$$\mathrm{At}z=L_f:-D_i\frac{\partial C_i}{\partial z}=-D_i\frac{\partial C_i}{\partial r}=0$$

(10)

$$\mathrm{At}r=0:\frac{\partial C_i}{\partial r}=0$$

(11)

We assumed all the fouling particles are completely rejected at the membrane wall.

$$\mathrm{At}r=R_f:u_r{C}_f-D_f\frac{\partial C_f}{\partial r}=0$$

(12)

In contrast, we set the mass flux continuity boundary condition for permeating particles.

$$\mathrm{At}r=R_f:u_r{C}_p-D_p\frac{\partial C_p}{\partial r}=J_p$$

(13)

$$J_p=U_{perm}{C}_{p,perm}$$

(14)

We established a relationship between the permeating solute concentration at the permeate side ($C_{p,perm}$) and the concentration on the lumen side ($C_{p,r=R_f}$) through the definition of the membrane partition coefficient ($K_{mem}$) for the permeating material.

$$K_{mem}=\frac{C_{p,perm}}{C_{p,r=R_f}}$$

(15)

However, we postulated that the selective partitioning of permeating particles across the membrane vary in the presence of fouling. To account for this, we developed an equation for the membrane partition coefficient ($K_{mem}$) that takes into consideration the fouling effects. The details of this fouling-dependent equation will be elaborated in the subsequent sections.

2.2. Resistance parameters

In the constant permeate flux setup of cross-flow filtrations of macromolecules, fouling often results in an increase in transmembrane pressure (TMP). This increment of TMP is necessary to overcome the growing resistances caused by fouling. However, the fouling of graphite led to a different TMP behavior. By the visual inspection of the TMP profile over time for the filtration of graphene from the non-exfoliated graphite (Figure S2.d of reference [20]), we detected three distinct regimes: an initial rise in TMP, a slight decline in TMP, and a stabilization of TMP. We hypothesize the initial abrupt increment in the TMP corresponds to the formation of the polarization layer, and the final stage corresponds to pore blockage and cake formation, which are related to the TMP through Darcy relation (Equation (16)) in our formulation.

$$U_{perm}=\frac{\mathrm{TMP}}{\mu R_{total}}$$

(16)

To characterize these two regimes, we adopted the the formulations developed by Marcos et al. [24] to represent the resistances associated with polarization ($R_{pol}$), and cake formation ($R_{cake}$).

$${\tau }_{pol}\frac{d{R}_{pol}}{dt}=(R_{pol,ss}-R_{pol})$$

(17)

$$R_{pol}(t=0)=0$$

(18)

The parameter $R_{pol,ss}$ corresponds to the steady state value of polarization resistance, indicating the onset of the pore blockage and cake formation regime. The ${\tau }_{pol}$ is the polarization time constant and relates to the concentration of the foulants on the inner wall of the membrane through the following equation.

$${\tau }_{pol}=\frac{k_{pol}}{C_{f,r=R_f}}$$

(19)

For the sake of more clarity, the $R^{*}$${}_{block}$ term in the Marcos et al. [24] report, which is the combined effects of pore blockage and cake formation resistances, is denoted as $R_{cake}$ in this work. The $R_{shear}$ accounts for shear stress in the flow, which helps wash away foulants and decrease the fouling resistance, and $R_{ads}$ corresponds to the initial adsorption of the foulants to the inner wall membrane and blockage of the pores.

$${\tau }_{cake}\frac{d{R}_{cake}}{dt}=(R_{cake}-R_{shear})$$

(20)

$$R_{cake}(t=0)=R_{ads}$$

(21)

However, the observed slight reduction in TMP during the filtration of 2D graphenes from the ethanol dispersion [20] is unique to the fouling of graphitic layered structures and has not been extensively studied before. To gain insights into this behavior, we turned to studies on graphene-coated membranes [33,34,35,36,37,38,39] , which have yielded enhanced filtraton performance. For example, the coating of the polyamide [37] and polysulfone [38] membranes with graphene oxide (GO) and modification of the polyvinylidene fluoride membrane [39] with sulfonated graphene oxide (SGO) has been reported to increase the flux rate of the permeating water. The higher flux rate of water for the coated membrane cases has been attributed to the greater hydrophilicity of GO and SGO compared to the polymeric surfaces, which leads to better spreading of the water droplets on the membrane surface, and accordingly higher permeation rate of water molecules through the membrane pores. Moreover, it has also been reported that ethanol-containing droplets spread on the surface of pristine graphene flakes [40], and the wettability of the surface increases (contact angle decreases) with higher ethanol content [41] . Therefore, we speculate when the graphite compounds foul on the inner wall of the membrane, smoother and more wettable surfaces are provided for the ethanol flow. This leads to a decrease in the contact angle between the ethanol droplets and the membrane surface, and they spread on the surface, promoting the rate of permeation of ethanol molecules through the pores. Therefore, to maintain the desired permeate flux in the constant permeate flux setup, the driving force (TMP) decreases slightly until the cake formation regime becomes dominant, at which point the TMP levels off. Thereby, to capture the effects of increased wettability of membrane surface upon fouling on the TMP behavior, we defined another resistance term, as below:

$${\tau }_{wet}\frac{d{R}_{wet}^{\prime }}{dt}=(R_{wet,ss}^{\prime }-R_{wet}^{\prime })$$

(22)

As this resistance parameter operates in an opposite direction of the previously defined resistances, that is it explains the reduction of the TMP, we used the prime symbol (’) to differentiate it from the other terms. The $R_{wet,ss}^{\prime }$ is the limit of fouling effect on the enhanced wettability before the cake formation regime takes over the control of the TMP behavior. Also, as the wettability resistance term depends on the concentration of the foulents on the membrane, we used a foulant’s concentration-dependent formulation for the time constant of the wettability effects.

$${\tau }_{wet}=\frac{k_{wet}}{C_{f,r=R_f}}$$

(23)

The total resistances used in the Darcy relation (16) to calculate the TMP is expressed in Equation (24). The $R_m$ term represents the clean membrane resistance, which is intrinsic to the membrane’s characteristics. Since the $R_{wet}^{\prime }$ term reduces the magnitude of the resistances in our formulations, it is associated with a negative sign in Equation (24).

$$R_{total}=R_m+R_{pol}+R_{block}-R_{wet}^{\prime }$$

(24)

This equation represents the local resistance of the membrane along the z-direction, which gives rise to local TMP values according to Equation (16). However, the experimentally obtained TMP values are determined by the following equation.

$$\mathrm{TMP}=\frac{P_{in}+P_{ret}}{2}-P_{perm}$$

(25)

Thus, we averaged the TMP values calculated using Equation (16) to (24) along the z-direction to compare them with the experimental data.

2.3. Permeate flux concentration

To determine the permeate concentration of the permeating solute (exfoliated graphene), we utilized a partition coefficient (Equation (15)). This coefficient depends on various factors, including membrane properties, solute characteristics, and thermodynamic conditions. However, in our case, the inner surface of the membrane undergoes changes during operation and fouling. As a result, we developed a fouling-dependent expression for the partition coefficient (Equation (26)), which satisfies three conditions:

a) Accounting for the size of the permeating particles. b) Maintaining a constant partition coefficient for very small solute particles molecules, such as ions. c) Reflecting a significant decline in the partition coefficient for large particles that are completely rejected with partial pore blockage.

$$K_{mem}\left(t\right)=\frac{R_{pol}^0+R_{cake}^0}{R_{pol}^0+\alpha \left(R_{pol}-R_{pol}^0\right)+R_{cake}^0+\beta \left(R_{cake}-R_{cake}^0\right)}$$

(26)

The formulation introduces two parameters, $\alpha$ and $\beta$, which represent the ratios of rejected permeating particles to rejected ethanol molecules in the polarization and cake formation stages, respectively. When dealing with small particles like ions where $\alpha$ = $\beta$ = 0 (Figure 2.a), the partition coefficient remains uniform, resembling the bulk concentration. Conversely, for larger molecules with $\alpha$ , $\beta$$>>$ 0 (Figure 2.b), as the solvent flows towards the membrane, the particles may be rejected by the foulants during the concentration polarization or cake formation regimes, leading to a decrease in the partition coefficient. Therefore, the fouling-dependent partition coefficient formulation (Equation (26)) accounts for changes in the membrane’s inner surface, differentiates between solute sizes, and ensures the accurate calculation of the permeate concentration.

3. Results

We carried out the optimization of model parameters through two consecutive steps. First, the resistance parameters were determined by fitting the calculated TMP values (Equation (16)) to the experimental data. Then, the $\alpha$ and $\beta$ parameters were obtained by fitting the calculated permeate concentration against experimental results. The experimental permeate concentration and TMP values were inherited from Figure S2.c and S2.d of the reference [20].

3.1. TMP dynamics

Table 1 summarizes the resistance parameters used to reproduce experimental TMP data. The $R_m$ value was estimated based on the previously measured hydraulic resistances for the microfiltration hollow fiber membranes (with pore sizes greater than 0.1 $\mu$m) [42,43,44]. These hydraulic resistance values consistently fall around the scale of ${10}^{11}$$\frac{1}{\mathrm{m}}$. The $R_{cake}^0$ value was deduced from the TMP value measured at the initial time step. Furthermore, we designated a value of 2 × 10${}^{-3}$ Pa·s for the feed viscosity to position it between the viscosity of pure ethanol (1.1 × 10${}^{-3}$ Pa·s) [45] and the measured viscosity of the thickened permeate stream (approximately 2.5 × 10${}^{-3}$ Pa·s) [20]. To ensure the consistency of the experimental and simulation data, we excluded the initial 3 minutes of the experiment and set t = 3 min as the initial time step. By optimizing the remaining resistance parameters, we successfully captured the three stages of fouling. The optimized parameters yielded a TMP profile that closely matched the experimental results (Figure 3), with a Mean Absolute Error (MAE) of 0.007 psi (Figure 3).

The TMP variation with time reflects the fouling of the foulant particles (non-exfoliated graphites) near the membrane’s inner wall as depicted in Figure 4. The higher concentration of foulants in the proximity of the inner wall ($\frac{r}{R_f}$≃ 1 ) indicates the reduction of permeating area for the permeating particles (exfolaited graphene flakes) and, thereby, reduced membrane partition coefficient for them. The effects of different stages of fouling on the filtration of the permeating particles are coupled to $\alpha$ and $\beta$ parameters, which will be elaborated in the following section.

3.2. Permeating particle dynamics

Analyzing Figure 3 reveals an initial abrupt increment of the TMP value until it reaches its maximum at t = 500 s, indicating the concentration polarization regime. Subsequently, the fouling particles reduce the pore size and form a cake layer, which begins to dominate the TMP behavior at the end of the membrane operation (t > 6000 s). The effects of the polarization stage on the rejection of permeating solute are encapsulated in the $\alpha$ parameter, while the $\beta$ parameter represents the influence of the cake layer growth on the disruption of the solute isolation. To investigate our model’s capability to capture the effects of $\alpha$ and $\beta$ parameters on the behavior of the permeating solute at the end of the polarization regime, we examined the concentration gradient in the neighborhood of the inner wall of the membrane, i.e., $\frac{r}{R_f}$≃ 1 (Figure 5), at t = 500 s. As Figure 5.a displays, for the very small solute particles ($\alpha$ = $\beta$ = 0), the uniformity of the inlet concentration (0.075 g/L) is preserved across the domain. For the particles with larger $\alpha$ and $\beta$ values, i.e., $\alpha$ = 10, $\beta$ = 10 (Figure 5.b), partial rejection of the particles by the foulant compounds occurs, leading to the build-up of the permeating solute particles and elevation of the average concentration value to 0.21 g/L on the membrane wall. Further increment of the $\alpha$ value ($\alpha$ = 100, $\beta$ = 10) leads to broader concentration gradient and promotes the average concentration on the wall to 0.32 g/L (Figure 5.c). However, we detected minimal impact by increasing the $\beta$ value, i.e., $\alpha$ = 10, $\beta$ = 100, (comparing Figure 5.d and Figure 5.b) on the behavior of permeating solute rejection. Hence, during the concentration polarization stage, solute filtration is primarily controlled by the $\alpha$ parameter, with negligible influence from $\beta$ values.

The effects of variation in $\alpha$ and $\beta$ parameters on the rejection behavior of permeating solute after more than 2 hours, i.e., t = 7500 s, where the cake formation contribution dominates, are displayed in Figure 6. For the small particles ($\alpha$ = $\beta$ = 0), the solute distribution remains uniform until the end of the simulation (Figure 6.a). With $\alpha$ = 10, and $\beta$ = 10, there is an insignificant increment of the wall concentration of permeating solute from 0.21 g/L at t = 500 s (Figure 5.a) to 0.23 g/L at t = 7500 s (Figure 6.b). Increasing the $\alpha$ value to 100 yields identical concentration gradient profile of permeating solute near the wall at both t = 500 s (Figure 5.c), and t = 7500 s (Figure 6.c). The slow dynamics of the permeating solute behavior near the wall can be attributed to the fact that for the large $\alpha$ values, the membrane partition coefficient converges to its lowest values at the end of the concentration polarization regime, maximizing the solute rejection. Therefore, a stable concentration boundary layer of the permeating solute is formed after the initial stage. Finally, in contrast to minimal effects of $\beta$ values on the average value of permeating solute concentration on the membrane wall at t = 500 s (comparing Figure 5.d and Figure 5.b), raising the $\beta$ value to 100 promoted the wall concentration of permeating solute to 0.27 g/L (Figure 6.d) compared to the 0.21 g/L obtained for the $\beta$ = 10 case at the end of the simulation (Figure 6.b). Therefore, our results suggest that the $\alpha$ parameter determines the level of permeating solute rejection during the concentration polarization regime, while the $\beta$ values control the rejection during the cake formation stage.

To further understand the effects of the $\alpha$ and $\beta$ parameters on membrane characteristics, we examined the variation of the relative membrane partition coefficient ($\frac{K_{perm}}{K_{perm}^0}$) over time (Figure 7.a), where $K_{perm}^0$ is the initial partition coefficient. As mentioned earlier, the particles with $\alpha$ = $\beta$ = 0 impose no impacts on the membrane characteristics, resulting in a constant partition coefficient (red plot). For the case with the highest $\alpha$ value ($\alpha$ = 100), the non-dimensional partition coefficient mitigates to 0.024 after 500 s and remains nearly constant (green plot). The two simulations for the $\alpha$ = 10 case initially exhibited identical partition coefficient behaviors (until t = 500 s), but diverged over time, with a slight decline in the case with $\beta$ = 10 (blue plot) and a steeper descent in the case with $\beta$ = 100 (purple plot).

The permeating solute concentration profile at the permeate stream is driven from the membrane partition coefficient and the wall concentration of the permeating solute (on the lumen side) using Equation (15). Figure 7.b demonstrates that the non-dimensional permeate concentration ($\frac{C_{p,perm}}{C_{p,perm}^0}$) profile follows a similar pattern to the membrane partition coefficient. However, except for the $\alpha$ = $\beta$ = 0 case, where the non-dimensional permeate concentration remains unity throughout the operation time, the permeate concentration is shifted upward compared to the membrane partition coefficient. This difference can stem from two factors. First, as the membrane partition coefficient decreases, a denser build-up of the solute particles is formed near the wall, leading to an increase in $C_{p,r=R_f}$ in the denominator of Equation (15), providing more solute to permeate through the remaining pores. In addition, the permeating solute concentration on the inner wall ($C_{p,r=R_f}$) is not uniformly distributed through the z-dimension of the membrane, and therefore, the permeate concentration at the permeate side does not exactly resemble the membrane partition coefficient behavior.

Ultimately, to assess the predictive capability of our model for the permeate concentration profile, we conducted the optimization of $\alpha$ and $\beta$ parameters to fit the calculated permeate concentration profile to the experimental data obtained for three different inlet concentrations. The experimental data were drawn from the Figure S2 of reference [20] , where C${}_{feed}$ denotes the combined graphene and non-exfoliated graphite content of the feed. In our simulations, we assumed the graphene compounds and non-exfoliated graphites accounted for 25 wt% and 75 wt% of the C${}_{feed}$, respectively. As displayed in Figure 8, the simulation yielded a satisfactory correlation with the simulations with a MAE of 0.0002, 0.0003, and 0.0022 g/L for C${}_{feed}$ values of 0.05, 0.1, and 0.3 g/L, respectively. Additionally, we compared the total amount of filtered permeating solute predicted by our model to the experimental data. The comparison revealed strong agreement, with absolute error values of 0.0003, 0.0005, and 0.0003 g/L for C${}_{feed}$ values of 0.05, 0.1, and 0.3 g/L, respectively.

As elaborated in the Methods section, the $\alpha$ and $\beta$ parameters represent the relative number of rejected solute particles compared to the rejected solvent molecules. Therefore, these parameters are inherently size and concentration-dependent. We speculate both the $\alpha$ and $\beta$ values ascend with increasing concentration. However, the optimized parameters for different concentrations failed to exhibit a meaningful correlation (Table 2). This lack of pattern may be attributed to our approach to obtain the $\alpha$ and $\beta$ parameters. As Equation (15) suggests, the partition coefficient is impacted by the interplay between the resistance parameters and the $\alpha$ and $\beta$ parameters. For the three different feed concentration studied systems, we adopted similar values for the resistance parameters, such as R${}_{pol,ss}$, k${}_{pol}$, R${}_{ads}$, etc., while these fouling-dependent parameters might vary upon change of the foulant content of the system. Therefore, although our model demonstrated the capability to accurately predict the permeate concentration profile of 2D materials, in order to determine concentration-dependent $\alpha$ and $\beta$ parameters, it is essential to target multiple experimental TMP and permeate concentration profiles over time. These profiles will serve as the basis for deriving the relationship between solute rejection behavior and varying concentrations. This approach can be explored in future studies to enhance our understanding of the concentration-dependent effects on fouling parameters.

4. Conclusion

In this study, we presented a CFD model for the filtration of 2D materials using hollow fiber membrane filtration. Navier-stokes and mass conservation equations have been integrated into the model to numerically solve the flow and concentration behaviors in the system. The fouling phenomena were coupled to the model through the Darcy relation and utilization of the resistance-in-series technique. The resistance contributions were assumed to arise from concentration polarization and combined effects of pore reduction and cake formation on the inner surface of the membrane. Additionally, we proposed that fouling of graphitic layered structures enhances the membrane surface wettability, leading to an increment in solvent permeation through the pores. The wettability effect was incorporated into the model by defining the $R_{wet}^{\prime }$ term and combined with the other resistance terms used in the Darcy relation. The model yielded an excellent fit with the experimental transmembrane pressure (TMP) data with the Mean Absolute Error of 0.007 psi.

Furthuremore, the model was also armed with a fouling-dependent membrane partition coefficient for the permeating solute, which formulated the effects of fouling regimes on the rejection of the solute by the definition of $\alpha$ and $\beta$ parameters. Our results revealed that the $\alpha$ parameter determines the level of the rejection during the concentration polarization regime, while the $\beta$ parameter controls the rejection of the permeating solute during the cake formation stage. The optimized $\alpha$ and $\beta$ parameters enabled the model to predict the experimental permeate concentration profile with Mean Absolute Error of 0.0002, 0.0003, and 0.0022 g/L for three different feed concentrations ($C_{feed}$) of 0.05, 0.1, and 0.3 g/L, respectively. The current model initiates the modeling of isolation of 2D materials using hollow fiber membranes, paving the way for further advancements in the eco-friendly production of 2D materials through filtration technique.