Computational modeling of drug separation from aqueous solutions using octanol organic solution in membranes

Continuous membrane separation of pharmaceuticals from an aqueous feed was studied theoretically by development of high-performance mechanistic model. The model was developed based on mass and momentum transfer to predict separation and removal of ibuprofen (IP) and its metabolite compound, i.e. 4-isobutylacetophenone (4-IBAP) from aqueous solution. The modeling study was carried out for a membrane contactor considering mass transport of solute from feed to organic solvent (octanol solution). The solute experiences different mass transfer resistances during the removal in membrane system which were all taken into account in the modeling. The model’s equations were solved using computational fluid dynamic technique, and the simulations were carried out to understand the effect of process parameters, flow pattern, and membrane properties on the removal of both solutes. The simulation results indicated that IP and 4-IBAP can be effectively removed from aqueous feed by adjusting the process parameters and flow pattern. More removal was obtained when the feed flows in the shell side of membrane system due to improving mass transfer. Also, feed flow rate was indicated to be the most affecting process parameter, and the highest solute removal was obtained at the lowest feed flow rate.


Scientific Reports
| (2020) 10:19133 | https://doi.org/10.1038/s41598-020-76189-w www.nature.com/scientificreports/ at industrial scales. The most widely used treatment method is adsorption, which is an energy-consuming process in order to regenerate the applied sorbent and can cause processing problems and pollution. On the other hand, solid-phase extraction approach allows low limits of detection, while consumes considerable time and has been found to be insufficient for application in sewage treatment plants 12 . In this sense, liquid membranebased separation strategies can be the appropriate substitute for the traditional techniques for the elimination of pharmaceutical pollutants from water. In current years, an innovative technology has been reported to separate trace concentrations of pharmaceutical molecules from water and wastewater streams. This membrane-based technology, e.g. hollow-fiber membrane systems, possesses tremendous privileges such as high analyte capacity, low consumption of organic solvent, easy handling, low analysis cost, membrane high surface area, and independently controllable flow rates [13][14][15][16][17][18][19] . Williams et al. used a hollow-fiber membrane contactor with immobilized canola oil for extraction of 4-IBAP from contaminated water 20 . With simultaneous regeneration of the immobilized solvent, 4-IBAP was extracted from the aqueous phase passing through the shell side of membrane contactor. This phenomenon made this system a green process through which 80% of 4-IBAP was removed by utilizing only 70-100 mL of canola oil. In order to extract three pharmaceutical compounds (salicylic acid, ibuprofen, and diclofenac) from the wastewater, a polypropylene hollow-fiber membrane was used by Payán et al. 21 . Detection limits were determined 20, 100, and 300 ng lit −1 for salicylic acid, diclofenac, and ibuprofen, respectively and this process had a very low organic solvent consumption.
In addition to experimental investigation of pharmaceutical contaminated water treatment by membranebased approaches, a mathematical model which can assess the effects of various processing parameters on the removal efficiency is required for future applications of these separation procedures. This research study aims to understand the removal of ibuprofen and its metabolite compound (4-IBAP) from an aqueous feed by contacting with an organic solvent (octanol) in a hollow-fiber membrane contactor system. In order to achieve this objective, a mechanistic modeling along with computational fluid dynamics (CFD) simulation is employed. Besides, the influence of different parameters and the flow rate of aqueous solution on the removal efficiency of IP and 4-IBAP will be evaluated.

Development of mathematical model
Here, the molecular removal of IP and 4-IBAP inside a membrane contactor is investigated by development of mechanistic model and CFD simulations. Figure 1 depicts the ibuprofen and 4-IBAP mass transfer and the geometrical structure of the membrane system. As seen, the tubular membrane contactor module contains many hollow fibers, however only a single fiber is considered for the modeling and numerical simulations. The aqueous feed containing pharmaceuticals flows inside the hollow fiber, whereas the octanol solution flows in the shell side, and in the opposite direction (counter current mode). However, it should be pointed out that the change of configuration (i.e. feed in the shell) will be studied using the model to evaluate the effect of flow configuration on the mass transfer efficiency. Owing to the hydrophobic nature of porous fiber, the organic phase (octanol) can penetrate into the membrane pores and reaches the other side (feed). Furthermore, as Williams et al. 8 have reported in an experimental study, by keeping a pressure difference which is higher on the aqueous feed, the organic phase was cannot penetrate the aqueous phase and mix with the feed, as the operation must be run as non-dispersive process. As such, the extraction takes place at the interface between the organic and aqueous phases (r 1 in Fig. 1, right side).  www.nature.com/scientificreports/ In this research, for estimating the effective hypothetical radius of shell encircling each membrane inside the membrane contactor, Happel's free surface model is applied [22][23][24] . Several simplifying assumptions have been made to develop the mechanistic model as follows: 1. Isothermal and steady state condition; 2. An axisymmetrical geometry of fiber is assumed; 3. Fluid velocity profile inside the tube side is in the fully developed condition; 4. The pores of hollow fiber are only filled with the organic phase; 5. Organic and aqueous solutions flow under a laminar regime.
The governing equation for describing IP and 4-IBAP mass transport from feed to the organic phase (octanol) is the continuity equation. This equation can be obtained as below 25,26 : where the concentration of solute, velocity vector, diffusive flux for component i, and term of reaction are respectively denoted as C i , V, J i , and R i . It should be noted that the term R i , is omitted in the model's equations because there is no chemical reaction inside the system. In order to estimate the diffusive fluxes for component i, Fick's diffusion law has been applied. Moreover, the velocity vector can be determined analytically or through combination of momentum and continuity equations. The convective and diffusional mass transports are respectively denoted by (∇C i V ) and (∇J i ) , which have been obtained as below 25,27 : In the above-mentioned equation (Eq. 3), mass flux vector is defined by N i . The value of diffusivity can be calculated from Stokes-Einstein equation by using the solute-solvent interaction 28,29 . Based on Happel's free surface model, and assuming laminar flow regime in both shell and tube sides, velocity distribution can be derived from Navier-Stokes equation as below 25,30 : In Eq. (4), V z , η, ρ , p, and F represent the velocity vector in the z direction, dynamic viscosity, fluid density, pressure, and body force term, respectively. The effective hypothetical shell radius surrounding each porous fiber (r 3 ) can be computed by Eq. where, the number of fibers and the module radius are described by n and R , respectively. The amount of r 3 is found to be 3.158 × 10 -4 m by combination of equations of (5) and (6).
The boundary conditions for solution of equations are listed below: Inside fiber: @inlet, z = 0: C i = C 0 @outlet, z = L: convective flux @r = 0, axial symmetry. @r = r 1 , Membrane: Scientific Reports | (2020) 10:19133 | https://doi.org/10.1038/s41598-020-76189-w www.nature.com/scientificreports/ Shell: @z = L: C s = 0 @z = 0: convective flux @r = r 2 , C s = C m @r = r 3 , Insulation/symmetry where C s , C m , and C i refer to the solute concentration in the shell, membrane, and tube side, respectively. Also, m refers to the partition coefficient which is defined as the ratio of concentration of solute in the organic phase to that of aqueous phase (m = C org /C Aq ).

Numerical simulations
The partition coefficient values of IP and 4-IBAP represented by m, and their values are summarized in Table 1.
For solving the Navier-Stokes and continuity equations and their relevant boundary conditions in the membrane system, COMSOL package version 5.2 has been employed. The package operates based on finite element method (FEM), which is a numerical technique used for solving partial differential equations (PDEs) [35][36][37][38] . Using FEM to solve the PDEs possesses significant advantages such as easy performance of non-uniform meshes and easy handling the complicated geometries [39][40][41] . A 64-bit operating system with an Intel core™ i5-6200U CPU at 2.30 GHz and an 8 Gigabyte RAM was utilized for the simulations. For solution of the model's equations, the required time has been computed to be about 1 min. Table 1 presents the physicochemical parameters and specifications of employed membrane system for developing the model and CFD-based simulations. Mapped meshing technique is aimed to be executed in this study in order to discretize the shell, tube and membrane domains of the contactor for assessing the change of important functional/design parameters at each domain point. A significant enhancement in the computational accuracy can be observed when the number of mapped meshes inside each domain of contactor increases, and thereby the calculations' error is reduced. The meshes in three domains have been shown separately in Fig. 2.
To assess the convergence of numerical scheme used for the simulations, the convergence status of the solution as function of iteration was plotted in Fig. 3. It is seen that the system reaches the acceptable error after 11 iterations, and the solution stops at this point. This can confirm the selection of efficient solver for this set of equations developed in this study. Mesh independence test was carried out to find the optimum number of grids for the numerical solution, and also to ensure that the solution is independent of number of meshes. Indeed, increasing the number of grids in the CFD simulations would decrease computational errors and, consequently enhance the simulations accuracy. However, it would significantly increase the computational expenses. Therefore, the optimum number of meshes need to be precisely determined to avoid any computational expenses and instability. Figure 4 represents the impact of mesh numbers on the IP/4-IBAP outlet concentration. It is seen that after the 212th mesh, no meaningful variations in the outlet concentration of IP/4-IBAP concentration occurs, which corroborates the convergence of the computational simulation results. Therefore, 212 can be selected as the optimum number of meshes for the simulations, and any simulations with grids more than 212 is independent of the mesh.

Results and discussion
Concentration distribution of pharmaceuticals. Figures 5 and 6 illustrate concentration distribution of IP and 4-IBAP in three distinct sections of the membrane contactor which have been obtained by solving continuity equation. It should be pointed out that the simulations have been conducted for two cases, i.e. when   Figure 7 shows the concentration profile when the feed flows in the tube side, and Fig. 8 shows the concentration profile when the feed passes through the shell side. As seen, the outlet solute concentration is lower for the case when feed flows in the shell side, indicating higher removal. This could be attributed to  www.nature.com/scientificreports/ the higher surface area for mass transfer when feed flows in the shell side. Moreover, it is seen that the removal of 4-IBAP is higher than IP which is due to higher partition coefficient of 4-IBAP in octanol.

Effect of aqueous flow rate on separation.
Influence of feed flow rate on the solutes removal is shown in Fig. 9 for both cases when feed flows in the tube and shell. The percentage of solute separation is determined as 42 : As seen in Fig. 9, the outlet concnetration of solutes is increased with increasing the feed flow rate in the membrane contactor. Indeed, separation percentages decline as the feed flow rate rises up, which is due to reduction of feed residence time in the contactor. As such, short residence time can decrease the mass transfer efficiency significantly, and in order to obtain high removal efficiency, the feed flow rate should be kept at the lowest possible value. When the residence time is long, the feed solution has more time to allow the solutes transfer to the organic phase which in turn increases the separation percentage.   tion of solutes. It should be pointed out that in membrane contactor systems, fibers do not provide selectivity for particular components, and the role of fiber is as physical barrier for contacting two phases. However, the membrane porosity can alter the mass transfer mechanism and resistance. Using the developed model in this study, the solutes outlet concentrations at different fiber porosities were evaluated, and the results are shown in Fig. 10. As seen, the simulations for both solutes and cases have been provided. It is indicated that the porosity has positive influence on the separation efficiency of pharmaceutical components, and higher separation has been achieved at high fiber porosity values. This behavior can be attributed to the facilitation of solute mass transport by increasing porosity which decreases the mass transfer resistance. In fact, at higher fiber porosity, the collision and interaction between the solute molecules and the fiber wall would be minimized, leading to higher mass transfer rate through the pores of fiber.   Fig. 11. It is shown that the outlet concentration of both solutes is decreased with increasing the number of fibers, which is favorable in mass transfer and removal of solute. In fact, variations in the number of fibers alters the mass transfer area between the feed and solvent, thereby higher mass transfer occurs with increasing the number of fibers 43,44 . It is clearly seen that for both cases, the effect of increasing the number of fibers on mass transfer is positive and can be considered as an important parameter in design and optimization of membrane-based separation of pharmaceuticals.
Effect of the aquoues/organic flows' arrangmement on IP/4-IBAP separation. Solutes removal efficiency as function of feed flow rate is illustrated in Fig. 12 for all operational cases. Comparisons have been made between the separation efficiency of ibuprofen and 4-IBAP when feed flows in the sell side and tube side.  www.nature.com/scientificreports/ As indicated, the highest removal percentage has been obtained for 4-IBAP when the feed flows in the shell side. The reason for higher removal when feed flows in the shell side could be due to higher mass transfer area, and mass transfer coefficient as well. Indeed, as the membrane is hydrophobic it woud be filled with the organic phase and the interface between two phases would be formed at the fiber outer surface, when feed flows in the shell side. This will provide higher mass transfer area for removal of solutes. Furthermore, it can be seen that aqueous flow rate has major effect on the removal of solutes using the membrane systems.

Conclusions
In this research, computational modeling drug separation in membrane contactors were carried out using computational fluid dynamics method. The mass transfer equations along with Navier-Stokes equations were derived in tubular coordinate and solved to determine the outlet concentration of solutes. Ibuprofen (IP) and 4-isobutylacetophenone (4-IBAP) which is IP's metabolite compound were considered as solutes in the simulations. The simulations were performed for counter-current mode, but the flow pattern was changed to investigate its effect on the removal efficiency. It was revealed that higher efficiency was obtained when the aqueous feed flows in the  www.nature.com/scientificreports/ shell side of membrane contactor which is due to higher mass transfer area. Furthermore, the simulation results indicated that lowering feed flow rate, increasing fiber porosity and number, would have positive effect on the removal efficiency of both solutes. Also, the removal of 4-IBAP was higher than IP in the membrane contactor due to its higher affinity towards the organic solvent. License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.