A Lyapunov-Sylvester numerical method for solving a reverse osmosis model

1.   Introduction

Water scarcity is a pressing global challenge that has prompted innovative solutions to ensure a sustainable and reliable freshwater supply. As the world's population continues to grow and traditional water sources face increasing stress due to climate change, population density, and industrialization, the need for alternative water resources is imperative. One such solution at the forefront of addressing this issue is water desalination. Water desalination is the process of removing salt and other impurities from seawater or brackish water to make it suitable for human consumption, agriculture, and industrial purposes.

Reverse osmosis (RO) is currently the most advanced technique for desalinating brackish water and seawater, due to its excellent versatility, ease of use, and high efficiency. In particular, the RO process is distinguished by its ability to handle a variety of plant designs and capacities, as well as its high salt rejection rate of 99% and recovery rate of up to 40%.

The process of osmosis occurs by applying pulsating pressure through the feed water that penetrates the fine pores of semi-permeable membranes. In practice, the reverse osmosis semi-permeable membranes are competitive, especially for seawater desalination. Recently, with the need to recycle many materials and wastes, reverse osmosis process has been implemented in municipal wastewater treatment.

Such a technique relies on mathematical models based on advection-diffusion PDEs to model the reverse osmosis system. Recall, that the concentration distribution process is always known as a time-dependent function depending on two-dimensional space variables, leading to a parabolic partial differential equation known as the advection-diffusion equation. The advection process describes heat transfer resulting from flow, while the diffusion phenomenon is due to the propagation of heat as a wave from higher-temperature water to lower-temperature water.

In the literature, many numerical methods have been applied to approximate solutions of the advection-diffusion or the reverse osmosis, such as one-dimensional discrete models based on Riemann sum numerical integration methods combined with the stretching transformations (^{[18,26,27,28]}), or fractional differential equations combined with wavelets (^[2,3,4,6,21]). Hadadian et al. ^[20] developed a mathematical simulation model for the RO process based on the mass, material, and energy balances considering the concentration polarization. Chen and Qin ^[10] developed a mathematical model to separate glucose from water through a reverse osmosis membrane. They optimized the model parameters, including the hydraulic permeability constant, reflection coefficient, and solute transport coefficient. Elnour et al. ^[15] developed an algorithm for designing RO systems. The algorithm is validated using operational data from a local plant and was found to be capable of simulating the actual plant operation with an average error of less than 5% for the majority of the system variables. Maure and Mungkasi ^[28] studied a mathematical equation in a reverse osmosis system of salt water at a given position and time. They converted the partial differential equation into ordinary differential equations with a transformation method. Noeiaghdam et al. ^[29] presented a novel mathematical method to find the optimal solution for the RO process. They applied the single exponential and double exponential decay to find the approximate solution of the RO.

In the present paper, numerical methods are developed to approximate the solutions of the continuous problem governing a reverse osmosis model by replacing the spacial partial derivatives with 2-dimensional finite-difference approximators to transform the continuous boundary-value problem into a linear algebraic system. The resulting method is analyzed for error, stability, solvability, and convergence. The method is shown to be more performant relative to tridiagonal systems. It transforms the continuous problem to the so-called Lyapunov-Sylvester equation of the form $AX+XB = C$. Thus, the problem then involves the study of the linear operator $\mathcal{L}(X) = AX+XB$ defined on an appropriate matrix space, for fixed matrices $A$, $B$, and $C$ depending on the discretization method and the problem parameters. The authors refer to the most recent work due to ^[23] for the inversion of such operators.

To check the approximate solution, we always adapt a small space step-size, or equivalently, a large discrete system, which hurts the increase in the algorithms' time of finding the solution. In this case, many methods such as Krylov's are applied to reduce the time of reckoning. In ^[30], for example, two ways have been investigated based on Krylov's method for solving a large-scale system of differential equations, by applying in one way an exponential projection technique, and in the other by dropping the original large-size system to a small-size one. In ^[31], a close problem and technique to the present paper are used related, especially, to the Sylvester matrix. More precisely, large-scale differential Sylvester matrix equations characterized by low-rank right-hand sides are investigated via an extended global Arnoldi process. The authors first approximated the exponential matrix in the exact solution via a global extended Krylov method, and secondly, they applied an extended Arnoldi algorithm to get a low-rank approximation of the Sylvester equation solution.

In the same direction, the authors in ^[32,33] investigated a large-scale generalized Sylvester equation by projecting the initial problem onto a block Krylov subspace to get a small-scale problem, which is then solved via an orthogonality Galerkin condition and Rosenbrock or the so-called backward differential formula method. The same authors applied both a modified global Arnoldi algorithm and an extended global Arnoldi algorithm for solving a large-scale Sylvester system of ordinary differential equations. The basic idea also involves a projection technique on the extended global Krylov or global Krylov subspaces. The main property of the systems investigated is the possibility of a non-invertible corresponding matrix, which requires the projection of a matrix exponential function to get solvable low-dimensional systems.

Again using Krylov's method in one of its variants, a numerical approach was developed in ^[34] to solve a large differential Stein matrix equation in a nonsymmetric case. The authors used an extended block Arnoldi algorithm and a backward differentiation formula. In ^[35], a Lyapunov differential equation was studied by applying an extended nonsymmetric block Lanczos method, which allows for the transformation of the original large-scale problem into a small-scale system using the extended block Krylov subspace projection and the extended nonsymmetric block Lanczos algorithm. Next, backward differentiation formula is applied to solve the low-dimensional Lyapunov equation. We conclude from these studies that the strange structures in dynamical systems need compatible tools to be well investigated, such as fractional derivatives, which, although being old concepts, have recently seen renewed interest in many fields.

Many problems and models involving PDEs are re-considered in the fractional calculus framework, where more adequate solutions to real-world problems in nature are provided accordingly. In ^[36], an exponential-cotangent derivative was involved in the well-known Riemann-Liouville approach to study the solutions of some fractional differential equations using the Laplace transform. The theoretical findings are applied to an SIR-type model. In ^[37], the problem of stability was investigated for a class of conformable linear systems in an infinite-dimensional form. $\alpha$-exponential stability, asymptotic stability, and uniform stability variants have been studied via conformable fractional derivatives and semigroup theory. Sadek et al, ^[38] studied a Lyapunov-Sylvester type differential equation using the Bernstein polynomials collocation method. As in the present paper, this technique allowed for the transformation of the original problem to a linear system of algebraic equations, which is then solved by using an iterative method.

In ^[39], a fascinating link between the frameworks of fractional calculus and the Lyapunov-Sylvester equation was given. The authors considered a class of fractional differential matrix equations and applied some fractional backward differentiation formulas techniques to numerically approximate the solutions of the original problem. An important problem when investigating PDEs concerns the controllability of the provided scheme or solution, and especially the observation of strange behaviors such as eventual chaotic or fractal behavior. In this way, in ^[40], the authors proposed a control technique for a non-homogeneous continuous-time fractal dynamical system derived from the Lyapunov equation.

The paper is organized as follows. In Section 2, some basic ideas from ^[23] on the invertibility of operators of the form $\mathcal{L}(X) = AX+XB$ are recalled. In Section 3, the discrete scheme is introduced. The continuous model describing the reverse osmosis process is converted into a discrete algebraic system involving the Lyapunov-Sylvester equation in a discrete form. The solvability of the discrete problem obtained in Section 3 is investigated in Section 4. Section 5 is devoted to the convergence, consistency, and stability of the method. Section 6 is concerned with the experimental results and their discussion. Section 7 is an appendix devoted to some preliminary results, some comments on existing models, and some mathematical proofs. Section 8 is a conclusion.

2.   On the Lyapunov-Sylvester operator

The resolution of the so-called Lyapunov equation

where $A$ and $C$ are given matrices and $X$ is an unknown one, is often of interest, and has been the object of some studies until the 60's ^[22]. However, these studies have been characterized by restrictive aspects. For example, in ^[22], a Kronecker product-based method was developed to resolve Eq (2.1). The method consisted of a recursive scheme using the well-known Caylay result, and it gives explicit solutions in 2 and 3-dimensional cases. Such cases are not complicated and a verification or a computation can be proved with simpler and more direct methods than the one in ^[22]. The equation has been re-considered by many authors due to its relation and application in PDEs, the theory of stability, and dynamical systems. Kohaupt has focused on general forms of the cited equation. The author in ^[23] considered the matrix eigenvalue problem

A large study of this problem has been provided, and special cases, such as $B = A^T$, $B = A^*$, $A$, and $B$ diagonalizable, etc., have been discussed. The main result of ^[23] is presented in the following theorem.

Theorem 2.1. ^[23] Let $A\in\mathbb{C}^{n\times n}$ and $B\in\mathbb{C}^{m\times m}$ both be diagonalizable. Let $\alpha_i$ be the eigenvalues of $A$ and $u_i$ be the associated right eigenvectors for $i = 1, ..., n$. Let $\beta_i$ be the eigenvalues of $B$ and $v_i$ be the associated left eigenvectors for $i = 1, ..., m$, i.e.,

Then,

are the $mn$ eigenvalues of (2.2) and

are the associated eigenmatrices.

3.   The mathematical model and its discretization

Osmosis is essentially a natural process that occurs when water molecules migrate from a solution with low solute concentration (low osmotic pressure) to a solution with high solute concentration (high osmotic pressure) via a semipermeable membrane (Figure 1(a)). The osmosis process is said to be at an equilibrium condition when the chemical potentials across the membrane are equal (Figure 1(b)). The flow of water molecules can be halted or reversed by applying external pressure to the solution (the feed solution) that has a higher concentration. If the applied pressure differential is greater than the osmotic pressure difference across the membrane, water molecules are compelled to flow against the direction of the natural osmosis phenomenon. Reverse osmosis is the process occurring in this case, s illustrated by Figure 1(c). The symbol $P$ designates the process or external pressure, while the symbol $\pi$ is the osmotic pressure.

The feed water stream separates into two streams in a continuous RO process, as represented in Figure 2. Water molecules that have passed through the membrane constitute the first stream. Permeate or product water is the name of this low-solute stream. A smaller number of water molecules and the rejected solutes are the constituents of the second stream. This stream is known as brine, concentrate, or reject because it contains more solutes than feed. The feed channel is characterized by three parameters: feed flow rate $F_d$, concentration $C_d$, and pressure $P_d$ (^{[24,25,41,42,43,44,45]}).

The permeate channel and the feed channel are located on opposite sides of the membrane. The feed channel is filled with feed solution at one end, while the rejected solution emerges at the other. The fluid flow in the feed channel is parallel to the $x$-coordinate, while in the permeate channel, it is parallel to the $y$-coordinate. Therefore, a two-dimensional advection-diffusion equation will be investigated as a governing model for the reverse osmosis phenomenon, which is expressed by the following time-dependent 2D-model,

where $(x, y)$ is the pair of $\mathbb{R}^2$-space variables, and $C = C(t, x, y)$ is the concentration of salt solution in a semipermeable membrane at point $(x, y)$ and the instant $t$. The quantity $\dfrac{\partial C}{\partial t}$ is the time derivative of $C(t, x, y)$, $\dfrac{\partial C}{\partial x}$ is the derivative of $C(t, x, y)$ relative to the first space variable $x$, and $\Delta C = \dfrac{\partial^2C}{\partial y^2}+\dfrac{\partial^2C}{\partial x^2}$ is the Laplacian of $C(t, x, y)$ resulting from the Laplace operator applied to it. The real number parameter $q$ represents the water flow rates in semipermeable distribution, $D$ is the salt diffusivity in water, $\delta$ is the distance from the semipermeable boundary to the center of the channel, $c_d$ is the concentration away from semipermeable membranes, $\nu_0$ is the horizontal velocity measured at a distance $\delta$ from the semipermeable boundary, $\alpha = \dfrac{D\delta}{\nu_0}$. The time-dependent two-dimensional advection-diffusion equation describes the concentration of a salt solution. The schematic illustration given in Figure 3 below explains the phenomenon, where $v(y)$ is the speed of the fluid flowing through the duct, parallel to the $x$-axis. The domain $\Omega = \Omega_1\times\Omega_2$ is a rectangular strip as in Figure 3. Here, we recall that the components $\Omega_1$ and $\Omega_2$ are in general different, (for example, $\Omega_1 = [0, L_x]$ and $\Omega_2 = [0, L_y]$). Numerically, we may take two different steps $h_x = \frac{L_x}{N}$ in $x$ and $h_y = \frac{L_y}{N}$ in $y$ with the same number of discretization points $N$, or we may also take two different sizes $N_x$ and $N_y$ leading to $h_x = \frac{L_x}{N_x}$ in $x$ and $h_y = \frac{L_y}{N_y}$ in $y$. However, it is well known in numerical analysis that these facts do not have a great effect on the approximations. For this, we restrict ourselves to an 'artificial' square region for the discretization as well as the numerical developments.

From a physical or natural interpretation point of view, we take a semi-infinite region in the $\mathbb{R}^2_+$ quadrant. At $t = 0$, the system is assumed to be still at rest, which means that the concentration of salts is still the same as the seawater to be filtered, which in turn leads to a natural initial condition $C(0, x, y) = C_0$. The notation or condition $x = \infty$ means that we are 'very' far from the semi-permeable membrane at $y = 0$, and the 'artificial' condition $y = \infty$ means that we are 'very' far from the membrane. In our case, we assume that the concentration away from the membrane is the same as that before flowing in the filter. In the discrete and numeric validation, we symbolize this as

The numerical method is based on a simple classical 2-dimensional finite difference scheme, where the differential operators are estimated as follows. Let $N\in\mathbb{N}$ and $s$ satisfying $Ns = L_1-L_0$ be a space step. For the reasons raised above, without loss of generality, we take in our numerical solutions $L_0 = 0$ and $L_1 = L$. Next, consider the net space

Let $u = \left(u_{j, m}\right)$ be a mesh function defined on $\Omega_s$. Next, we introduce the operators

and

The first-order derivatives in $x$ and $y$ are approximated by

respectively, while the second-order operators in $x$ and $y$ are

and, finally, the Laplace operator is

Let $l = \Delta\, t$ be the time step and $t_n = t_0+nl$, $n\in\mathbb{N}$ be the discrete time grid. For $(j, m)\in{I} = \{0, 1, \dots, N\}^2$ and $n\geq0$, $C_{j, m}^n$ will be the net function $C(t_n, x_j, y_m)$ (also used to designate the numerical solution). The following discrete approximations will be applied for the time derivatives,

Our numerical scheme is based on a calibrating method which consists of a barycenter ponderations according to the time for the final space derivative to be applied in the approximation model. Precisely, we set

where the weight parameter $\mu\in[0, 1]$.

The calibration method is usually used to allow for the best numerical approximation of the real solution. As this solution is, in fact, unknown, we did not know a priori which time position is the closer to it. The variation of the weights allows us to investigate this problem. See for example ^{[7,8,9,11,12]}

As a result, the approximate solution $C_{j, m}^n$ is defined as a solution to the discrete problem

For the rest of the paper, for $n\geq0$, denote by $C^n = \Bigl(C_{j, m}^n\Bigr)_{0\leq j, m\leq N}$ the matrix with coefficients $C_{j, m}^n$.

4.   Solvability of the numerical model

The first result in this paper concerns the solvability of the system (3.2) above, and is stated as follows.

Theorem 4.1. System (3.2) is uniquely solvable.

Proof. Denote

From system (3.2), it is obtained that

Now, denote

and, similarly,

and regrouping the terms in $n+1$, $n$, and $n-1$, we obtain

Otherwise, we may write this system as the inner product of vectors as

where the upper script $[.]^T$ designates the transpose. Next, let $A_m$, $B_m$, and $R$ be the tridiagonal matrices of size $(N+1, N+1)$ whom coefficients are, respectively,

and

System (4.2) above may be transformed into a Lyapunov-Sylvester form as

where $I = I_{N+1}$ is the identity matrix of size $N+1$. Now, if we denote for two matrices $A, B\in M_{N+1}(\mathbb{C})$, the Lyapunov-Sylvester operator

we obtain a dynamic Lyapunov-Sylvester equation of the form

Notice that the matrices $(I-A_m)$ and $\lambda R$ have different eigenvalues (see Appendix 7.1, Lemma 7.1). As a consequence, due to Theorem 2.1 (see also ^[23]), the linear Lyapunov-Sylvester operator $\mathcal{L}_{I-A_m, -\lambda R}$ is invertible. So, the result follows.

5.   Analysis of the consistency, stability, and convergence

5.1.   Consistency of the numerical model

To show the consistency of the numerical scheme, we use the local truncation error due to (3.2). The result is given in the following lemma.

Lemma 5.1. The numerical scheme is consistent at the order $(l^2+s^2)$.

Proof. To show this lemma, we will compute the principal part due to the truncation error applied to (3.2). To do this, we start by evaluating the necessary elements of the discrete scheme. So, by assuming that $C$ is sufficiently regular, standard computations yield that

Similarly, we have

By canceling the terms in $s$ (taking $s = 0$), we obtain

Now, by replacing $x$ with $y$ in Eqs (5.1) and (5.2) above, we deduce that

Similarly, we have

Now, by canceling $l$ (taking $l = 0$) in all the equations above, we get

Now, as previously, by replacing $l$ with $-l$ in the equations above, we get

In the same way, we get

Similarly, by canceling $s$ in the last equation, we get

Now, by interchanging the role of $j$ and $m$ in (5.7), we obtain

Now, by replacing $l$ with $-l$ in (5.10), we obtain

All the equations above permit us to deduce the following principal part of the truncation error as

which tends to 0 as $(l, s)$ tends to 0. This means that the method is consistent.

5.2.   Stability of the numerical model

To examine the stability of the numerical method, we will apply the Von-Neumann technique by investigating the impact of the scheme on an isolated Fourier mode.

Denote $C_{j, m}^n$ as the solution of the numerical scheme, and $\check{C}_{j, m}^n$ as the numerical value obtained subject to computer round-off errors. Next, denote

where $\psi = \psi_1+i\psi_2$ is a complex number, $\varphi_x$ and $\varphi_y$ are real, and $i^2 = -1$. By applying the Von-Neumann criterion for stability, we shall show that

Indeed, from (3.2) or its developed form (4.1), after canceling $X^{n-1}Y^{j-1}Z^{m-1}$, and regrouping the terms in $X^2$ and $X$, it holds that

Now, notice that

Equation (5.13) may be written in a simplified form

where

To establish the Von-Neumann criterion for stability, it suffices to show that

Let us now examine the difference $|B(Y, Z)|^2-|A(Y, Z)|^2$. Careful computations yield that

For $0\leq\mu < \dfrac{1}{3}$, from (5.17), we get

By denoting

we obtain

The left-hand term reads

Next, for $\dfrac{1}{3} < \mu < 1$, from (5.17), we get

By denoting

we obtain

which, similar to the previous case, gives

Finally, for $\mu\in\{\dfrac{1}{3}, 1\}$, Eq (5.17) yields that

This leads to

Thus, the stability follows as a result of (5.18), (5.19), and (5.20).

5.3.   Convergence of the numerical model

To show the convergence of the numerical scheme, here we apply the Fourier mode concept combined with the Vonn-Neuman criterion to establish a control on the time and space steps. Indeed, denote

The intersection of (5.18), (5.19), and (5.20) gives

Taking $l = s^{\eta+2}$ for some $\eta > 0$ small enough (which is always possible), we get

Thus, convergence holds.

6.   Numerical results and simulations

In this section, we will test the mathematical model above for salt rejection, and mainly for inorganic salts such as NaCl or Na2SO4. We take the parameter values $h = 10^{-3}$, $D = 10^{-9}$, $\nu_0 = 10^{-3}$, and $q = 10^{-3}$. The numerical tests and simulations are conducted using MATLAB R2018a on a Dell Latitude 5530, Windows 10 64-bit Operating System, x64-based processor, 12th Gen Intel(R) Core(TM) i7-1255U 1.70 GHz, and 16.0 GB RAM.

In existing models based on the advection-diffusion equation, many constraints always applied, such as the dependence of the liquid (as well as the concentration of salts) flux speed on the perpendicular direction of the membrane (the $x$-axis here), and thus the concentration acceleration in the same direction. In this case, the model obtained (reviewed next in Appendix 7.2) may be converted via the well-known stretching or self-similar transformations into an ODE, which may be solved even for exact solutions.

The idea may be explained according to the speed of the liquid, the proportions of the different pollutants, such as salts, and due to many factors related to the filtration/osmosis system (the nature or the homogeneity of the membrane, and so on), which affects the result(s) of the osmosis. This task resembles the water in valleys and rivers, as it brings with it many types of pollutants and materials. Relative to the speed of feed water and according to the quantity, the percentages, and the quality of materials in the seawater, several sediments are formed on both sides of the membrane. Over time, these sediments cause the stream itself to narrow. Thus, we are faced with a phenomenon of interconnected parameters and dimensions, where the parameters, factors, and dimensions, that are initially neglected, excluded, and considered ineffective, come back to influence the model and the entire phenomenon.

Figure 4 illustrates the variation of the concentration $C$ of the salt on the semi-permeable boundary with the distance $x$ along the plate. Notice easily that this curve joins well with existing works, such as ^[18] or its variants ^[26,27,28].

We notice from Figure 4 that, when relaxing the parameter $t$, and thus getting a stationary solution $C(t, x, y) = C(x, y)$, we come back to the classical equilibrium discussed in ^[18]. In this case, at the boundary of the membrane $y = 0$, the phenomenon of the precipitation of the salts may persist, even being small, but as time passes, we notice effectively that this small amount of salts may have an impact on the whole process. Moreover, model (7.1) may lead to ambiguities relative to the boundary conditions assumed to be satisfied there, such as $C(0, y) = C(x, \infty) = C_0$, which cannot be fulfilled via solution (7.4). In such a case, simple computations permit us to get $C(x, y)\sim C_0y$, whenever $x\rightarrow0$ or $y\rightarrow \infty$. The last approximations may contradict the physical reality of the whole process. These remarks mean that, for model (7.1) to approach reality, some artificial conditions should be added, which may again affect the compatibility with the real process. Moreover, this again raises the impact of neglecting some parameters or some parts in the real 2D time-dependent model, and thus motivates our approach in considering the whole 2D time-dependent model (3.1). Examinations of the complete model have led to Figures 5, 6, and 7.

In Figure 5, the time-wise behavior of the salt concentration is illustrated at some interior points $(x, y)$, and compared to the concentration at the boundary point $(0, 0)$, where it is initially $c_0$. We notice that the RO process works well at interior points (especially far from the boundary $y = 0$), as the concentration of the salt decreases naturally until reaching a value less than half in a reduced time interval.

In Figures 6 and 7, the profile of the concentration of salts is investigated for a fixed time instants when the process is active. Notice from these figures that the concentration is somehow perturbed at the boundaries and takes a quasi-linear asymptotic profile at the interior of the channel. A logical explanation of this behavior is that, effectively, it is not completely true that the concentration of salts is omitted from the liquid, but, in contrast, the effect of the perpendicular direction to the flow of the osmosis is real, and it induces precipitation of the salts along the boundary of the membrane, which becomes considerable over time, and thus has to be taken into consideration in any model. Taking it into consideration in our opinion needs a whole revision of the system, such as the membranes, and their components, homogeneity, heterogeneity, construction materials, temperature, and also the permeability flow, which may not be constants, and may depend on the salts' precipitation along the process.

Besides, with the discovery and rapid development of modern technologies, such as nano-scales and nanomaterials, the variations in scales and levels that were not observable in classical technology are now well detected by modern instruments. This means that the ignorance of many parameters in the equations or models governing physical and/or chemical phenomena such as reverse osmosis no longer has any motivation or justification, as a maybe major part of these phenomena takes place in macro and sometimes nano levels. This motivates our involvement of both the variations according to the time and space variables in the model, especially the second-order space derivatives which explain, even if small, the acceleration of the process propagation in the space, and the time derivative which reflects the rate of the process in time.

7.   Appendix

7.1.   Appendix A: LS equation invertibility

Lemma 7.1. The matrices $A$ and $\lambda R$ above have different eigenvalues.

Proof. For the sake of simplicity and clearness, we will develop the first two special cases of $2\times2$- and $3\times3$-matrices. In the first, we have to investigate the intersection of the spectra of the matrices

It is easy to see that the eigenvalues of $\lambda R$ are $\lambda$ and $3\lambda$. But, notice that both $det((1-\lambda)I-A_m)\not = 0$, and $det((1-3\lambda)I-A_m)\not = 0$. Hence, $\lambda$ and $3\lambda$ are not eigenvalues of $(I-A_m)$. Therefore, Theorem 2.1 (see also ^[23]) implies that the linear Lyapunov-Sylvester operator $\mathcal{L}_{I-A_m, -\lambda R}$ is invertible.

For the $3\times3$-case, we get

It is easy to see that $\lambda R$ has the only eignevalues $2\lambda$, $(2-\sqrt2)\lambda$ and $(2+\sqrt2)\lambda$. It is easy to check that $det(A-xI)\not = 0$, for $x\in\{2\lambda, (2-\sqrt2)\lambda, (2+\sqrt2)\lambda\}$. Thus, the lemma is proved.

7.2.   Appendix B: Critics about the existing model

In the existing works, as criticized in the introduction, approximate models have been investigated, characterized by two major drawbacks. First, such models ignore or neglect the role of some variables, based on assumptions on the minerals' concentration, such as, salts, which is assumed to be rapidly varying according to the $y$-direction, compared to the $x$-one, allowing one to neglect the diffusion in the $x$-direction. Second, such models assume that an equilibrium in time exists, and thus a stationary variant is investigated. In this case, we get the following governing equation to the reverse osmosis model:

However, this model may be solved for a general solution by utilizing the classical change of variables $\zeta=\frac{y}{\sqrt[3]{x}}$, which, by denoting $C(x, y) = \sqrt[3]{x}\varphi(\zeta)$, permits one to transform Eq (7.1) into

Notice that this last equation has a solution $\varphi_0(\zeta) = \zeta$. To find a general solution, denoting $\varphi(\zeta) = \zeta\phi(\zeta)$, and substituting in (7.2), we get

which is a classical ordinary differential equation whose solution is

The general solutions of (7.2) reads

where $K_0$ and $K_1$ are constants.

8.   Conclusion and some future directions

The aim of the present paper is to develop a numerical method to solve a reverse osmosis model, where the governing (1+2)-dimensional diffusion-convection PDE due to the salt concentration in the osmosis membrane is converted into a Sylvester equation that is proved to be uniquely solvable, stable, consistent and convergent. The numerical scheme developed was validated with experimental data from the literature, and via some numerical simulations. The objective of the present paper is therefore twofold. As a result, many future directions may be raised.

Related to the Lyapunov-Sylvester algebraic system, an interesting direction that interested us is to link fractional derivatives to convection-diffusion evolutionary problems by considering fractional derivatives operators for the time and space variables in the PDE. With the LS equation, in algebra and especially when the operators $\mathcal{L}_{A, B}$ are defined for $A, B$ belongs to some ring $R$ by

we often call them generalized derivations on the ring $R$. These have been studied in the context of algebras in some normed spaces. Indeed, a generalized derivation of an algebra $R$ is any map of the form $\mathcal{L}_{A, B}$ where $A$ and $B$ are fixed elements in $R$. In the theory of operator algebras, they are considered as an important class of so-called elementary operators $\mathcal{L}(X) = \sum\limits_{i = 1}^nA_iXB_i$, where, as previously, the $A_i$'s and the $B_i$'s are elements of the algebra $R$. Remark that the operator $\mathcal{L}_{A, B}$ satisfies

where $\mathcal{D}_{B}(Y) = YB-BY$ is the so-called inner derivation. Remark that

which looks like the well-known product derivation rule on function spaces $(fg)' = f'g+fg'$. Hence, the nomination of derivation on algebraic structures. So, it is questionable that an integration operator $\mathcal{I}_{A, B}$ could or could not be defined so that one has for any $X$ and $Y$ in $R$,

Taking the simple example, $A = \left(\begin{array}{cc} 0 & 1 \\ 2 & 1 \\ \end{array} \right)$ and $B = \left(\begin{array}{cc} 1 & 1 \\ 0 & 2 \\ \end{array} \right)$, the associated Lyapunov-Sylvester operator $\mathcal{L}_{A, B}$ is invertible for the ring $M_2(\mathbb{R})$, but non-invertible on $M_2(\mathbb{Z}/3\mathbb{Z})$. A second example, for $A = \left(\begin{array}{cc} 1 & 2 \\ 2 & 1 \\ \end{array} \right)$ and $B = \left(\begin{array}{cc} 1 & 1 \\ 2 & 1 \\ \end{array} \right)$, leads to an invertible operator on $M_{2}(\mathbb{R})$, but its inverse has no longer the same form. A general form that may be expected is

So, here also, what could be the minimum values of the parameters $m$ and $n$ given the operator $\mathcal{L}^{-1}$? Could we expect a signature-type $(n, m)$ for these operators?

We now return to the physical phenomenon subject of the application, which deals with the RO model. In this context, water represents an essential factor for the existence of all creatures in this universe. Even wastewater or brackish water constitute suitable environments for the life of several types of organisms, such as bacteria, fishes, marine plants, and so on. This is in addition to its basic role in environmental balance, such as maintaining temperatures and climate.

Clean, fresh, and/or potable water is the most important type of water and liquids overall, due to its essential role in the continuation of the life of the most important type of beings in the universe, the human, and such a role is irreplaceable by other liquids.

With the increase in demand for potable water due to frequent use, accompanied by a shortage of natural sources such as rain, water napes, and many other reasons, humans have resorted to inventing several ways to benefit from other water sources, such as rivers and seas, thus confronting the problem of desalinating water and extracting impurities from it to make it suitable for drinking in particular.

The present paper investigates one of the methods applied to extract fresh water from seawater using the reverse osmosis process. A numerical method is developed to solve a two-dimensional time-dependent original mathematical model of a reverse osmosis process. The governing PDE is transformed to a Lyapunov-Sylvester algebraic system, (which has been solved in its original geometric space 3D-dimension Euclidean space), without transforming it into a 1-dimensional space as in the majority of existing works dealing with higher-dimensional PDEs. The system is investigated for solvability, stability, consistency, and convergence. Besides, experimental results have been provided based on simulation results, and eventually compared with existing studies. The findings showed that the RO may be suitable for the aim of fresh water extraction, whenever it is accompanied by many improvements, related to the space parameters for the used models, such as the dimensions, the time factor, the materials used for the physical RO system, the permeability of the membrane walls, and also the way of controlling the system during its progress.

On another side, RO, or more generally desalination tools, may also be linked to energy consumption, which in turn needs optimal models to reduce such consumption. The membrane module configuration, for example, could balance the flux and, hence, minimize energy consumption. New technologies using nanofiltration for pressure-driven membrane processes are linked to both the diffusion and the convection through the pores and/or charged membrane. To understand more and improve these facts, the scientific community needs to involve related parameters into the mathematical models such as the present convection-diffusion one. Optimality based on the membrane properties, such as porosity and geometry, is still limited in comparison to the use of substances' properties (^{[1,5,13,14,16,17,19]}.

Author contributions

S. Helali: Supervision, Conceptualization, Investigation, Writing–original draft, Writing–review and editing, Funding acquisition; A. Ben Mabrouk: Conceptualization, Methodology, Formal analysis, Investigation, Writing–original draft, Writing–review and editing, Visualization; M. Rashad: Conceptualization, Formal Analysis; N. Bel Haj Ali: Conceptualization, Formal analysis; M. A. Alsharif: Formal analysis; M. A. Alanazi: Writing–review and editing, E. M. Al-Ali: Methodology, Writing–review and editing; L. A. Alharbi: Writing–review and editing; M. S. Mustafa: Writing–review and editing. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank and extend their appreciation to the deanship of scientific research at the University of Tabuk, KSA, and the Deputyship for Research & Innovation, Ministry of Education, Saudi Arabia for funding this research work through the project number (0097-1443-S).

Conflict of interest

The author declares here that no conflicts of interest in this paper.