remarkRemark hypothesisHypothesis claimClaim A finite element model for
membrane channel Nicolás Carro, David Mora and Jesus Vellojin

# A finite element model for concentration polarization and osmotic effects in
a membrane channel ††thanks: Submitted to the editors DATE. The authors were
partially supported by ANID-Chile through project Anillo of Computational
Mathematics for Desalination Processes ACT210087. The second author was
partially supported by DICREA through project 2120173 GI/C Universidad del
Bío-Bío, by the National Agency for Research and Development, ANID-Chile
through FONDECYT project 1220881, and by project Centro de Modelamiento
Matemático (CMM), ACE210010 and FB210005, BASAL funds for centers of
excellence.

Nicolás Carro GIMNAP-Departamento de Matemática, Universidad del Bío - Bío,
Casilla 5-C, Concepción, Chile<EMAIL_ADDRESS>David Mora GIMNAP-
Departamento de Matemática, Universidad del Bío - Bío, Casilla 5-C,
Concepción, Chile<EMAIL_ADDRESS>Jesus Vellojin Corresponding author.
GIMNAP-Departamento de Matemática, Universidad del Bío - Bío, Casilla 5-C,
Concepción, Chile<EMAIL_ADDRESS>
###### Abstract

In this paper we study a mathematical model that represents the concentration
polarization and osmosis effects in a reverse osmosis cross-flow channel with
porous membranes at some of its boundaries. The fluid is modeled using the
Navier-Stokes equations and Darcy’s law is used to impose the momentum balance
on the membrane. The scheme consist of a conforming finite element method with
the velocity–pressure formulation for the Navier-Stokes equations, together
with a primal scheme for the convection–diffusion equations. The Nitsche
method is used to impose the permeability condition across the membrane.
Several numerical experiments are performed to show the robustness of the
method. The resulting model accurately replicates the analytical models and
predicts similar results to previous works. It is found that the submerged
configuration has the highest permeate production, but also has the greatest
pressure loss of all three configurations studied.

###### keywords:

Navier-Stokes equations, Reverse osmosis, Finite elements, Nitsche method

##  1 Introduction

Membrane-based desalination methods have been widely used for tackling the
water scarcity world problem, due to its relatively low energy consumption
compared to thermal-based methods such as multi-stage flash [39], and its
capability to use renewable or low-grade energy resources [1, 45]. Currently,
reverse osmosis leads this trend, being used in $69\%$ of industrial
desalination plants around the world [17], as it is the most energy-efficient
desalination process and its membrane technology is the most investigated,
accounting for more than $60\%$ of articles concerning the most relevant
membrane-based desalination techniques in 2019 [28].

Models associated with concentration polarization in reverse osmosis (RO) are
usually based on the Navier-Stokes and convection-diffusion equations,
although the Brinkman equations for porous media can also be considered.
Variations in the approach of these equations for the formulation of
concentration polarization equations have been widely discussed due to the
difficulty in obtaining a numerical solution. Geometric parameters (such as
the length and thickness of the channel or the existence of spacers) and
physical parameters (low density and high diffusivity) are associated with
unstable behavior of the numerical methods proposed. This is usually improved
by considering sufficiently fine meshes, at the cost of increasing the memory
space required to run the algorithm, or using the SUPG method, which is
capable of stabilize the scheme for convective dominated problems [10].
Although naturally different, the equations modeling RO have a clear
similarity with the Boussinesq equations, and that is the coupling between
fluid and convection-diffusion models. Behind the Boussinesq equations there
is an extensive bibliography, analysis and modeling, where different
configurations and environmental conditions are assumed in order to analyze
their influence in the development of numerical methods for their solution
(see [3, 7, 14, 15] and the references therein). However, the consideration of
a nonlinear boundary condition to model the water flow in the membrane makes
the difference between both models.

A major difficulty in these systems is the boundary conditions, which can
significantly complicate the study of a problem, especially if they are
boundary conditions on vector functions. In the models associated with RO we
can observe that there are boundary conditions, such as the permeability
associated with the membrane in terms of the osmotic pressure, which present a
nonlinear behavior. The numerical implementation always represents a challenge
at the computational level because the numerical scheme selected to solve the
problem must be adjusted to consider this condition. Now, several techniques
have been proposed over the years to deal with boundary conditions on vector
fieds, such as the slip condition in the Navier-Stokes equations, which has
been implemented using penalty methods, Lagrange multipliers, or the Nitsche
technique [16, 21, 40, 41, 42].

The Nitsche technique was introduced in 1971 [35] as an idea to impose
Dirichlet conditions without the use of Lagrange multipliers. The scheme is
similar to a mesh-size-dependent penalty method, but consistency terms are
added for optimal convergence. Over the years, its use has been extended to
interface problems, contact, and boundary conditions in general.

At the discrete level, the majority of RO models is simulated with Finite
Difference (FD) or Finite Volumes methods (FV), thus little has been
researched in finite elements (FE) numerical algorithms applied to this kind
of problems outside the preexisting commercial software, e.g. COMSOL [38, 29,
23, 18, 37, 36] or ANSYS CFX [30, 43, 26], where boundary implementation is
done through a user interface or user defined macros. The coupling between
velocity and concentration in RO models are treated in a iterative way by many
CFD softwares using FE or FV. For example, we have the works of Wiley and
Fletcher [44, 2, 20], where they implemented a CFD model for fluid and
transport processes. Here, it is observed that CFD software are versatile when
complex geometries are considered. There is also the work of Salcedo-Diaz et
al. [44], where the use of Comsol Multiphysics were use to study the behavior
in a three dimensional model of the coupled system, and a ﬂuid -dynamics and
solute transfer numerical simulation that was studied through the Ansys
software [32].

Some of the ideas behind this CFD softwares uses the approach made in [33],
where we observe one of the first stabilized finite element approximations on
a RO channel. Here, the membrane dimensions are taken into account, together
with a SUPG stabilization for coercivity conservation in the convective
dominated model of the coupling between Navier-Stokes and convection-diffusion
equations. An alternative approach was studied in [6] a numerical method is
developed based on the Prandtl equations that couples laminar fluid and mass
transport for the study of concentration polarization and osmotic effects in a
membrane. The validation of the model is done by comparison with classical
analytical solutions.

In our study, the goal is to propose a finite element method using the Nitsche
technique on an RO model. The difference with respect to other studies in the
literature is that the Nitsche trick allow to impose the Darcy-Starling law by
means of a consistent finite element method. This proposition is inspired by
[25], along with several studies on Nitsche’s method on interface and contact
problems [11, 12, 13], where it is shown that Nitsche’s technique can be
applied to problems with general boundary conditions. The main contribution of
this work is that the resulting discrete model is easy to implement because
linearizing the model depends only on linearizing the flow. Furthermore, to
the best of the authors’ knowledge, this is the first time that the Nitsche
technique is used to impose a nonlinear boundary condition on a cross-flow
model that depends on the unknowns. The model is implemented with the open
source software FEniCS [4], and the results are validated with analytical
solutions of simplified problems.

The structure of this article is as follows: in section 2 the assumptions for
the studied system are stated, and the transport equations with the respective
boundary conditions. The weak formulation for the continuous model is given in
section 2.1. The discretization of the model using the Nitsche method is given
in section 2.3. Then, the numerical results for various operating conditions
and the discussion of the obtained results are presented in section 3.

##  2 Model problem

Consider the domain $\Omega=[0,L]\times[0,d]$, which consists of a channel
with semi-permeable walls at $\Gamma_{m}$ and circular spacers at
$\Gamma_{w}$, with flow input in $\Gamma_{\text{in}}$ and output in
$\Gamma_{\text{out}}$ (see Figure 1). The assumptions for this model are:

  1. 1.

The system is in steady state.

  2. 2.

The system is isothermal (therefore, the heat equation in the system is
omitted).

  3. 3.

The fluid in the channel representing seawater is only composed of water and
salt (sodium chloride). Therefore, the fouling effects on the membrane’s
performance due to scaling of less soluble salts or biofouling caused by
bacteria proliferation, are ignored.

  4. 4.

The fluid is considered to be Newtonian and incompressible with density
$\rho_{0}$.

  5. 5.

The fluid has uniform viscosity $\mu_{0}$ and solute diffusivity through the
solvent $D_{0}$.

  6. 6.

The channel membranes are porous in nature, and its behavior is described by
Darcy-Starling Law.

  7. 7.

The osmotic pressure effect of the dissolved solute is represented by Van’t
Hoff’s Law.

  8. 8.

The effect of pressure drop inside the channel due to viscous effects on the
permeate flux (Darcy-Starling Law) is negligible. Therefore, $\Delta P$ in
equation 10 can be considered constant. As an approximation, we will consider
this valid for a channel whose length $L$ is such that the pressure loss due
to viscous effects is less than $1\%$ the value of $\Delta P$.

  9. 9.

The membrane’s performance does not become affected by wearing, so $I_{0}$ is
constant.

Figure 1: Model of cross-flow membrane filtration and concentration
polarization.

For the described model, we will use a Navier-Stokes (based on assumptions 4
and 5) and a convection-diffussion model (based on assumptions 3 and 5), whose
coupling is well-known as the Boussinesq equations. To make matters precise,
we have

The continuity condition:

(1) $\mathop{\mathrm{div}}\nolimits\,\boldsymbol{u}=0,\text{ in }\Omega.$

The Navier-Stokes equation

(2) $\rho_{0}\boldsymbol{u}\cdot\nabla\boldsymbol{u}=-\nabla
p+\mu_{0}\Delta\boldsymbol{u},\text{ in }\Omega.\\\ $

The Convection-Diffusion equation:

(3) $\boldsymbol{u}\cdot\nabla\theta=D_{0}\Delta\theta,\text{ in }\Omega.\\\ $

With respect to the boundary conditions, we assume an inlet parabolic velocity
profile with a fixed concentration to simulate a fully developed flow of
seawater entering the channel, a normal momentum balance to set the pressure
and a zero diffusive salt flux boundary condition in the outlet. Additionally,
a zero fluid velocity, together with no salt penetration are assumed in the
spacers or impermeable walls. This is summarized in the following equations

(4) $\displaystyle\boldsymbol{u}\cdot\boldsymbol{n}$
$\displaystyle=6u_{0}\frac{y}{d}(1-\frac{y}{d}),\;\theta=\theta_{0},$
$\displaystyle\text{ on }\Gamma_{\text{in}},$ (5)
$\displaystyle\mu_{0}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}-p\boldsymbol{n}$
$\displaystyle=\boldsymbol{0},$ $\displaystyle\text{ on }\Gamma_{\text{out}},$
(6) $\displaystyle-D_{0}(\nabla\theta)\cdot\boldsymbol{n}$ $\displaystyle=0,$
$\displaystyle\text{ on }\Gamma_{\text{out}},$ (7)
$\displaystyle\boldsymbol{u}$ $\displaystyle=\boldsymbol{0},$
$\displaystyle\text{ on }\Gamma_{s},$ (8)
$\displaystyle\theta(\boldsymbol{u}\cdot\boldsymbol{n})-D_{0}\nabla\theta\cdot\boldsymbol{n}$
$\displaystyle=0,$ $\displaystyle\text{ on }\Gamma_{m},\Gamma_{s},$ (9)
$\displaystyle\boldsymbol{u}\cdot\boldsymbol{t}$ $\displaystyle=0,$
$\displaystyle\text{ on }\Gamma_{m},$ (10)
$\displaystyle\boldsymbol{u}\cdot\boldsymbol{n}$ $\displaystyle=\frac{\Delta
P-\kappa\theta}{I_{0}},$ $\displaystyle\text{ on }\Gamma_{m}.$

In the above, $\boldsymbol{u}=(u_{x},u_{y})$, $p$, and $\theta$ represent the
fluid velocity, pressure and molar concentration profile, respectively. Its
worth noting that the pressure $p$ is the gauge pressure and not the absolute
pressure $P$, as the assumption $4$ allows to write the Navier-Stokes
equations in terms of $p$. The constant $D_{0}$ represents the solute
diffusivity trough the solvent, $I_{0}$ is the membrane resistance, whereas
$\kappa$ is a constant given by the Darcy-Starling law. $\Delta P$ represents
the transmembrane pressure between the membrane and the “outside” (the
permeate channel).

###  2.1 Formulation

Assuming free flow conditions in $x=L$, we define the following spaces.

$\displaystyle\mathrm{H}:=\\{\boldsymbol{v}\in
H^{1}(\Omega)^{2}\;:\;\boldsymbol{v}\cdot\boldsymbol{n}=u_{0}\text{ on
}\Gamma_{\text{in}}\\},$ $\displaystyle Y:=\\{\boldsymbol{v}\in
H^{1}(\Omega)^{2}\;:\;\boldsymbol{v}=\boldsymbol{0}\text{ on
}\Gamma_{\text{in}}\cup\Gamma_{s},\;\boldsymbol{v}\cdot\boldsymbol{t}=\boldsymbol{0}\text{
on }\Gamma_{m}\\},\qquad Q=L^{2}(\Omega),$ $\displaystyle Z:=\\{\tau\in
H^{1}(\Omega)\;:\;\tau=\theta_{0}\text{ on }\Gamma_{\text{in}}\\},\qquad
X:=\\{\tau\in H^{1}(\Omega)\;:\;\tau=0\text{ on }\Gamma_{\text{in}}\\}.$

Testing (2) with $\boldsymbol{v}\in Y$, integrating by parts and using the
boundary conditions we obtain

$\displaystyle\int_{\Omega}\rho_{0}\left(\boldsymbol{u}\cdot\nabla\boldsymbol{u}\right)\cdot\boldsymbol{v}$
$\displaystyle=\left[\int_{\Omega}p\mathop{\mathrm{div}}\nolimits\boldsymbol{v}-\int_{\partial\Omega}(p\boldsymbol{n})\cdot\boldsymbol{v}\right]+\mu_{0}\left[-\int_{\Omega}\nabla\boldsymbol{u}:\nabla\boldsymbol{v}+\int_{\partial\Omega}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}\cdot\boldsymbol{v}\right]$
$\displaystyle=-\int_{\Omega}\mu_{0}\nabla\boldsymbol{u}:\nabla\boldsymbol{v}+\int_{\Omega}p\mathop{\mathrm{div}}\nolimits\boldsymbol{v}+\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}-p\boldsymbol{n}\right)(\boldsymbol{v}\cdot\boldsymbol{n})$

Testing (3) with $\tau\in X$ and integrating by parts we obtain

$\displaystyle\int_{\Omega}(\boldsymbol{u}\cdot\nabla\theta)\tau$
$\displaystyle=D_{0}\left[-\int_{\Omega}\nabla\theta\cdot\nabla\tau+\int_{\partial\Omega}\frac{\partial\theta}{\partial
n}\tau\right]$
$\displaystyle=-\int_{\Omega}D_{0}\nabla\theta\cdot\nabla\tau+\int_{\Gamma_{m}\cup\Gamma_{s}}\theta(\boldsymbol{u}\cdot\boldsymbol{n})\tau.$

Lets introduce the following forms:

(11) $\displaystyle
a(\boldsymbol{u},\boldsymbol{v}):=\mu_{0}\int_{\Omega}\nabla\boldsymbol{u}:\nabla\boldsymbol{v},\qquad\widetilde{a}(\boldsymbol{w},(\boldsymbol{u},\boldsymbol{v})):=\rho_{0}\int_{\Omega}(\boldsymbol{w}\cdot\nabla\boldsymbol{u})\cdot\boldsymbol{v},$
$\displaystyle
d(\theta,\tau):=D_{0}\int_{\Omega}\nabla\theta\cdot\nabla\tau,\qquad
b(\boldsymbol{v},q):=-\int_{\Omega}q\mathop{\mathrm{div}}\nolimits\boldsymbol{v}$
$\displaystyle\widetilde{c}(\boldsymbol{w},(\theta,\tau)):=\int_{\Omega}(\boldsymbol{w}\cdot\nabla\theta)\tau-\int_{\Gamma_{m}\cup\Gamma_{s}}\theta\left(\boldsymbol{w}\cdot\boldsymbol{n}\right)\tau.$

Hence, we have the following formulation for the proposed cross-flow model:
Find $(\boldsymbol{u},p,\theta)\in\mathrm{H}\times Q\times Z$ such that (10)
holds, and

(12) $\displaystyle
a(\boldsymbol{u},\boldsymbol{v})+b(\boldsymbol{v},p)+\widetilde{a}(\boldsymbol{u},(\boldsymbol{u},\boldsymbol{v}))-\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}-p\boldsymbol{n}\right)(\boldsymbol{v}\cdot\boldsymbol{n})$
$\displaystyle=0,$ $\displaystyle\forall\boldsymbol{v}\in Y,$ $\displaystyle
b(\boldsymbol{u},q)$ $\displaystyle=0,$ $\displaystyle\forall q\in Q,$
$\displaystyle d(\theta,\tau)+\widetilde{c}(\boldsymbol{u},(\theta,\tau))$
$\displaystyle=0,$ $\displaystyle\forall\tau\in X.$

Note that if homogeneous boundary conditions for both $\boldsymbol{u}$ and
$\theta$ are considered in $\Gamma_{m}$, then we fall into a typical
formulation of Boussinesq equations. For an analysis on viscous flows around
porous with non-homogeneous boundary data, we refer to [24]. On the other
hand, note that the nonlinearity inherited by the boundary condition (8) is
present in the trilinear form $\widetilde{c}$. For the bilinear form
$b(\cdot,\cdot)$ there holds [22, Theorem 3.7]

(13) $\sup_{\boldsymbol{v}\in
Y,\boldsymbol{v}\neq\boldsymbol{0}}\frac{b(\boldsymbol{v},q)}{\|\boldsymbol{v}\|_{1,\Omega}}\geq\beta\|q\|_{0,\Omega},\qquad
q\in W.$

###  2.2 Nitsche method

Let us consider a shape-regular family of partitions of $\Omega$, denoted by
$\\{\mathcal{T}_{h}\\}_{h>0}$. Let $h_{T}$ be the diameter of a triangle
$T\in{\mathcal{T}}_{h}$ and let us define
$h:=\max\\{h_{T}\,:\,T\in{\mathcal{T}}_{h}\\}$.

One of the main difficulties in the discrete solution of problems involving
the Navier-Stokes equations is the choice of elements satisfying the discrete
version of (13). Of all the possible existing options, we choose the Taylor-
Hood elements [8, 9], namely, given $k\geq 1$, the finite element spaces to
develop the numerical scheme for the velocity and pressure are the following

$\displaystyle\mathbf{H}_{h}$ $\displaystyle:=\\{\boldsymbol{v}_{h}\in
C(\Omega)^{2}\,:\,\boldsymbol{v}_{h}|_{T}\in\mathbb{P}_{k+1}(T)^{2}\quad\forall
T\in\mathcal{T}_{h}\\}\cap Y,$ $\displaystyle Q_{h}$
$\displaystyle:=\\{q_{h}\in
C(\Omega)\,:\,q_{h}|_{T}\in\mathbb{P}_{k}(T)\quad\forall
T\in\mathcal{T}_{h}\\}.$

Another widely used and computationally very economical option is the lowest
order mini-element pair $(\mathbf{H}_{h},Q_{h})$, where $Q_{h}$ is as above
with $k=1$, and $\mathbf{H}_{h}$ is given by

$\mathbf{H}_{h}:=\\{\boldsymbol{v}_{h}\in
C(\Omega)^{2}\,:\,\boldsymbol{v}_{h}|_{T}\in\mathbb{P}_{1,b}(T)^{2}\quad\forall
T\in\mathcal{T}_{h}\\}\cap Y,$

where $\mathbb{P}_{1,b}$ are the bubble-enriched piecewise continuous elements
defined with the barycentric coordinates of each element $T$ (see, e.g. [19,
22]).

In turn, for $\ell\geq 1$, the finite element space for the discrete
concentration is just

$M_{h}:=\\{\tau_{h}\in
C(\Omega)\,:\,\tau_{h}|_{T}\in\mathbb{P}_{\ell}(T)\quad\forall
T\in\mathcal{T}_{h}\\}\cap X.$

It is important to mention that condition (10) is not imposed in throughout
the integration by parts. Hence, we use the Nitsche’s method, whose goal is to
relax and penalize Dirichlet boundary conditions such that they are satisfied
in the limit of small mesh size. In this study, we will use this technique in
order to impose the boundary condition (10), which characterizes the pressure
difference between the interior of the channel (feed) and the outside
(permeate). In order to propose our numerical scheme using the Nitsche method,
we introduce the following bilinear forms

(14)
$A((\boldsymbol{u}_{h},p_{h}),(\boldsymbol{v}_{h},q_{h})):=a(\boldsymbol{u}_{h},\boldsymbol{v}_{h})+b(\boldsymbol{v}_{h},p_{h})+b(\boldsymbol{u}_{h},q_{h})\\\
-\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{u}_{h}}{\partial\boldsymbol{n}}-p_{h}\boldsymbol{n}\right)(\boldsymbol{v}_{h}\cdot\boldsymbol{n})-\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{v}_{h}}{\partial\boldsymbol{n}}-q_{h}\boldsymbol{n}\right)\left(\boldsymbol{u}_{h}\cdot\boldsymbol{n}\right)\\\
+\alpha\int_{\Gamma_{m}}h^{-1}(\boldsymbol{u}_{h}\cdot\boldsymbol{n})(\boldsymbol{v}_{h}\cdot\boldsymbol{n}),\\\
F(\boldsymbol{v}_{h},q_{h}):=\alpha\int_{\Gamma_{m}}h^{-1}\left(\frac{\Delta
P}{I_{0}}\right)(\boldsymbol{v}_{h}\cdot\boldsymbol{n})-\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{v}_{h}}{\partial\boldsymbol{n}}-q_{h}\boldsymbol{n}\right)\left(\frac{\Delta
P}{I_{0}}\right)\\\
B(\theta_{h},(\boldsymbol{v}_{h},q_{h})):=\alpha\int_{\Gamma_{m}}h^{-1}\left(\frac{\kappa\theta_{h}}{I_{0}}\right)(\boldsymbol{v}_{h}\cdot\boldsymbol{n})-\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{v}_{h}}{\partial\boldsymbol{n}}-q_{h}\boldsymbol{n}\right)\left(\frac{\kappa\theta_{h}}{I_{0}}\right),$

where $h$ represents the local mesh size. The parameter $\alpha>0$ is a
constant that has to be chosen sufficiently large. With this forms at hand, we
propose the following Nitsche method: Find
$(\boldsymbol{u}_{h},p_{h},\theta_{h})\in\mathbf{H}_{h}+u_{0,I}\times
Q_{h}\times M_{h}+\theta_{0,I}$ such that

(15) $\displaystyle
A((\boldsymbol{u}_{h},p_{h}),(\boldsymbol{v}_{h},q_{h}))+\widetilde{a}(\boldsymbol{u}_{h},(\boldsymbol{u}_{h},\boldsymbol{v}_{h}))+B(\theta_{h},(\boldsymbol{v}_{h},q_{h}))$
$\displaystyle=F(\boldsymbol{v}_{h},q_{h}),$ $\displaystyle
d(\theta_{h},\tau_{h})+\widetilde{c}(\tau_{h},(\boldsymbol{u}_{h},\theta_{h}))$
$\displaystyle=0,$

for all $(\boldsymbol{v}_{h},q_{h})\in\mathbf{H}_{h}\times Q_{h}$ and for all
$\forall\tau_{h}\in M_{h}$, where $u_{0,I}$ and $\theta_{0,I}$ are suitable
interpolants of $u_{0}$ and $\theta_{0}$.

Below we prove that (15) contains all the elements to be consistent.

######  Lemma 2.1.

If $(\boldsymbol{u},p,\theta)$ is the solution of (1)-(3) together with
boundary conditions (6)-(10), then $(\boldsymbol{u},p,\theta)$ solves (15).

######  Proof 2.2.

We start by multiplying (10) by $\boldsymbol{v}_{h}\cdot\boldsymbol{n}$, with
$\boldsymbol{v}_{h}\in\mathbf{H}_{h}$, integrate over $\Gamma_{m}$, and
multiply by $\alpha h^{-1}$ in order to have

(16)
$\alpha\int_{\Gamma_{m}}h^{-1}(\boldsymbol{u}\cdot\boldsymbol{n})(\boldsymbol{v}_{h}\cdot\boldsymbol{n})=\alpha\int_{\Gamma_{m}}h^{-1}\left(\frac{\Delta
P-\kappa\theta}{I_{0}}\right)(\boldsymbol{v}_{h}\cdot\boldsymbol{n}).$

Similarly, we now multiply $\eqref{model-eq9}$ by
$\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{v}_{h}}{\partial\boldsymbol{n}}-q_{h}\boldsymbol{n}\right)$,
where $q_{h}\in Q_{h}$, to obtain

(17)
$\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{v}_{h}}{\partial\boldsymbol{n}}-q_{h}\boldsymbol{n}\right)(\boldsymbol{u}\cdot\boldsymbol{n})=\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{v}_{h}}{\partial\boldsymbol{n}}-q_{h}\boldsymbol{n}\right)\left(\frac{\Delta
P-\kappa\theta}{I_{0}}\right).$

Note that multiplying (1)-(3) with
$(\boldsymbol{v}_{h},q_{h},\tau_{h})\in\mathbf{H}_{h}\times Q_{h}\times
M_{h}$, integrating by parts and using boundary conditions (6)-(9) we have
that

(18) $\displaystyle
a(\boldsymbol{u},\boldsymbol{v}_{h})+b(\boldsymbol{v}_{h},p)+\widetilde{a}(\boldsymbol{u},(\boldsymbol{u},\boldsymbol{v}_{h}))-\int_{\Gamma_{m}}\boldsymbol{n}^{t}\left(\mu_{0}\frac{\partial\boldsymbol{u}}{\partial\boldsymbol{n}}-p\boldsymbol{n}\right)(\boldsymbol{v}_{h}\cdot\boldsymbol{n})$
$\displaystyle=0,$ $\displaystyle b(\boldsymbol{u},q_{h})$ $\displaystyle=0,$
$\displaystyle
d(\theta,\tau_{h})+\widetilde{c}(\boldsymbol{u},(\theta,\tau_{h}))$
$\displaystyle=0.$

Then, by adding (16), (17) and (18) we obtain that

$\displaystyle
A((\boldsymbol{u},p),(\boldsymbol{v}_{h},q_{h}))+\widetilde{a}(\boldsymbol{u},(\boldsymbol{u},\boldsymbol{v}_{h}))+B(\theta,(\boldsymbol{v}_{h},q_{h}))$
$\displaystyle=F(\boldsymbol{v}_{h},q_{h}),$ $\displaystyle
d(\theta,\tau_{h})+\widetilde{c}(\tau_{h},(\boldsymbol{u},\theta))$
$\displaystyle=0,$

for all $(\boldsymbol{v}_{h},q_{h})\in\mathbf{H}_{h}\times Q_{h}$ and for all
$\forall\tau_{h}\in M_{h}$. The proof is complete.

The above lemma gives the modified variational formulation using the Nitsche
method, and its consistency. Note that the nonlinearity is present in the
convective terms and the boundary term
$\int_{\Gamma_{m}\cup\Gamma_{s}}\theta(\boldsymbol{u}\cdot\boldsymbol{n})\tau$,
where the nonlinearity in the equations is caused only by the velocity. Hence,
the system (15) can be rearranged such that a fixed point iteration technique
is straightforward.

##  3 Numerical experiments: validation and discussion

In this section we will show some numerical results in different channel
configurations. The implementation of the scheme is perfomed in FEniCS, and
the meshes were created using Gmsh. The numerical experiments have been
developed with Taylor-Hood elements. However, similar results were obtained
using the mini-element. For the concentration case, it was sufficient to
consider the lowest case with $\ell=1$. All the cases consider the Nitsche
constant $\alpha=1$.

We consider the geometrical parameters of a common membrane channel “unit”
whose length is defined by a subsection of the channel that allows a fully
developed flow [31]. The channel thickness and spacer diameter is that of a
module with standard $28\,mil$ gap spacer mesh [27]. For the channel width $W$
(perpendicular to the simulation domain) we will consider $W=L$. More
precisely, we have

$L=1.5\cdot 10^{-2}\,m,\quad d=7.4\cdot 10^{-4}\,m,\qquad W=1.5\cdot
10^{-2}\,m,\qquad d_{S}=3.6\cdot 10^{-4}\,m,$

while the physical parameters [6, 34] are taken as

$\displaystyle\kappa=4955.144\,Pa,\quad I_{0}=8.41\cdot 10^{10}Pa\,\cdot
s/m,\quad\mu_{0}=8.9\cdot 10^{-4}kg\,m^{-1}s^{-1},$
$\displaystyle\rho_{0}=1027.2\,kg/m^{3},\quad D_{0}=1.5\cdot 10^{-9}m^{2}/s.$

With respect to the boundary data, we will consider combinations among the
following:

$\displaystyle u_{0}=1.29\times k\cdot 10^{-1}m/s,\quad k=0.5,1,2.$
$\displaystyle C_{0}=600\text{ mol }m^{3}.$ $\displaystyle\Delta
P=4053000\,Pa\text{ or }5572875\,Pa.$

As suggested by several studies, all the experiments use a highly refined mesh
close to the membrane $\Gamma_{m}$.

###  3.1 Validation test

This experiment aims to test the convergence of our method. To this end, we
consider a channel with no spacers, i.e., we set $\Gamma_{s}=\emptyset$. The
domain $\Omega=(0,L)\times(0,d)$ is meshed with uniform elements such that
$n_{x}=(L/d)n_{y}$, where $n_{x}$ and $n_{y}$ denote the number of cells in
the $x$ and $y$ directions, respectively. We divide this test in two cases,
that we detail below.

####  3.1.1 Case 1

The first case consists in validating the method by comparing the pressure
drop with those of classical analytical models of momentum, namely, Poiseuille
and Berman flow models. These two models have the advantage that the pressure
drop, denoted by $\Delta p(x,d/2):=p(0,d/2)-p(x,d/2)$, can be obtained by
solving the equations of motion with boundary conditions

$\boldsymbol{u}\cdot\boldsymbol{n}=\frac{\Delta
P}{I_{0}},\quad\boldsymbol{u}\cdot\boldsymbol{n}=0,$

on $\Gamma_{m}$, for Berman and Poiseuille flow, respectively. Note that, for
the Poiseuille flow, we have a no-slip and no-penetration boundary condition
(Dirichlet boundary condition) that can be imposed strongly on the space.
Hence, we have the following pressure drop equation for a constant permeable
wall

$\Delta
p(x,d/2):=\left(\frac{1}{2}\rho_{0}u_{0}^{2}\right)\left(\frac{24}{Re}-\frac{648}{35}\frac{Re_{\boldsymbol{n}}}{Re}\right)\left(1-\frac{2Re_{\boldsymbol{n}}}{Re}\frac{x}{\widetilde{d}}\right)\left(\frac{x}{\widetilde{d}}\right),$

where a constant permeate velocity $\boldsymbol{u}\cdot\boldsymbol{n}$ is
assumed. On the other hand, the case of impermeable walls assumes that
$Re_{\boldsymbol{n}}=0$, so the pressure drop equation is reduced to

$\Delta
p(x,d/2):=\left(\frac{1}{2}\rho_{0}u_{0}^{2}\right)\left(\frac{24}{Re}\right)\left(\frac{x}{\widetilde{d}}\right).$

In the above equations, $Re:=\frac{4\rho_{0}\widetilde{d}u_{0}}{\mu_{0}}$ is
the cross flow Reynolds number, whereas
$Re_{\boldsymbol{n}}=\frac{\rho_{0}\tilde{d}(\boldsymbol{u}\cdot\boldsymbol{n})}{\mu_{0}}$
is the Reynolds number at the channel walls. The parameter $\widetilde{d}$
denotes half of the channel height ($\widetilde{d}=d/2$).

In Figure 2 we observe the exact pressure drop compared with our computed
pressure drop, denoted by $\Delta p_{h}(x,d/2)$. To compare with our
approximate solution, we have selected $10$ equally spaced pressure drop
values in the channel. We note that the proposed method agrees perfectly with
the analytical axial pressure drop for both, the permeable and impermeable
channel. Also, a similar pressure drop was observed for $\Delta P=4053000$ Pa
and $\Delta P=5572875$ Pa since they correspond to
$\boldsymbol{u}\cdot\boldsymbol{n}\approx 4.819\cdot 10^{-5}$ m/s and
$\boldsymbol{u}\cdot\boldsymbol{n}\approx 6.626\cdot 10^{-5}$ m/s,
respectively. Therefore, $\Delta P$ has little impact on the pressure drop on
the channel.

Figure 2: Test 1. Comparison of the exact and approximate pressure drop.

####  3.1.2 Case 2

We now aim to test the method when the non-linear boundary condition (10) is
considered. To this end, we compute the total flow (see equation (19)) over
the membrane channel in order to observe the stability of the result for
different number of elements:

(19) $\dot{m}_{tot}=\rho W\int_{\Gamma_{m}}\lvert\boldsymbol{u}_{y}(s)\rvert
ds.$

The importance of this test lies in the determination of the number of
elements necessary to obtain accurate results. Some studies suggest that
refinement towards the membrane are required for better results [29, 31]. A
total of four refinement near $\Gamma_{m}$ have been used. We set
$N=2j,\,j=1,\ldots,19,$ for uniform refinement, whereas $N=5j$ is used for the
meshing towards $\Gamma_{m}$.

In Figure 3 we observe the average flux for several inlet velocities and
osmotic pressure when using uniform refinements. In both cases, it is observed
that from $5\cdot 10^{3}$ elements the total mass flow is stable. All the
simulations in this work were tested in a similar manner, preferring the
meshing towards $\Gamma_{m}$ with a variable amount of elements, but always
above $10^{4}$ cells.

Figure 3: Test 1. Comparison of the stability in the total mass flow
computation stability when using uniform refinement and meshes refined towards
$\Gamma_{m}$.

###  3.2 Test 2. Channel with different spacer configurations

In this section we will depict the velocity fields, pressure distribution and
concentration polarization along a channel with different spacer
configurations. Discussion of the simulations are presented in Section 3.3
below.

####  3.2.1 Spacers in cavity configuration

The numerical results for the spacers in cavity configuration are presented in
Figures 4-6.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=5572875Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=5572875Pa$

Figure 4: Test 2. Comparison of concentration levels using a cavity spacers
configuration.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$ $u_{0}=2.58\cdot
10^{-1}m/s,\Delta P=4053000Pa$

Figure 5: Test 2. Comparison of velocity fields using a cavity spacers
configuration.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=5572875Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=5572875Pa$

Figure 6: Test 2. Comparison of relative pressure using a cavity spacers
configuration.

####  3.2.2 Zig-zag spacers configuration

The numerical results for the spacers in zig-zag configuration are presented
in Figures 7-9.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=5572875Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=5572875Pa$

Figure 7: Test 2. Comparison of concentration levels using a zig-zag spacers
configuration.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

Figure 8: Test 2. Comparison of velocity fields using a zig-zag spacers
configuration.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=5572875Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=5572875Pa$

Figure 9: Test 2. Comparison of relative pressure using a zig-zag spacers
configuration.

####  3.2.3 Submerged spacers configuration

The numerical results for the spacers in submerged configuration are presented
in Figures 10-12.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=5572875Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=5572875Pa$

Figure 10: Test 2. Comparison of concentration levels using a submerged
spacers configuration.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

Figure 11: Test 2. Comparison of velocity fields using a submerged spacers
configuration.

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=4053000Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=4053000Pa$

$u_{0}=6.45\cdot 10^{-2}m/s,\Delta P=5572875Pa$

$u_{0}=2.58\cdot 10^{-1}m/s,\Delta P=5572875Pa$

Figure 12: Test 2. Comparison of relative pressure using a submerged spacers
configuration.

###  3.3 Discussion

This study focuses on the effect of three main parameters on the behavior of
the membrane channel: the transmembrane pressure, the inlet velocity (i.e.,
the inlet mass flow) and the spacer configuration.

Analysis in Figures 4, 5 and 6 show that the increase in mass flow at the
inlet (increase in inlet velocity) increases the length of the boundary layer
and diminishes its concentration. Also we can establish that $\Delta P$ does
not have any appreciable effect on the hydrodynamics of the channel for equal
inlet mass flows. Indeed, the permeate flow difference between different
$\Delta P$, although significative in terms of the obtained permeate flow, is
two to three orders of magnitude smaller than the inlet velocity. These
behaviors are further confirmed by Figure 13.

In all of the performed cases, the main cause of pressure loss is the presence
of a bottleneck caused by the spacers. As seen in Figure 5, the spacers cause
a sudden augmentation of velocity and velocity gradients in those zones, which
in turn derives into major pressure losses. The increment in inlet flow gets
bigger as this gradient changes, and the pressure losses become greater.
However, for increasing $\Delta P$ there is a subtle decrease in pressure
loss, because the increment in permeate flux at higher $\Delta P$ decreases
the amount of mass flow that has to cross the bottleneck, effectively reducing
the pressure loss in each stage.

Both opposing behaviors leads to conclude that there is an optimum operating
condition for the system from the energy efficiency point of view, in which
both the increasing pressure loss due to increasing pressure flow and the
reducing pressure loss from the increment in transmembrane pressure leads to
an overall minimum pressure loss in the channel. This optimization has well
established boundaries, as $\Delta P$ must be lower than the liquid entry
pressure (LEP) of the membrane which would cause the membrane’s rupture, but
must be higher than the osmotic pressure of the feed solution, as otherwise
the process would turn into Forward osmosis (FO). To determine these optimum
values, further numerical essays need to be performed, which will be addressed
in a subsequent article.

The previous analysis is the same for a zig-zag and submerged configurations
upon comparing Figures 7–9 and 10–12, respectively. However, for the case of a
submerged configuration, the boundary layer of salt is not developed as seen
in Figure 10, mainly due to the absence of recirculation zones close to the
membrane walls, depicted in Figure 11. The recirculating fluid is located in
the middle of the channel behind the spacers, away from the influence of the
membrane wall concentration boundary layer (see Figure 10). Therefore, the
effects of accumulation of salt in the recirculation zone will be greater
along the channel, provided that the boundary layer is large enough.

In the case of pressure losses in submerged configurations, in contrast to the
cavity and zig-zag configurations, has two bottlenecks in each spacer that are
smaller, causing more pressure loss inside the channel. As the bottlenecks are
reduced to half the size, both the velocity gradients and the convective flux
at the bottenecks increase, therefore increasing the pressure losses by
viscous friction and convective losses.

The main variables studied in this paper for all conditions are shown in
Figures 13 through 16. From inspection of the permeate velocity curves in
Figures 13–14, we can say that incrementing the inlet mass flow helps to
diminish the concentration in the whole membrane, therefore reducing the
osmotic pressure effect in the permeate flux. Also, the increase in pressure
leads to an increase in concentration at the membrane, as more water is being
separated out, leaving behind more salt. Also, as expected, the permeate flux
increases when $\Delta P$ is increased.

Pressure profiles in Figures 13–14 show the pressure loss along the membrane
channel for the cavity configuration for all different systems. The two main
conclusions are that $\Delta P$ has a negligible effect on the pressure loss
of the system, as stated earlier in the discussion, and the increment in the
inlet velocity has a positive contribution in the pressure loss. Indeed, for
twice the inlet velocity, the pressure losses are increased by a factor
greater than two. The negative pressure values in all Figures are valid due to
the incompressible formulation of the problem, which allows to write the
variable $p$ in terms of the gauge pressure.

Figure 13: Test 1. Comparison between permeate velocity, concentration level
and relative pressure drop on the membrane walls along the line $y=0$ and
$x\in[0,L]$ using the cavity spacers configuration. Figure 14: Test 1.
Comparison between permeate velocity, concentration level and relative
pressure drop on the membrane walls along the line $y=0$ and $x\in[0,L]$ using
the submerged spacers configuration. Figure 15: Test 1. Comparison between
permeate velocity, concentration level and relative pressure drop on the
membrane walls along the line $y=0$ and $x\in[0,L]$ using the zig-zag spacers
configuration. Figure 16: Comparison between absolute permeate of different
configurations along the line $y=0$ and $x\in[0,L]$.

Additionally, the effect of spacer configuration on the membrane’s performance
needs to be analyzed. For this, we fix the inlet velocity and transmembrane
pressure at $2.58\cdot 10^{-1}\,m/s$ and $\Delta P=5572875\,Pa$, respectively.
From Figure 16, it can be concluded that spacers in both cavity and zig-zag
configurations have a negative impact on permeate velocity on their vicinity,
as shown by the sharp decrease in the permeate velocity at the positions where
the spacers are located. The opposite effect occurs for the submerged
configuration, where a local maximum can be located at each spacer position.
This is due to the increase in flow in the channel bottlenecks, which decrease
the salt concentration and consequently the osmotic resistance to permeate
flux. However, the mixing effects are less effective than the other two
configurations, as the concentration at the end of the submerged configuration
channel is higher and the permeate flux is the lowest of all three
configurations. We can also infer that the high inlet velocity in the zig-zag
configuration as seen in Figure 8 produces a dominant fluid path, which
diminishes the concentration between the first two spacers and therefore
increases the overall permeate flux. Lastly, the total volumetric flow per
unit width for all membrane channels in different configurations, for the same
conditions used in Figure 16, is given in Table 1. This quantity can be
obtained from direct integration of the definition of volumetric flow from the
flux:

(20) $\frac{\dot{V}}{W}=\int_{\Gamma_{m}}\lvert u_{y}(s)\rvert ds,$

where $\dot{V}$ is the total volumetric flow of the permeate from both
membrane walls, $W$ is the width of the channel, and $u_{y}\left(s\right)$ is
the permeate flux at a distance from the inlet $s$. The absolute value is used
because on the lower side of the membrane the velocity vector is negative with
respect to the coordinate origin.

Table 1: Total permeate flow per unit length $\dot{V}/W$ for all membrane channel simulations. Configuration | $u_{0}\,[m/s]$ | $\Delta P\,[Pa]$ | $\dot{V}/W\,[m^{3}/(s\cdot m)]$
---|---|---|---
Cavity | 0.258 | 4053000 | $2.39347\cdot 10^{-7}$
Cavity | 0.129 | 4053000 | $2.16131\cdot 10^{-7}$
Cavity | 0.064 | 4053000 | $1.93002\cdot 10^{-7}$
Cavity | 0.258 | 5572875 | $5.59461\cdot 10^{-7}$
Cavity | 0.129 | 5572875 | $5.00098\cdot 10^{-7}$
Cavity | 0.064 | 5572875 | $4.41913\cdot 10^{-7}$
Submerged | 0.258 | 4053000 | $2.51089\cdot 10^{-7}$
Submerged | 0.129 | 4053000 | $2.28837\cdot 10^{-7}$
Submerged | 0.064 | 4053000 | $2.03092\cdot 10^{-7}$
Submerged | 0.258 | 5572875 | $5.88266\cdot 10^{-7}$
Submerged | 0.129 | 5572875 | $5.31309\cdot 10^{-7}$
Submerged | 0.064 | 5572875 | $4.66112\cdot 10^{-7}$
Zig-zag | 0.258 | 4053000 | $2.38891\cdot 10^{-7}$
Zig-zag | 0.129 | 4053000 | $2.16246\cdot 10^{-7}$
Zig-zag | 0.064 | 4053000 | $1.93055\cdot 10^{-7}$
Zig-zag | 0.258 | 5572875 | $5.57996\cdot 10^{-7}$
Zig-zag | 0.129 | 5572875 | $5.01593\cdot 10^{-7}$
Zig-zag | 0.064 | 5572875 | $4.42017\cdot 10^{-7}$

From the values of Table 1 it follows that the submerged configuration has the
most production of permeate of all three configurations for all of the
considered variations of $u_{0}$ and $\Delta P$. Also, zig-zag and cavity
configurations show similar results, with zig-zag being slightly more
efficient than its alternate counterpart. However, for longer channels this
tendency may be inverted, as the submerged configuration concentration profile
may rise well above the other configurations, therefore diminishing its
permeate flux. The results also show that quadrupling the value of $u_{0}$
while mantaining $\Delta P$ can generate an increase in $\dot{V}/W$ of up to
$24\%$ for the lower pressure and $26\%$ for the higher pressure.
Additionally, increasing the pressure from $4053000\,Pa$ to $5572875\,Pa$
drastically rises the value up to about $130\%$ for the lower inlet velocity
and up to $134\%$ for the higher value.

##  4 Conclusions

In this work we have presented the performance of the numerical method for a
RO model using the Nitsche technique. The method has been successfully
validated by testing the accuracy with respect to the pressure drop, as well
as the number of cells required to have independence with respect to the mesh
size. However, for the highest speed of operation, it was required to
implement the SUPG scheme to stabilize the algorithm.

Also, we can follow the same approach in this study to obtain the classical
penalty formulation, which also serves to impose the permeability condition.
However, it is well known that, for very large values of $\alpha$, the system
will be ill conditioned, and consequently, the convergence of the method will
be affected [5].

The obtained tendencies for permeate flux against transmembrane pressure and
inlet mass flow are in accordance with previous works. It was also determined
that the bottlenecks caused by the spacers were the main contributors to
pressure loss (and therefore, energy loss), being the biggest drop associated
with the submerged configuration. Considering all process parameters constant,
the greatest permeate flux was obtained for the submerged configuration,
although the higher concentration at the exit in comparison with its
counterparts suggest that this advantage may not hold for bigger channels, as
the submerged configuration is less effective at disrupting the boundary layer
at the membrane walls.
For future development of the proposed model, further studies are needed to
establish the influence of the Nitsche technique on the stability of the
method. Additionally, other factors must be considered to accurately capture
the complex behavior of a membrane module in operational conditions, e.g., the
inclusion of fouling agents and the effects of three-dimensional spacer
mixing.

## References

  * [1] F. E. Ahmed, A. Khalil, and N. Hilal, Emerging desalination technologies: Current status, challenges and future trends, Desalination, 517 (2021), p. 115183.
  * [2] A. Alexiadis, D. Wiley, A. Vishnoi, R. Lee, D. Fletcher, and J. Bao, CFD modelling of reverse osmosis membrane flow and validation with experimental results, Desalination, 217 (2007), pp. 242–250.
  * [3] K. Allali, A priori and a posteriori error estimates for Boussinesq equations, Int. J. Numer. Anal. Model, 2 (2005), pp. 179–196.
  * [4] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, The FEniCS project version 1.5, Archive of Numerical Software, 3 (2015).
  * [5] J. W. Barrett and C. M. Elliott, Finite Element Approximation of the Dirichlet problem using the boundary penalty method, Numerische Mathematik, 49 (1986), pp. 343–366.
  * [6] B. Bernales, P. Haldenwang, P. Guichardon, and N. Ibaseta, Prandtl model for concentration polarization and osmotic counter-effects in a 2-d membrane channel, Desalination, 404 (2017), pp. 341–359.
  * [7] C. Bernardi, B. Métivet, and B. Pernaud-Thomas, Couplage des équations de Navier-Stokes et de la chaleur: le modèle et son approximation par éléments finis, ESAIM: Mathematical Modelling and Numerical Analysis, 29 (1995), pp. 871–921.
  * [8] D. Boffi, Stability of higher order triangular Hood-Taylor methods for the stationary Stokes equations, Mathematical Models and Methods in Applied Sciences, 4 (1994), pp. 223–235.
  * [9] , Three-dimensional Finite Element methods for the Stokes problem, SIAM journal on numerical analysis, 34 (1997), pp. 664–670.
  * [10] A. N. Brooks and T. J. Hughes, Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer methods in applied mechanics and engineering, 32 (1982), pp. 199–259.
  * [11] F. Chouly and P. Hild, On convergence of the penalty method for unilateral contact problems, Applied Numerical Mathematics, 65 (2013), pp. 27–40.
  * [12] F. Chouly, P. Hild, and Y. Renard, A Nitsche Finite Element method for dynamic contact: 1. space semi-discretization and time-marching schemes, ESAIM: Mathematical Modelling and Numerical Analysis, 49 (2015), pp. 481–502.
  * [13] , A Nitsche Finite Element method for dynamic contact: 2. stability of the schemes and numerical experiments, ESAIM: Mathematical Modelling and Numerical Analysis, 49 (2015), pp. 503–528.
  * [14] E. Colmenares, G. N. Gatica, and R. Oyarzúa, Analysis of an augmented mixed-primal formulation for the stationary Boussinesq problem, Numerical Methods for Partial Differential Equations, 32 (2016), pp. 445–478.
  * [15] , An augmented fully-mixed Finite Element method for the stationary Boussinesq problem, Calcolo, 54 (2017), pp. 167–205.
  * [16] I. Dione and J. M. Urquiza, Penalty: finite Element Approximation of Stokes equations with slip boundary conditions, Numerische Mathematik, 129 (2015), pp. 587–610.
  * [17] J. Eke, A. Yusuf, A. Giwa, and A. Sodiq, The global status of desalination: An assessment of current desalination technologies, plants and capacity, Desalination, 495 (2020), p. 114633.
  * [18] K. El Kadi, I. Janajreh, and R. Hashaikeh, Numerical simulation and evaluation of spacer-filled direct contact membrane distillation module, Applied Water Science, 10 (2020), pp. 1–17.
  * [19] A. Ern and J.-L. Guermond, Applied mathematical sciences, Theory and Practice of Finite Elements, 159 (2004).
  * [20] D. Fletcher and D. Wiley, A computational fluids dynamics study of buoyancy effects in reverse osmosis, Journal of Membrane Science, 245 (2004), pp. 175–181.
  * [21] J. Freund and R. Stenberg, On weakly imposed boundary conditions for second order problems, in Proceedings of the Ninth Int. Conf. Finite Elements in Fluids, Venice, 1995, pp. 327–336.
  * [22] V. Girault and P.-A. Raviart, Finite Element Approximation of the Navier-Stokes equations, vol. 749, Springer Berlin, 1979.
  * [23] B. Gu, C. S. Adjiman, and X. Y. Xu, The effect of feed spacer geometry on membrane performance and concentration polarisation based on 3d CFD simulations, Journal of Membrane Science, 527 (2017), pp. 78–91.
  * [24] R. Ingram, Finite Element Approximation of nonsolenoidal, viscous flows around porous and solid obstacles, SIAM journal on numerical analysis, 49 (2011), pp. 491–520.
  * [25] M. Juntunen and R. Stenberg, Nitsche’s method for general boundary conditions, Mathematics of Computation, 78 (2009), pp. 1353–1374.
  * [26] M. Keshtkar, M. Eslami, and K. Jafarpur, A novel procedure for transient CFD modeling of basin solar stills: Coupling of species and energy equations, Desalination, 481 (2020), p. 114350.
  * [27] J. Kucera, Reverse osmosis: industrial processes and applications, John Wiley & Sons, 2015.
  * [28] W. J. Lee, Z. C. Ng, S. K. Hubadillah, P. S. Goh, W. J. Lau, M. Othman, A. F. Ismail, and N. Hilal, Fouling mitigation in forward osmosis and membrane distillation for desalination, Desalination, 480 (2020), p. 114338.
  * [29] M. Li, T. Bui, and S. Chao, Three-dimensional CFD analysis of hydrodynamics and concentration polarization in an industrial ro feed channel, Desalination, 397 (2016), pp. 194–204.
  * [30] Y. Liang, M. Chapman, G. F. Weihs, and n. D. Wiley, CFD modelling of electro-osmotic permeate flux enhancement on the feed side of a membrane module, Journal of Membrane Science, 470 (2014), pp. 378–388.
  * [31] J. Luo, M. Li, and Y. Heng, A hybrid modeling approach for optimal design of non-woven membrane channels in brackish water reverse osmosis process with high-throughput computation, Desalination, 489 (2020), p. 114463.
  * [32] E. Lyster, J. Au, R. Rallo, F. Giralt, and Y. Cohen, Coupled 3-d hydrodynamics and mass transfer analysis of mineral scaling-induced flux decline in a laboratory plate-and-frame reverse osmosis membrane module, Journal of Membrane Science, 339 (2009), pp. 39–48.
  * [33] S. Ma, L. Song, S. L. Ong, and W. J. Ng, A 2-d streamline upwind petrov/galerkin Finite Element model for concentration polarization in spiral wound reverse osmosis modules, Journal of Membrane Science, 244 (2004), pp. 129–139.
  * [34] K. G. Nayar, M. H. Sharqawy, L. D. Banchik, et al., Thermophysical properties of seawater: A review and new correlations that include pressure dependence, Desalination, 390 (2016), pp. 1–24.
  * [35] J. Nitsche, Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind, in Abhandlungen aus dem mathematischen Seminar der Universität Hamburg, vol. 36, Springer, 1971, pp. 9–15.
  * [36] A. Parasyris, M. Discacciati, and D. B. Das, Mathematical and numerical modelling of a circular cross-flow filtration module, Applied Mathematical Modelling, 80 (2020), pp. 84–98.
  * [37] V. Perfilov, V. Fila, and J. S. Marcano, A general predictive model for sweeping gas membrane distillation, Desalination, 443 (2018), pp. 285–306.
  * [38] C. Picioreanu, J. Vrouwenvelder, and M. Van Loosdrecht, Three-dimensional modeling of biofouling and fluid dynamics in feed spacer channels of membrane devices, Journal of Membrane Science, 345 (2009), pp. 340–354.
  * [39] C. Skuse, A. Gallego-Schmid, A. Azapagic, and P. Gorgojo, Can emerging membrane-based desalination technologies replace reverse osmosis?, Desalination, 500 (2021), p. 114844.
  * [40] J. M. Urquiza, A. Garon, and M.-I. Farinas, Weak imposition of the slip boundary condition on curved boundaries for Stokes flow, Journal of Computational Physics, 256 (2014), pp. 748–767.
  * [41] R. Verfürth, Finite Element Approximation on incompressible Navier-Stokes equations with slip boundary condition, Numerische Mathematik, 50 (1986), pp. 697–721.
  * [42] R. Verfürth, Finite Element Approximation of incompressible Navier-Stokes equations with slip boundary condition II, Numerische Mathematik, 59 (1991), pp. 615–636.
  * [43] G. F. Weihs and D. Wiley, CFD analysis of tracer response technique under cake-enhanced osmotic pressure, Journal of Membrane Science, 449 (2014), pp. 38–49.
  * [44] D. E. Wiley and D. F. Fletcher, Techniques for computational fluid dynamics modelling of flow in membrane channels, Journal of Membrane Science, 211 (2003), pp. 127–137.
  * [45] L. Xu, S. Xu, X. Wu, P. Wang, D. Jin, J. Hu, L. Li, L. Chen, Q. Leng, and D. Wu, Heat and mass transfer evaluation of air-gap diffusion distillation by $\varepsilon$-ntu method, Desalination, 478 (2020), p. 114281.