# Mathematical Modeling and Simulation of SWRO Process Based on Simultaneous Method

Authors: Aipeng Jiang; Qiang Ding; Jian Wang; Shu Jiangzhou; Wen Cheng; Changxin Xing
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Reverse osmosis (RO) technique is one of the most efficient ways for seawater desalination to solve the shortage of freshwater. For prediction and analysis of the performance of seawater reverse osmosis (SWRO) process, an accurate and detailed model based on the solution-diffusion and mass transfer theory is established. Since the accurate formulation of the model includes many differential equations and strong nonlinear equations (differential and algebraic equations, DAEs), to solve the problem efficiently, the simultaneous method through orthogonal collocation on finite elements and large scale solver were used to obtain the solutions. The model was fully discretized into NLP (nonlinear programming) with large scale variables and equations, and then the NLP was solved by large scale solver of IPOPT. Validation of the formulated model and solution method is verified by case study on a SWRO plant. Then simulation and analysis are carried out to demonstrate the performance of reverse osmosis process; operational conditions such as feed pressure and feed flow rate as well as feed temperature are also analyzed. This work is of significant meaning for the detailed understanding of RO process and future energy saving through operational optimization.


---

## Body

## 1. Introduction

China, which has the largest population in the world and is called the world’s factory, has a much more serious problem of freshwater shortage. This problem is expected to worsen in the whole world with the growth of industrialization, as well as climate change [1–3]. Seawater desalination based on RO (reverse osmosis) membrane technique is one of the most promising ways to obtain freshwater, especially in the coastal regions and islands. More than 40% market shares of seawater desalination based on RO technique was occupied worldwide in the past decade [4]. Reverse osmosis (RO) is also a membrane technology driven by high pressure and is used to separate salt and water of the same order of molecular size. High pressure is applied on the feed side of the membrane to overcome the osmotic pressure and cause transport of the solvent from feed to permeate side. Among the membrane modules, spiral-wound RO occupies the largest market share because of its relative ease of cleaning, fabrication technology, and very large surface area per unit volume [5, 6]. The general process of SWRO can be found in the specialized literature [7, 8]. To achieve a high recovery in a single stage, several spiral-wound elements in series, held in a cylindrical pressure vessel, are usually used in seawater desalination applications. Generally a reverse osmosis desalination unit is composed of dozens of cylindrical pressure vessels, and several units with a permeate water storage tank were used as the key part of a seawater desalination plant.The feed pressure as well as other parameters such as operational temperature and properties of the membrane has significant effect on the performance of SWRO process. It is quite important to understand the mechanism of SWRO process and improve the performance and reduce the energy consumption. For this reason, based on solution-diffusion model and film theory, the model of RO process was researched, found, analyzed, and used to optimize the operation of SWRO plant [9, 10].Libotean and Khayet established the performance prediction model based on data-driven method such as artificial neural networks [11, 12]. Since modeling based on rigid first principle can reveal more information, Senthilmurugan et al. [9] applied the solution-diffusion model modified with the concentration polarization theory for analyzing the operation and performance of RO systems. Abbas analyzed and simulated an industrial medium-scale brackish water reverse osmosis plant with semirigorous model, and then he used it to analyze the optimal operation of RO system [13]. Geraldes considered both the investment cost and operation cost of SWRO system and provided the differential equations of the RO process with distributed method [14]. Then, to avoid computing difficulties, he simplified the simulation process and solved his optimal problem with simple difference equations. Sundaramoorthy et al. build more complex analytical model of spiral wound reverse osmosis modules, based on complex model parameters estimation; the model is computed based on sequential method [15]. To get more accurate results, more modeling and simulation of reverse osmosis system are carried out with CFD technique [16–18]. But since these models are not founded with simultaneous method, computing efficiency for simulation and optimization cannot be guaranteed for online application.In this paper, to achieve accurate prediction of RO process performance and to accelerate the computing efficiency for real time simulation, the mathematical modeling is carried out based on rigid first principle, then simultaneous method within which the differential variables are fully discretized by finite element collocation. The mathematical model is verified and then simulation under different conditions will be discussed.


## 2. Mathematical Modeling of Spiral-Wound RO Process

The performance of SWRO plants is quite sensitive to the quality of the feed water and plant operating conditions. Based on solution-diffusion model and film theory and according to the schematic diagram of the SWRO process shown in Figure2, more accurate membrane transport equations of steady-state in the form of distributed parameters can be derived.For the RO membrane, the overall fluid and solute mass balance equations are(1)
Q
p
=
Q
f
-
Q
r
Q
f
C
f
=
Q
r
C
r
+
Q
p
C
p
Q
p
=
n
l
W
∫
0
L
J
v
d
z
.Here subscriptsf, r, and p refer to the feed, reject (brine), and permeate (product) streams. Q and C refer to the flow rate and salt concentration, respectively. The local water flux and salt flux can be calculated from the Kimura-Sourirajan solution-diffusion mass transport relations. n
l, W, and L express the number of leaves and the width and length of the RO module, respectively. The solution-diffusion model is assumed to be valid for the transport of solvent and solute through the membrane. According to this model, solvent flux J
v and solute flux J
s through membrane are expressed by the following equations:
(2)
J
v
=
A
w
(
P
f
-
P
d
-
P
p
-
Δ
π
)
J
s
=
B
s
(
C
m
-
C
b
)
.
Let
(3)
P
b
=
P
f
-
P
d
,
Δ
P
=
(
P
b
-
P
p
)
.
Then
(4)
J
v
=
A
w
(
Δ
P
-
Δ
π
)
,
where A
w is the solvent transport parameter, P
f is the feed pressure, P
d is the pressure drop along a RO spiral-wound module, P
b is the pressure along the channel of spiral-wound module, and P
p is the pressure of permeate side, which in general is assumed as environment pressure. Δ
π is the pressure loss of osmosis pressure. B
s is the solute transport parameter, C
m and C
s
p are solute concentration at the membrane surface on the feed side and solute concentration on the permeate side, respectively, and C
p is the value of C
b at the end of module. That means C
p
=
C
b
(
L
). A
w and B
s are sensitive with operational temperature and related factors; the relationship is expressed as follows:
(5)
A
w
=
A
w
0
exp
⁡
(
α
1
T
-
273
273
-
α
2
(
P
f
-
P
d
)
)
B
s
=
B
s
0
exp
⁡
(
β
1
T
-
273
273
)
,
where A
w
0 and B
s
0 are intrinsic transport parameter in standard condition and α
1, α
2, and β
1 are constant parameters for transport. T represents operational temperature with unit of kelvin degree. The osmotic pressure is nearly linearly related to concentration by the equation
(6)
Δ
π
=
R
T
(
C
m
-
C
p
)
,
where R is the gas law constant.Solution of the above equations requires knowledge of the RO process specification and parameters as well as the solute concentration ofC
m at the membrane wall, which is quite different from the bulk concentration C
b due to the CP (concentration polarization) phenomenon. Through steady-state material balance around the boundary layer and CP theory, the following simple expression is developed:
(7)
ϕ
=
C
m
-
C
p
C
b
-
C
p
=
exp
⁡
(
J
v
k
c
)
.The bulk concentrationC
b and solvent transports J
v vary along the membrane channel. The computation of the mass transfer coefficient k
c is as follows [19]:
(8)
Sh
=
k
c
d
e
D
A
B
=
0.065
Re
0.875
Sc
0.25
,
where
(9)
Re
=
ρ
V
d
e
μ
Sc
=
μ
(
ρ
D
A
B
)
.ρ is the density of permeate water, d
e is hydraulic diameter of the feed spacer channel, μ is kinematic viscosity, and D
A
B is dynamic viscosity; it can be calculated by the regression equation
(10)
D
A
B
=
6.725
×
1
0
-
6
·
exp
⁡
(
0.1546
×
1
0
-
3
C
b
-
2513
273.15
+
T
)
.
The relationship between J
v and J
s is
(11)
J
s
=
J
v
C
p
.
According to Figure 2, the pressure loss along the RO channel can be formulated as
(12)
d
P
d
d
z
=
-
λ
ρ
d
e
V
2
2
,
where
(13)
λ
=
6.23
K
λ
Re
-
0.3
.
K
λ is the empirical parameter. Since the pressure along RO P
b
=
P
f
-
P
d, so
(14)
d
P
b
d
z
=
-
d
P
d
d
z
=
λ
ρ
d
e
V
2
2
,
at z
=
0, P
b
=
P
f, at z
=
L, P
b
=
P
r.V is axial velocity in feed channel and satisfies
(15)
d
V
d
z
=
-
2
J
v
h
sp
at
z
=
0
,
V
=
V
f
=
Q
f
n
e
W
h
sp
at
z
=
L
,
V
=
V
r
=
Q
r
n
e
W
h
sp
.
h
sp is height of the feed spacer channel.The bulk concentrationC
b varies along the membrane channel and can be given as
(16)
d
C
b
d
z
=
2
J
v
h
sp
V
(
C
b
-
C
p
)
,
and at z
=
0, C
b
=
C
f; at z
=
L, C
b
=
C
r.From the solution of the above equations,Q
p and C
p at given operational conditions and specification of membrane can be obtained, from which the water recovery rate Rec and specific energy consumption SEC can be calculated by the equations
(17)
Rec
=
Q
p
Q
f
(18)
SEC
=
(
P
f
Q
f
/
ε
p
-
P
r
Q
r
ε
e
f
)
Q
p
.Salt passage and salt rejection coefficient are also two important parameters reflecting the performance of the RO process. They are formulated as the following equations:(19)
Sp
=
C
p
C
f
×
100
%
(20)
Ry
=
(
C
f
-
C
p
)
C
f
×
100
%
,
where ε
p and ε
p
f are the mechanical efficiency and energy recovery efficiency, respectively.Equations (1)–(20) explain the detailed reverse osmosis process with spiral-wound membrane module. If the operational variables, such as feed flow rate and feed pressure, are fixed, the permeate quality and related performance can be obtained by the solution of the equations. Otherwise, the operational variables and related parameters can be calculated by (1)–(20) under certain constraints. The constraints listed in (21) include requirement of equipment safety, water quality, and parameter of CP:
(21)
P
f
l
o
≤
P
f
≤
P
f
u
p
Q
f
l
o
≤
Q
f
≤
Q
f
u
p
C
p
l
o
≤
C
p
≤
C
p
u
p
T
l
o
≤
T
≤
T
u
p
ϕ
≤
ϕ
max
⁡
H
t
min
⁡
≤
H
t
≤
H
t
max
⁡
.


## 3. Discretization and the Solution Method

Since (1)–(20) include a set of strong nonlinear equations and even DAEs (differential and algebraic equations) and since there are more constraints than those of the equation mentioned above, it is fairly tough to solve this kind of problem efficiently with high accuracy [20–22]. Here the DAEs (1)–(21) are converted into the form of NLP (nonlinear programming) through simultaneous approach by approximating state and control profiles with a family of polynomials on finite elements (shown as Figure 3) [20, 23]. These polynomials can be represented as power series, such as orthogonal polynomials, or in Lagrange form. Here the following monomial basis representation for the differential profile was selected, which is popular for Runge-Kutta discretizations [23]:
(22)
z
(
t
)
=
z
i
-
1
+
h
i
∑
q
=
1
K
Ω
q
(
t
-
t
i
-
1
h
i
)
d
z
d
t
i
,
q
.Herez
i
-
1 is the value of the differential variable at the beginning of element i, h
i is the length of element i, d
z
/
d
t
i
,
q denotes the value of its first derivative in element i at the collocation point q, and Ω
q is a polynomial of order K, satisfying
(23)
Ω
q
(
0
)
=
0
,
q
=
1
,
…
,
K
Ω
q
′
(
ρ
r
)
=
δ
q
,
r
,
q
=
1
,
…
,
K
,
where ρ
r is the location of the rth collocation point within each element. Continuity of the differential profile is enforced by
(24)
z
i
=
z
i
-
1
+
h
i
∑
q
=
1
K
Ω
q
(
1
)
d
z
d
t
i
,
q
.Based on a number of studies of Larry group [20, 24–26], Radau collocation points were selected in this study because they allow constraints to be set at the end of each element and to stabilize the system more efficiently if high index DAEs are present. In addition, the control and algebraic profile are approximated using a Lagrange basis representation which takes the following form:
(25)
y
(
t
)
=
∑
q
=
1
K
ψ
q
(
t
-
t
i
-
1
h
i
)
y
i
,
q
u
(
t
)
=
∑
q
=
1
K
ψ
q
(
t
-
t
i
-
1
h
i
)
u
i
,
q
.
Here y
i
,
q and u
i
,
q represent the values of the algebraic and control variables, respectively, in element i at collocation point q, t is the value satisfying t
i
-
1
⩽
t
⩽
t
i, and Ψ
q is the Lagrange polynomial of degree K satisfying
(26)
ψ
q
(
t
)
=
∏
r
=
1
K
(
t
-
t
r
)
(
t
q
-
t
r
)
=
∏
r
=
1
K
(
t
-
t
r
)
/
h
i
(
t
q
-
t
r
)
/
h
i
=
∏
r
=
1
K
(
t
-
t
r
)
/
h
i
(
ρ
q
-
ρ
r
)
ψ
q
(
ρ
r
)
=
δ
q
,
r
,
q
,
r
=
1
,
…
,
K
.
And, at collocation points,
(27)
ψ
q
(
t
r
)
=
δ
q
,
r
=
{
1
,
q
=
r
0
,
q
≠
r
(
q
,
r
=
1
,
…
,
K
)
.
The differential variables are required to be continuous throughout the time (or space) horizon, while the control and algebraic variables are allowed to have discontinuities at the boundaries of the elements.Through the simultaneous method DAEs of (1)–(21) were transformed into NLP problem. Some state variables in (1)–(21) can be discretized as flow variables: P
d
i
,
j, V
i
,
j, C
b
i
,
j, Z
i
,
j, C
b
i
,
0, V
i
,
0, P
d
i
,
0,  Z
i
,
0, and J
v
i
,
j, J
s
i
,
j, λ
i
,
j
,
Sc
i
,
j, Re
i
,
j, k
c
i
,
j, C
m
i
,
j, Δ
π
i
,
j. Operational variables Q
f, P
f and other environmental variables such as T
f and C
f are not discretized because they are operational variables and feed condition variables. The equations are expressed in the following form:
(28)
A
w
(
i
,
j
)
=
A
w
0
exp
⁡
⁡
(
α
1
T
-
273
273
-
α
2
(
P
f
-
P
d
(
i
,
j
)
)
)
B
s
=
B
s
0
exp
⁡
(
α
2
T
-
273
273
)
J
v
i
,
j
=
A
w
(
P
f
-
P
d
i
,
j
-
Δ
π
i
,
j
)
J
s
i
,
j
=
B
w
(
C
m
i
,
j
-
C
p
i
,
j
)
C
p
i
,
j
=
J
s
i
,
j
J
v
i
,
j
Δ
π
i
,
j
=
R
T
(
C
m
i
,
j
-
C
p
i
,
j
)
C
m
i
,
j
=
[
(
C
b
i
,
j
-
C
p
i
,
j
)
exp
⁡
(
J
v
i
,
j
k
c
i
,
j
)
]
+
C
p
i
,
j
k
c
i
,
j
=
0.065
Re
i
,
j
0.875
Sc
i
,
j
0.25
d
e
D
A
B
Sc
i
,
j
=
μ
i
,
j
(
ρ
D
A
B
)
λ
i
,
j
=
6.23
K
λ
Re
i
,
j
-
0.3
.
ODE formulations are
(29)
Re
i
,
j
=
ρ
V
i
,
j
d
e
μ
i
,
j
d
V
i
,
j
d
x
=
-
2
J
v
i
,
j
h
s
p
d
C
b
i
,
j
d
x
=
2
J
v
i
,
j
(
C
b
i
,
j
-
C
p
i
,
j
)
(
h
s
p
×
V
i
,
j
)
d
P
d
i
,
j
d
x
=
λ
i
,
j
ρ
d
e
×
(
V
i
,
j
2
)
2
.
Collocation equations are
(30)
C
b
i
,
j
=
C
b
i
,
0
+
h
i
∑
k
=
1
n
c
p
Ω
k
(
ρ
j
)
d
C
b
d
x
i
,
k
V
i
,
j
=
V
b
i
,
0
+
h
i
∑
k
=
1
n
c
p
Ω
k
(
ρ
j
)
d
V
d
x
i
,
j
P
d
i
,
j
=
P
d
i
,
0
+
h
i
∑
k
=
1
n
c
p
Ω
k
(
ρ
j
)
d
P
d
d
x
i
,
j
.
Continuity equations are
(31)
C
b
i
,
0
=
C
b
i
-
1,0
+
h
i
-
1
∑
k
=
1
n
c
p
Ω
k
(
1
)
d
C
b
d
x
i
-
1
,
k
V
i
,
0
=
V
i
-
1,0
+
h
i
-
1
∑
k
=
1
n
c
p
Ω
k
(
1
)
d
V
d
d
x
i
-
1
,
k
P
d
i
,
0
=
P
d
i
-
1,0
+
h
i
-
1
∑
k
=
1
n
c
p
Ω
k
(
1
)
d
P
d
d
x
i
-
1
,
k
.
Initial conditions are
(32)
C
b
1,0
=
C
f
V
1,0
=
V
f
P
d
1,0
=
0
.
Here 0
≤
x
i
j
≤
L,  i
=
1
⋯
100
j
=
1
⋯
K.Figure3 gives the detailed discretization process and distribution of discretized variables. Since the NLP problem is formulated with large scale equations and variables through discretization of the original model, it can be solved by large scale NLP solvers such as IPOPT. If all the operational variables (feed pressure and feed flow rate) are fixed, the solution of the problem results in the simulation of the RO process; profiles of key performance parameters as well as detailed status variables can be obtained, which is helpful for us to understand the mechanism and the whole process.


## 4. Model Simulation and Analysis Based on Seawater Desalination Plant

Based on the mathematical model and collocation method on finite element, the solution of the problem will yield the profiles of key performance parameters as well as the distribution of state variables. To verify the validity of the established model, data of a real seawater desalination plant is studied and then used for the simulation and analysis. The main process of the plant is the same as Figure1. SW30HR-380 membrane modules were selected in series in the pressure vessel, and 55 pressure vessels were used to form two RO process units. Some conditional parameters are listed in Table 1.Table 1
Feed condition of the seawater plant.
Feed concentration (kg/m3)
30
Feed temperature (°C)
25
Feed pressure (Bar)
59
Feed PH
5–8
Densityρ (kg/m3)
1000 (at 0°C)
Kinematic viscosity (kg/m/s)
1.02
e
-
6 (at 0°C)Detailed schematic diagram of RO process and water storage process.
(a)
A schematic representation of the process diagram of desalination
(b)
Detailed SWRO process and storage tank information for the desalination plantFigure 2
Scheme of the rectangular channel model of spiral wound module feed channel.Figure 3
Discretization by collocation on finite element.Field data of a RO process are used for model verification and parameter identification. The key performance parameters such as water recovery rate and salt rejection are compared to indicate the accuracy of the model. Comparison of the field data with those obtained from established models and ROSA9.0 is listed in Table2. ROSA is a powerful software package developed by Dow Chemical Company to help design the RO process system [27]. As can be seen from Table 2, the overall results obtained from the established model are in good agreement with those obtained from both ROSA and the field data. Though the water recovery ratio and salt reject are lower than our model and that of ROSA9.0, the relative error is quite small. So the established model can be well used for the simulation and analysis of the seawater RO process. Since the established model was affected by feed temperature obviously, the properties and correction factor are also compared with the real data or the data from literature [13].Table 2
Comparison of the actual performance with those from ROSA9.0 and our model.
Item
Rec (%)
Ry (%)
Field data
41.6
99.37
ROSA9.0
42.7
99.58
Model
42.1
99.52Figures4, 5, and 6 are the comparison results between model value and real value; it can be seen that these property parameters agree with the real data with relatively high accuracy. The RO process model is fairly robust at the temperature range from 0°C to 60°C.Figure 4
Value of seawater viscosity along temperature.Figure 5
Comparison of DAB with literature.Figure 6
Comparison of TCF curve along temperature.Based on the collocation method on finite element, simulation results are yielded through the solution of the discrete model. Simulation results are shown in Figures7, 8, 9, 10, 11, 12, 13, 14, 15, and 16. From these figures it can be seen that the velocity of solvent, permeate flux, and concentration polarization in modules decrease nonlinearly along the length of membrane. The reason is that the feed flow is separated as two parts, and the pressure difference between the two sides decreases quickly. The osmosis pressure increase is the main factor for the decrease of pressure difference, which is shown in Figure 9. It also can be seen that the decrease of CP along the length of membrane is relatively small, and the rejection coefficient is higher than 99.4% along the length of membrane; the decrease is fairly small and the simulation value is almost the same as the real value. Figure 11 shows the relationship of feed TDS and permeate TDS; as the feed TDS increase linearly, the permeate TDS increase nonlinearly.Figure 7
Flow velocity along the length of membrane.Figure 8
Permeate flux along the length of membrane.Figure 9
Osmosis pressure along the membrane.Figure 10
Concentration polarization along membrane.Figure 11
Values of feed TDS and permeate TDS along membrane.Figure 12
Rec and Ry as the function of feed pressure.Figure 13
Rec and Ry as the function of feed flow rate.Figure 14
Rec and Ry as the function of feed concentration.Figure 15
Rec and Ry as the function feed temperature.Figure 16
Sec value change curve along the feed pressure.From (18), the value SEC is affected by the energy recovery efficiency ε
p
f and the mechanical efficiency ε
p. Here they are set as 0.85 and 0.65, respectively. At the given feed condition of Table 1, calculation result shows that about 4.373 kw·h electricity energy should be consumed to produce one cubic meter freshwater. From Figure 16 it can be found that the value of SEC decreases at initial time when the feed pressure increases; then the SEC increases again as the pressure increases continually; there is a minimizing value as the feed pressure increases from about 30 bars to 90 bars. This means that through improving the operational condition, significant energy saving can be achieved.Based on the well-established model, profiles of the relationship between feed conditions and key performance parameters can be obtained, which are listed in Figures12–16. From Figure 12 it can be seen that Rec (water recovery ratio) increases quickly and almost linearly with the increase of feed pressure, because the increase of feed pressure can cause the increase of permeate flux J
v. The other parameter of salt reject increases nonlinearly and slowly with the increase of feed pressure; at the feed pressure of 59 bar, the salt reject is more than 99.5%. From Figure 13 it can be seen that Rec reduces quickly as the feed flow rate increases, while the Ry has the positive relationship with the feed flow rate; this is because of the fact that velocity in the feed channel has positive linear relationship with feed flow rate. And when the feed flow rate increases, the velocity of feed channel increases correspondingly, but the residence time of seawater in the channel will be largely decreased; there is also not enough time for salt to be transported to the permeate water side, which results in the reduction of Rec. From Figure 14, it can be seen that, as the feed concentration C
f increases, the water recovery and salt reject decrease substantially; this is explained by the needing of much more transmembrane osmotic pressure and more difficulty to resist the diffusion of brine. The effect of temperature on the RO process performance is shown in Figure 15, from which it can be seen that, with the increase of feed temperature, Rec increases substantially, but the Ry decreases accordingly. This can be explained by the increase of transport parameters A
w and B
s and the corresponding increase of permeate flux and solute flux.


## 5. Conclusion

In this paper, based on the solution-diffusion and mass transfer theory, a spiral-wound SWRO process model is established. The model is expressed by DAEs (differential and algebraic equations) of (1)–(21) with some inequality and equality constraints of equipment and water quality. To solve the tough problem with high efficiency and accuracy, orthogonal collocation on finite element is proposed to transform the problem into NLP. All the state variables and intermediate variables are discretized and (1)–(21) are transferred into (28)–(32). Then IPOPT solver under GAMS platform is used for the efficient solution. And the efficient solution yields detailed simulation results related to operational variables and key performance. Profiles of feed pressure, feed concentration, and feed velocity along feed channel are investigated as well as relationship of water recovery with temperature and feed pressure. SEC, which is the key performance parameter of energy consumption, is also calculated. The minimizing value in the middle curve shows that, through improving the operational condition, significant energy saving can be achieved. Factors including the feed conditions and operational temperature are also analyzed. All the work is helpful to further understand the mechanism of SWRO process, and shows the significant potential for energy saving through optimal operation.



---

*Source: 908569-2014-06-23.xml*