Skim milk concentration is an important filtration process in the food industry in which milk proteins are separated from other components of skim milk, such as salt, lactose, microsorganisms, and water. This process is widely used in cheese making [1, 2] and production of milk protein concentrates, such as whey, which is a vital ingredient used in infants’ food and also valuable for preparing barley sourdough [3, 4, 5].
Many techniques can be used for skim milk concentration, such as spray drying, freeze drying, evaporation, salt, pH and heat treatment. However, for large-scale milk concentration, these methods are neither economical nor suitable because they degrade the nutritional properties of the final product and their overall productivity is also very low. Deterioration of the natural state of the milk due to these traditional techniques is associated with the exposure to heat treatment and changing the salt and pH content of the feed [6, 7].
Instead, ultrafiltration membrane contactors are widely applied in milk concentration and dehydration processes. These filtration systems have attracted great attention over the past few decades not only for their improved capacity and efficiency but also for their considerable decrease in power consumption [8-12]. In fact, ultrafiltration processes have several advantages over traditional concentration techniques such as evaporation. Firstly, since the separation is generally based on the size differences between milk components and membrane pores, much more control over the final composition of the product is achieved. Secondly, energy consumption of the separation process decreases significantly since there is no need to heat the feed. Last but not least the final quality and nutritional properties of the product is better preserved because the long time heat exposure is eliminated [2, 13].
Unfortunately, the performance of ultrafiltration membrane processes is severely limited by a negative phenomenon called fouling, which is defined as the deposition of solutes on the surface of the membrane or the inner wall of the membrane’s pores . This negative phenomenon results in the permeate flux decline over a longer period of time and reduces the efficiency of the filtration process. Major fouling materials in milk solutions are proteins and fats. The size of these molecules is larger than the pores of the membrane, so they are rejected by the membrane and consequently accumulate on the surface of the membrane in the form of a cake layer. Removal of this layer is of vital importance to maintain the performance of the membrane at a desirable level. Therefore, accurate models are needed to predict the fouling at various operational conditions, in order to consider and plan suitable cleaning processes for removing that.
Fouling is heavily dependent on the size of the feed solutes, intrinsic properties of the membrane itself and operational conditions, such as trans-membrane pressure (TMP), temperature and velocity of the feed flow [15, 16]. Thus, working under the optimum operating conditions and also the treatment of the membrane make it possible to remarkably slow the rate of the fouling down. For instance, Ahmad et al. used the sol-gel technique to prepare improved antifouling membranes by placing a mesoporous TiO2 layer on a porous alumina support . They used the treated membrane to analyze the fouling behavior of the model dye compounds in membrane reactors and found that the fouling rate dropped sharply after membrane treatment by this technique.
So far, extensive studies have been conducted to model fouling in membrane contactors using numerical or mathematical methods. These methods either use macroscopic or microscopic models in their simulations [17, 19-22]. Microscopic models focus on solving Navier-Stokes equations along with appropriate boundary and initial conditions. Based on this technique, Reza Kazemi et al. proposed a model for concentration distribution of water in ammonia/water solution in a hollow-fiber membrane system. They used COMSOL Multiphysics software to solve Navier-Stokes equations of mass and momentum transfer numerically by Finite Element Method (FEM). Their reported results were in good agreement with experimental data .
In macroscopic models, however, a mathematical equation for flux decline is derived from major fouling mechanisms. These fouling mechanisms include complete blocking mechanism, intermediate blocking mechanism, cake formation mechanism and standard blocking mechanism . According to intermediate or complete pore blocking mechanisms, the available membrane area shrinks with filtered volume, as a consequence of the blockage of the membrane’s pores with the feed solutes. These two models are similar to each other, but the assumption behind the complete blocking model is stricter compared with the intermediate blocking model. According to the complete blocking model, the particles seal off the pores without superimposition upon each other. However, the intermediate blocking model alleviates complete blocking mechanism assumption by assuming that some of the solutes block pores and the others superimpose on one another. Standard blocking mechanism assumes that small particles deposit on the inner wall of the membrane’s pores, deteriorating the membrane structure and adding additional resistance to the feed flow. This fouling mechanism is considered the most severe fouling, which is very difficult or in many cases impossible to be removed by the usual cleaning processes. According to the cake formation mechanism, particles accumulate on the surface of the membrane in the form of a permeable cake layer, adding an excess resistance to flow . Previous studies on ultrafiltration separation processes showed that this fouling mechanism is the major mechanism responsible for ultrafiltration fouling in wastewater treatment, milk concentration and nitrogen removal [13, 14]. However, there are more than one individual fouling mechanisms involved in the flux decline. Therefore, great attention has been paid to the combination of these mechanisms and their simultaneous effects on flux loss. In this regard, Ho and Zydney  developed a mathematical model for flux decline during the purification of bovine serum albumin solution, which accounted for initial fouling by complete blocking and consequently by cake formation mechanism. Their proposed mathematical model, which explicitly described the inhomogeneity in the cake layer due to the complete blocking mechanism was the first model presented for the combined effects of fouling mechanisms and was in good harmony with the empirical data.
Following Ho and Zydney modeling work, a method was utilized to develop combined models of fouling with two fitted parameters . The method was to insert explicit equations of resistance and available membrane area as a function of time or filtered volume into Darcy’s equation. The equations then were integrated to obtain explicit equations of permeate volume as a function of time in constant pressure operation or pressure as a function of time in constant flow operation. The two fouling mechanisms were assumed to occur simultaneously. Based on this method, five models were developed that explained the combined effects of cake filtration-complete blocking, cake filtration-intermediate blocking, complete-standard blocking, intermediate-standard blocking, and cake filtration-standard blocking. The models then were compared with experimental data during microfiltration and ultrafiltration of bovine serum albumin and human IgG. The authors found that these models provided better predictions of flux decline compared to individual fouling mechanisms. Their findings also showed that the combined cake filtration-complete blocking model provided the best data fit for fouling of biological fluids.
Briao and Tavares also used macroscopic modeling in their analysis and found that the data for the fouling of dairy wastewater in tubular ultrafiltration membrane could be fitted initially by the complete blocking and subsequently by the cake formation mechanism . They also found that the cake filtration mechanism was the main mechanism responsible for fouling in spiral wound membranes.
Another fouling modeling was conducted by Corbaton-Baguena et al.  who fitted an exponential model to the data for the fouling of polyethylene glycol aqueous solution. By taking concentration polarization, particle accumulation on the surface of the membrane and long-term fouling into consideration, they found that this model gave accurate predictions of the fouling in harsh operational conditions.
As understood, cake filtration model and combined cake formation-complete blocking model are the major models in ultrafiltration membrane’s fouling. In this study, these two models were used to analyze the flux decline during skim milk concentration at various TMPs and temperatures. These models were compared with each other to provide a deeper insight into the fouling mechanisms involved in this process. Simulations were carried out by CurveExpert Professional software version 1.6 that is a powerful tool for fitting mathematical equations to the experimental data. Permeate volume equations proposed by Bolton et al. were differentiated to obtain flux equations as a function of separation time. Then, the equations were fitted to the empirical data using CurveExpert Professional software. Two different sets of skim milk ultrafiltration data were used here. The first set of data was obtained from testing under constant temperature and flow rate operations with reconstituted skim milk solution. Reconstituted skim milk was filtered through a spiral wound ultrafiltration membrane module. The effect of varying trans-membrane pressures on flux and total hydraulic resistance was then studied. The simulation results revealed that the combined cake filtration-complete blocking model was in excellent agreement with experimental data for skim milk concentration by spiral wound ultrafiltration membrane module, providing a better prediction of flux decline compared to cake filtration model. The second set of experimental data was obtained from partially skimmed milk concentration testing through a Pellicon cassette ultrafiltration membrane module. Unlike the first experiment, the combined cake filtration-complete blocking model was not capable of predicting flux decline and failed catastrophically by CurveExpert software (so there is no curve for this model in the final results) while the cake formation model was in excellent agreement with data sets.
The resistance data was also modeled using Bolton et al. proposed resistance equations of standard blocking model and cake filtration model. The results indicated that standard blocking model predictions were incompatible with experimental data and thus was not appropriate for predicting resistance. However, cake filtration model was in good accordance with the resistance empirical data and can be used to solutions whose solutes have a similar behavior to the milk components.
The flow rate can be calculated according to the Darcy’s law
where P is the trans-membrane pressure (Pa), μ is the viscosity of feed solution (Pa.s), A is the membrane area (m2) and R is the resistance (m-1).
Complete blocking model
In this model, it is assumed that membrane consists of parallel pores with constant radius and length and that each solid particle arriving at the membrane blocks a portion of the pores, without superimposition of particles on one another. The following equation shows the relationship between available membrane area and filtered volume 
Where Ccb, A0 and J0 denote fitted parameter constant for complete blocking model (s-1), initial membrane frontal area (m2) and initial flux (m/s), respectively. The equation can be used both in constant flow rate and constant trans-membrane pressure conditions.
At constant flow rate conditions, Eq. (2) can be written as a function of time 
By inserting Eq. (2) into Darcy’s law and integrating, the equation of permeate volume as a function of time can be obtained 
This equation can then be differentiated to obtain the equation for flux in terms of processing time
Cake filtration model
In the cake filtration model, it is assumed that particles accumulate on the surface of the membrane and superimpose upon each other in the form of a permeable cake. As the thickness of the cake increases with time so does the resistance to flow. Therefore, the resistance of the cake layer along with the membrane intrinsic resistance contributes to the total resistance. The total resistance increases with permeate volume according to Eq. (6) and with time according to Eq. (7) 
where Ccf denotes cake filtration constant parameter which has units of s.m-2. The filtered volume can be obtained from the following equation 
This equation can then be differentiated to obtain equation for flux as a function of time
Standard blocking model
In this model, it is assumed that membrane consists of straight cylindrical pores whose radius decreases with time due to the accumulation of solid particles on the pore walls of the membrane [18, 26]. The following equations show the relationship between resistance and volume (Eq. (10)) or time (Eq. (11)) 
In these equations, constant Csb denotes fitted parameter for standard blocking model with units of m-1. The equation of permeate volume as a function of time thus can be obtained 
Combined cake filtration-complete blocking model
Bolton et al. developed a combined model based on the effects of pore blockage and cake formation mechanisms. Filtration area loss predicted by the complete blocking mechanism was combined with the resistance from caking. By inserting Eq. (2) and Eq. (7) into Darcy’s equation, the equation for volume filtered as a function of time can be obtained 
The equation can then be differentiated to determine the equation for flux in terms of time
In case of caking, the equations of resistance (R) versus time, Eq. (7), was used. The equation of R as a function of V cannot be used here since volume filtered is defined relative to the available membrane area which is declining during the experiment.
In case of complete blocking, the equation of area in terms of volume filtered was used. The equation of area loss as a function of time is not valid here since the rate of complete blocking with respect to time is slower than cake formation. A detailed description of these models is provided by Bolton et al. .
A summary of these models is provided in Table 1.
Skim milk ultrafiltration using a spiral wound module
The experimental data sets used in this paper were obtained from previous studies of skim milk ultrafiltration done by Razavi et al. . They used a spiral wound membrane ultrafiltration unit (Biocon company, Moscow, Russia). The membrane material was polysulfone amide with a 20-kDa molecular weight cut off.
The membrane system was equipped with a tubular heat exchanger and a temperature sensor was applied to keep the feed at a constant temperature of 40 °C. Membrane unit was also installed to a feed tank (20 L) and a flow meter measured the feed flow rate which was held constant at 15 L.min-1 .
The membrane’s length, inner radius, and surface area were 0.47 m, 0.11 m and 0.33 m2, respectively.
Two pressure gauges measured inlet pressure (Pin) and outlet pressure (Pout). The TMP was obtained by the following equation 
The feed solution was prepared by mixing skim milk powder with water at a temperature of about 50 °C in a blender. The final pH of the feed was obtained to be 6.54. The average content of the total solid (TS) and water in produced solution was 8.443 and 91.557 %, respectively. Experiments were carried out at different TMPs (50, 100, 150, 200 and 250 kPa) to investigate the effect of pressure on flux decline. All experimental runs repeated twice. For each filtration run, the water was first contacted with the membrane to measure water flux. Membrane intrinsic resistance then was calculated according to Darcy’s law 
where μpermeate denotes permeate dynamic viscosity (water dynamic viscosity) (Pa.s) and Ri,m denotes membrane intrinsic resistance (m-1). In the next step, the water was replaced with reconstituted skim milk at 40 °C. The permeate flux was measured and recorded every 30 s.
After 30 min filtration, the operation was stopped and the membrane system was cleaned by distilled water and NaOH solution according to the protocol described by the manufacturer.
Skim milk ultrafiltration using a Pellicon cassette module
Rinaldoni et al. used a Pellicon cassette membrane ultrafiltration unit (Millipore, USA) . The membrane material was modified polysulfone with 10-kDa molecular weight cut off and total membrane surface area of 0.5 m2.
The feed solution (partially skimmed milk supplied by MILKAUT) was first heated in a water bath. It was then pumped through the membrane system with a constant feed flow rate of 29 ± 0.05 L.min-1.
The effect of different pressures (0.5, 1 and 1.5 bar) and different temperatures (20, 30 and 40 °C) on flux was investigated. All experiments were carried out twice.
Skim milk concentration with spiral wound ultrafiltration module
The combined cake filtration-complete blocking model was compared with individual cake filtration model and experimental data for the concentration of skim milk. Experiments were performed using a polysulfone amide ultrafiltration membrane with twelve kilograms of reconstituted skim milk at constant temperature and flow rate (40 °C and 15 L.min-1, respectively) . To investigate the effect of pressure on flux decline, experiments were carried out at 50, 100, 150, 200 and 250 kPa trans-membrane pressures. The permeate flux and total resistance were measured as a function of time and recorded every 30 s.
The flux data and the model predictions for skim milk samples at different pressures are shown in Fig. 1. The combined cake filtration-complete blocking model (solid green lines) provided the best fits of the flux data at moderate pressures (100 and 150 kPa) with error values remarkably lower than the individual cake filtration model (Fig. 1(b) and (c), Table 2). At high and low pressures, however, the error values of the combined cake filtration-complete blocking model were slightly higher than the individual cake formation model. This behavior can be attributed to the prevalent effect of cake formation mechanism at higher TMPs and negligible complete blocking mechanism effects. This could be accurate, since high TMPs will result in a more compact cake layer, with pore inlets much smaller to be blocked by particles, which tend to seal off membrane pores in a complete blocking manner.
Values of Ccb and Ccfwere applied to evaluate the contributions of complete blocking and cake formation to the combined cake formation-complete blocking model. The small value of which was approximately 0, indicated that the contributions of the two individual mechanisms were similar.
Since milk solution has distinct species with various size and shape, probably the fouling is caused by caking of large molecules such as proteins and fats arriving at the membrane earlier in the run. Fouling then will be followed by complete blocking of smaller particles like microorganisms. These particles can go through the cake, reach the membrane surface and seal off membrane pores. Thus complete blockage occurs at the end of the experiment. As it can be seen from Table 2, with increasing trans-membrane pressure the differences between experimental data and models generally increased. This behavior can be attributed to other mechanisms involved in fouling at higher pressures, like diffusion of small particles through the cake layer. Thus a fifth mechanism must be considered at higher pressures, which is expected to be investigated at further works.
The results of total resistance data modeling are shown in Fig. 2. Since the resistance data at 150 kPa and 100 kPa were very close to each other, only data at 100 kPa are shown and modeled. As it can be seen from the figure, standard blocking model was poor in predicting resistance to flow (dashed green lines), while cake formation resistance (solid red lines) was in agreeable agreement with the resistance data for skim milk at lower pressures with low error values, as seen in Table 3.
In ultrafiltration of skim milk, the fouling is begun with an accumulation of large particles like proteins and fats in the form of a cake layer on the surface of the membrane. Smaller particles then will settle down on top of this layer. Therefore, the cake filtration mechanism seems to be the main mechanism in this process which appropriately predicts fouling behavior. This was verified by the poor predictions of the standard blocking model and good predictions of total resistance due to the cake formation mechanism.
The calculated fitted parameters versus iteration using CurveExpert are shown in Fig. 3. As it can be seen in the figure, the number of iterations used for calculating fitted parameters was more than 40 and 15 for combined cake filtration-complete blocking model (Fig. 3(a)) and cake filtration model (Fig. 3(b) and (c)), respectively. This means that CurveExpert is converged after more than 40 and 15 iterations for these models. According to Fig. 3(d), after approximately 2 to 3 iterations, a fitted parameter for the standard blocking model is not changed. This means that CurveExpert is converged after approximately 2 to 3 iterations for standard blockage mechanism.
Fig. 4 and 5 show the effect of different trans-membrane pressures on flux and total resistance, respectively. Averaged values of flux and resistance were used and the data was modeled using Hoerl and Exponential Association regressions. For Hoerl regression (blue solid line), the following equation was used
In which Y denotes either flux or resistance, x denotes TMP and a, b and c are constant parameters. These parameters are obtained by regression and their values are shown in Table 4 and 5 for flux and resistance data, respectively. For Exponential Association regression, the following equation was used
As it can be seen from Fig. 4, flux increased by TMP. Both regressions were in excellent agreement with actual values. However, predictions of Hoerl regression were better for flux prediction with error values slightly lower than Exponential Association regression (Table 4). For resistance prediction (Fig. 5), two regressions were approximately the same.
Skim milk concentration with Pellicon cassette ultrafiltration module
Experiments were performed using a modified polysulfone ultrafiltration membrane with partially skimmed milk at constant flow rate (29 ± 0.05 L.min-1) . To investigate the effect of pressure on flux decline, experiments were carried out at 0.5, 1 and 1.5 bar trans-membrane pressures. The effect of different temperatures (20, 30 and 40 °C) on flux decline was also investigated.
The combined cake filtration-complete blocking model failed catastrophically by CurveExpert Professional, so there is no curve for this model in the final results. The flux data at different pressures and different temperatures thus was modeled using only cake filtration model. Applying cake formation mechanism resulted in a very good data fit. The flux data and the model predictions for skim milk samples at different pressures and different temperatures are shown in Fig. 6 and 7, respectively.
As it can be seen from Table 6, standard errors for the cake filtration model were almost the same at different trans-membrane pressures. This indicates that TMP had no significant effect on simulation of flux decline for ultrafiltration of partially skimmed milk.
The calculated cake filtration fitted parameter versus iteration using CurveExpert is shown in Fig. 8. As it can be seen in this figure, the number of iterations used for calculating fitted parameters in cake filtration model was more than 15. This means that, after more than 15 iterations, cake filtration fitted parameter has not changed and CurveExpert is converged.
The effect of different trans-membrane pressures and temperatures on flux decline is shown in Fig. 9 and 10, respectively. Averaged values of flux were used and the data of flux versus TMP and temperature was modeled using Modified Power and Geometric regressions. For Modified Power regression, the following equation was used
For Geometric regression, the following equation was used
in which x denotes either TMP or temperature. As it can be seen from the figures, flux grew by TMP and temperature. However, flux growth was higher with a 0.5 bar increase in TMP than a 10 °C increase in temperature. Geometric regression was better for predicting flux versus temperature, with error values slightly lower than Modified Power regression, as shown in Table 8. However, predictions of Modified Power regression were better for predicting flux vs. TMP, with error values considerably lower than Geometric regression error values (Table 7).
To model fouling in milk concentration process, the combined cake filtration-complete blocking model proposed by Bolton et al. and typical cake filtration model were applied to ultrafiltration of skim milk at different TMPs. The combined model used two fitted parameters and simplified to the equations for the typical mechanisms when the effects of second fouling mechanism were negligible.
The applicability of the model to the data for the skim milk concentration at different trans-membrane pressures was tested. The combined cake filtration-complete blocking model was in excellent accordance with experimental flux data. It was able to predict fouling of skim milk concentration through spiral wound ultrafiltration module and provided good fits compared to individual cake formation mechanism. Thus, the combined cake filtration-complete blocking model can be regarded as an effective model for predicting flux decline in solutions where the volume filtered decreases in a manner between the extremes of cake formation and complete blocking. The results also showed that cake formation model was in good harmony with experimental data for flux decline in ultrafiltration of partially skimmed milk while combined cake formation-complete blocking model failed catastrophically. The results of resistance modeling also showed that the cake formation model provided good fits of the data sets while the standard blocking predictions were poor.
CONFLICT OF INTEREST
The authors declare that there are no conflicts of interest regarding the publication of this manuscript.
A available membrane area (m2)
A0 initial membrane area (m2)
a regression constant parameter
b regression constant parameter
c regression constant parameter
Ccb complete blocking constant (s-1)
Ccf cake filtration constant (s.m-2)
J flux (m.s-1)
J0 initial flux (m.s-1)
Jpermeate permeate flux (m.s-1)
Pin inlet pressure (Pa)
Pout outlet pressure (Pa)
Ppermeate permeate pressure (Pa)
Q flow rate (m3.s-1)
R resistance to flow (m-1)
R0 initial resistance to flow (m-1)
Ri,m intrinsic membrane resistance (m-1)
t time (s)
TMP trans-membrane pressure (Pa)
V permeate volume (m3.m-2)
µ feed viscosity (kg.m-1.s-1)
µpermeate permeate dynamic viscosity (kg.m-1.s-1)
If you have any problem in viewing images please download the PDF version.