Ind. Eng. Chem. Res. 2007, 46, 5673-5685
5673
Application of Response Surface Methodology and Experimental Design in Direct Contact Membrane Distillation Mohamed Khayet,*,† Cornel Cojocaru,‡ and Carmen Garcı´a-Payo† Department of Applied Physics I, Faculty of Physics, UniVersity Complutense of Madrid, Spain, and Department of EnVironmental Engineering and Management, Faculty of Chemical Engineering, Technical UniVersity of Iasi, Romania
The experimental design and response surface methodology (RSM) is applied to a direct contact membrane distillation process. The factors considered for experimental design were the feed and permeate flow rate, the mean temperature, and the initial feed concentration of salt (NaCl) aqueous solution. The significant factors were optimized using a central composite design of orthogonal type. The quadratic models between the responses (permeate fluxes) and the independent parameters were built for both commercial and laboratory made membranes of different characteristics. The response surface models were tested with analysis of variance (ANOVA). For optimization purposes, the canonical analysis was employed. An algorithm was developed following the gradient method and using step adjustment in order to explore the response surface only inside the region of experimentation. The obtained optimal points were located in the valid region. The predicted permeate fluxes were compared with the experimental ones. In general, there is a good agreement between the experimental and the predicted permeate fluxes by RSM. 1. Introduction Membrane distillation (MD) is a thermally driven process mainly suited for applications in which water is the major component present in the feed solution to be treated. It is a thermally driven transport of vapor through nonwetted porous hydrophobic membranes. The driving force for the partial pressure difference between each side of the membrane pores. The liquid feed to be treated by MD must be maintained in direct contact with one side of the membrane without penetrating its dry pores unless a transmembrane pressure higher than the membrane liquid entry pressure (i.e., breakthrough pressure, LEP) is applied. In fact, the hydrophobic nature of the membrane prevents liquid solutions from entering its pores due to the surface tension forces. As a result, liquid/vapor interfaces are formed at the entrances of the membrane pores. Various MD modes differing in the technology applied to establish the driving force can be used. The differences between them are localized only in the permeate side. MD can be configured as direct contact membrane distillation (DCMD), air gap membrane distillation (AGMD), vacuum membrane distillation (VMD), or sweeping gas membrane distillation (SGMD).1-3 More details and extended explanations of the MD process can be found in the two reviews reported by Lawson and Lloyd2 and El-Bourawi et al.3 In the MD process, when aqueous solutions with nonvolatile components such as NaCl are considered, only water vapor flows through a membrane and the achieved solute rejection is very close to 100%. In fact, the potential applications of MD for production of high-purity water, concentration of ionic, colloid, or other nonvolatile aqueous solutions, and removal of trace volatile organic compounds (VOCs) from wastewater are well recognized.2 * To whom correspondence should be addressed. E-mail: khayetm@ fis.ucm.es. Tel.: +34(91)3944454. Fax: +34(91)3945191. † University Complutense of Madrid. ‡ Technical University of Iasi.
It must be pointed out here that MD offers various advantages in comparison to the traditional distillation and pressure-driven membrane processes.2,3 Nonetheless, from a commercial standpoint, MD is not implemented yet in industry and is still at the laboratory research and development stage. This is attributed to various reasons. Among them are the design of adequate membranes and modules for MD, the energy inefficiency, and the lack of optimization studies in MD systems. In general, the reported MD studies deal with conventional methods of experimentation in which one of the MD parameters is varied maintaining the others constant at different levels. Such conventional or classical methods of experimentation usually involve many experimental runs, which are timeconsuming, ignore interaction effects between the considered parameters of the process, and lead to a low efficiency in optimization issues. These limitations of the classical method can be avoided by applying the response surface methodology (RSM) that involves statistical design of experiments in which all factors are varied together over a set of experimental runs. In fact, RSM is a combination of mathematical and statistical techniques useful for developing, improving, and optimizing processes and can be used to evaluate the relative significance of several affecting factors even in the presence of complex interactions.4 Frequently, experimenters faced the same problem in many technical fields, where in general they needed to find the mathematical model in order to predict the response variable of interest as a function of a set of system variables. Frequently, the elaboration of an analytical mathematical model describing a physical and/or chemical process is a laborioustask owing to the complexity of the phenomena involved. Alternatively, an empirical model can be developed obtaining the response from the experiment for various input signals without involving the corresponding governing equations. Such a type of model, which can be developed by applying RSM, sets the correlation between input and output signals of a process experiment over the range of independent variables (factors). It must be stated here that one of the objectives of RSM is to
10.1021/ie070446p CCC: $37.00 © 2007 American Chemical Society Published on Web 07/24/2007
5674
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
determine the optimum operational conditions for systems or to determine a region that satisfies the operating specification. Identifying and fitting from experimental data an appropriate response surface model requires first some use of statistical design fundamentals, i.e., design of experiments (DoE), as well as regression modeling techniques and optimization methods. All these three topics are usually combined into RSM. RSM has been applied successfully in various scientific and technical fields such as applied chemistry and physics, biochemistry and biology, chemical engineering, environmental protection, and so on. In the adsorption field, Ravikumar et al.4 applied a 24 full factorial central composite for the optimization of batch process parameters for dye removal. Lacin et al.5 applied RSM for modeling of adsorption and ultrasonic desorption of cadmium(II) and zinc(II) on local bentonite by developing both linear and interaction empirical models based on the 26-1 fractional design for adsorption and 25-1 fractional factorial design for desorption. RSM has also been applied in biochemistry and biological fields.6,7 For example, RSM based on a five-variable central composite rotatable design was used for the synthesis of isobutyl isobutyrate by direct esterification of isobutyric acid and isobutyl alcohol.7 In chemical engineering, Zhao et al.8 determined which factors had significant effects on hydrogen cyanide (HCN) oxidation over a 0.5% Pt/Al2O3 catalyst using a fractional factorial design. Lemoine et al.9 used a central composite design (CCD) and analysis technique in order to study the effect of operating conditions on the hydrodynamic and mass transfer parameters for the toluene oxidation process. In the wastewater treatment field, Ahmadi et al.10 applied CCD and RSM to the advanced treatment of olive oil processing wastewater using Fenton’s peroxidation. Also, Ramirez et al.11 applied RSM for experimental design and optimization of the degradation process of the synthetic dye Orange II using Fenton’s reagent as catalyst. Under the optimum conditions, performances of 99.7 and 70.7% for color and total organic content (TOC) removal, respectively, were experimentally reached. Rana-Madaria et al.12 designed batch experiments by RSM using Box-Behnken design to optimize the key parameters for maximizing the removal percentage of chromium from aqueous solutions by treatment with carbon aerogel electrodes. In the engineering and applied physics area, the robust design for torque optimization using RSM was reported by Gao et al.13 Robust design optimization was performed on torque performance to investigate the electromagnetic behavior of the motor with changes of design parameters. In order to reduce the experimental runs, a mixed resolution rotatable central composite design was introduced. In membrane science and technology, Idris et al.14 applied RSM for the development of a thin film composite (TFC) membrane by an interfacial polymerization reaction between aqueous and organic phases on a polysulfone support. The experimental plan was based on CCD. The considered factors were the composition of the aqueous phase, which includes the ratio of m-phenyldiamine (MPDA) to hydroquinone (HQ) as the monomer ratio, percent amount of tetrabutylammonium bromide (TBAB) as catalyst, and percent amount of sodium hydroxide (NaOH) as acid acceptor. The performance of the TFC membranes in terms of the rejection and flux were the response variables investigated.
Table 1. Membrane Characteristics: Membrane Thickness, δM; Liquid Entry Pressure of Water, LEPw; Void Volume, E; Mean Pore Size, µp; Effective Porosity, E/Lp membrane
δM (µm)
LEPw (bar)
(%)
µp d (nm)
/Lp d (1/m)
TF200a TF450a GVHPb M12c
54.8 60.0 117.7 50.9
2.76 1.49 2.04 3.41
68.7 64.3 70.1
198.96 418.82 283.15 22.86
7878.14 7439.02 2781.93 3773.12
a Membrane supplied by Gelman (pore size 0.2 µm, porosity 80%). The measured total thickness was 165.2 µm for TF200 and 170.4 µm for TF450. b Membrane supplied by Millipore with the trade name Durapore (pore size 0.22 µm, porosity 75%). c Laboratory-made membrane.15 d The terms µp and /Lp were determined from the gas permeation test.16
Table 2. Independent Variables, Their Coded Levels, and Actual Values Used for DCMD Experimental Design real values of coded levels variable
symbol
stirring velocity, W (rpm) feed temperature, T1 (°C) NaCl concentration, C (M)
x1 x2 x3
-R
-1
0
+1
+Ra
88.2 150 437.5 725 786.8 22.3 25 37.5 50 52.7 0.007 0.2 1.1 2 2.193
a R ) 1.215 (star or axial point for orthogonal CCD in the case of three independent variables).
In the present paper, RSM, statistical design of MD experiments, MD variable factor interactions, and optimization of MD processes together with prediction of permeate fluxes were applied in the DCMD of salt (NaCl) aqueous solutions. 2. Experimental 2.1. Membranes and Direct Contact Membrane Distillation Experiments. Three commercial microporous hydrophobic flat-sheet membranes of polytetrafluoroethylene (TF200 and TF450, Gelman) supported by a polypropylene net and polyvinylidene fluoride (GVHP, Millipore) and a laboratory-made membrane (M12)15 were used. Their principal characteristics are summarized in Table 1. DCMD experiments were conducted for salt (NaCl) aqueous solutions using the laboratory system described elsewhere.16 The central part of the system is a stainless steel cell composed of two cylindrical chambers. One of the chambers is connected to a heating system through its jacket to control the temperature of the liquid feed. The other chamber is connected to a cooling system to control the temperature of the permeate. The membrane was placed between the two chambers (feed side and permeate side). The effective membrane area is 2.75 × 10-3 m2. The bulk feed and permeate temperatures were measured, after steady state was reached, inside each chamber by a pair of sensors connected to a digital meter with an accuracy of (0.1 °C. Both the feed and permeate liquids were stirred inside the cell by graduated magnetic stirrers. The DCMD flux was calculated in every case by measuring the condensate collected in the permeate chamber for a predetermined period. In each experimental run, the feed concentration was maintained constant. The salt concentration in the feed and permeate aqueous solutions were measured by a conductivimeter 712 Ω Metrohm. The solute rejection factor, RS, was calculated using the following expression:
(
RS ) 1 -
)
Cp 100 C
(1)
where Cp and C are the solute concentration in the permeate and in the bulk feed solutions, respectively.
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5675 Table 3. Central Composite Design and DCMD Experimental Results responses (permeate fluxes) (10-6 m/s)
input variables stirring rate run number (N) and typea 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 a
TF200
TF450
M12
GVHP
W (rpm)
levelb x1
T1 (°C)
levelb x2
[NaCl] (M)
levelb x3
YTF200
YTF450
YM12
YGVHP
725 150 725 150 725 150 725 150 786.8 88.2 437.5 437.5 437.5 437.5 437.5 437.5
1 -1 1 -1 1 -1 1 -1 R -R 0 0 0 0 0 0
50 50 25 25 50 50 25 25 37.5 37.5 52.7 22.3 37.5 37.5 37.5 37.5
1 1 -1 -1 1 1 -1 -1 0 0 R -R 0 0 0 0
2 2 2 2 0.2 0.2 0.2 0.2 1.1 1.1 1.1 1.1 2.193 0.007 1.1 1.1
1 1 1 1 -1 -1 -1 -1 0 0 0 0 R -R 0 0
2.141 1.083 0.802 0.407 2.734 1.418 1.037 0.541 1.548 0.704 2.227 0.637 1.094 1.649 1.368 1.338
2.455 1.387 0.993 0.577 3.358 1.887 1.409 0.897 1.964 1.108 2.641 0.906 1.312 2.034 1.669 1.659
2.625 1.165 1.123 0.589 3.912 2.135 1.787 1.141 2.263 1.006 3.034 0.958 1.899 2.311 1.892 1.859
0.314 0.226 0.187 0.108 0.764 0.532 0.255 0.158 0.398 0.206 0.828 0.193 0.328 0.435 0.362 0.353
O1 O2 O3 O4 O5 O6 O7 O8 S1 S2 S3 S4 S5 S6 C1 C2
feed temperature
feed concentration
O ) orthogonal design points, C ) center points, S ) star or axial points. b -1 ) low value, 0 ) center value, +1 ) high value, (R ) star point value.
2.2. Experimental Design. The central composite design (CCD) is a method that can be efficiently applied to develop the second-order response model with few numbers of factors n (2 e n e 6). The CCD consists of three distinct portions: (1) full (2n) or fractional (2n-z g n + 1) design where the factor levels are coded to the usual low (-1) and high (+1) values; (2) axial points or “star” points on the axis of each variable at a distance R from the designed center; and (3) center points that can be replicated to provide an estimate of experimental error variance. Based on CCD, a quadratic approximation is commonly used for developing a response surface model using experimental design. In general, the response model of second-order is written as17 n
Yˆ ) b0 +
∑ i)1
n
bixi +
∑ i)1
n
biixi2 +
bijxixj + ξ ∑ i 600 rpm. This is attributed, as it is well-known, to the temperature polarization effects in both the feed and permeate sides of the membrane as well as to the concentration polarization effects in the feed
membrane side.21,22 On the other hand, the increasing of feed temperature (mean temperature) leads to an “exponential” increase of the DCMD flux indicating the exponential variation of the vapor pressure of the feed solution as the driving force of the transmembrane vapor pressure. Moreover, the decrease of the DCMD flux with increasing the feed salt concentration was also observed, as it was expected. This is attributed to the decrease of the feed vapor pressure with the NaCl concentration and also to the concentration polarization effects. More explanations on the effects of these individual DCMD parameters may be found in other DCMD papers.1,2,16,23-26 The same effects of the operating variables upon the response (predicted flux) were observed by graphical response surface analysis of the other membranes, namely, TF450 and M12. For this reason, the response surfaces plot of these membranes are not included in the paper. It is worth quoting that the response surface plots in the case of the membrane GVHP (Figure 3) exhibit higher curvature effects comparing with other membranes. At higher stirring rates, the DCMD flux decreases slightly according to the approximated response surface. This can be explained by all quadratic coefficients that appear in the approximation equation
5678
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
Figure 3. Response surface plots and a contour-line plot of predicted DCMD flux for the membrane GVHP (a) as a function of stirring velocity (W) and temperature of the feed solution (T1) at solute concentration C ) 1.1 M, (b) as a function of T1 and C at W ) 437.5 rpm, and (c and d) as a function of W and C at T1 ) 37.5 °C.
for the GVHP membrane as well as by the highest negative value of the quadratic coefficient attributed to the solute concentration factor (Table 4). In fact, the membrane GVHP exhibits the lowest DCMD fluxes compared to the other three membranes. According to graphical response surface analysis carried out for all membranes, the most significant individual effect upon the response (permeate flux) is attributed to the feed temperature rather than to the stirring rate effect. The feed solute concentration has the lowest individual effect upon the DCMD permeate flux in comparison to the feed temperature and stirring rate. In general for feed temperature, the positive effect upon the response was observed for all membranes, i.e., increasing T1 involves increasing the response (predicted flux). For the stirring rate, the positive effect on the response persists for all membranes up to 500-600 rpm, after that an asymptotic effect has been observed for TF200, TF450, and M12 membranes and slightly a negative effect has been observed for the membrane GVHP. The solute concentration factor in general exhibits a negative effect upon the response, i.e., increasing C leads to lower values of predicted fluxes for all membranes. A particular case for the membrane GVHP should be noted where the response surface plot indicates a slightly positive effect of the
C factor at the lower feed temperature (25-30 °C). This could be attributed to the highest positive linear empirical coefficient corresponding to the C factor for this membrane and to the low value of the interaction effect between the temperature and concentration, γ23, (Table 4). The surface response models cited previously (eqs 4 and 5) were used to predict the permeate flux at any particular values of stirring velocity, temperature, and initial solute concentration within the limits of the experimentation region. Therefore, other experimental runs than those given in Table 3 were carried out. It must be stated here that these last experimental runs were not involved for the empirical model building. More than 20 confirmations runs (additional points) were performed for each membrane in order to check experimentally the adequacy of the developed empirical models. Comparisons between the predicted permeate fluxes and the experimental ones are summarized in Figure 4. ANOVA was also conducted for the DCMD fluxes of each membrane and the obtained coefficients (R2 and Radj2) are shown in each case. It should be mentioned that in the case of membranes TF200, TF450, and M12, the calculated R2 values are lower than those given in Table 5 but still remain greater than 0.9. It is interesting to note that for the GVHP membrane the R2
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5679
Figure 4. Predicted DCMD flux vs experimental one for each membrane.
Table 6. Canonical Analysis Performed for Each Objective Function of the Tested Membranes in DCMD stationary points
canonical coefficients
membrane
x1S
x2S
x3S
TF200 M12 TF450 GVHP
-0.58 -8.809 -1.066 0.338
-2.204 -18.228 -2.505 -1.629
3.601 -11.991 2.132 0.428
a
Yˆ S
(10-6 m/s) 0.272 -4.672 0.327 0.277
λ1
λ2
λ3
nature of stationary point
-0.19 -0.294 -0.155 -0.099
-0.015 0.01 -0.02 0.058
0.105 0.132 0.124 -0.06
saddle point (OVRa) saddle point (OVR) saddle point (OVR) saddle point (OVR)
OVR ) outside of the valid region.
value was increased compared to the original one given in Table 5 indicating that the model for the GVHP membrane revealed a better prediction than was expected. In all cases, the predicted R2 values are in reasonable agreement with the adjusted statistics Radj2 values. Thus, RSM can be used for prediction of DCMD performance decreasing the number of experimental runs.
1 c -1 xjS ) - B bhL 2
(7)
The nature of the stationary point is defined by eigenvalues c by solving the characteristic equation18,19 λh of B
3. Optimization of the DCMD Process 3.1. Canonical Analysis. The canonical analysis is a compulsory technique in RSM that tells us if the objective function has a global optimal point (stationary point) and if this point is located inside the region of experimentation (valid region). Since the second-order models were found to be adequate for predictions, the canonical analysis was performed to determine the locations and the natures of the stationary points of the second-order models. The estimated second-order approximations can be written in matrix form as follows
c xj Yˆ ) b0 + xjTbhL + xjTB
c are the estimates of the intercept, linear, where b0, bhL, and B and second-order coefficients, respectively. The stationary point xjS of the second-order polynomial is determined as17
(6)
c - λhE c) ) 0 det(B
(8)
c is identity matrix. where E c are negative (positive), then the If all eigenvalues of B quadratic surface has a maximum (minimum) at the stationary point xjS. If the eigenvalues have mixed signs, then the stationary point xjS is a saddle point.17-19 The canonical analysis was carried out using MathCAD software, and the results of the canonical analysis performed for each objective function (i.e., empirical model developed for each membrane) are shown in Table 6.
5680
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
Table 7. Search for the Optimum Point of the Membrane TF200, as an Example, by Means of the Step Adjusting Gradient Method coded values of stirring rate
coded values of feed temperature
coded values of solute concentration
k
x1
m1
s1
x2
m2
s2
x3
m3
s3
Yˆ (k) (10-6 m/s)
e
0 1 2 3 4 5 6 7 8 9 10
0 0.467 0.712 0.869 0.972 1.041 1.091 1.126 1.151 1.169 1.182
0.533 0.401 0.345 0.316 0.299 0.289 0.281 0.277 0.273 0.271 0.269
0.877 0.611 0.455 0.327 0.233 0.172 0.126 0.092 0.067 0.048 0.035
0 0.419 0.646 0.809 0.925 1.008 1.067 1.109 1.139 1.161 1.176
0.81 0.877 0.898 0.908 0.913 0.916 0.917 0.919 0.919 0.92 0.92
0.517 0.259 0.181 0.128 0.091 0.065 0.046 0.033 0.024 0.017 0.012
0 -0.528 -0.765 -0.906 -0.998 -1.06 -1.104 -1.136 -1.159 -1.175 -1.186
-0.245 -0.264 -0.272 -0.276 -0.279 -0.28 -0.281 -0.282 -0.283 -0.283 -0.283
2.157 0.896 0.518 0.334 0.223 0.158 0.113 0.08 0.057 0.04 0.029
1.366 1.932 2.256 2.487 2.653 2.771 2.858 2.92 2.965 2.997 3.021
0.415 0.168 0.102 0.067 0.045 0.031 0.022 0.015 0.011 0.0077
Table 8. Optimal Points in Terms of Actual Operating Variables Determined by the Step Adjusting Gradient Method and Located inside of the Region of Experimentation: 88.2 e W e 786.8 rpm; 22.3 e T1 e 52.7 °C; 0.0065 e C e 2.1935 M membrane
W* (rpm)
T/1 (°C)
C* (M)
calc Yˆ * (10-6 m/s)
exp Y* (10-6 m/s)
TF200 M12 TF450 GVHP
777.4 777.4 777.2 616.0
52.21 52.21 52.21 51.88
0.032 0.036 0.034 0.051
3.021 4.193 3.687 0.839
2.897 4.343 3.642 0.889
It was found that all stationary points corresponding to each objective function were located outside of the valid region (i.e., region of experimentation). The canonical coefficients indicate that the shapes of response surfaces are hyperbolic paraboloids and the stationary points, in fact, are saddle points. This means that stationary points cannot be accepted as optimal solutions of the objective functions. In order to determine the feasible (i.e., local) optimal points inside of the valid region, a numerical optimization technique was applied using the gradient method.27,28 3.2. Numerical Optimization: Step Adjusting Gradient Method. The gradient basic method of steepest ascent is a cyclic search technique, which does not follow the continuous line of steepest ascent, but approximates it by a succession of straight lines. Each line corresponds to a stage in the search in which the local function, gradient, is initially determined. A set of moves is then made along the resulting direction of steepest ascent, searching one-dimensionally for the approximate location of the restricted optimum. This point serves as the new base point for the next stage, and the same procedure is repeated.27 The algorithm of the gradient method is thoroughly shown in the Appendix. The calculation procedure proposed for step adjustment in the gradient method is more important from an RSM standpoint where searching for the optimal point should be performed in a small region of interest. More details of the proposed procedure are indicated in the Appendix by a flowchart algorithm. The results of searching by means of the gradient method using the step adjusting procedure are shown in Table 7 for the case of the TF200 membrane as an example. The regression model in terms of coded variables is considered as the objective function (eq 4). The search starts from the base point xj(k)0) ) {0, 0, 0}T, and the direction of the steepest ascent is determined, i.e., m j (k)0) ) {0.533, 0.810, -0.245}T. The corresponding cj (k)0) will be cj(k)0) ) {0.533, 0.810, 0.245}T. This means that the following conditions are satisfied for R ) 1.215: Rβ5 e c(k)0) 1 < Rβ6; Rβ7 e c(k)0) < Rβ8; Rβ3 e c(k)0) < R β4 hence δ h (k)0) ) 2 3 { 2.6, 2.9, 2.3 }T. According to eq A.4 (in the Appendix), the adjusting searching step will be js (k)0) ) {0.877, 0.517, 2.157}T. Knowing the direction of steepest ascent and the adjusted
searching step, the next point is established by eq A.3, i.e., xj(k)1) ) {0.467, 0.419, -0.528}T. A new local exploration is then performed starting from xj(k)1), and a new direction of steepest ascent is determined. The calculation is repeated adjusting the step length at each searching stage according to eq A.4 in order to keep the values of coded variables in the valid region. The optimum searching is finally stopped when the convergence criteria e e e0 (e0 ) 0.01) is reached.
e)
|
|
Yˆ (k+1) - Yˆ (k) Yˆ (k)
e e0
(9)
Table 7 shows that the whole search for the optimal point using the proposed adjusting step procedure was performed inside the valid region. The feasible optimal point given in terms of coded variables is the maximal one and is located in the region of interest, i.e., valid interval -1.215 e xj e 1.215 (j ) 1, 2, 3), where the values (1.215 corresponds to the axial (star) point R in orthogonal CCD for n ) 3 factors. A similar procedure was carried out for all membranes. It must be mentioned here that, in the case of objective functions with more than two independent variables, the evolution profile of optimum searching is difficult to represent by means of response surface contour plots. In this respect, the radar diagrams are more suitable for suggestive representation of the optimum searching. Such diagrams illustrate the searching evolution from the base point to the optimum point outlining the number of points (i.e., iterations) in the searching profile and the corresponding values of the objective function in these points (permeate flux in our case). The value of the objective function in the center of the radar diagram is minimal, and at the edge, it is maximal. For all tested membranes, the radar diagrams of optimal searching are illustrated in Figure 5. As can be seen, by applying the proposed step adjusting gradient method, the values of objective function were calculated in 10 iterations (i.e., 10 points) for the TF200, TF450, and M12 membranes and in 9 iterations for the GVHP membrane. Thus, the convergence of the method is almost the same for all objective functions. The calculated feasible optimal point corresponds to the maximal values of predicted permeate fluxes inside the valid region. The founded optimal points in terms of actual operating variables were verified experimentally, and the results are shown in Table 8. It can be observed from Table 8 that for the TF200, TF450, and M12 membranes the optimal operational conditions are practically the same and are located near the boundary of the valid region. In this case, the average values of optimal points can be considered as optimal operational conditions for these membranes: W* ) 777.3 rpm; T/1 ) 52.21 °C; C* ) 0.034 M.
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5681
the gradient method, was developed in order to determine the feasible optimal points in the region of experimentation of the DCMD process for each membrane. The optimum conditions of the DCMD process were found to depend on the type of membrane. It was found that the membranes having higher permeate fluxes (TF200, TF450, and M12) exhibit different optimal conditions (W* ) 777.3 rpm; T/1 ) 52.21 °C; C* ) 0.034 M) than the membrane GVHP having lower flux (W* ) 616.0 rpm; T/1 ) 51.88 °C; C* ) 0.051 M). The existence of the highest values of DCMD fluxes using the determined optimum DCMD conditions was confirmed experimentally. Appendix Algorithm of the Step Adjusting Gradient Method. According to the gradient algorithm, the optimal point was searched by means of eq A.1.
jx(k+1) ) xj(k) + m j (k)‚s(k)
(A.1)
(k) (k) T where xj(k) ) {x(k) 1 x2 ... xn } is the vector of coded variables, (k) (k) (k) m j (k) ) {m1 m2 ... mn }T is the vector indicating the direction of steepest ascent, s(k) is the step of searching, and k is the stage of searching. The directions of steepest ascent were calculated by means of the gradient equation, i.e., eq A.2, as follows27,28
m j (k) )
Figure 5. Progress profile of optimum searching by means of step adjusting gradient method.
For the membrane GVHP, the optimal point is slightly far from the boundary of the valid region compared to the other three membranes exhibiting higher DCMD fluxes. It must be pointed out that in all cases the optimal points have been checked experimentally and were found to be the best experimental points giving maximal fluxes in the studied region of experimentation. 4. Conclusions The response surface methodology and experimental design were applied in direct contact membrane distillation when using salt (NaCl) aqueous solutions. It was proved that it is adequate in assessing the permeate fluxes for both commercial- and laboratory-made membranes. The existence of interactions between the operating factors, namely, the stirring rate, feed temperature, and solute concentration, were observed and graphically presented. As it was expected, the operational DCMD parameters affect the DCMD permeate flux in the same way as observed from the classical method of experimental testing. An algorithm, taking into consideration step adjustment by exploring the response surface only in the valid region using
∇Yˆ (xj(k)) ||∇Yˆ (xj(k))||
(A.2)
In the gradient basic method of steepest ascent, the step of searching s(k) is usually considered to be constant at each stage (k) of searching, being the same for each independent variable.27 In response surface methodology (RSM), the objective function is valid over the small range of experimentation defined by the constraints of type -R e xi e +R, being i ) 1, 2, ..., n and R the axial or star point depending on the type of CCD used and the number of factors. If the searching step s(k) in the gradient method is selected to be constant and high, it is possible to readily obtain a point outside of the valid region in the next iterations. The using of a normal (average) constant step for all variables does not control the searching near the boundary of the experimentation region well, and in this case, the risk of valid region violation may occur. If the step of searching will be taken as a constant and too small, a lot of calculation will be needed for optimum localization. For this reason, in this paper, an alternative procedure is proposed for adjusting the searching step in the gradient method in order to determine the feasible optimal point by exploring the response surface only in the valid region. To apply the proposed algorithm, one should use the objective function in terms of coded variables. The start or base point is recommended to be the center point of the experimental design, i.e., xj(0) ) {0 0 ... 0}T. The basic idea of the proposed step adjusting procedure relys on the fact that the step of searching will be higher near the center point of the experimental design and smaller near the limit of the valid region. In the proposed procedure (see Figure A.1), the step of searching depends on the value of the star point R, the values of directions of steepest ascent m j (k), and the values of vector components of coded variables xj(k). In this case, the step of searching could be different for each variable at each stage of searching, and for this reason,
5682
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5683
Figure A.1. Flowchart algorithm of the step adjusting gradient method.
the step will be considered as a vector js(k) for each stage of searching. Thus, eq A.1 becomes
jx(k+1) ) xj(k) + m j (k)‚sj(k)
(A.3)
(k) (k) T where js(k) ) {s(k) 1 s2 ... sn } and
s(k) j
(k) 1 R - |xj | , j ) 1, 2, ..., n ) (k) δj |m(k) j |
(A.4)
where δ(k) is the step diminishing coefficient that should j be at least higher than 2 in order to ensure a displacement
from the point (xj(k)) up to the next point (xj(k+1)) lower than 50% of the distance between xj(k) and the boundary of the valid region (R). Therefore, the interval 2 e δ(k) j e 3.5 was adopted for the step diminishing coefficient. In this case, the δ(k) j (k) (k) coefficient is minimal, i.e., δ(k) j ) 2, if the value cj ) |xj + m(k) j | is near 0 (i.e., near to the center of the experimental (k) (k) design) and is maximal, i.e., δ(k) j ) 3.5, if the value cj ) |xj (k) + mj | is near R (i.e., near to the boundary of the experimental design). The c(k) value gives information about the position (in j absolute value) of the “next point” (fictive point) considering
5684
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007
Table A1. Subdomain Coefficients for Optimum Searching Using the Step Adjusting Gradient Method h
1
2
3
4
5
6
7
8
9
10
11
βh rh
0 2
0.1 2.15
0.2 2.30
0.3 2.45
0.4 2.60
0.5 2.75
0.6 2.90
0.7 3.05
0.8 3.20
0.9 3.35
1 3.5
the step of searching equal to unity (i.e., js(k) ) 1). The variation (k) interval for the δ(k) j coefficient depends on the cj value inside the valid region interval (positive interval, i.e., 0 ( R). In this paper, it is also proposed to divide the valid region interval (0 ( R) into 10 subdomains from the center point up to the limit of the valid region defined by the axial point R. The limits of each subdomain are given by the relationship (R· βh), where the range of the βh coefficient is 0 e βh e 1 for h ) 1, 2, ..., q + 1, q being the number of subdomains (q ) 10). (k) (k) If the c(k) j value is in the range Rβh e cj < Rβh+1, then δj ) rh, where 1 e rh e 3.5 according to the values given in the Table A.1. Nomenclature b0, bi, bii, bij ) regression coefficients bh (u × 1) ) vector of regression coefficients b0 ) intercept regression coefficient bhL ) vector of linear regression coefficients c ) symmetric matrix of second-order regression coefficients B c(k) j ) step adjusting coefficient Cp ) solute concentration in the permeate solution C ) solute concentration (NaCl) in the feed aqueous solution, M C* ) optimal value of solute concentration in the feed aqueous solution, M c ) identity matrix E e ) convergence criteria chosen for optimization e0 ) imposed (stated) convergence F ) ratio of variances k ) stage (or iteration) of searching for the optimum in the gradient method LEPw ) liquid entry pressure of water, bar m j (k) ) vector indicating the direction of steepest ascent of the objective function N ) the number of experimental runs n ) number of factors (independent variables) q ) the number of subdomains considered into the valid region R2 ) coefficient of multiple determination Radj2 ) adjusted statistic coefficient rh ) subdomain coefficient in step adjusting gradient method SSRegression ) sum of squares of the regression model SSResidual ) sum of squares of the residual s(k) and js(k) ) thestep of searching for the optimal point in the gradient method T1 ) feed temperature, (°C) T/1 ) optimal value of feed temperature, °C ∆T ) temperature difference between the bulk feed and permeate, °C u ) the number of significant regression coefficients in the response surface model W ) stirring rate, rpm W* ) optimal value of stirring rate, rpm c (N × u) ) matrix of the independent variables levels X xjS ) vector of coordinates of the stationary point in coded terms xj(k) ) vector of coded variables x1, x2, x3, xI ) the coded levels of the factors (independent or control variables)
Y h (N × 1) ) vector of the experimental observations (permeate flux) Y ) permeate flux (experimental value) Yˆ ) permeate flux (predicted value by response surface model) Y* ) experimental maximal value of permeate flux according to optimal parameters Yˆ * ) maximal predicted value of permeate flux according to optimal parameters ∇Yˆ (xj(k)) ) gradient of the objective function z ) integer number used to determine the number of runs in fractional factorial design Greek Letters R ) axial points or “star” points in CCD RS ) the solute rejection factor, percent βh ) subdomain coefficient in step adjusting gradient method δ(k) j ) step diminishing coefficient in step adjusting gradient method δM ) membrane thickness, µm ) void volume, percent p/Lp ) effective porosity, 1/m γi ) empirical coefficients related to variables in actual scale λh ) vector of eigenvalues of B h matrix µp ) mean pore size, nm ξ ) statistical error in response surface model AbbreViations ANOVA ) analysis of variances CCD ) central composite design DCMD ) direct contact membrane distillation MLR ) multilinear regression method Literature Cited (1) Mengual, J. I.; Pen˜a, L. Membrane distillation. Colloid Surf. Sci. 1997, 1, 17-29. (2) Lawson, K. W.; Lloyd, D. R. Membrane distillation: review. J. Membr. Sci. 1997, 124, 1-25. (3) El-Bourawi, M. S.; Ding, Z.; Ma, R.; Khayet, M. A framework for better understanding membrane distillation separation process: review. J. Membr. Sci. 2006, 285, 4-29. (4) Ravikumar, K.; Pakshirajan, K.; Swaminathan, T.; Balu, K. Optimization of batch process parameters using response surface methodology for dye removal by a novel adsorbent. Chem. Eng. J. 2005, 105, 131138. (5) Lacin, O.; Bayrak, B.; Korkut, O.; Sayan, E. Modeling of adsorption and ultrasonic desorption of cadmium(II) and zinc(II) on local bentonite. J. Colloid Interface Sci. 2005, 292, 330-335. (6) Hamsaveni, D. R.; Prapulla, S. G.; Divakar, S. Response surface methodological approach for the synthesis of isobutyl isobutyrate. Process Biochem. 2001, 36, 1103-1109. (7) Palamakula, A.; Nutan, M. T. H.; Khan, M. A. Response Surface Methodology for Optimization and Characterization of Limonene-based Coenzyme Q10 Self-Nanoemulsified Capsule Dosage Form AAPS. Pharm. Sci. Tech. 2004, 5, 1-8. (8) Zhao, H.; Tonkyn, R. G.; Barlow, S. E.; Peden, C. H. F.; Koel, B. E. Fractional factorial study of HCN removal over a 0.5% Pt/Al2O3 catalyst: effects of temperature, gas flow rate and reactant partial pressure. Ind. Eng. Chem. Res. 2006, 45, 934-939. (9) Lemoine, R.; Behkish, A.; Morsi, B. I. Hydrodynamic and masstransfer characteristics in organic liquid mixtures in a large-scale bubble column reactor for toluene oxidation process. Ind. Eng. Chem. Res. 2004, 43, 6195-9212. (10) Ahmadi, M.; Vahabzadeh, F.; Bonakdarpour, B.; Mofarrah, E.; Mehranian, M. Application of the central composite design and response surface methodology to the advanced treatment of olive oil processing wastewater using Fenton’s peroxidation. J. Hazard. Mater. B 2005, 123, 187-195. (11) Ramirez, J. H.; Costa, C. A.; Madeira, L. M. Experimental design to optimize the degradation of the synthetic dye Orange II using Fenton’s reagent. Catal. Today 2005, 107-108, 68-76.
Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5685 (12) Rana-Madaria, P.; Nagarajan, M.; Rajagopal, C.; Garg, B. S. Removal of Chromium from aqueous solutions by treatment with carbon aerogel electrodes using response surface methodology. Ind. Eng. Chem. Res. 2005, 44, 6549-6559. (13) Gao, X. K.; Low, T. S.; Liu, Z. J.; Chen, S. X. Robust design for torque optimization using response surface methodology. IEEE Trans. Magn. 2002, 32, 1141-1144. (14) Idris, A.; Kormin, F.; Noordin, M. Y. Application of response surface methodology in describing the performance of thin film composite membrane. Sep. Purif. Technol. 2006, 49, 271-280. (15) Khayet, M.; Matsuura, T. Preparation and characterization of polyvinylidene fluoride membranes for membrane distillation. Ind. Eng. Chem. Res. 2001, 40, 5710-5718. (16) Khayet, M.; Mengual, J. I.; Matsuura, T. Porous hydrophobic/ hydrophilic composite membranes: application in desalination using direct contact membrane distillation. J. Membr. Sci. 2005, 252, 101-113. (17) Carley, K. M.; Kamneva, N. Y.; Reminga, J. Response Surface Methodology; Center for Computational Analysis of Social Organizational Systems, CASOS Technical Report; CMU-ISRI-04-136, Carnegie Mellon University, School of Computer Science, Pittsburgh, PA, 2004. (18) Akhnazarova, S.; Kafarov, V. Experiment Optimization in Chemistry and Chemical Engineering; Mir Publishers: Moscow, 1982. (19) Taloi, D. Optimization of Technological Processes (in Romanian); Ed. Acad.: Bucharest, 1987. (20) Mahat, M. K.; Illias, R. M.; Rahman, R. A.; Rashid, N. A. A.; Mahmood, N. A. N.; Hassan, O.; Aziz, S. A.; Kamaruddin, K. Production of cyclodextrin glucanotransferase (CGTase) from alkalophilic Bacillus sp.
TS1-1: media optimization using experimental design. Enzyme Microb. Technol. 2004, 35, 467-473. (21) Martı´nez-Die´z, L.; Va´zquez-Gonza´lez, M. I. Temperature and Concentration polarization in membrane distillation of aqueous salt solutions. J. Membr. Sci. 1999, 156, 265-273. (22) Phattaranawik, J.; Jiraraananon, R.; Fane, A. G. Heat transport and membrane distillation coeficients in direct contact membrane distillation. J. Membr. Sci. 2003, 212, 177-193. (23) Schofield, R. W.; Fane, A. G.; Fell, C. J. D.; Macoun, R. Factors affecting flux in membrane distillation. Desalination 1990, 77, 279-294. (24) Alklaibi, A. M.; Lior, N. Comparative study of direct contact and air gap membrane distillation processes. Ind. Eng. Chem. Res. 2007, 46, 584-590. (25) Khayet, M.; Godino, M. P.; Mengual, J. I. Modelling transport mechanism through a porous partition. J. Non-Equilib. Thermodyn. 2001, 26 (1), 1-14. (26) Li, B.; Sirkar, K. K. Novel membrane and device for direct contact membrane distillation-based desalination process. Ind. Eng. Chem. Res. 2004, 43, 5300-5309. (27) Beveridge S. G., Schechter, R. S. Optimization: Theory and Practice; McGraw-Hill: NewYork, 1970. (28) Curievici, I. Optimization in Chemical Engineering (in Romanian); Ed. Did. Ped.: Bucharest, 1980.
ReceiVed for reView March 27, 2007 ReVised manuscript receiVed May 31, 2007 Accepted June 11, 2007 IE070446P
