Analysis of the Dynamic Response of a Reverse Osmosis Membrane

Ind. Eng. Chem. Res. 2005, 44, 7823-7834

7823

Analysis of the Dynamic Response of a Reverse Osmosis Membrane to Time-Dependent Transmembrane Pressure Variation A. Alexiadis,*,† J. Bao,† D. F. Fletcher,‡ D. E. Wiley,† and D. J. Clements§ School of Chemical Engineering and Industrial Chemistry and School of Electrical Engineering and Telecommunications, The University of New South Wales, High Street, Sydney NSW 2052, Australia, and Department of Chemical Engineering, The University of Sydney, Paramatta Road, Sydney NSW 2006, Australia

In this paper, the dynamic behavior of reverse osmosis (RO), high-pressure (60-80 bar) membranes is investigated using an approach that combines both computational fluid dynamics (CFD) and system identification methods. Detailed CFD models of membrane systems allow predictions of the response of such systems to operating parameter changes. These models, however, are complex and require significant computing resources, which make them impractical for use in membrane system operation. In the present work, system identification theory is applied to data generated from a CFD model in order to obtain approximate transfer functions, which describe the dynamics of the system. The transfer function parameters are then related to the dynamic changes in momentum and concentration. The results show that the local concentration plays a fundamental role in the membrane performance. A feedback model, which quantifies the effect of the concentration boundary layer on the permeate flux, is finally proposed. 1. Introduction In a reverse osmosis (RO) process, two channels (called the feed and permeate channels) are separated by a membrane, which acts as a barrier to the flow, allowing selective passage of a particular species (usually the solvent) while other species (solutes) are retained partially or completely. If a transmembrane pressure is applied between the two channels, the solution is driven through the membrane and a phase rich in solvent and poor in solutes is collected in the permeate channel. Separation occurs at the membrane, and consequently, a layer with a higher solute concentration with respect to the bulk feed solution is formed at the membrane surface on the feed side. This layer is called the concentration boundary layer, while the whole phenomenon is known as concentration polarization and has a negative consequence on the membrane performance because the solvent flux decreases due to the increase in osmotic pressure. This phenomenon is, however, reversible, and the original situation can be restored if the boundary layer is destroyed. Nevertheless, unless a strategy against it is adopted, the concentration polarization can lead to an even worse scenario. In the case of a saltwater solution, for example, the concentration in the boundary layer can exceed the solubility limit of the salt, causing it to precipitate on the membrane surface, which can in turn increase the resistance to flow and reduce the solute flux permanently. A common strategy to counteract the phenomenon of concentration polarization involves the use of flow disturbances, which can temporarily destroy the * To whom correspondence should be addressed. Tel.: +61 (2) 9385 4382. Fax: +61 (2) 9385 5966. E-mail: a.alexiadis@ unsw.edu.au. † School of Chemical Engineering and Industrial Chemistry, The University of New South Wales. ‡ Department of Chemical Engineering, The University of Sydney. § School of Electrical Engineering and Telecommunications, The University of New South Wales.

boundary layer by means of a sudden modification of the hydrodynamics. Furthermore, these disturbances also enhance the shear on the membrane surface and can sweep away deposited particles. The importance of these flow disturbances in increasing the membrane performance has already been demonstrated in the literature (e.g., refs 1-3); however, the works published until now provide experimental observations with limited theoretical insights. In this paper, we investigate the physics of the mechanisms involved in the destruction and reformation of the boundary layer; a correct knowledge of these mechanisms would be useful for improving module performance using strategies that rely on flow disturbances. To properly understand the mechanisms of boundary layer formation, a large amount of data on the hydrodynamics, such as the local concentrations and velocities in the membrane channels, is required. The type of detailed information required is very extensive and cannot be obtained simply through experiments. Computational fluid dynamics (CFD), however, is the ideal tool to use to obtain such information, because it solves the equations of motion at each time step and at each point in the system. Although the results obtained by a CFD code provide a complete solution for all the variables involved, they are, nevertheless, unique to the specific period, amplitude, and shape of the disturbance used in each numerical simulation. A change in any of the characteristics of the disturbance would necessitate a new simulation. For this reason, a theoretical method that can go beyond the specific situation and provide the intrinsic response of the system regardless of the shape of the disturbance is desirable. “Transfer functions”, often used by control or electrical engineers, provide a good technique for this method because they work “from function to function” rather than “from value to value”. Mathematically, they are operators rather than functions. This means that, given the shape of the transmembrane pressure (TMP) disturbance, the transfer function automatically provides the shape of the permeate flux disturbance. If the

10.1021/ie050290y CCC: $30.25 © 2005 American Chemical Society Published on Web 09/02/2005

7824

Ind. Eng. Chem. Res., Vol. 44, No. 20, 2005

Figure 1. Schematic representation of a membrane channel. Table 1. Channel Geometry and Membrane Parameters feed channel length (L1) feed channel height (H1) permeate channel length (L2) permeate channel height (H2) membrane permeation coefficient (K) membrane rejection coefficient (R)

274 mm 2.5 mm 260 mm 2.0 mm 1.81 × 10-6 m s-1 bar-1 0.995

input shape changes, the same transfer function can be straightforwardly used without any new simulation (as is required in the CFD case). Initially, CFD computations are required in order to calibrate the transfer function, but once this stage is accomplished, the use of transfer functions is simpler and faster. In this paper, we use some modified system identification techniques to develop the transfer functions that relate pressure and concentration disturbances to permeate flux for a simple RO system filtering sodium chloride. The data for the RO system is evaluated using a model developed previously.4 1.1. CFD Modeling of Flow in Membrane Channels. Modeling of fluid flow in membrane systems has been undertaken since the early 1960s (e.g., refs 5, 6, and 7), but until the development of CFD software, the results achieved were limited to steady-state conditions and were based on major simplifications of, for instance, the velocity field and flux through the membrane. Recently (e.g., refs 4 and 8), complete CFD models have been developed and validated against classical cases available in the literature. The model used in this work is adopted from ref 8, where the permeate flux through the membrane is computed using a law analogous to the Darcy law for porous media and the mass and momentum transfer across the membrane is modeled using appropriate source terms. Viscosity, density, and diffusion coefficients are considered to vary with the solute concentration. The study is carried out using CFX 4.4 for a saltwater system and the geometrical configuration reported in Figure 1 and Table 1 (details are reported in ref 8). The salt concentration ranges between 25 and 35 g L-1. The transmembrane pressure varies between 60 and 80 bar. These conditions are consistent with highpressure reverse-osmosis seawater desalination.9 The time step used for the numerical computations was set to 2 × 10-5 s for t e tω (the time for which the perturbation was applied) and, after the disturbance (t > tω), was gradually increased to 5 × 10-3 s (nonuniform sampling). Issues connected with pH variation were not considered. 1.2. Variational Variables. A pulse can be defined as a sudden variation in the input variable(s) x from

the starting point x1 to a new value x2 and a sudden restoration to the initial conditions. A variational (or perturbation) variable expresses the variation (x - x1) of the variable x from a starting steady-state condition x1. Dividing it by the overall variation of the impulse between x2 and x1 results in a dimensionless variational variable.

x*(t) )

x(t) - x1 x2 - x1

(1)

In the case of permeate flux J, for example, J* is the dimensionless variational permeate flux, J1 is the flux value at TMP1 (steady state 1), and J2 is the flux value at TMP2 (steady state 2). Using dimensionless variational variables as in eq 1, the two steady-state conditions are always x*1 ) 0 and x*2 ) 1. 2. Transfer Functions In this work, CFD data were analyzed numerically in order to identify the transfer functions for a simulated reverse osmosis process. Transfer functions describe the relationship between the input and output variables in the Laplace domain (see Appendix for an overview or ref 10 for a complete discussion). Often, they are represented as a rational function of the type

G(s) )

(τas + 1)(τbs + 1)(τcs + 1) ‚‚‚ (τ1s + 1)(τ2s + 1)(τ3s + 1) ‚‚‚

(2)

where the parameters -1/τ1,2,3 are called poles and the parameters -1/τa,b,c are zeros. The number of poles and zeros determines the model family with a typical dynamic behavior (a model with n poles and m zeros is often called an n/m model). The numerical values of the parameters τ1,2,3,a,b,c identify a single member of that specific model family, and usually, they are determined by fitting experimental (or, as in this paper, CFD) data. The first step consists of choosing the most significant variables of the process to be used as the input(s) and the output(s) of the transfer function. The easiest way to generate flow disturbances in an RO module is a rapid variation in transmembrane pressure; for this reason, the TMP is chosen as the input variable. The membrane performance is measured by the permeate flux, and consequently, this variable is chosen as the output variable. The permeate flux response to a 20 bar transmembrane pressure pulse (∆P ) TMP1 - TMP2) was computed using the detailed membrane model of

Ind. Eng. Chem. Res., Vol. 44, No. 20, 2005 7825

Figure 2. Amplitude ratio AR (output ) permeate flux) for a TMP displaced cosine impulse with duration tω ) 0.001 s (c0 ) 25 g L-1 and TMP1 ) 80 bar).

Figure 3. Angle shift φ (output ) permeate flux) for a TMP displaced cosine impulse with duration tω ) 0.001 s (c0 ) 25 g L-1 and TMP1 ) 80 bar).

ref 4. The results (expressed in terms of dimensionless variational permeate flux) were found to be independent of the disturbance amplitude (in the range of 5 bar e ∆P e 30 bar). Consequently, linear dynamic models can be adopted. This does not mean that this phenomenon is linear but only that it can be sufficiently approximated by a linear dynamic model, within the range of the parameters and the amplitude of the perturbations considered. Under these conditions, the nonlinear partial differential equations describing the relationship between input and output parameters can be simplified to a series of linear ordinary differential equations and the method proposed in this paper can be used. As shown in the Bode diagrams (i.e., amplitude ratio AR and phase shift φ) in Figures 2 and 3, a displaced cosine impulse with tω ) 0.001 s was used at TMP ) 80 bar and c0 ) 25 g L-1. The independence of the Bode diagrams to tω or impulse shapes was verified. The “without osmotic pressure” points in Figures 2 and 3 will be discussed in a following section.

The analysis of the CFD data for the permeate flux as a function of time (the triangles in Figures 2 and 3) suggests a model with an equal number of poles and zeros. For the different models considered (1-pole/1-zero, 2-poles/2-zeros, and 3-poles/3-zeros in Figures 4 and 5), the model parameters (τ1, τa, etc.) were determined in order to minimize the error between the CFD data and the fitted model (Table 2). A 3/3 model (see Figures 4 and 5) provides the most adequate fit to the dynamic behavior of the system. More complex models (4/4 or 5/5) increase the number of parameters without significant improvements in the fitting. The same analysis was performed for different initial conditions in the range of 60-80 bar and salt concentration of 25-35 g L-1. The optimal fitting parameters change only slightly from the ones reported in Table 2, and consequently, they can be assumed constant over the operating range considered. 2.1. Characteristic Times. The time constants, τ1, τ2, and τ3 reported in Table 2 (-1/τ defines the pole,

7826


Figure 4. Fit to the data from Figure 2, with different types of models (1/1, 2/2, and 3/3).

Figure 5. Fit to the data from Figure 3, with different types of models (1/1, 2/2, and 3/3). Table 2. Fitting Parameters for 1/1, 2/2, and 3/3 Models 1/1 model

2/2 model

3/3 model

τ1 τ2 τ3 τa τb τc

0.23 ( 0.01

0.42 ( 0.01 0.02 ( 0.001

0.72 ( 0.04

1.07 ( 0.03 0.03 ( 0.001

0.6 ( 0.005 0.06 ( 0.001 0.006 ( 0.0001 1.27 ( 0.009 0.09 ( 0.002 0.007 ( 0.002

χ2

9.41

1.09

0.06

while τ is the time constant), deserve particular attention because they represent the characteristic times of the process. This can be seen from the following example. The response of the dynamic 3-pole/3-zero model to a step input condition (∆P is constant, M) can be represented analytically as

(

J*(t) ) M 1 +

)

R1 -t/τ1 R2 -t/τ2 R3 -t/τ3 e + e + e τ1 τ2 τ3

(3)

where R1, R2, and R3 are functions of τ1,2,3,a,b,c. The magnitudes of τ1, τ2, and τ3 are more significant than

those of τa, τb, and τc, because they appear in the exponential terms of eq 3 as characteristic times of decay. The dynamic component e-t/τ, which is associated with the time constant of τ, reaches, for instance, 86% and 99% of the new steady state in t ) 2τ and t ) 5τ, respectively. In Figure 6, we show the variation in dimensionless concentration c* response to a step transmembrane pressure input at a point near the membrane surface (10 µm above the surface) and in the middle of the membrane (15 cm from the inlet). Figure 7 shows the variation in dimensionless fluid velocity (v/w) normal to the membrane at the same point (since the point is very close to the membrane, this value can be considered approximatively equal to the local velocity v/w through the membrane). In the simulations, the dynamics of v/w(t) and J(t) are equivalent. When and why this assumption is correct is discussed in the next section. The dynamics shown in Figures 6 and 7 can be associated with the characteristic times in Table 2 and can be understood mechanistically as depicted in Figure 8.


Figure 6. c*(t) response to a 20 bar TMP step at a point 10 µm from the membrane and 150 mm from the inlet (TMP1 ) 80 bar and c0 ) 25 g L-1).

Figure 7. v/w(t) response to a 20 bar TMP step at a point 10 µm from the membrane and 150 mm from the inlet (TMP1 ) 80 bar and c0 ) 25 g L-1).

Figure 8. Consequences of a ∆P pulse on the local concentration.

Figure 8a shows the situation at steady state 1: the concentration boundary layer δ0 defines the zone where the concentration is higher than that of the bulk solution, c(δ0) ) 0.01(cw - cbulk). The velocity v0w reacts immediately to ∆P and increases to v′w (see Figure 7),

but other variables, such as the concentration (which is proportional to ∆π), have slower dynamics (Figures 6 and 8b). Only after 0.001 s is there a reduction in the osmotic pressure (in Figure 6 there is a slight decrease of c* between 0.001 and 0.01 s) because the higher

7828


Figure 9. Amplitude ratio AR (output ) local concentration) for a TMP displaced cosine impulse with duration tω ) 0.001 s (c0 ) 25 g L-1 and TMP1 ) 80 bar).

velocity “thins” the boundary layer (Figure 8c). As time increases (t > 0.015 s), new salt is accumulated at the membrane surface. The osmotic pressure consequently increases and the velocity v′w decreases to vFw according to the equation

vw(l,t) ) K[P - ∆π(l,t)]

(4)

where vw is the local transmembrane velocity, ∆π is the local osmotic pressure, and l is a particular location along the membrane. (This equation is nonlinear because the osmotic pressure depends on the concentration, which in turn depends on the local flow field, which is determined from the nonlinear (in velocity) NavierStokes equations.) The thickness δ′ keeps on increasing until steady state 2 is reached (Figure 8d). In summary, the dynamics of the membrane process reflects the characteristic times given in Table 2: (i) t ) 0 s: initial steady state (1). (ii) 0 s < t < 0.001 s (Phase I): At the time of t ) 0.001 s (τ3/6), the dynamics associated with the characteristic time τ3 starts to be visible, reaching 15% of its final value. This time is insufficient for the surface concentration to change (δ0 does not change) but sufficient for the velocity to increase (v0w increases to v′w). (iii) 0.001 s < t < 0.015 s (Phase II): At the time of t ) 0.015 s (which is 2τ3 and τ2/4), the dynamics associated with τ3 is ending (reaching 92% of its final value), while the effect of τ2 starts to be perceptible (reaching 22% of its final value). In this time period, the concentration decreases slightly because of the higher flux through the membrane (δ0 decreases to δ′). (iv) 0.015 s < t < 2 s (Phase III): Over this time, the effects of τ2 and τ1 overlap. Finally, at the time of 2 s, the slowest dynamics associated with τ1 almost ends (reaching 97% of its final value). During this time frame, new salt is transported to the surface by convection, and as it accumulates, the osmotic pressure increases, thus reducing the flux (v′w decreases to vFw as δ′ increases to δF). (v) t > 2 s: final steady state (2).

2.2. Relationship between Permeate Flux and Concentration. The time domain response of the concentration (see Phases I, II, and III in Figure 6) displays characteristics in common with the vw response (see Phases I, II, and III in Figure 7), and it is evident that the TMP effect on the permeate flux can be explained only through changes in the concentration boundary layer. This suggest the need to study the dynamic responses of the boundary layer concentration (whose transfer function is denoted GC(s)) to a TMP disturbance. However, the concentration in the boundary layer is not an integral variable like the permeate flux. It is a local variable since it changes along the membrane surface. But, under certain conditions (which will be investigated later), the Bode diagrams of the dimensionless variational local concentration are qualitatively similar at all points on the membrane surface. Therefore, it is possible to obtain the average concentration dynamics. In Figures 9 and 10, the amplitude ratio (AR) and the phase shift (φ) diagrams for the local concentration at the monitoring point (10 µm from the membrane and 150 mm from the inlet) are reported. The analysis shows that a suitable model for GC(s) ) C(s)/∆P(s) (where C(s) and ∆P(s) are the Laplace transform of c(t) and ∆P(t)) has 2 zeros and 3 poles with τ1 ) 0.5 s, τ2 ) 0.05 s, τ3 ) 0.008 s, τa ) -0.005 s, and τb ) 0.18 s. The most significative results regarding G(s) and GC(s) are compared in Table 3. Although the dynamic models are different (2/3 for the concentration and 3/3 for the permeate flux), the number of poles are the same and the characteristic times are very similar for the two cases. 3. Effect of the Inlet Velocity The analysis of the errors associated with mathematical fitting is very useful because it can give an idea of the range within which the present approximation is valid. In Tables 2 and 3, G(s) and GC(s) were studied over a particular range of TMP and inlet concentration. A third important input parameter is the inlet velocity U0 (or the time of residence L1/U0) In Table 4 the fitting parameters and χ2 error at different inlet velocities are reported for the variable


Figure 10. Angle shift φ (output ) local concentration) for a TMP displaced cosine impulse with duration tω ) 0.001 s (c0 ) 25 g L-1 and TMP1 ) 80 bar). Table 3. Comparison between G(s) and GC(s) transfer function

G(s)

Gc(s)

definition model characteristic times (s) Bode diagrams

J(s)/∆P(s) 3/3 τ1 ) 0.6, τ2 ) 0.06, τ3 ) 0.006 Figures 9 and 10

C(s)/∆P(s) 2/3 τ1 ) 0.5, τ2 ) 0.05, τ3 ) 0.08 Figures 2 and 3

Table 4. Fitting Parameters and χ2 Error at Different U0 (TMP1 ) 80 bar, c0 ) 25 g L-1) U0 (m/s) τ1 (s) τ2 (s) 0.005 0.5 1.5 4.5 15

5.8 0.7 0.6 0.3 0.1

1.22 0.22 0.06 0.03 0.02

τ3 (s)

τa (s)

-0.021 0.018 0.006 0.0043 0.0034

5.9 0.9 1.3 0.6 0.3

τb (s)

τc (s)

χ2

-0.04 -0.567 0.705 0.38 0.030 1.007 0.09 0.007 0.056 0.05 0.005 0.007 0.02 0.004 0.006

permeate flux. The χ2 error increases significantly at low velocities (U0 e 0.5 m s- 1), and the resulting characteristic times lose their physical meaning. For instance, τ3 at U0 ) 0.005 m s- 1 is negative. Consequently, when velocities are low, the approach described in the previous section cannot be used unless a very high-order transfer function is employed. The error increment can be explained through the fact that permeate flux is the integral of the local fluid velocities through the membrane (Figure 11). In theory, the membrane process is a distributed parameter system, and its dynamics is the result of all the local vw dynamics computed at each point along the membrane surface. The models proposed here, however, are of the lumped-parameter type, meaning that any dependent variable can be assumed to be a function only of time and not of the spatial position. The local dynamics of the velocity at different locations are similar when U0 is high (Figure 12) but very different when U0 is low (Figure 13), as confirmed by the χ2 error in Table 4. In Figure 12, the time constants only change slightly from point to point, but the shape remains the same. Therefore, a model representing the average behavior can be obtained. When U0 is low, however, the local responses change significantly and display completely different shapes (as shown in Figure 13). The presence of multiple peaks, particularly at distances of 24 and 25 cm from

Figure 11. Schematic representation of the local velocities and concentrations on the membrane.

the inlet, is a clue to the more complicated nature of the system at low inlet velocities. In mathematical terms, this can be explained by the fact that the permeate flux is the integral of vw through the membrane surface (S),

J(t) )

∫∫S vw(t,x,y) dx dy S

(5)

and, in terms of transfer functions,

G(s) )

J(s) ∆P(s)

)

∫∫S

VW(s) dx dy S∆P(s)

∞

≈

Gi(s)∆S ∑ i)0

(6)

where ∆S ) dx dy, VW is the Laplace transform of vw, and Gi is the ith local transfer function. If all the Gi are similar along the membrane (as in the case of high U0), then it is possible to use an overall transfer function for the entire process. When the Gi are very different,

7830


membrane (see Figure 16). Except for the section near the entrance, most of the boundary layer lies in the fully developed region and thus, for a membrane with the same or greater length than the one reported in this paper, can be approximated by a constant value. Therefore, the integration over the length of membrane can be approximated by a linear equation. 5. Feedback Model

Figure 12. v/w local AR-Bode diagram at different points on the membrane surface (U0 ) 1.5 m/s, TMP1 ) 80 bar, and c0 ) 25 g L-1).

the overall dynamics is very complicated and, thus, cannot be represented using simple lumped-parameter models. 4. Linearity and Steady States In the paper, we used the term “linearity” with two different meanings. In the first case (which considers the dynamic behavior), the term linearity is used in the sense that the resulting ordinary differential equations are linear (which usually implies exponential-like solutions as functions of time). In ref 4 and, consequently, in eq 4 of this paper, the model used to express the permeation flux is nonlinear because the osmotic pressure depends on the concentration, which, in turn, depends on the flow field obtained from the NavierStokes equations. One of the key results of the paper is to find the conditions where the complex relationship implied in the cited equation (which usually leads to nonlinear partial differential equations) can be simplified to linear ordinary differential equations. We found that, within the range investigated, the behavior of the above nonlinear system can be sufficiently approximated by linear dynamical systems, as evident from the similar steady-state gains and time constants when ∆P varies in the range of 5-30 bar. In the second case, the term linearity is used in the sense that the relationship between TMP and J (at steady state) can be approximated by a linear static equation. The importance of this statement lies in the fact that the transfer function must be independent of the perturbation amplitude (at least within a certain range), and this can be assured if TMP and J have a linear dependence over a certain range. The CFD results show that this assumption is true for 60 bar < TMP < 100 bar and 25 g L-1 < c0 < 35 g L-1 (see Figure 14). The reason for the system’s almost linear behavior in the above range can be found, again, from eq 4. The permeate flux is the integral of vw over the entire membrane (eq 5). This equation can be approximated by a linear equation only if the following two circumstances occur. First, the local concentration must be almost linear with the pressure. Second, the value of the integral over the membrane must preserve this linearity. The first point is illustrated in Figure 15 computed from the CFD results. The second point depends on the shape of the boundary layer along the

A further investigation can be performed to study the relationship between the transfer function G(s) ) J(s)/ ∆P(s), which associates the transmembrane pressure input to the permeate flux output, and the transfer function GC(s) ) C(s)/∆P(s), which associates a transmembrane pressure input to the (average) concentration C(s). In Table 3, we have shown that G(s) and GC(s) have different dynamics (3/3 and 2/3) but the same number of poles and similar characteristic times. Was this chance or is there a theoretical explanation behind this similarity? In other words, does an overall model exist, which includes both G(s) and GC(s)? The answer is yes, and we show here that this model has the features of a feedback model. The block diagram of this model is shown in Figure 17, where the interactions between the variables ∆P, J, vw (which, for simplicity, is denoted as v from now on), c (in the boundary layer), and ∆π are illustrated. This feedback model is developed on the basis of the previous explanation of the dynamics (as shown in Figure 8): ∆P affects the velocity v, which affects the local concentration c, which in turn again affects v. It is worth noting that the feedback model is not just a trivial consequence of eq 4. In eq 4, in fact, the local transmembrane velocity and osmotic pressure are interconnected through the Navier-Stokes equations and depend on the momentum balance of the entire channel. Dividing the membrane into N parts and noting that the applied TMP is the same at all points along the membrane, the local variables and transfer functions of Figure 17 can be described through the following vectors, v ) vi, c ) ci, and gv ) Gvi, and the matrix Gc-v ) Gci-vj (i ) 1, ..., N and j ) 1, ..., N). The permeate flux J and the integral osmotic pressure ∆Π are expressed respectively by

J)

1

N

∑ vi

Ni)1

(7)

and

∆Π )

A1 N

∑ci

N i)1

(8)

(Here it is assumed, for simplicity, that ∆Π ) A1c, but the following discussion is still valid for more complex formulations such as ∆Π ) A1c + A2c2 + A3c3). According to Figure 17 and eq 8,

v ) gv(∆P - A1c) ) gv(∆P - A1Gc-vv)

(9)

By combining eqs 7 and 9,

J)

A1 ∆P gvi - gv[Gc-vv] N N

where i is a row vector of ones [1, 1, ..., 1].

(10)


Figure 13. v/w local AR-Bode diagram at different points on the membrane surface (U0 ) 0.5 m/s, TMP1 ) 80 bar, and c0 ) 25 g L-1).

Figure 14. Steady-state results.

Figure 15. Dimensionless concentration versus TMP at a point 10 µm from the membrane and 150 mm from the inlet.

At high velocities (as shown in the previous section), it is possible to approximate local values with a global average and consequently v ≈ iv, c ≈ ic, gv ≈ iGv, and Gc-v ≈ IGc-v, where I is the identity matrix; eq 10 becomes

J ) ∆PSGv - A1JGvGc-v

(11)

and, finally, the overall transfer function G ) J∆P in Figure 17 can be expressed by

G)

∆PSGv 1 + RTGvGc-v

(12)

Equation 12 is the mathematical form of a simple feedback model, as illustrated in Figure 17 with G1(s) ) ∆PSGv(s) and G2(s) ) A1Gc-v(s)/∆PS. Consequently, the global system shown in Figure 17 can be simplified to the simpler feedback model of Figure 18. The transfer function G1 can be found by opening the feedback loop. This can be achieved by setting the osmotic pressure to zero in the CFD computations. This treatment, obviously, is not physically feasible. However, mathematically, it allows us to express G1 as J/∆P and perform the usual pulse test. This is another advantage of using the CFD approach for analysis, because in experiments the osmotic pressure cannot be eliminated. The results were shown in Figures 2 and 3. According to these

7832


Figure 16. Dimensionless concentration along the membrane at a point 10 µm from the membrane surface.

Simplifying eq 15 gives

(τ′as + 1)(τ′bs + 1)(τ′cs + 1) G(s) ) K2 (τ1s + 1)(τ2s + 1)(τ3s + 1)

Figure 17. Complete feedback model block-scheme of global and local transfer function.

(16)

where K2, τ′a, τ′b, and τ′c, are functions of τa, τb, τ1, τ2, and τ3 (K2 also contains the factors β and A1). What needs to be stressed is that both GC(s) and G(s) have the same denominator. Previously, it was found empirically that the transfer functions G(s) and GC(s) have the same number of poles and very similar characteristic times (see Table 3). Now, using the feedback model, we arrive at the mathematical conclusion that they must have exactly the same poles. The small differences found before can be explained by the approximation made in averaging c(x,y) and vw(x,y) at different locations. 6. Conclusions

Figure 18. Simplified feedback model block-scheme of Figure 17.

results, it is possible to assume G1(s) ) β, where β is a constant. This leads to

G(s) )

β 1 + βG2(s)

(13)

The final transfer function to be determined is G2. Previously, the transfer function GC(s) ) C(s)/∆P(s) was investigated (see Table 3). The two functions G2(s) and GC(s) are related because G2(s) ) A1C(s)/J(s). Consequently, G2(s) ) A1GC(s)/G(s). By substituting this equation into eq 13 we obtain

G(s) ) β[1 - RTGC(s)]

(14)

Is eq 14 consistent with the results previously found? It was found that G(s) had 3 poles and 3 zeros and GC(s) had 3 poles and 2 zeros. By substituting GC(s) in eq 14 with a 2-zero/3-pole model, we obtain

[

]

(τas + 1)(τbs + 1) G(s) ) β 1 - A1 (τ1s + 1)(τ2s + 1)(τ3s + 1)

(15)

In RO membranes, the difference in concentration between the permeate and the retentate side produces an osmotic pressure difference, which opposes the TMP and reduces the permeate flow. Flow disturbances, however, can maintain high membrane performance, because they interfere with the boundary layer where the salt is accumulated. Such methods are already used in industrial applications, but their implementation is mostly based on empirical experience rather than systematic study. In this work, the dynamic responses of an RO membrane system to time-dependent flow disturbances have been studied by using transfer functions. The first transfer function we have investigated is G(s), which relates the TMP with the permeate flux. This is the most logical choice, because flow disturbances are often generated by transmembrane pressure changes and the permeate flux can be used as a measure of the membrane performance. The study of the physical characteristics connected with G(s), however, highlights the importance of the concentration boundary layer. For this reason, a second transfer function GC(s), which relates the TMP with the concentration at the membrane surface, was investigated. The two transfer functions G(s) and GC(s) show some common characteristics, which in turn seem to point to the existence of an overall model. A feedback model is thus proposed in order to clarify the connection between transmembrane pressure, boundary layer concentration, and permeate flux. This model leads to a better understanding of the fundamental membrane dynamics and, in particular, of the role of the boundary layer. This feedback model may be useful for designing better dynamic operations and,


in particular, flow disturbances aimed at improving the membrane performance through the control of the concentration polarization.

the input pulse duration Tx. Then, the frequency response of the process can be calculated as follows:

∫0T y(t) e-iωt dt G(iω) ) T ∫0 x(t) e-ωt dt y

Acknowledgment

x

The authors gratefully acknowledge the financial support of the Australian Research Council (Discovery Project No. DP0343073). Appendix: Transfer Functions, Pulse Tests, and Bode Diagrams In this paper, transfer functions are used to model the dynamics of membrane systems. A transfer function describes the relationship between the input (the known or assumed conditions) and output (the unknown) variables in the Laplace domain. If a process has an input of x(t) and an output of y(t), the process dynamics can be represented by the following equation

G(s) )

Y(s) X(s)

(17)

where Y(s) and X(s) are the Laplace transform of the output and the input, respectively. The advantage of this approach is that the resulting equations are independent of the initial conditions or the type of input functions (see ref 10). With the transfer function, the dynamics of the output in response to any dynamic input can be obtained by multiplying the transfer function by the Laplace transform of the input. Generally speaking, system identification methods based on time series analysis are preferred for obtaining process models from input/output data with uniform time sampling periods (in this case, a discrete time transfer function will be developed). However, in this paper, nonuniform time steps were used in our CFD simulation to avoid an excessive amount of data being generated. Therefore, a different method based on frequency response analysis is adopted to obtain the process transfer function. The dynamic behavior of the membrane system can be studied using frequency response analysis. According to control theory, the steady-state response of a linear system to a sinusoidal input x(t) ) A sin(ωt) would be of the same frequency but with different magnitude and phase shift, given by y(t) ) A′ sin(ωt + φ). The amplitude ratio (AR ) A′/A) and phase shift (φ) can be computed at different frequencies and graphically represented by Bode diagrams. The frequency response (denoted as G(iω) ) AR(ω)∠φ(ω)) tells how the process would respond to input variationsshow large and how fast the response is. The frequency response of a process can be obtained from input and output data in the time domain, which are, in turn, acquired from experiments or simulation studies. One possible approach is the “pulse test”, where a pulse input x(t) is used and the corresponding process output y(t) is collected. The process transfer function is given as follows:

Y(s) ) G(s) ) X(s)

∫0∞y(t) e-st dt ∫0∞x(t) e-st dt

A pulse input contains many frequency components and, thus, excites the system at different frequencies. Therefore, it is possible to obtain the frequency response at different frequencies from a single pulse test. Not all the pulse shapes have the same frequency spectra, and certain shapes are more useful than others. Suitable pulse shapes to use for input forcing are those that have a high magnitude over a wide range of frequencies. Hougen11 analyzed a number of pulse shapes and found that the displaced cosine pulse

x(t) ) 1 - cos(2πt/tω), 0 < t e tω

Because y(t) and x(t) are deviation variables, the numerator and denominator integrals in eq 18 need only to be evaluated over the output pulse duration Ty and

(20)

has a more favorable frequency spectra compared with the rectangular pulse. The displaced cosine pulse has an effective frequency range up to ωtω ) 8, while the frequency range for the rectangular pulse is only up to ωtω ) 5. A specific code was written and used to integrate eq 19 to obtain the frequency response G(iω) ) AR(ω)∠φ(ω) at different frequencies ω. The results can then be represented by Bode diagrams. From the slopes of the AR curve at different frequencies and the corresponding phase shift, we can determine the number of poles and zeros that should be included in the transfer function and their approximate locations (details can be found in ref 10). This provides the structure of the process transfer function as well as the initial values of the time constants in the denominator and numerator. For example, if the Bode plot suggests 3 poles and 3 zeros in the transfer function, we will have the following estimated transfer function with unknown parameters:

(τ′as + 1)(τ′bs + 1)(τ′cs + 1) G ˆ (s) ) K (τ1s + 1)(τ2s + 1)(τ3s + 1)

(21)

A nonlinear optimization can be performed such that the estimated transfer function gives the frequency response at different frequencies ω as close as possible to that given in eq 19 (obtained from the CFD simulation):

min

K,τ1,τ2,τ3,τa,τb,τc

|G(iω) - G ˆ (iω)|

(22)

In some cases, because of the complexity of the dynamic process, the number of poles and zeros cannot be easily estimated from the Bode plots. An iterative optimization procedure may be used, which starts with a simple transfer function structure. More poles and zeros are added, and the optimization in eq 22 is repeated until the model error cannot be further significantly reduced. Symbols

(18)

(19)

A1 ) osmotic pressure constant, bar m3 kg-1 AR) amplitude ratio c ) salt concentration, kg m-3 K ) membrane permeation coefficient, m s-1 bar-1 R ) membrane rejection coefficient P ) total pressure (TMP + ∆P), bar

7834


t ) time, s S ) membrane surface, m2 TMP ) transmembrane pressure, bar tω ) impulse duration, s U0 ) lean inlet velocity, m s-1 u ) local fluid velocity in the x-direction, m s-1 vw ) local fluid velocity through the membrane, m s-1 ∠ ) phase shift Greek Letters ∆P ) transmembrane pressure variation, bar ∆π ) local osmotic pressure, bar ∆Π ) integral osmotic pressure, bar J ) permeate flux, m s-1 φ ) angle shift, rad χ2 ) sum of square error ω ) frequency, s-1 Subscripts

(3) Wenten, I. G.; Koenhen, D. M.; Roesnik, D. W.; Rasmussen, A.; Jonsson, G. The Blackshock ProcesssA Novel Backflush Technique in Microfiltration. Presented at Process Engineering Membrane Process II, Environmental Application, Il Ciocco, Tuscany, Italy, 1994. (4) Fletcher, D. F.; Wiley, D. E. A computational fluid dynamics study of buoyancy effects in reverse osmosis. J. Membr. Sci. 2004, 245, 175. (5) Brian, P. L. T. Concentration polarization in reverse osmosis desalination with variable flux and incomplete salt rejection. Ind. Eng. Chem. Fundam. 1965, 4, 439. (6) Gill, W. N.; Tien C.; Zeh, D. W. Concentration polarization effects in reverse osmosis system. Ind. Eng. Chem. Fundam. 1965, 4, 433. (7) Sherwood, T. K.; Brian P. L. T.; Fisher, D. R. Salt concentration and phase boundaries in desalination by reverse osmosis. Ind. Eng. Chem. Fundam. 1965, 4, 113. (8) Wiley, D. E.; Fletcher, D. F. Computational fluid dynamics modelling of flow and permeation for pressure-driven membrane processes. Desalination 2002, 145, 183.

0 ) inlet 1 ) steady state 1 2 ) steady state 2

(9) Ho, W. S. W.; Sirkar, K. K. Membrane Handbook; Van Nostrand Reinhold: New York, 1992.

Superscripts

(10) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control; Wiley: New York, 1989.

* ) variational variable defined in eq 1

Literature Cited (1) Li, H.; Bertram, C. B.; Wiley D. E. Mechanisms by which pulsatile flow affects cross-flow microfiltration. AIChE J. 1998, 44, 1950. (2) Madaeni, S. S.; Fane, A. G.; Wiley, D. E. Factors Influencing Critical Flux of Membrane Filtration of Activated Sludge. J. Chem. Technol. Biotechnol. 1999, 74, 539.

(11) Hougen, J. O. Experiences and experiments with process dynamics. Chem. Eng. Prog., Monogr. Ser. 1964, 60, 4. (Cited in ref 10.)

Received for review March 1, 2005 Revised manuscript received August 8, 2005 Accepted August 9, 2005 IE050290Y

Analysis of the Dynamic Response of a Reverse Osmosis Membrane

Recommend Documents