Soft sensing of system parameters in membrane distillation

Abstract

Various examples of methods and systems are provided for soft sensing of system parameters in membrane distillation (MD). In one example, a system includes a MD module comprising a feed side and a permeate side separated by a membrane boundary layer; and processing circuitry configured to estimate feed solution temperatures and permeate solution temperatures of the MD module using monitored outlet temperatures of the feed side and the permeate side. In another example, a method includes monitoring outlet temperatures of a feed side and a permeate side of a MD module to determine a current feed outlet temperature and a current permeate outlet temperature; and determining a plurality of estimated temperature states of a membrane boundary layer separating the feed side and the permeate side of the MD module using the current feed outlet temperature and the current permeate outlet temperature.

Description

BACKGROUND

The rapid growth of population around the world has increased the demand on freshwater and put the scarce natural freshwater resources under great pressure. Many water stressed countries rely on desalination plants to meet the ever increasing freshwater demand. The most common water desalination techniques are reverse osmosis (RO), multi-stage flash distillation (MSF), and multi-effect distillation (MED). However, most of these conventional desalination methods rely on fossil fuels to operate, which results in pollution of the environment by emissions and concentrated rejected brine. Membrane distillation (MD) offers a more sustainable and environment-friendly alternative. MD is a thermal separation process that utilizes a specific type of membranes for different applications, such as water desalination. Membranes are selected to be microporous, non-wetted, and highly hydrophobic. These properties allow the transfer of only water vapor or other volatile molecules through the membrane dry pores and so produce clean water. The separation is achieved under the validation of vapor-liquid equilibrium principle for molecules separation. MD is an attractive technology, with promising potentials to take the lead over the other existing water desalination methods. These potentials contribute toward promoting the MD process to be a cost-effective desalination process that is capable of being powered with renewable energy sources, such as geothermal, solar energy or low grade waste heat.

SUMMARY

Embodiments of the present disclosure are related to soft sensing of system parameters in membrane distillation (MD).

In one embodiment, among others, a system, comprises a MD module comprising a feed side and a permeate side separated by a membrane boundary layer; and processing circuitry configured to estimate feed solution temperatures and permeate solution temperatures of the MD module based at least in part upon monitored outlet temperatures of the feed side and the permeate side. In one or more aspects of these embodiments, the processing circuitry can be configured to estimate a membrane mass transfer coefficient of the MD module. The feed solution temperatures and permeate solution temperatures can include boundary temperatures along the membrane boundary layer.

In one or more aspects of these embodiments, the estimation of the feed solution temperatures and permeate solution temperatures can be based on an Advection-Diffusion Equation (ADE) model of feed and permeate solutions in the MD module. The ADE model can be a two-dimensional ADE model including dimensions of the membrane. The ADE model can comprise a plurality of grid points extending across the feed side and the permeate side of the MD module. The feed solution temperatures and permeate solution temperatures can include temperatures corresponding to at least a portion of the plurality of grid points.

In one or more aspects of these embodiments, the system can comprise temperature sensors configured to monitor the outlet temperatures of the feed side and the permeate side of the MD module. The MD module can be a counter current direct-contact MD module. The processing circuitry can be configured to control operation of the MD module based at least in part upon the feed solution temperatures and permeate solution temperatures.

In another embodiment, a method comprises monitoring outlet temperatures of a feed side and a permeate side of a MD module to determine a current feed outlet temperature and a current permeate outlet temperature; and determining a plurality of estimated temperature states of a membrane boundary layer separating the feed side and the permeate side of the MD module based at least in part upon the current feed outlet temperature and the current permeate outlet temperature. In one or more aspects of these embodiments, the method can comprise estimating a membrane mass transfer coefficient of the MD module.

In one or more aspects of these embodiments, the method can comprise adjusting characteristics of the feed solution or the permeate solution supplied to the MD module based at least in part upon the current feed outlet temperature and the current permeate outlet temperature. The estimated temperature states can be based on an Advection-Diffusion Equation (ADE) model of feed and permeate solutions in the MD module. The estimated temperature states can include temperatures corresponding to grid points extending across the feed side and the permeate side of the MD module.

Other systems, methods, features, and advantages of the present disclosure will be or become apparent to one with skill in the art upon examination of the following drawings and detailed description. It is intended that all such additional systems, methods, features, and advantages be included within this description, be within the scope of the present disclosure, and be protected by the accompanying claims. In addition, all optional and preferred features and modifications of the described embodiments are usable in all aspects of the disclosure taught herein. Furthermore, the individual features of the dependent claims, as well as all optional and preferred features and modifications of the described embodiments are combinable and interchangeable with one another.

BRIEF DESCRIPTION OF THE DRAWINGS

Many aspects of the present disclosure can be better understood with reference to the following drawings. The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views.

FIG. 1 is a schematic diagram of an example of a counter current direct-contact membrane distillation (DCMD) module in accordance with various embodiments of the present disclosure.

FIG. 2 is a schematic diagram illustrating an example of a discretization grid for the feed and permeate sides of the DCMD module of FIG. 1 in accordance with various embodiments of the present disclosure.

FIG. 3A is a schematic diagram illustrating an example of DCMD experimental set up in accordance with various embodiments of the present disclosure.

FIGS. 3B and 3C include tables of membrane and feed water characteristics for the DCMD experimental set up of FIG. 3A in accordance with various embodiments of the present disclosure.

FIGS. 4 and 5 are plots comparing experimental and modeled results in accordance with various embodiments of the present disclosure.

FIG. 6 illustrates the effect of membrane pore size variation on the membrane flux of the DCMD module of FIG. 1 in accordance with various embodiments of the present disclosure.

FIG. 7 includes a table of membrane distillation parameters used in the model analysis in accordance with various embodiments of the present disclosure.

FIGS. 8-13 are plots illustrating the performance of soft sensing estimation using the model for the DCMD module of FIG. 1 in accordance with various embodiments of the present disclosure.

FIGS. 14 and 15 are schematic diagrams of examples of DCMD modules in accordance with various embodiments of the present disclosure.

FIG. 16 is a table illustrating electrical thermal analogue elements in accordance with various embodiments of the present disclosure.

FIG. 17 is a schematic diagram illustrating an electrical thermal analogue network of the n^th cell of the DCMD module of FIG. 15 in accordance with various embodiments of the present disclosure.

FIG. 18 is a table of elements of the electrical thermal analogue network of FIG. 16 in accordance with various embodiments of the present disclosure.

FIG. 19 is a schematic diagram illustrating an example of a completed electrical thermal analogue network of the DCMD module in accordance with various embodiments of the present disclosure.

FIG. 20 is a graphical representation of an experimental lab scale setup in accordance with various embodiments of the present disclosure.

FIG. 21 is a table of DCMD module parameters of the experimental lab scale setup of FIG. 20 in accordance with various embodiments of the present disclosure.

FIGS. 22-27 are plots illustrating the performance of soft sensing estimation using the model for the DCMD module of FIG. 14 in accordance with various embodiments of the present disclosure.

FIGS. 28A-28E illustrate simulation results of a descriptor adaptive observer in accordance with various embodiments of the present disclosure.

FIGS. 29A and 29B are schematic diagrams illustrating membrane fouling and an electrical thermal analogue network modeling the membrane fouling in accordance with various embodiments of the present disclosure.

FIGS. 30A-30C illustrate simulation results of membrane fouling detection in accordance with various embodiments of the present disclosure.

FIGS. 31A and 31B illustrates examples of convergence of the estimated states and error using the nonlinear Lyapunov-based observer in accordance with various embodiments of the present disclosure.

FIG. 32 is a schematic diagram illustrating an example of a nonlinear observer-based Lyapunov controller in accordance with various embodiments of the present disclosure.

FIG. 33A is a table of membrane parameters used in the controller analysis in accordance with various embodiments of the present disclosure.

FIGS. 33B-33F illustrate simulation results of the controller in accordance with various embodiments of the present disclosure.

FIG. 34 is a schematic block diagram that illustrates an example of processing circuitry employed by the DCMD module of FIGS. 1 and 14 in accordance with various embodiments of the present disclosure.

DETAILED DESCRIPTION

Disclosed herein are various embodiments of methods related to soft sensing of system parameters in membrane distillation (MD), which can include estimation of membrane characteristics such as temperature, mass coefficient transfer, etc. Reference will now be made in detail to the description of the embodiments as illustrated in the drawings, wherein like reference numbers indicate like parts throughout the several views.

Membrane distillation (MD) is a hybrid process in which both thermal and membrane separation techniques are combined. In this process, hot salty stream passes along one side (the feed side) of a microporous hydrophobic membrane, while the other side (the permeate side) of the membrane is kept at a lower temperature (coolant). Due to this temperature difference, water evaporates at the feed-membrane interface, and only water vapor passes through the dry pores of a hydrophobic membrane, driven by the differential saturation vapor pressure across the feed and the permeate sides of the membrane. Water vapor then condenses on the permeate-membrane interface. Since only water vapor is allowed through the hydrophobic membrane, MD has a 100% theoretical rejection rate of ions. Several configurations have been proposed for MD, each with its own pros and cons. One configuration is direct contact membrane distillation (DCMD), where both the hot feed and cold permeate streams are in direct contact with the membrane. Another configuration is the air gap membrane distillation (AGMD), in which the permeate stream is separated from the membrane by a stagnant air gap. A third setup is vacuum membrane distillation (VMD), where a vacuum is applied at the permeate side in order to enhance the flux transfer. Another configuration is material gap membrane distillation (MGMD), which combines DCMD and AGMD.

Unlike thermal-based conventional desalination methods, MD does not need intensive heat which means it can be integrated with solar thermal energy and waste heat sources. Since MD operates at a lower hydraulic pressure, it is less susceptible to scaling and fouling which is a general drawback of pressure-driven membrane based desalination techniques, like RO. Moreover, the general operational concept of MD has less demanding membrane characteristics and properties. All these features make MD ideal for sustainable water desalination in remote areas.

To capitalize on the MD advantages, the MD water desalination process has to be well understood, modeled and then optimized for maximum distilled water flux. The membrane distillation (MD) process starts when the solution in the feed side, after being heated to a desired temperature, evaporates along a membrane boundary layer. The vapor passes through the membrane dry pores, driven by the difference of pressure between the two sides of the membrane maintained by the difference of temperature, and gets condensed in the permeate side providing clean water. The flux of the MD process depends on the vapor gradient pressure between the two sides of the membrane, in addition to the characteristics of the module and the membrane.

The four common configurations include air-gap membrane distillation (AGMD), vacuum membrane distillation (VMD), sweeping-gas membrane distillation (SGMD), and direct-contact membrane distillation (DCMD), although other configurations have also been developed such as water gap MD. Condensation which occurs in the permeate side is the main difference among these configurations. Several applications have targeted the use of DCMD technology, where dewatering like desalination, waste water treatment and concentration of liquid solutions are among the top applications.

Many theoretical studies have been carried out to model and predict the performance of the MD process. Such models contribute to the deep understanding of the involved parameters and thus enhance the performance of the process. The models are mainly empirical and hold the advantages of relating the involved process parameters with each other to show the effect of the overall performance. The majority of the models consider the steady-state performance more important than the transient performance. The kinetic theory of gases through porous media was adopted in the models and the process performance depending on the selected MD configuration was studied. These models lacked the compatibility to be powered by intermittent energy supply, such as the use of solar energy source as they could not accommodate their unsteady nature. Notwithstanding, the use of heating storage tanks was proposed as a remedy to the intermittent unsteadiness effect. However, they introduced additional problems when they missed enough attention to some challenging drawbacks of the heating storage tanks, such as the incompatibility to the space requirement of some applications, and the extra expense for thermal insulation and maintenance, and the relatively low storage capacity they provide, to list a few.

In recent years there has been some effort to model the MD process and understand the effect of its operational parameters. The steady-state space-independent analysis has set the foundation of understanding the process. To capture the steady-state spatial effects on the process, several approaches have been proposed. Energy and mass conservation laws may be used to compare the energy efficiency of three MD configurations, namely DCMD, AGMD, and VMD. Simulations of the laminar steady-state Navier-Stokes equation have been compared to experimental measurements of DCMD. A more recent work, also based on energy and mass balances, was compared to experimental measurements of a DCMD experimental setup.

However, as solar radiation or heat sources vary with time, the inlet temperatures of the MD process can be unsteady throughout the operation. Therefore, due to the spatial and temporal variations of the temperatures distribution inside the MD module and other effects, it is important to learn more about the dynamic response of the process. In addition, when coupling the MD process with renewable energy, the additional challenge of the unsteady nature of the solar energy has to be considered as it evolves in time. All of these reasons have motivated the dynamical modeling of the MD process. Dynamic partial-differential equations (PDEs) have been considered to model AGMD and DCMD, while a black-box model based on neural-networks and a dynamical model for DCMD based on a lumped-parameter method have also been proposed.

There are two main limitations for most of the reported models. They are either only applicable for steady-state dimensionless analysis or computationally costly, specifically for large module sizes where model reduction approaches may be needed. Unlike black-box models, physical models offer more insight into the process and its operational parameters (i.e., the feed and permeate inlet temperatures and flow rates as well as the geometry of the module and the membrane properties). Alternative approaches are based on lumped-parameter models that offer the benefits of physical models, while maintaining a reasonable computational cost. These models are therefore suitable for real-time control and optimization applications.

This disclosure is directed to enhancing the productivity of the MD process accompanied with deep understanding of its behavior, and a dynamic model for the DCMD process is presented with experimental validation in both transient and steady state phases. The model is based on the two-dimensional (2D) Advection-Diffusion Equation (ADE), and it is considered to be in two dimensions including the dimensions of the membrane as well as the feed and permeate solutions. The model appears can contribute toward improving the performance of the process and reducing energy consumption. The model predicts the process to be powered with intermittent energy supply considering its unsteady nature, e.g., renewable energy sources or feeding the module with variable water temperatures. Moreover, an accurate description of the time evolution of the different heat transfer mechanisms inside the MD module is provided, which can account for sudden changes that might happen to the process, besides the steady state phase. This leads eventually to making the process more suitable for automation and control systems, and therefore, resulting in better performance and efficiency. A model is described for the DCMD configuration due to its convenience for laboratory tests and simplicity, but the model can be extended to include other MD configurations. The different heat and mass transfer mechanisms which drive the MD process for improving the performance will now be examined.

Mass Transfer Mechanism in MD Process.

Maintaining the permeate solution to be cooler than the feed solution forces vapor molecules to transfer through the membrane dry pores, from feed side to permeate side where it condenses. Mass transfer starts when the feed solution reaches the evaporation limit, which lets the vapor form in a small region of the membrane boundary layer. The vapor then will be ready to pass through the dry pores of the membrane toward the permeate side. The transfer is driven by the vapor pressure gradient between the feed solution and the permeate solution. Mass transfer is directly proportional to the permeability of the membrane as well as the vapor pressure gradient through:J(T)=C(T)(P₁(T)−P₀(T)),  (1)where J(T) is the trans-membrane heat flux, P(T) is the pressure along the membrane boundary layers and C(T) is the membrane mass transfer coefficient of the system, and which are functions of the induced temperature T. The mass transfer coefficient can be modeled with a Knudsen model, a Poiseuille model or a combination of both in a molecular model.

The Knudsen model is the most suitable to describe the mass transfer in the DCMD process, since its criteria of having a mean free path of the transferred molecules to be larger than the membrane pore size is achieved, and the effect of the molecular diffusion for the air trapped within each membrane pores is extremely small relative to the water flux, and therefore it can be neglected. It can be expressed as:

$\begin{array}{cc}J\left(T\right)=1.064\frac{rϵ}{{\mathrm{\chi \delta }}_m}{\left(\frac{M}{{\mathrm{RT}}_{\mathrm{mean}}}\right)}^{0.5}\left(P_1\left(T\right)-P_0\left(T\right)\right),& \left(2\right)\end{array}$ where r is the pore size, ϵ is the porosity, χ is the tortuosity, δ_m is the membrane thickness, R is the gas universal constant, M is the water molecule weight, T_mean is the mean temperature along the membrane interfaces, and P is the pressure. It is worth mentioning that Knudsen model utilizes knowledge of some membrane parameters, such as pore size, porosity and thickness. These parameters can be estimated by applying the models to experimental gas fluxes through a given membrane.

Heat Transfer Mechanism in MD Process.

Heat transfers with the vapor and volatile molecules from the feed solution to the permeate solution. This transfer follows the second law of thermodynamics. The transferred heat takes a form of sensible and latent heat inside the membrane, and conduction and convection heat throughout the bulk flow in the feed and the permeate sides.

The majority of the heat transfer models are set according to factors, such as heat balance equations, and conservation of energy and mass laws in thermodynamics. As a result, the involved amount of transferred heat should be equal in both process sides as well as through the membrane dry pores. The involved amount of heat to be generated in the process can be shown as:

$\begin{array}{cc}{Q}_f={\alpha }_f\left(T_f-T_1\right),& \left(3\right)\\ Q_p={\alpha }_p\left(T_2-T_p\right),& \left(4\right)\\ Q_m=\frac{k_m}{{\delta }_m}A\left(T_1-T_2\right)+J\left(T\right)H\left(T\right),& \left(5\right)\end{array}$ where Q is the total heat transferred from the feed through the membrane to the permeate, k_m is the thermal conductivity of the membrane. A is the membrane area, α is the convective heat transfer coefficient, and H(T) is the latent heat of vaporization.

DCMD Models.

Many MD models have been developed to form mathematical formulas for the membrane boundary layers. However, the majority of these models are built upon empirical relations and laboratory experiments. These models are either valid for steady-state phase, or consider only zero-order performance with no dynamic behavior. This means that they are valid only after a sufficient period of time has passed, where the process reaches the steady-state phase and the performance becomes stationary. A group of mathematical relations for the temperature of the membrane boundary layer region have been proposed and these relationships derived from semi-empirical relations and laboratory experiments. Moreover, these relations are based on the fact that it is difficult to measure interfacial temperatures and concentrations directly. Another model has been proposed to estimate those parameters considering a full knowledge of the heat and mass transfer coefficients throughout the MD module and bulk flow temperatures and concentrations.

In another implementation, a group of differential equations with respect to space can be assumed. The model gives a formula for the membrane interface temperature by relating parameters of the process, however it doesn't account for the time evolution of the process. This model is similar to the model to estimate parameters considering the heat and mass transfer coefficients, but with additional complexity. More elaboration is spent on the same empirical relations that were used. The contribution is through development of an optimized iterative method that reduces the error between initial guess for the membrane interface temperature with the assumed values, so they succeeded to know the temperature in each experiment. An intuition to consider a dynamic model for time evolution of the heat in the MD process was addressed, where the interest was focused on introducing a numerical simulation for the Navier-Stokes equations. Interest in dynamic modeling for the MD process has continued, with many studies tackling the dynamic analysis of heat and mass transfer using Computational Fluid Dynamics (CFD) techniques.

The literature models offer a predictive estimation for the temperature of the boundary layer region of the membrane, but they are unable to provide information for the temperature of the bulk solutions inside. Moreover, the CFD modeling requires high computational effort and cost, and the steady-state models do not account for the time evolution of the process, and the needs of the intermittent sources. For example, if the use of solar energy source to drive the MD process is considered, the model should account for the variable output intensity of the source throughout the day.

Dynamic Modeling of DCMD.

A dynamic mathematical model for the DCMD process is disclosed, which considers the time evolution of the heat and mass transfer mechanisms throughout the feed and the permeate solutions. In addition, the dynamic model offers a complete overview of the process performance starting with the transient phase until the steady-state phase is reached. Heat inside feed and permeate solutions transfers in a form of a sensible and latent heat, and a convective heat. Therefore, mechanisms of heat conduction and heat convection are considered within the evolution of heat inside the DCMD module.

Referring now to FIG. 1, shown is a schematic diagram of an example of a counter current DCMD module 100a with the coordinates illustrated for each process component. The DCMD module 100a includes a feed side 103 and a permeate side 106 separated by a membrane 109. The feed side 103 includes a feed inlet 112 and feed outlet 115 of hot/cold water and the permeate side 106 includes a permeate inlet 118 and permeate outlet 121 of cold/hot water. Due to the very small thickness of the membranes 109 that are used, heat transfer mechanisms inside the feed side 103 and the permeate side 106 are affected by each other, and coupled together. This attachment is clear via the heat flux that flows from the feed solution to the permeate solution through the dry pores of the membrane 109. A 2D ADE can be used to model the dynamic heat transfer mechanisms. The interesting properties of the ADE model make it possible to describe the convection and conduction mechanisms of heat transfer inside the feed and permeate solutions. Moreover, the model accounts for the evolution with time and considers the transient behavior of the MD process.

The DCMD mathematical model can be given by:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial t}+v_2\frac{\partial T\left(x,z,t\right)}{\partial z}+v_x\frac{\partial T\left(x,z,t\right)}{\partial x}={\alpha }_x\frac{{\partial }^{2}T\left(x,z,t\right)}{\partial x^2}+{\alpha }_2\frac{{\partial }^{2}T\left(x,z,t\right)}{\partial z^2},& \left(6\right)\end{array}$ with0<t<τ_final,x∈=[0,x_fm]∪[x_pm,X],0<z<Z, where x∈[0,x_fm] is in the feed and x∈[x_pm,X] is in the permeate. ν_x and ν_z are the hot and the cold water stream flow rates along the x, z axes, respectively. They can be assumed to be constant along the operation of the process. α_x and α_z are constants that depend on physical properties of the used water such as: thermal conductivity κ_e, specific heat c_p, and the density ρ in:

$\begin{array}{cc}\alpha =\frac{k_e}{\rho c_p}.& \left(7\right)\end{array}$ The initial conditions can be set to be constant throughout the model domain as:T(x,z,0)=T_finitial for x∈[0,x_fm],  (8)T(x,z,0)=T_pinitial for x∈[x_pm,X].  (9)

The associated boundary conditions with the model (e.g., the feed solution) can be assumed to be a Dirichlet condition for the inlet of each subsystem container like in:T(x,0,t)=T_0f for x∈[0,x_fm],  (10)T(x,Z,t)=T_0p, for x∈[x_pm,X],  (11)and Neumann conditions for the rest of the boundaries. The process can be isolated from the left side boundary of the feed side as in:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=0}=0,& \left(12\right)\end{array}$ and from the right boundary side for the permeate side as in:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=X}=0,& \left(13\right)\end{array}$ and can have a trans-membrane heat transfer formula for the heat flux and the latent heat of evaporation with the other subsystem (e.g. the permeate subsystem) such as in:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=x_{\mathrm{fm}}}=\left[J\left(T\right)H\left(T\right)+\frac{k_m}{{\delta }_m}\left(T\left(x_{\mathrm{fm}},z,t\right)-T\left(x_{\mathrm{pm}},z,t\right)\right)\right]\text{/}{k}_f& \left(14\right)\end{array}$ for the feed side and as in:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=x_{\mathrm{pm}}}=\left[J\left(T\right)H\left(T\right)+\frac{k_m}{{\delta }_m}\left(T\left(x_{\mathrm{fm}},z,t\right)-T\left(x_{\mathrm{pm}},z,t\right)\right)\right]\text{/}{k}_p& \left(15\right)\end{array}$ for the permeate side, and finally the process has a heat flux release at the bottom boundary for the feed side like in:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial z}{|}_{z=Z}=\varphi ,& \left(16\right)\end{array}$ and at the top boundary for the permeate side as in:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial z}{|}_{z=0}=\varphi ,& \left(17\right)\end{array}$ where ϕ∈ is a positive constant to represent the amount of flux that is released outside the process containers. It is considered to be zero in the following simulations. H(T) is the latent heat of evaporation that is responsible of the released sensible heat. The empirical formula for the latent heat of water can be shown as:H(T)=(2500.8−2.362500.8−2.36T+0.0016T²−0.00006T³)J/g.  (18)

In the ADE model, it can be assumed that the heat convection transfer is represented by the transport terms

$\left(\frac{\partial T\left(x,z,t\right)}{\partial z},\frac{\partial T\left(x,z,t\right)}{\partial x}\right).$ and the heat conduction transfer is represented by the second derivative terms

$\left(\frac{{\partial }^{2}T\left(x,z,t\right)}{\partial z^2},\frac{{\partial }^{2}T\left(x,z,t\right)}{\partial x^2}\right).$

Semi-Discretized Model.

It is usually difficult to obtain analytic solutions for Partial Differential Equation (PDE) models, hence numerical methods are an appropriate alternative to get an acceptable approximation for these solutions. Numerical approaches commence by discretizing the equations and forming the discretization grid. This process is fulfilled through dividing the desired domain into specific numbers of points denoted by (n), each point in the grid represents the temperature at that location. In order to indicate the temperature location on the grid, each temperature is indexed with two subscripts T_i,j to indicate the rows (i) and columns (j). The distance between two consecutive points in the same row in the grid is equal to the horizontal axis step size Δx, where it is equal to Δz between two points below each other in the vertical axis. The discretization grid specifications can be expressed as:x_i=iΔx,  (19)z_j=jΔz, and  (20)T(x_i,z_j)=T_i,j.  (21)

The ADE model is semi-discretized (space only) using finite difference discretization schemes. The convective and diffusive terms are solved with first order upwind scheme, and central difference scheme, respectively. The semi-discretization transforms the PDE model from infinite dimensional problem into a finite dimensional group of Ordinary Differential Equations (ODEs). The relations of the grid points in the domain are significant to solve the model. FIG. 2 shows an example of a grid 203 of discretization for the feed side 103 and permeate side 106, where the points T_i,j that are located between x∈[0, x_fm] represent the feed solution temperatures and the points T_i,j that are located between x∈[x_pm,X] represent the permeate solution temperatures. In DCMD module 100a, the permeate solution starts right after the thickness of the membrane (δ_m) and the width of the feed solution (x_fm), as in:x_pm=x_fm+δ_m.  (22)

The discrete version of the 2D ADE is a nonlinear high dimension system of ODEs, given by:

$\begin{array}{cc}\underset{\_}{\stackrel{.}{T}}\left(x,z,t\right)=\frac{\partial \underset{\_}{T}\left(x,z,t\right)}{\partial t}=A\underset{\_}{T}\left(x,z,t\right)+B\underset{\_}{U}\left(t\right)+\underset{\_}{b}\left(\underset{\_}{T}\left(x,z,t\right)\right),& \left(23\right)\end{array}$ where vector T∈^2n2×1 is the state vector which represents the temperature distribution at each grid point inside the feed domain (T_f), and the permeate domain (T_p) for the time evolution, such as in

$\begin{array}{cc}\underset{\_}{T}=\left[\begin{array}{c}{\underset{\_}{T}}_f\\ {\underset{\_}{T}}_p\end{array}\right],& \left(24\right)\end{array}$ and thus T_f,T_p∈^2n2×1 with n the number of discretization points in x and z. Matrix A∈^2n2×2n2 represents the operator matrix that has the dynamics of the feed and the permeate subsystems with their couplings on the boundaries as shown in:

$\begin{array}{cc}A=\left[\begin{array}{cc}{A}_{f_{i,j}}& C_{\mathrm{fp}}\\ C_{\mathrm{pf}}& A_{p_{i,j}}\end{array}\right],& \left(25\right)\end{array}$ where i,j are the discretization grid indices, A_f∈^n2×n2, A_p∈^n2×n2, C_fp∈^n2×n2, and C_pf∈^n2×n2 are in:

$\begin{array}{cc}{A}_{f_{i,j}}=\left(-\frac{v_z}{\Delta z}-\frac{v_x}{\Delta x}-\frac{2{\alpha }_x}{\Delta x^2}-\frac{2{\alpha }_z}{\Delta z^2}\right)T_{i,j}+\left(\frac{v_z}{\Delta z}+\frac{{\alpha }_z}{\Delta z^2}\right)T_{i,j-1}+\left(\frac{v_x}{\Delta x}+\frac{{\alpha }_x}{\Delta x^2}\right)T_{i-1,j}+\left(\frac{{\alpha }_z}{\Delta z^2}\right)T_{i,j+1}+\left(\frac{{\alpha }_x}{\Delta x^2}\right)T_{i+1,j},& \left(26\right)\\ A_{p_{i,j}}=\left(\frac{v_z}{\Delta z}+\frac{v_x}{\Delta x}-\frac{2{\alpha }_x}{\Delta x^2}-\frac{2{\alpha }_z}{\Delta z^2}\right)T_{i,j}+\left(-\frac{v_z}{\Delta z}+\frac{{\alpha }_z}{\Delta z^2}\right)T_{i,j+1}+\left(\frac{{\alpha }_x}{\Delta x^2}\right)T_{i-1,j}+\left(\frac{{\alpha }_z}{\Delta z^2}\right)T_{i,j-1}+\left(-\frac{v_x}{\Delta x}+\frac{{\alpha }_x}{\Delta x^2}\right)T_{i+1,j},& \left(27\right)\\ C_1=\frac{{\alpha }_x}{\Delta x^2}\frac{\Delta {\mathrm{xk}}_m}{{\delta }_m{k}_l},l=f,p.& \left(28\right)\end{array}$

Matrix B∈^2n2×2 is the input matrix as shown in:

$\begin{array}{cc}B=\left[\begin{array}{cc}{\left(\frac{v_z}{\Delta z}+\frac{{\alpha }_z}{\Delta z^2}\right)}_1& 0\\ ⋮& ⋮\\ {\left(\frac{v_z}{\Delta z}+\frac{{\alpha }_z}{\Delta z^2}\right)}_n& 0\\ 0& 0\\ ⋮& ⋮\\ 0& 0\\ 0& {\left(-\frac{v_z}{\Delta z}+\frac{{\alpha }_z}{\Delta z^2}\right)}_1\\ 0& ⋮\\ 0& {\left(-\frac{v_z}{\Delta z}+\frac{{\alpha }_z}{\Delta z^2}\right)}_n\end{array}\right],& \left(29\right)\end{array}$ and U(t) is the input vector that contains feed and permeate temperature inlets as in:

$\begin{array}{cc}\underset{\_}{U}=\left[\begin{array}{c}{T}_{{\mathrm{fN}}_f}\\ T_{{\mathrm{fN}}_p}\end{array}\right].& \left(30\right)\end{array}$ Vector b(T) contains the boundary conditions associated with the membrane boundary layer cells for the feed and the permeate subsystems, such that the trans-membrane heat flux and the latent heat of evaporation are both functions of the membrane interfaces temperature difference, such as in:

$\begin{array}{cc}\underset{\_}{b}=\left[\begin{array}{c}{\zeta }_1\\ ⋮\\ {\zeta }_n\\ {\eta }_1\\ ⋮\\ {\eta }_n\end{array}\right].& \left(31\right)\end{array}$ The polynomial nonlinearity structures in the boundary vector are represented by the vectors:

$\begin{array}{cc}\underset{\_}{\zeta }=\left[\begin{array}{c}0\\ ⋮\\ \frac{{\alpha }_x}{\Delta x}\frac{J\left(T\right)H\left(T\right)}{k_f}\end{array}\right]& \left(32\right)\\ \underset{\_}{\eta }=\left[\begin{array}{c}-\frac{{\alpha }_x}{\Delta x}\frac{J\left(T\right)H\left(T\right)}{k_p}\\ 0\\ ⋮\\ 0\end{array}\right]& \left(33\right)\end{array}$ The measured output can be given by:y=CT  (34)with the output matrix C∈^1×2n2 being C=[0 . . . 0 1 1 0 . . . 0].

Referring next to FIG. 3A, shown is a schematic diagram of a fully automated DCMD experimental set up used to evaluate the dynamic model. A flat sheet membrane module with an active surface area of 0.005 m², with flow channel dimensions of 0.1 m×0.05 m×0.002 m, was designed and fabricated using poly methyl methacrylate (PMMA) material. A composite membrane with a poly tetra fluoroethylene (PTFE) active layer and a nonwoven poly propylene support layer was used for MD testing. Some of the membrane parameters are presented in the Table of FIG. 3B, while the Table of FIG. 3C presents the feed water characteristics.

As shown in FIG. 3A, Red Sea water was preheated to desired temperatures and circulated through the feed side 103 of the membrane 109, while de-ionized (DI) water was circulated through the premeate side 106 of the membrane 109 simultaneously in a countercurrent mode. Feed and permeate temperatures were controlled using thermo-regulators. Fresh seawater was used for each experiment. The permeate generated during the MD process resulted in an increase in the volume of coolant and overflow through the outlet of the permeate/coolant tank 303. The overflown permeate was collected in a separate container 306 placed on a weighting balance. The increase in weight of the container 306 was continuously monitored and recorded. Pressure, temperature and flow rates were continuously monitored at the inlet and outlet of the feed and permeate sides 103/106 using respective sensors connected to a data acquisition unit equipped with a Lab View Software. Conductivities of both the permeate and feed solutions were continuously monitored and measured using conductivity meters (Oakton Eutech Instruments, Malaysia).

The influence of feed water temperature ranging from 40° C. to 70° C. on water vapor flux was studied by maintaining a constant temperature on the permeate side at 20° C. All the data were recorded after reaching the steady state temperature. Similarly, the influence of feed water temperature ranging from 30° C. to 75° C. on trans-membrane flux was also studied by ramping the feed solution temperature at a rate of 0.05° C. per minute. During ramping experiments, the coolant temperature was also maintained at 20° C. Feed and coolant flow rates were kept constant at 90 Liters/hour and 60 Liters/hour, respectively, for all the experiments.

The salt rejection was calculated by: SR=(1−C_p/C_f)*100, where C_p and C_f are the salt concentrations of permeate and feed solutions, respectively. A salt rejection of 99:99% has been observed in all experiments. The water vapor flux can be determined by J_v=mw/At, where mw is the weight of collected permeate at a particular time interval t, and A is the effective membrane area.

Model Validation.

The ADE model has been simulated with a developed MATLAB® software and then validated with experimental data through two sets of experiments. For validation purposes, the same DCMD parameter values and membrane specifications as well as environment conditions were applied.

Steady-State Validation.

Referring to FIG. 4, shown is a plot of the validation of the dynamic model at the steady-state phase between the experimental data and the simulations for different feed inlet temperatures. The absolute error is also presented in FIG. 4. The model acts best when compared to experiments at high feed inlet temperatures varying from 55° C. to 70° C., while the absolute error is a bit larger for lower temperatures below 40° C. This may be attributed to the water vapor flux being higher at high feed temperatures, which reduces the errors in measurements compared to the very low flux obtained at low feed inlet temperatures.

Transient Validation.

In the other set of experiments, the feed inlet temperature was increased in a ramp fashion from 30° C. to 75° C. The increment was 0.1° C. every 2 minutes, with the outlet feed temperature and permeate flux being recorded. The comparison between the experimental results and the modeled outlet feed temperature and permeate flux is shown in FIG. 5 for the transient phase. In this observation, the transient response of the outlet feed temperature and the permeate flux showed high consistency with the experimental data. This contributes towards the validity of the developed dynamic model to be extended to the transient phase.

The MD process is highly sensitive to considerable variations in the membrane parameter values. For example, FIG. 6 illustrates the effect of membrane pore size variation on the membrane flux, where a higher pore size leads to higher permeate flux. At the same time, the Liquid Entry Pressure (LEP) should not be exceeded to ensure no pore wetting occurs inside the membrane pores. The LEP is equal to 15 PSI for the used membrane, while the actual pressure difference at the liquid-vapor interface ΔP is lower or equal to (7.64-8.33 PSI) for a pore size equal to 0.5 μm and along temperature variation between (20° C.-60° C.), respectively. The variation of the membrane parameter values is very likely to happen when manufacturing the membranes, and so it can introduce some divergence in some cases.

Model Analysis.

The ADE model was applied inside the feed solution as well as the permeate solution. For simulations and up-scaling purposes a software (in a MATLAB® environment) was developed to solve the dynamic model and obtain behavior plots for analysis and improvement aspects. In the simulations, real numeric values for the membrane parameters are used as listed in the Tables in FIGS. 3B and 7. The simulations estimated the response in transient and steady-state phases. The feed and coolant temperatures are set at 60° C. and 20° C., respectively.

Referring to FIG. 8, shown is an example of the evolution of the feed and permeate solution temperatures at specific instant of time during the transient phase. FIG. 8 shows the feed bulk temperature to be higher than the temperature at the membrane boundary layer from the feed side, and in a same manner the temperature of the membrane boundary layer from the permeate has higher temperature than permeate bulk. It also shows that the difference of temperature between the two sides of the membrane (driving force) is higher at the top side (inlet) of the feed side 103 and decreases gradually towards the bottom side (outlet) of the feed side 103, which is in agreement with the experimental data.

Referring next to FIG. 9, shown is an example of the temperature of the membrane boundary layers at a specific point along the membrane length. FIG. 9 illustrates how close is the amount of temperature drop in the feed side and gained in the permeate side. The feed interlace temperature ranges between 55° C. and 60° C., and the permeate interface temperature ranges between 20° C. and 25° C. along the membrane length.

The flow rate of the feed and the permeate streams have a significant role in managing the heat transfer inside the DCMD module. The flow rate controls the amount of the convective heat that transfers throughout the bulk solution, in which increasing its value makes the heat transfers via convection in the bulk more than the conduction through membrane pores. It is useful to have high convection effect, however continuous increasing in the flow rate is unprofitable, because the heat flux would saturate eventually. A smart compromise should take place between convection and conduction mechanisms. The relation between the feed flow rate and the temperature for different flow rates is shown in FIG. 10. Temperature of water molecules near the inlet increases rapidly while increasing flow rate, and therefore, temperature polarization and residence time decrease, while energy consumed by the flow pumps increases.

The temperature polarization was simulated to check the time evolution of its coefficient (TPC). The formula of the TPC is:

$\begin{array}{cc}\mathrm{TPC}=\frac{T_{\mathrm{mf}}-T_{mp}}{T_{\mathrm{bf}}-T_{\mathrm{bp}}}.& \left(35\right)\end{array}$ The value of TPC varies from zero for systems that suffer from large boundary layer resistances to unity for well designed systems. FIG. 11 shows examples of the evolution of TPC with time. During the transient phase of the process the amount of heat transfer through the membrane is low due to the starting of heat accumulation, and then it magnifies and increases after reaching steady-state phase.

The dynamic model can be efficiently utilized for DCMD process applications with an intermittent energy supply, such as a solar energy source. The model predictions are strongly correlated with the experimental measurements in both steady state and transient phases, with model predictions being within ±5% of the experimentally obtained values.

The DCMD dynamic model is nonlinear, and it is originated from the functions J(T) and H(T). The present nonlinearity is Lipshitz continuous, which can be shown by examining each function separately. H(T) is Lipshitz continuous following the fact that it is smooth and infinitely differentiable polynomial function. On the other hand, define J(T) to be

$\begin{array}{cc}𝒥\left(T\right)={c_1\left(\frac{1}{T}\right)}^{\frac{1}{2}}{10}^{\frac{-c_2}{T}},& \left(36\right)\end{array}$ where c₁,c₂ are constants, and T∈[a,b], b>a>0. Therefore, J(T) is continuous and differentiable on this interval [a,b], and bounded like in∥J(T)∥≤.  (37)As a result, there exists a Lipshitz constant λ that satisfies:∥f(T₁)−f(T₂)∥≤λ∥T₁−T₂∥,  (38)and the nonlinearity in the DCMD dynamic model is Lipshitz continuous.

Nonlinear Observer.

Nonlinear processes are challenging in stability and control, due to the high sensitivity to initial and forcing conditions. The operation of heat transfer, inside the MD module, is nonlinear and needs sufficient temperature difference along the membrane boundaries to operate efficiently. The temperature distribution along the membrane boundaries can be used in the control of the temperature difference to manage the water production inside the DCMD process. While measuring the temperature distribution is not viable in practice, the process model is observable and the present nonlinearity is Lipshitz continuous, therefore, an observer can be used to estimate the temperature states and feed them back to the controller.

Consider the nonlinear Lyapunov-based observer in the following form:=+U+(⁻¹)+L(y−).  (39)where ∈^2n2×11 is the estimated state vector for the temperature difference vector ({tilde over (T)}=T) along the membrane boundaries, and L is the observer gain matrix.In addition, ∈^n×2n2 is a non-square transformation matrix that form the difference of each correspondent boundaries, such as:

$\begin{array}{cc}\mathscr{H}=\left[\begin{array}{ccccccccc}0& \dots & 1& 0& -1& \dots & 0& \dots & 0\\ 0& \dots & 0& 1& \dots & -1& 0& \dots & 0\\ 0& \dots & 0& 0& \ddots & 0& \ddots & 0& 0\\ 0& \dots & 0& \dots & 0& 1& \dots & -1& 0\end{array}\right],& \left(40\right)\end{array}$ and ⁻¹=^T, =A⁻¹ ∈^n×n represents the dynamics of the membrane boundaries with representing the dynamics matrix which includes all dynamics of the feed and the permeate with their correspondent couplings on the boundaries as in Eq. (25), =B∈^n×2 is the new input matrix, (⁻¹{tilde over (T)})=f(T)∈^n×1 is the nonlinear effect on the membrane boundaries, =C⁻¹ ∈^1×n is the output matrix, and T=T. The observer error dynamic equation can be formulated as:{dot over (T)}=(A−LC)T+(⁻¹{tilde over (T)})+(⁻¹{circumflex over ({tilde over (T)})})  (41)where T={tilde over (T)}−{circumflex over ({tilde over (T)})} is the state estimation error. The following theorem provides sufficient conditions for the exponential stability of the state observer error of Eq. (41).

Theorem.

There exists an exponential stable observer of the form of Eq. (39), if and only if, there exists a symmetric positive definite matrix, P, such that the following standard algebraic Ricatti equation is satisfied:P+P^T+PR₁P+Q₁=0,  (42)where

$\begin{array}{cc}{R}_1=\left({\gamma }^{2}I-\frac{C^{T}C}{ϵ}\right),Q_1=I\left(ϵ+1\right),ϵ>0,ϵ\in {ℝ}^{+},& \left(43\right)\end{array}$ and the observer gain matrix

$L=\frac{{\mathrm{PC}}^T}{2ϵ}.$

The performance of the nonlinear observer was tested with the convergence time and error. FIG. 12 shows the fast convergence of the estimated states using the nonlinear Lyapunov-based observer for different locations along the membrane length. These locations are marked as a percentage of the total length of the membrane 109 (FIGS. 1 and 2), starting from the top of the membrane to the bottom. As such, FIG. 12 shows the observer convergence on the top (6%), the middle (34%) and the bottom (67%) of the whole membrane length. The fast convergence of the observer was guaranteed by taking smaller sampling times. The error analysis shows the exponential convergence of the error to zero as depicted in FIG. 13.

Soft sensing of the temperature inside a MD module 100a (FIG. 1) along with the estimation of the membrane characteristics is possible using the dynamic model described above. The membrane temperatures are not accessible and are not provided by standard temperature's sensors. The distributed model, which includes a two-dimensional advection diffusion equation, can be used to describe the heat transfer in the MD module 100a. A soil sensing system (or algorithm) based on the observer concept can be used to estimate the temperature and the membrane characteristics of the dynamic system, using measured outlet temperature at both feed and permeate sides.

Estimating the temperature, at the membrane level, for MD plants can provide more information on the process, which can be used for monitoring the system in case something wrong happens. The estimated temperature can be provided to a controller, which controls the performance of the system, especially when connected to a variable source such as a solar energy source. The membrane distillation is driven by the difference of temperature between the two sides of the membrane. With the solar energy source being non-uniform and unstable, soft sensing of the temperature enables the whole process to be controlled in order to insure a desired fresh water production. Along with the temperature measurements, the soft sensing system can provide estimations of other membrane characteristics that cannot be directly measured, such as the mass transfer coefficient.

The present disclosure also builds upon a model based on an electrical analogy to thermal systems. The performance of this model is studied by analyzing the DCMD and to study its predictive capacity. The disclosed model accounts for the spatial and temporal simultaneous heat and mass transfer phenomena in the DCMD. This model enables dynamic predictions of distilled water flux and the temperature distribution along the flow direction by directly solving a system of differential algebraic equations (DAEs). It also offers fixable platform to carryout scale up and design studies, and enable real-time estimation of some intrinsic variables that are inaccessible in physical setups, like the membrane mass transfer coefficient and the temperature at the membrane interfaces. Accordingly, this leads to optimal operating conditions for the MD process. Moreover, the same foundation can be easily extended to other MD configurations.

Dynamic Model of DCMD Formulation.

Referring to FIG. 14, shown is a schematic diagram of a flat-sheet DCMD module 100b. In this configuration, hot water is passed along a hydrophobic membrane 109 on one side (the feed side 103) and cold water flows in the counter direction along the other side (the permeate side 106). Water vapor is driven from the feed side 103 across the membrane 109 and into the permeate side 106 by the induced partial vapor pressure difference. Both heat and mass transfer processes occur simultaneously as water evaporates at the feed-membrane interface and condenses at the permeate-membrane interface. As a result, the temperature at the membrane boundary layers differs from the bulk temperature of the feed and permeate streams, this is known as the temperature polarization effect. In order to quantify this effect, the temperature polarization coefficient (TPC) is given as:

$\begin{array}{cc}\mathrm{TPC}=\frac{T_{\mathrm{mf}}-T_{mp}}{T_{\mathrm{bf}}-T_{\mathrm{bp}}}.& \left(44\right)\end{array}$ This effect reduces the mass transfer driving force, and as a result lowers the production rate of the MD desalination system. Moreover, because of the variation of the feed/permeate temperature throughout the module 100b, the TPC varies along the membrane module. Therefore, by dividing the DCMD module 100b into small control volumes of length dz, the temperature changes along the membrane can be accounted for in both the feed and permeate sides. The following subsections detail the mass and heat transfer processes.

Mass Transfer in DCMD.

The transport phenomena can be described by the classic gas permeation and heat transfer theories. The mass flux (J) in the DCMD 100b is related to the saturated vapor pressure difference ΔP across the membrane 109 through the membrane mass transfer coefficient B_m, as follows:J=B_mΔP=B_m(P_mf−P_mp).  (45)The mechanism dominating the mass transfer through the porous membranes 109 depends on the pore radius (r) and the mean free path of the vapor molecules (λ). For membranes with pore radius in the range of 0.5λ<r<50λ, the membrane mass transfer coefficient is expressed as a parallel combination of Knudsen diffusion (B_Kn) and molecular diffusion (B_D) coefficients given by:

$\begin{array}{cc}{B}_m=\frac{1}{1/B_{\mathrm{Kn}}+1/B_D},& \left(46\right)\end{array}$ where:

$B_{\mathrm{Kn}}=\frac{4}{3}\frac{\varepsilon r}{{\chi }^{\delta }}\sqrt{\frac{2m_w}{\pi \stackrel{\_}{R}T}},\mathrm{and}{B}_D=\frac{\varepsilon }{{\chi }^{\delta }}\frac{\mathrm{PD}}{P_{\alpha }}\frac{m_w}{\stackrel{\_}{R}T},$ where χ is the tortuosity of the membrane pores, δ is the membrane thickness, ε is the membrane porosity, m_w is the molecular mass of water, R is the gas constant, P is the vapor pressure, D is the diffusivity of the water vapor and air mixture, P_α is the air pressure, and T is the temperature.

The saturated vapor pressure of pure water (P_w^sat[T]) as a function of temperature is given by the Antoine equation:

$\begin{array}{cc}{P}_w^{\mathrm{sat}}\left[T\right]=\exp \left(23.1964-\frac{3816.44}{T+227.02}\right).& \left(47\right)\end{array}$ Dissolved salt in the feed stream reduces the saturated vapor pressure. Therefore, to compensate for this, the following relation can be used:P_mf=(1−x_NaCl)(1−0.5x_NaCl−10x_NaCl²)P_w^sat[T_mf],  (48)where x_NaCl is the mole fraction of NaCl in the feed stream. However, the permeate is pure and the saturated vapor at the membrane-permeate interface is P_mp=P_w^sat[T_mp].

Heat Transfer in DCMD.

To consider the spacial variations on the temperature along the feed and permeate flow directions, the DCMD module 100b can be divided into control-volume cells. Then, based on the lumped-capacitance method, a dynamical model for heat transfer can be developed using the energy conservation law. Conservation laws can be applied on membrane distillation systems in order to develop a mathematical model which describes the dynamical relations between the process variables. This approach is based on lumped-capacitance method, which reduces the DCMD system to discrete control-volume cells. FIG. 15 depicts the n^th DCMD cell, where the bulk temperatures (T_bfn, T_bpn) are uniform throughout the cell except at the membrane interfaces due to the temperature polarization effect. Therefore, heat transfer takes place in three stages, by conduction and due to mass transfer. In the first stage, heat is transferred from the hot bulk feed stream to the boundary layer at the feed-membrane interface, the heat transfer rate is expressed as Q_mf by:Q_mf=A_m(h_f(T_bfn−T_mfn)+J_nc_pT_bfn).  (48)The rate of change of the bulk feed stream energy in the n^th cell can now be expressed as:

$\begin{array}{cc}{C}_{\mathrm{bfg}}\frac{{\mathrm{dT}}_{{\mathrm{bf}}_n}}{\mathrm{dt}}=Q_{f_n}-Q_{f_{n+1}}-A_m\left(h_f\left(T_{{\mathrm{bf}}_n}-T_{{\mathrm{mf}}_n}\right)+J_n{c}_p{T}_{{\mathrm{bf}}_n}\right),& \left(49\right)\end{array}$ where Q_fn and Q_fn+1 are the heat transfer rates into and out of the n^th feed cell, respectively.

At the second stage, heat is transferred through the membrane via three mechanisms: the first mechanism (Q_m1) is the latent heat of vaporization (H_v) transported by the mass flux (J_n) through the n^th cell, expressed as:Q_m1=A_mJ_nH_v[T_mf]=B_m(P_mfn−P_mpn)H_v[T_mf],  (50)where the latent heat of vaporization H_v in (KJ/Kg) is expressed as a function of temperature:H_v[T]=−2.42T+2503.  (51)The second and third mechanisms are heat conduction through the membrane material and air trapped in the membrane pores which are combined in (Q_m2) as:Q_m2=A_mh_m(T_mf−T^mp),  (52)where the membrane heat transfer coefficient h_m is given as:

$h_m=\frac{k_g\varepsilon +k_m\left(1-\varepsilon \right)}{\delta }.$

Combining these mechanisms to write the energy balance at the membrane interfaces gives the following equation:Q_mf=Q_mp  (53)where the heat transfer rate at the permeate-membrane interface Q_mp is expressed as:Q_mp=A_m(h_p(T_mpn−T_bpn)+J_nc_pT_mpn)=Q_m1+Q_m2.  (54)

Finally, the third stage of heat transfer is where the water-vapor condenses at the permeate-membrane interface and heat is transferred to the bulk permeate mass. The rate of change of energy for the bulk permeate stream is given by:

$\begin{array}{cc}{C}_{\mathrm{bp}}\frac{{\mathrm{dT}}_{{\mathrm{bp}}_n}}{\mathrm{dt}}=Q_{p_n}-Q_{p_{n+1}}-A_m\left(h_f\left(T_{m{p}_n}-T_{{\mathrm{bp}}_n}\right)+J_n{c}_p{T}_{m{p}_n}\right),& \left(55\right)\end{array}$ where Q_pn and Q_pn+1 are also the heat transfer rate into and out of the n^th permeate cell, respectively.

The two heat transfer coefficients at the membrane interfaces (h₁, h_p) can be calculated from empirical correlations. These correlations depend on the flow characteristic (laminar or turbulent) and vary accordingly. The following relation can be used for both heat transfer coefficients:

$\begin{array}{cc}h=0.13{\mathrm{Re}}^{0.64}{\mathrm{Pr}}^{1/3}\frac{k_w}{D_h},& \left(56\right)\end{array}$ where Re and Pr are the Reynolds and Prandtl numbers, respectively.

The analysis done so far has not quantified the coupling terms between neighboring cells (Q_fn, Q_fn+1, Q_pn, and Q_pn+1) Equations (49)-(55) can be represented by an electrical analogue circuit. This representation has another advantage, besides facilitating simulation, it allows to couple neighboring cells by modeling the thermal inertia of the system, as discussed in the next subsection.

Dynamic model formulation. The dynamical model is based on electrical analogy to thermal systems, which can be derived from the basic laws of each system. Appendix A details the derivation process and a summary of the electrical analogs of the thermal system is shown in the table of FIG. 16. Indeed, electrical-analogy based methods can be used to describe the dynamical behavior of many industrial and biological systems such as heat exchangers, and the human cardiovascular system. Moreover, the transient diffusion phenomena and the heat transfer due to non-steady fluid flow can be described by an electrical analogue.

Based on the equations derived for the n^th DCMD cell, an electrical analogue can be constructed to simulate heat and mass transfer processes. A schematic representation of the electrical analogue of the n^th cell of the DCMD module 100b is shown in FIG. 17. The thermal capacity of the feed and permeate bulk sides is represented by C_bf and C_bp respectively. In each of the three stages of heat transfer discussed above, the heat transfer rate by conduction is proportional to the temperature difference across the thermal resistances R_f, R_m, and R_p, whereas the heat transfer rate due to mass transfer is modeled by the current sources Q_mfⁿ, Q_m1ⁿ, and Q_mpⁿ. This completes the analogy of heat transfer within the same cell, and in order to couple neighboring cells, the series impedances (opposition to the heat transfer rate) Z_fⁿ and Z_pⁿ are introduced. Apart from the series impedances, the table of FIG. 18 details the expression of each element in the electrical analogue circuit (or electrical thermal network) of FIG. 17.

Another aspect of the DCMD electrical analogy is consideration of the heat transfer by the feed and permeate inlet mass flow rates. Therefore, the electrical analog network of FIG. 17 should be fed and terminated properly to account for the heat transfer rates into and out of the MD module 100b. This, as well as the series impedances, are discussed in Appendix B, and the final results are presented as follows. The temperature drop from one cell to the next in the feed side is proportional to the series impedance Z_fⁿ, which can be given as:Z_fⁿ=R_fzⁿ+jωL_fⁿ,  (57)where:

$\begin{array}{cc}{R}_{\mathrm{fz}}^n=\frac{1}{M_{f_n}^2{c}_p^2\left(R_f+R_m+R_p\right)},\mathrm{and}& \left(58\right)\\ L_f^n=\frac{R_{f_z}^2{C}_{\mathrm{bf}}}{4}.& \left(59\right)\end{array}$ The symbol jω is used to indicate complex impedance. In this case the thermal inertia of the n^th cell feed side is modeled by the inductor L_fⁿ, which accounts for the dynamic response that resist sudden changes to temperature by changing the inlet mass flow rate. On the other hand, the steady-state response is given by the thermal resistance R_fzⁿ. Similarly, the permeate side series impedance (Z_pⁿ) is given as:Z_pⁿ=R_pzⁿ+jωL_pⁿ,  (60)where:

$\begin{array}{cc}{R}_{\mathrm{pz}}^n=\frac{1}{M_{p_n}^2{c}_p^2\left(R_f+0.5R_m+R_p\right)},\mathrm{and}& \left(61\right)\\ L_p^n=\frac{R_{\mathrm{pz}}^2{C}_{\mathrm{bp}}}{4}.& \left(62\right)\end{array}$ Both of Z_fⁿ and Z_pⁿ parametrization was done based on experimental identification using real data.

The other aspect is consideration of the heat transfer rates into and out of the MD module 100b. The feed and permeate inlet temperatures can be simulated by voltage sources that can be set as needed. This allows for simulation of any desired inlet temperature profile. On the other hand, by conservation of energy, the feed outlet temperature cannot go below the permeate inlet temperature and the permeate outlet temperature cannot exceed the feed inlet temperature. Based on this, the feed and permeate outlet sides can be terminated across termination resistances (R_fterm and R_pterm receptively) as:0=T_fout−T_pin−R_ftermQ_fn+1, and  (63)0=T_pout−R_fin+R_ptermQ_p1.  (64)

Referring to FIG. 19, shown is a schematic representation of an example of a complete electrical analogue network including properly identified and parameterized elements. Based on this analogue circuit, electrical laws can be applied to derive the model of the DCMD 100b.

The coupling between neighboring cells can now be quantified by the current (in thermal analogy, current is the heat transfer rate) through the inductors. At the feed side, the rate of change of the heat transfer rate from the n−1 cell to the n^th cell is proportional to the temperature difference between them. Taking into consideration the series impedance Z_fⁿ, this can be expressed as:

$\begin{array}{cc}\frac{d{Q}_{f_n}}{\mathrm{dt}}=\frac{1}{L_f^n}{T}_{{\mathrm{bf}}_{n-1}}-\frac{R_{\mathrm{fz}}^n}{L_f^n}{Q}_{f_n}-\frac{1}{L_f^n}{T}_{{\mathrm{bf}}_n}.& \left(65\right)\end{array}$ Using Kirchoff's current law at the n^th feed node, it follows that the rate of change for the bulk feed temperature (T_bfn) is:

$\begin{array}{cc}\frac{d{T}_{{\mathrm{bf}}_n}}{\mathrm{dt}}=\frac{1}{C_{\mathrm{bf}}}{Q}_{f_n}-\frac{1}{C_{\mathrm{bf}}}\left(\frac{1}{R_f}+J_n{A}_m{c}_p\right)T_{{\mathrm{bf}}_n}-\frac{1}{C_{\mathrm{bf}}}{Q}_{f_{n+1}}+\frac{1}{C_{\mathrm{bf}}{R}_f}{T}_{{\mathrm{mf}}_n}.& \left(66\right)\end{array}$ Notice that Eq. (66) is equivalent to Eq. (49), but now Eq. (65) describes the dynamics of the heat transfer rates into and out of the n^th feed cell (Q_fn and Q_fn+1, respectively).

Similarly for the permeate side, the rate of change of the heat transfer rate (Q_pn) is:

$\begin{array}{cc}\frac{d{Q}_{p_n}}{\mathrm{dt}}=\frac{1}{L_p^n}{T}_{{\mathrm{bp}}_{n-1}}-\frac{R_{\mathrm{pz}}^n}{L_p^n}{Q}_{p_n}-\frac{1}{L_p^n}{T}_{{\mathrm{bp}}_n},& \left(67\right)\end{array}$ and the dynamics of the bulk permeate temperature (T_bpn) is

$\begin{array}{cc}\frac{d{T}_{{\mathrm{bp}}_n}}{\mathrm{dt}}=\frac{1}{C_{\mathrm{bp}}}{Q}_{p_n}-\frac{1}{C_{\mathrm{bp}}{R}_p}{T}_{{\mathrm{bp}}_n}-\frac{1}{C_{\mathrm{bp}}}{Q}_{p_{n+1}}+\frac{1}{C_{\mathrm{bp}}}\left(\frac{1}{R_p}+J_n{A}_m{c}_p\right)T_{{\mathrm{mp}}_n}.& \left(68\right)\end{array}$ The coupling between the feed and the permeate dynamics in the n^th cell is established through the algebraic constraints of Equations (53) and (54), which can be written in residue form as:

$\begin{array}{cc}0=\left(\frac{1}{R_f}+J_n{A}_m{c}_p\right)T_{{\mathrm{bf}}_n}-\frac{1}{R_f}{T}_{{\mathrm{mf}}_n}-\left(\frac{1}{R_p}+J_n{A}_m{c}_p\right)T_{{\mathrm{mp}}_n}+\frac{1}{R_p}{T}_{{\mathrm{bp}}_n},\mathrm{and}& \left(69\right)\\ 0=\left(\frac{1}{R_m}+\frac{1}{R_p}+J_n{A}_m{c}_p\right)T_{{\mathrm{mp}}_n}-\frac{1}{R_p}{T}_{{\mathrm{bp}}_n}-J_n{A}_m{H}_v\left[T_{{\mathrm{mf}}_n}\right]-\frac{1}{R_m}{T}_{{\mathrm{mf}}_n}.& \left(70\right)\end{array}$

The outlet temperatures at the terminal cells of the feed and permeate analogue are also given by the algebraic equations (63) and (64). The heat and mass transfer equations (63)-(70) represent a nonlinear differential-algebraic system.

The electrical thermal network (ETN) model can now be completed with all elements of the network analyzed and parameterized. In this model, the states are the temperatures in each cell and the heat transfer rates into and out of the cell, the manipulated variables are the inlet feed and permeate water temperatures and flow rates, the controlled variables are the water mass fluxes in each cell which when averaged together represent the overall water mass flux of the DCMD module 100b. In the next section, the equations for the DCMD electrical analogy will be driven based on the analysis that has been carried out in order to describe the mass and heat transfer processes.

Experimental Setup and Materials.

The DCMD lab scale set-up used in this experimental work was locally designed and fabricated. A schematic diagram of this setup is presented in FIG. 20. All experimental runs were conducted in a 0.1 m×0.05 m×0.002 m flat sheet membrane module 109 (FIG. 14) made of polymethylmethacrylate (PMMA) material. A composite membrane with a polytetrafluoroethylene (PTFE) active layer and a nonwoven polypropylene support layer was used in all the testing. The table of FIG. 21 presents parameters of the DCMD module 100b.

For the steady state experiments, Red Sea water was preheated to desired temperatures and circulated through the feed side 103 (FIG. 14) of the membrane, while deionized (DI) water was circulated through the permeate side 106 (FIG. 14) of the membrane simultaneously in a counter-current mode. Feed and permeate temperatures were controlled using thermo-regulators. Fresh seawater was used for each experiment. The vapor flux produced in the MD process increases the volume of water in the coolant tank 303 and, as a result, overflow occurs through the outlet of the coolant tank 303 and was collected in a separate container 306 placed on a weighing balance. The increase in weight of the container 306 was continuously monitored and recorded. Pressure (P), temperature (T) and flow rates (F) were continuously monitored at the inlets and outlets of the module 100b for both feed and permeate sides 103 and 106 using respective sensors.

For the transient response experiment, the influence of feed water temperature ranging from 30° C. to 68° C. on trans-membrane flux was also studied by ramping the feed solution temperature at a rate of 0.05° C. per minute. During ramping experiment, the coolant temperature was also maintained at 20° C. For all experiments, feed and coolant flow rates were kept constant at 90 l/h and 60 l/h, respectively, and all measurements were fed through a data acquisition unit into NI LabView software. The conductivity of both the permeate and feed solutions was continuously monitored and measured using conductivity meters (Oakton Eutech Instruments, Malaysia).

Model Validation Results.

This model was implemented using MATLAB® environment, which allows the user to easily adjust the model to simulate various experimental setups and conditions and accurately predict the steady-state and transient behaviors of the DCMD process. In addition, the same model enables scaled up simulations which facilitate module design and process optimization. By increasing the number of total cells (N), the desired level of performance can be achieved. For all simulations in this work, a total of 10 cells were used to simulate the DCMD flat-sheet setup detailed in the previous section and the model was solved using MATLAB® ode15s solver which gave accurate and fast results.

To validate the ETN model, experimental data was collected for both steady state and transient model analysis. The effect of the feed inlet temperature and mass flow rate were investigated and the results are presented and discussed in the following subsections.

DCMD Performance at Different Feed Inlet Temperatures.

One of the main factors contributing to the flux in DCMD is the feed inlet temperature. Therefore, it was demonstrated that the DCMD model adequately predicts the flux at different feed inlet temperatures. For this experiment, the feed and permeate flow rates were kept constant at 90 l/hr and 60 l/hr, respectively. The permeate inlet temperature was also maintained at 20° C. while the feed inlet temperature was increased from 40° C. to 70° C. in steps of 10° C. The results are plotted in FIG. 22, which shows very good agreement with the experimental measurements throughout the wide range of feed inlet temperatures. The maximum absolute error is 1.3 Kg/m² hr at 70° C., while an error of less than 1.0 Kg/m² hr was observed for the other temperatures.

Dynamic Response of the System.

To study the transient response of the DCMD process, this experiment was designed to ramp up the feed inlet temperature from 30° C. to 68° C. with increment of 0.1° C. per two minutes approximately, while maintaining the permeate inlet temperature constant at 20° C. Counter-current flow mode was used in this experiment, the feed flow rate was 90 l/h, and the permeate flow rate was 60 l/h. Measurements of the feed outlet temperature and water vapor flux were recorded once every minute during the course of approximately 14 hours experiment.

Referring to FIG. 23, shown is a plot of feed outlet temperature for a ramp feed inlet temperature from 30° C. to 68° C. The simulated (solid line) and the experimental measurements (dashed-line) of the feed outlet temperature are depicted. It can be seen that the feed outlet temperature linearly follows the inlet profile and it is clear that the simulation results closely match the experimental data. The simulation results of distilled water flux are illustrated in FIG. 24, which shows a plot of the vapor flux for a ramp feed inlet temperature from 30° C. to 68° C. The dots are the experimental measurements and the solid line is simulation results. Since the vapor flux is a nonlinear function of the trans-membrane temperature difference, an exponential increase in the flux with respect to the feed inlet temperature was noticed. Also, the figure shows good agreement between the measured and simulated flux.

DCMD Process Prediction.

Based on the satisfactory results obtained for the steady-state and dynamic response validations, the same model can be used to solve for intrinsic variables that are not accessible in the physical setups. In the next subsections, the effect of the feed and permeate inlet velocities on the TPC and the response of temperatures and vapor flux at a given cell to a step change in the inlet velocity are presented. In order to demonstrate the flexibility of the model, the DCMD module dimensions for all the following studies were: length 0.4 m, width 0.2 m, channel thickness 0.003 m and keeping the same membrane parameters shown in the table of FIG. 21.

Effect of Linear Velocity on the TPC.

To investigate the effect of feed/permeate inlet velocities on the TPC, the following simulation was designed. The feed and permeate inlet temperatures are kept constant at 60° C. and 20° C., respectively, while the inlet velocities were increased from 0.1 m/s to 0.65 m/s. FIG. 25 shows the TPC as a function of the feed inlet velocity. It is clear that increasing the inlet velocities has a desirable effect on the TPC. However, this effect begins to saturate starting from 0.5 m/s onwards. As the velocity inside the feed/permeate channels increases, the flow becomes more turbulent and, as a result, the thermal boundary layer thickness is reduced and therefore the temperatures at the membrane interfaces are brought closer towards that of the respective bulk temperature.

Dynamic Response to a Step Change of Inlet Velocities.

The other advantage of the ETN model is the ability to obtain a time response under dynamic changes in the process. The effect of increasing the feed inlet temperature was previously examined. In this simulation, the feed and permeate inlet temperatures were again kept constant at 60° C. and 20° C., respectively, while the feed and permeate inlet velocities were changed in steps from 0.2 m/s to 0.4 m/s at t=10 s.

Referring to FIGS. 26A and 26B illustrate the temperature response of the third cell under a step change in the inlet velocities. FIG. 26A shows the bulk feed temperature and FIG. 26B shows the membrane-feed temperature. First, note that the DCMD process response is fast, given sufficient actuation power. Second, note that the temperature increase at the membrane-feed interface is higher (about 1.5° C.) than that at the feed bulk (about 0.3° C.). This shows how influential the inlet velocities on the TPC and therefor on the distilled water flux, depicted in FIG. 27. The water vapor flux response for a step change in the feed/permeate inlet velocities vs. time is illustrated.

Membrane Mass Transfer Coefficient Identification.

The membrane mass transfer coefficient can be estimated to study the efficiency of the module design and monitor the condition of the membrane during operation. A simple test was devised to identify the membrane mass transfer coefficient in each cell which demonstrates the potential application for this model. Again, ten cells were used for this test, but now two cells had a membrane mass transfer coefficient of 4.1919e-07 (Kg/m² hr Pa) and the remaining eight cells had a mass transfer coefficient of 6.6919e-07 (Kg/m² hr Pa). Ten experiments were simulated with increasing linear velocities of the feed and permeate streams from 0.115 m/s to 0.55 m/s, with a feed inlet temperature of 60° C. and a permeate inlet temperature of 20° C., and the total flux was recorded for each experiment. The total flux of the module can be expressed as the sum of the individual cells fluxes:

$\begin{array}{cc}J=\sum _{n=1}^{10}{B}_{m_n}\Delta P_n.& \left(71\right)\end{array}$ A linear system of equations was set to identify the individual cell membrane mass transfer coefficients, in which the partial vapor pressure difference across the membrane in each cell was calculated from the temperature distribution along the flow direction. By solving the system of equations for B_mn using the truncated singular value decomposition method, the mass transfer coefficient of each cell was successfully identified and in the correct order.

The derivation of a dynamic model for DCMD based on electrical analogy to thermal systems and lumped parameter method was presented. This model was then used to predict and simulate the dynamic response of the water vapor flux under various conditions. The obtained results show great agreement with experimental measurements. The effect of inlet velocities on the process was then investigated. The TPC increased significantly with higher velocities reaching a maximum value of 0.88 at inlet velocities of 0.65 m/s, with asymptotic behavior at high velocities. It was also found that the increase of temperature at the membrane-feed interface is higher than that at the bulk feed. The model was then used to estimate the membrane mass transfer coefficient and identify faults that might be taking place at long operation.

Adaptive Observer for Descriptor Systems

Adaptive observers can be used for joint state and parameter estimation. Among the main motivations to study and design adaptive observers are adaptive control and fault detection and isolation (FDI). Under Lipschitz conditions for the nonlinear terms, one can derive the observer gain based on Lyapunov analysis by solving a set of linear matrix inequalities (LMI). Presented here is an adaptive descriptor observer design for a class of nonlinear descriptor systems with Lipschitz nonlinearities and unknown parameters. Unlike previous methods, the disclosed method does not require coordinate transformation. In addition the adaptive observer gain design can be facilitated by solving a set of LMI, which ensures sufficient conditions for asymptotic convergence of the states and parameters estimations to their true values.

Consider the following descriptor system:E{dot over (x)}=Ax+Bu+Dγ(t,x)θ+ψ(t,y,u),y=Cx,  (72)where x∈ⁿ is the system state, u∈^m is the system input (assumed known), y∈^r is the system output, θ∈^nθ is a vector of unknown time invariant parameters. The matrices E,A∈^n×n, B∈^n×m, D∈^n×q, and C∈^r×n are constants and rank(E)=s<n. The known nonlinear functions γ(t,x)∈^q×nθ and ψ(t,y,u)∈ⁿ are locally Lipschitz.

The following definitions are used for descriptor systems stability and observability:

• • Lemma 1. Let x and y be real vectors of the same dimension. Then, for any scalar ϵ>0 the following inequality holds:2x^Ty≤ϵx^Tx+ϵ⁻¹y^Ty.
  • Definition 1. If det(sE−A)≠0, s∈, then the pair (E, A) is said to be regular. A regular pair (E, A) is impulse-free if deg det(sE−A)=rank(E). If all the roots of the polynomial det(sE−A)=0, s∈ have negative real parts, then the pair (E, A) is stable. If the pair (E, A) is regular, impulse-free, and stable, then it is said to be admissible.
  • Lemma 4. The pair (E, A) is admissible if and only if there exists a matrix X∈^n×n that satisfied the following:E^TX=X^TE≥0,A^TX+X^TA<0.
  • Definition 2. The system of Eq. (72) is R-detectable if and only if its slow subsystem (differential equations) is detectable, i.e.:

$\mathrm{rank}\left[\begin{array}{c}sE-A\\ C\end{array}\right]=m,\forall s\in ,s\mathrm{finite}.$

• • Definition 3. The system of Eq. (72) is called impulse-observable if and only if its fast subsystem (algebraic equations) is impulse-observable, i.e.:

$\mathrm{rank}\left[\begin{array}{cc}E& A\\ 0& E\\ 0& C\end{array}\right]=n+\mathrm{rank}E.$ It is assumed that the system of Eq. (72) is R-detectable and impulse-observable, that for α₁>0, the function γ(t, x) satisfies:∥γ(t,x₁)−γ(t,x₂)∥≤α₁∥x₁−x₂∥  (73)And that the parameter vector θ is piecewise constant and bounded ∥θ∥≤α₂, α₂>0.

The full-order adaptive observer for the system of Eq. (72) is of the form:E{dot over ({circumflex over (x)})}=A{circumflex over (x)}+Bu+Dγ(t,{circumflex over (x)}){circumflex over (θ)}+ψ(t,y,u)+K(y−ŷ)ŷ=C{circumflex over (x)}  (74)where {circumflex over (x)}∈ⁿ and {circumflex over (θ)}∈^nθ are the state and parameter estimates, respectively, and K∈^n×r is the observer gain to be designed. Defining the state estimation error as e₁=x−{circumflex over (x)}, then from Eq. (72) and Eq. (74) the observer error dynamics are given by:Eê₁−(A−KC)e₁+(γ(t,x)θ−γ(t,{circumflex over (x)}){circumflex over (θ)}).  (75)The following theorem gives sufficient conditions for convergence of the adaptive observer error dynamics.

• • Theorem 1. Consider the following parameter estimation adaptation law:{dot over ({circumflex over (θ)})}=Σ⁻¹γ^T(t,{circumflex over (x)})ηC(x−{circumflex over (x)}),
  • where Σ⁻¹=Σ^−T>0 is an arbitrary constant matrix of appropriate dimension, which can be thought of as the adaptation rate, and η∈^q××r such that D^TP=ηC. Then, the system of Eq. (74) is an asymptotically stable adaptive observer for the system of Eq. (72), if there exists a non-singular matrix P∈^n×n, a matrix W∈^r×n such that K=(WP⁻¹)^T, and a scalar ϵ₁>0 for which the following matrix inequalities are solvable:

$\begin{array}{cc}{E}^{T}P=P^{T}E\ge 0,& \left(76\right)\\ \left[\begin{array}{cc}\Omega & P^T\\ P& -{ϵ}_{1}I\end{array}\right]<0,& \left(77\right)\end{array}$ where Ω=A^TP+P^TA−C^TW−W^TC+ϵ₁α₁²α₂²∥D∥²I.

• • Moreover, if the persistent excitation condition holds ∀t₀,∃ξ, δ>0 such that:∫_t0^t0+δDγ(τ,x)γ(τ,x)^TD^Tdτ>ξI,  (78)
  • then the parameter estimation error converges to zero, i.e {circumflex over (θ)}→θ as t→∞.

Proof.

If the system of Eq. (75) is shown to be asymptotically stable, then it follows that the system of Eq. (74) is a full-order adaptive observer for the system of Eq. (72). Let Ā=(A−KC) and define the Lyapunov candidate function for error dynamics as:V(t)=e₁^TE^TPe₁+e₂^TΣe₂,  (79)where e₂=θ−{circumflex over (θ)} is the parameter estimation error. The derivative of V(t) can be given as:

$\begin{array}{cc}\stackrel{.}{V}\left(t\right)=& {\left(E{\stackrel{.}{e}}_1\right)}^T{\mathrm{Pe}}_1+e_1^T{P}^{T}E{\stackrel{.}{e}}_1+2e_2^T\Sigma {\stackrel{.}{e}}_2,\\ =& {\left[\stackrel{\_}{A}{e}_1+D\gamma \left(t,x\right)\theta -D\gamma \left(t,\hat{x}\right)\hat{\theta }\right]}^T{\mathrm{Pe}}_1+\\ & e_1^T{P}^T\left[\stackrel{\_}{A}{e}_1+D\gamma \left(t,x\right)\theta -D\gamma \left(t,\hat{x}\right)\hat{\theta }\right]+\\ & 2e_2^T\Sigma {\stackrel{.}{e}}_2,\\ =& e_1^T\left[{\stackrel{\_}{A}}^{T}P+P^T\stackrel{\_}{A}\right]e_1+\\ & {2\left[D\gamma \left(t,x\right)\theta -D\gamma \left(t,\hat{x}\right)\hat{\theta }\right]}^T{\mathrm{Pe}}_1+2e_2^T\Sigma {\stackrel{.}{e}}_2.\end{array}$ Substituting {circumflex over (θ)}=θ−e₂ in the second term gives:

$\begin{array}{c}\stackrel{.}{V}\left(t\right)=e_1^T\left[{\stackrel{\_}{A}}^{T}P+P^T\stackrel{\_}{A}\right]e_1+\\ {2\left[D\gamma \left(t,x\right)\theta -D\gamma \left(t,\hat{x}\right)\theta \right]}^T{\mathrm{Pe}}_1+\\ {2\left[D\gamma \left(t,\hat{x}\right)e_2\right]}^T{\mathrm{Pe}}_1+2e_2^T\Sigma {\stackrel{.}{e}}_2.\end{array}$ Applying Lemma 1 (from above) to the second term gives:

$\begin{array}{c}\stackrel{.}{V}\left(t\right)⩽e_1^T\left[{\stackrel{\_}{A}}^{T}P+P^T\stackrel{\_}{A}\right]e_1+\\ {ϵ}_1||D\left(\gamma \left(t,x\right)-\gamma \left(t,\hat{x}\right)\right)\theta {||}^2+\\ {ϵ}_1^{-1}{e}_1^T{P}^T{\mathrm{Pe}}_1+\\ {2\left[D\gamma \left(t,\hat{x}\right)e_2\right]}^T{\mathrm{Pe}}_1+2e_2^T\Sigma {\stackrel{.}{e}}_2.\end{array}$ Using Eq. (73) on the second term gives:{dot over (V)}(t)≤e₁^T[Ωϵ₁⁻¹P^TP]e₁+2[Dγ(t,{circumflex over (x)})e₂]^TPe₁+2e₂^TΣė₂,  (80)where Ω=Ā^TP+P^TĀ+ϵ₁α₁²α₂²∥D∥²I.

Since θ is piecewise constant, {dot over (θ)}=0, thus ė₂=−{dot over ({circumflex over (θ)})}. Using the condition D^TP=ηC, it is possible to obtain:2[Dγ(t,{circumflex over (x)})e₂]^TPe₁+2e₂^TΣė₂=0,and the parameter adaptation law can be given by:{dot over ({circumflex over (θ)})}=Σ⁻¹γ^T(t,{circumflex over (x)})ηCe₁.Then, the inequality of Eq. (80) becomes:{dot over (V)}(t)≤e₁^T[Ω+ϵ₁⁻¹P^TP]e₁.  (81)Therefore, a sufficient condition for {dot over (V)}(t)<0 is:Ω+ϵ₁⁻¹P^TP<0.  (82)

The above inequality can be converted to an LMI, using Schur complement, as:

$\left[\begin{array}{cc}\Omega & P^T\\ P& -{ϵ}_{1}I\end{array}\right]<0.$ For some β>0 and from the inequality of Eq. (82), it follows that:Ω+ϵ₁⁻¹P^TP<−βI.  (83)Substituting Eq. (83) into Eq. (81) produces:{dot over (V)}(t)≤−βe₁^Te₁,  (84)which implies that V(t) is bounded, since V(t)≥0 and Eq. (84) shows that V(t) is non-increasing. It follows from Eq. (79) that e₁ and e₂ are also bounded. Integrating both sides of Eq. (84) from t=0 to t=t₁ yields:V(t₁)≤V(0)−β∫₀^t1e₁^T(τ)e₁(τ)d_τ.  (85)From Eq. (85), e₁∈L₂ can be obtained. In addition, the observer error dynamics of Eq. (75) implies that ė₁∈L_∞. Based on this and using Barbalat's lemma, it follows that lim_t→∞e₁=0 and, as a result, lim_t→∞ė₁=0. Therefore, from Eq. (75):

$\begin{array}{cc}\underset{t->\infty }{\lim }D\left(\gamma \left(t,x\right)\theta -\gamma \left(t,\hat{x}\right)\hat{\theta }\right)=0.& \left(86\right)\end{array}$ As lim_t→∞{circumflex over (x)}=x, Eq. (86) reduces to:

$\begin{array}{cc}\underset{t->\infty }{\lim }D\left(\gamma \left(t,x\right)\left(\theta -\hat{\theta }\right)\right)=0.& \left(87\right)\end{array}$ Moreover, the parameter estimation error converges to zero, i.e., {right arrow over (θ)}→θ if the persistent excitation condition of Eq. (78) holds, which is a typical requirement in system identification.

Robustness to Noise.

Consider the noise corrupted descriptor system:E{dot over (x)}=Ax+Bu+Dγ(t,x)θ+ψ(t,y,u)+w(t)y=Cx+v(t)  (88)where w(t)∈ⁿ is the system noise and v(t)∈^r is measurement noise.

• • Theorem 2. If the noises in the system of Eq. (88) w(t) and v(t) are bounded, then the state and parameter estimation errors, e₁(t) and e₂(t) respectively, of the adaptive descriptor observer of Eq. (74) for the system of Eq. (88) are also bounded.

In addition, if the noises w(t) and v(t) have zero mean for all t and are independent of γ(t,x), then as ←∞, the mean of the state and parameter estimation error e₁(t) and e₂(t) converge asymptotically to zero.

Proof.

In this case, the observer error dynamics can be given by:Eė₁=(A−KC)e₁+D(γ(t,x)θ−γ(t,{circumflex over (x)}){circumflex over (θ)})+w(t)−Kv(t),  (89)ė₂=Σ⁻¹γ^T(t,{circumflex over (x)})ηCe₁−Σ⁻¹γ^T(t,{circumflex over (x)})ηv(t).  (90)The homogeneous part of Eq. (89) is exponentially stable, thus the state estimation error e₁(t) can be bounded for γ(t,x) satisfying the local Lipschtiz condition of Eq. (73) and bounded noises w(t) and v(t). In the proof of Theorem 1 above, if was concluded that the asymptotic convergence of parameter estimation error e₂(t) to zero in the noise free case. Therefore, for bounded noises the e₂(t) is also bounded.

In addition, using the assumption that the noises w(t) and v(t) have zero mean for all time and are independent of y(t, x), then the expected value () of the error dynamics for the noise corrupted case (5.35) and (5.36) can be taken as:

$\begin{array}{cc}E\frac{d\left(e_1\right)}{\mathrm{dt}}=\left(A-DC\right)e_1+D\left(\gamma \left(t,x\right)\theta -\gamma \left(t,\hat{x}\right)\hat{\theta }\right),& \left(91\right)\\ \frac{d\left(e_2\right)}{\mathrm{dt}}=-{\sum }^{-1}{\gamma }^T\left(t,\hat{x}\right)\eta C{e}_1,& \left(92\right)\end{array}$ where the order of differentiation and expectation operations have been interchanged. Notice that the dynamics of E(e₁) and (e₂), of Eq. (91) and Eq. (92) respectively, are similar to the noise free case. Therefore, the state and parameter estimation errors converge to zero asymptotically in the mean sense.

The adaptive observer design is demonstrated through the following example. Consider the system of Eq. (72) with:

$E=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right],A=\left[\begin{array}{cccc}0& -4& 0& 3\\ 1& 2& -4& 0\\ 0& 2& 0& -5\\ -2& 0& 3& 0\end{array}\right]$ $B=\left[\begin{array}{cc}-1& 0\\ 0& 0\\ 0& 0\\ 0& 1\end{array}\right],D=\left[\begin{array}{cc}-1& 0\\ -1& 0\\ 0& 0\\ 0& 1\end{array}\right],C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]$ $\gamma =\left[\begin{array}{cc}{\sin }^2\left(0.125t\right)& x_1^2+\sqrt[3]{x_4}\\ x_3^3+x_2& {\cos }^2\left(0.125t\right)\end{array}\right]$ Setting the design parameters α₁=1 and α₂=4 and using the YALMIP toolbox to solve the LMI sufficient conditions of Eq. (76) and Eq. (77), the following results can be obtained:

$K=\left[\begin{array}{cc}-8.808& 6.250\\ 1.782& 2.788\\ 1.267& 1.650\\ 7.648& 0\end{array}\right],\eta =\left[\begin{array}{cc}0& -0.05937\\ 0.05791& 0.05845\end{array}\right]$ The initial states of the system are:x(0)=[−1.6624 0.8 2.8 0.2]^T,{circumflex over (x)}(0)=[−1.5017 1.0 2.25 0.5]^T,and the initial state of the adaptive observer is:{circumflex over (θ)}(0)=[0.4−3.5]^T.The adaptation rate was chosen as Σ⁻¹=diag([50,30]). For this simulation, the true values of the parameters θ₁ and θ₇ switch from 0.1 to 0.3 and from −4 to −2, respectively, at t=40 and back again at t=60.

The state and parameters estimation results for the noise-free systems are shown in FIGS. 28A and 28B, respectively. FIG. 28A illustrates examples of the observer response showing the two immeasurable states for the noise free system, with x₃ and {circumflex over (x)}₃ plotted in the top graph and with x₄ and {circumflex over (x)}₄ plotted in the bottom graph. FIG. 28B illustrates examples of the estimates of the unknown parameters, with θ₁ and {circumflex over (θ)}₁ plotted in the top graph and θ₂ and {circumflex over (θ)}₂ plotted in the bottom graph.

In addition, Gaussian noise with standard deviation of 0.05 was added to the outputs of the system y(t). Examples of the response of the system unknown states and their estimates are shown in FIG. 28C. The observer response illustrates the two immeasurable states for the noise corrupted system, with x₃ and {circumflex over (x)}₃ plotted in the top graph and with x₄ and {circumflex over (x)}₄ plotted in the bottom graph. Examples of the estimates of two unknown parameters are depicted in FIG. 28D. The estimates of the unknown parameters are illustrated, with θ₁ and {circumflex over (θ)}₁ plotted in the top graph and θ₂ and {circumflex over (θ)}₂ plotted in the bottom graph.

It can be observed that the estimates of the states and parameters are robust to abrupt changes to the system input signals. The estimate of the parameters reconverges to the new value in about 10 s. After each parameter change, the state estimation error e_i(t) converges to zero, in the mean sense, faster than the parameter error e₂(t). FIG. 28E illustrates examples of the state estimation error for the noise corrupted case.

The convergence of the descriptor adaptive observer was obtained by solving a set of sufficient linear matrix inequalities. The robustness under the noise corrupted system and measurements was established. For zero mean noises, the estimation errors of the states and parameters converge, in the mean, to zero. The performance of the method was demonstrated through a numerical example. This method can be used to systematically design adaptive observers for systems with a large number of states, which can be useful for fault detection and adaptive control applications.

Membrane Fouling Detection in DCMD

Like many membrane-based desalination processes, MD is prone to membrane fouling and scaling, which is a process of accumulated deposition of clay, salt, organic and/or biological materials onto the membrane surface or into the membrane pores. Membrane fouling affects the water quality and reduces the production rate and may even cause system shutdown for cleaning, which in turn, increases the operational cost of the MD system. In MD, the fouling layer deposited on the surface of the membrane introduces a thermal resistance to heat transfer between the bulk feed stream and the membrane surface. This can result in a lower temperature at the feed-membrane interface and hence a lower partial vapor pressure, which ultimately reduces the production rate of the system. FIG. 29A is a schematic diagram illustrating an example of the fouling layer is DCMD setups.

Several methods have been proposed for detecting membrane fouling in various membrane-based processes. Because incorporating sensors into the process can be challenging and expensive, using mathematical tools and readily available measurements to detect fouling is ideal for this application. Observers can be applied for fault detection and isolation (FDI). A general class of models are descriptor systems (known also as differential algebraic systems, DAE, or singular systems), which appear naturally when modeling electrical networks, constrained mechanical systems, and interconnected processes. Observers for linear descriptor systems can be designed using a system of ODE, which asymptotically converges to the true states by eigenvalue assignment. This approach can be extended for fault estimation for linear descriptor systems with unknown inputs and linear parameter varying systems. For Lipschitz nonlinear descriptor systems, a stabilizing gain for a Luenberger-like observer can be designed by solving a set of linear matrix inequalities (LMI).

The descriptor model for direct contact membrane distillation (DCMD) of FIG. 17 can be extended to account for membrane fouling. An adaptive observer can be used for the detection of membrane fouling. The disclosed method does not require coordinate transformation. In addition, the adaptive observer gain design can be facilitated by solving a set of LMI, which guarantees sufficient conditions for asymptotic convergence of the states and parameter estimations to their true values.

Membrane Fouling Modeling.

In MD, the fouling layer deposited on the surface of the membrane introduces a thermal resistance to heat transfer between the bulk feed stream and the membrane surface, as shown in FIG. 29A. This results in a lower temperature at the feed-membrane interface and hence a lower partial vapor pressure, which ultimately reduces the production rate of the system. The fouling thermal resistance (R_sc) can be incorporated into the model of FIG. 17. FIG. 29B shows an example of the electrical analogue of the nth cell of the DCMD module including the thermal resistance of the membrane fouling layer. Rewriting the heat transfer equations to account for the fouling resistance gives:

$\begin{array}{cc}\frac{{\mathrm{dT}}_{{\mathrm{bf}}_n}}{\mathrm{dt}}=\frac{1}{C_{\mathrm{bf}}}{Q}_{f_n}-\frac{1}{C_{\mathrm{bf}}}\left(\frac{1}{R_f}+J_n{A}_m{c}_p\right)T_{{\mathrm{bf}}_n}+\frac{1}{C_{\mathrm{bf}}}\left(\frac{R_{\mathrm{sc}}\left(1+R_f{G}_n{c}_p\right)}{R_f+R_{\mathrm{sc}}\left(1+R_f{G}_n{c}_p\right)}\right)T_{{\mathrm{bf}}_n}-\frac{1}{C_{\mathrm{bf}}}{Q}_{f_{n+1}}+\frac{1}{C_{\mathrm{bf}}\left(R_f+R_{\mathrm{sc}}\left(1+R_f{G}_n{c}_p\right)\right)}{T}_{{\mathrm{mf}}_n},& \left(93\right)\\ \left(\frac{1}{R_f}+J_n{A}_m{c}_p+\frac{R_{\mathrm{sc}}\left(1+R_f{G}_n{c}_p\right)}{R_f+R_{\mathrm{sc}}\left(1+R_f{G}_n{c}_p\right)}\right)T_{{\mathrm{bg}}_n}-\frac{1}{\left(R_f+R_{\mathrm{sc}}\left(1+R_f{G}_n{c}_p\right)\right)}{T}_{{\mathrm{mf}}_n}=\left(\frac{1}{R_p}+J_n{A}_m{c}_p\right)T_{m{p}_n}-\frac{1}{R_p}{T}_{{\mathrm{bp}}_u}.& \left(94\right)\end{array}$ The model can now be represented as:

$\begin{array}{cc}E\stackrel{.}{x}=\mathrm{Ax}+\mathrm{Bu}+D\gamma \left(t,x\right)B_m+H\left(x,R_{\mathrm{sc}}\right),y=\left[\begin{array}{c}\mathrm{Cx}\\ G\end{array}\right],& \left(95\right)\end{array}$ where H(x, R_sc) is a nonlinear function of the state x and the fouling thermal resistance R_sc.

Membrane Fouling Detection.

Expanding the function H(x, R_sc) around R_sc=0 and taking the linear terms, the system of Eq. (95) can be written as:

$\begin{array}{cc}E\stackrel{.}{x}=A\hat{x}+\mathrm{Bu}+D_1\gamma \left(t,\hat{x}\right)B_m+D_2\tilde{H}\left(x\right)\Delta {\hat{R}}_{\mathrm{sc}},\hat{y}=\left[\begin{array}{c}C\hat{x}\\ \hat{G}\end{array}\right],& \left(96\right)\end{array}$ where ΔR_sc is the deviation from the nominal value of zero fouling thermal resistance. Consider the following adaptive observer form for the system of Eq. (96):

$\begin{array}{cc}E\stackrel{\stackrel{.}{^}}{x}=A\hat{x}+\mathrm{Bu}+D_1\gamma \left(t,\hat{x}\right)B_m+D_2\tilde{H}\left(\hat{x}\right)\Delta {\hat{R}}_{\mathrm{sc}}+K\left(\mathrm{Cx}-C\hat{x}\right),\hat{y}=\left[\begin{array}{c}C\hat{x}\\ \hat{G}\end{array}\right]& \left(97\right)\end{array}$ Define the state estimation by e₁=x−{circumflex over (x)}, then from Eq. (96) and Eq. (74), the observer error dynamics can be given by:Eė₁=(A−KC)e₁+D(γ(t,x)θ−γ(t,{circumflex over (x)}){circumflex over (θ)}).  (98)

The following theorem gives sufficient conditions for convergence of the adaptive observer error dynamics.

• • Theorem 3. Consider the following adaptation law:{dot over ({circumflex over (R)})}=Σ⁻¹{tilde over (H)}({dot over (x)})ηC(x−{circumflex over (x)})+Σ⁻¹L(G−Ĝ)  (99)
  • where Σ⁻¹=Σ^−T>0 is an arbitrary constant matrix of appropriate dimension, which can be thought of as the adaptation rate. Then the system of Eq. (97) is an asymptotically stable adaptive observer for the system of Eq. (96), if there exists a non-singular matrix P∈^n×n, a matrix W∈^r×n such that K=(WP⁻¹)^T, and a scalar ϵ₁>0 for which the following matrix inequalities are solvable:

$E^{T}P=P^{T}E\ge 0,\left[\begin{array}{cc}\Omega & P^T\\ P& -{ϵ}_{1}I\end{array}\right]<0,$

• • where Ω=A^TP+P^TA−C^TW−W^TC+ϵ₁α₁²α₂²∥D₂∥²I.

Moreover, if the persistent excitation condition holds ∀t₀, ∃ξ, δ>0 such that:∫_t0^t0+δD₂{tilde over (H)}({circumflex over (x)})^TD₂^Tdτ>ξI,  (100)

• • then the parameter estimation error converges to zero, i.e Δ{circumflex over (R)}_sc→ΔR_sc as t→∞.

Proof.

If the system of Eq. (75) is shown to be asymptotically stable, then it follows that the system of Eq. (74) is a full-order adaptive observer for the system of Eq. (96). Let Ā=(A−KC) and define the Lyapunov candidate function for error dynamics as:V(t)=e₁^TE^TPe₁+e₂^TΣe₂,  (101)where e₂θ−{circumflex over (θ)} is the parameter estimation error. The derivative of V(t) can be given as:

$\begin{array}{c}\stackrel{.}{V}\left(t\right)={\left(E{\stackrel{.}{e}}_1\right)}^T{\mathrm{Pe}}_1+e_1^T{P}^{T}E{\stackrel{.}{e}}_1+2e_2^T\sum {\stackrel{.}{e}}_2,\\ ={\left[\stackrel{\_}{A}{e}_1+D_2\tilde{H}\left(x\right)\Delta R_{\mathrm{sc}}-D_2\tilde{H}\left(x\right)\Delta {\stackrel{.}{R}}_{\mathrm{sc}}\right]}^T{\mathrm{Pe}}_1+\\ e_1^T{P^T\left[\stackrel{\_}{A}{e}_1+D_2\tilde{H}\left(x\right)\Delta R_{\mathrm{sc}}-D_2\tilde{H}\left(\hat{x}\right)\Delta {\hat{R}}_{\mathrm{sc}}\right]}^T{\mathrm{Pe}}_1+\\ 2e_2^T\sum {\stackrel{.}{e}}_2,\\ =e_1^T\left[{\stackrel{\_}{A}}^{T}P+P^T\stackrel{\_}{A}\right]e_1+\\ {2\left[D_2\tilde{H}\left(x\right)\Delta R_{\mathrm{sc}}-D_2\tilde{H}\left(\hat{x}\right)\Delta {\hat{R}}_{\mathrm{sc}}\right]}^T{\mathrm{Pe}}_1+2e_2^T\sum {\stackrel{.}{e}}_2.\end{array}$ Substituting Δ{circumflex over (R)}_sc=ΔR_sc−e₂ in the second term gives:

$\stackrel{.}{V}\left(t\right)=e_1^T\left[{\stackrel{\_}{A}}^{T}P+P^T\stackrel{\_}{A}\right]e_1+{2\left[D_2\tilde{H}\left(x\right)\Delta R_{\mathrm{sc}}-D_2\tilde{H}\left(\hat{x}\right)\Delta R_{\mathrm{sc}}\right]}^T{\mathrm{Pe}}_1+{2\left[D_2\tilde{H}\left(\hat{x}\right)e_2\right]}^T{\mathrm{Pe}}_1+2e_2^T\sum {\stackrel{.}{e}}_2$ Applying Lemma 1 (from above) to the second term gives:

$\stackrel{.}{V}\left(t\right)⩽e_1^T\left[{\stackrel{\_}{A}}^{T}P+P^T\stackrel{\_}{A}\right]e_1+{ϵ}_1{D_2\left(\tilde{H}\left(x\right)-\tilde{H}\left(\hat{x}\right)\right)\Delta R_{\mathrm{sc}}}^2+{ϵ}_1^{-1}{e}_1^T{P}^T{\mathrm{Pe}}_1+{2\left[D_2\tilde{H}\left(\hat{x}\right)e_2\right]}^T{\mathrm{Pe}}_1+2e_2^T\sum {\stackrel{.}{e}}_2.$ Using Eq. (73) on the second term gives:{dot over (V)}(t)≤e₁^T[Ω+ϵ₁⁻¹P^TP]e₁+2[D₂{tilde over (H)}({circumflex over (x)})e₂]^TPe₁+2e₂^TΣė₂,  (102)where Ω=Ā^TP+P^TĀ+ϵ₁α₁²α₂²∥D₂∥²I.

Since ΔR_sc is piecewise constant, Δ{dot over (R)}_sc=0, thus ė₂=−_sc. Using the condition D^TP=ηC and substituting the adaptation law of Eq. (99), it is possible to obtain:2[D₂{tilde over (H)}({circumflex over (x)})e₂]^TPe₁+2e₂^TΣė₂=−e₂^TL(G−Ĝ),where L∈⁻ is a negative constant gain, then the term −e₂^TL(G−Ĝ) is always non-positive as it can be shown:

$\hspace{1em}\{\begin{array}{c}\mathrm{if}{e}_2<0,->\left(G-\hat{G}\right)>0,\mathrm{therefore}-e_2^{T}L\left(G-\hat{G}\right)<0\\ \mathrm{if}{e}_2>0,->\left(G-\hat{G}\right)<0,\mathrm{therefore}-e_2^{T}L\left(G-\hat{G}\right)<0\end{array}$ Then, the inequality of Eq. (102) becomes:{dot over (V)}(t)≤e₁^T[Ω+ϵ₁⁻¹P^TP]e₁.  (103)Therefore, a sufficient condition for {dot over (V)}(t)<0 is:Ω+ϵ₁⁻¹P^TP<0,  (104)

The above inequality can be converted to an LMI, using Schur complement, as:

$\left[\begin{array}{cc}\Omega & P^T\\ P& -{ϵ}_{1}I\end{array}\right]<0.$ For some β>0 and from the inequality of Eq. (104), it follows that:Ω+ϵ₁⁻¹P^TP<−I.  (105)Substituting Eq. (105) into Eq. (103) produces:{dot over (V)}(t)≤−e₁^Te₁,  (106)which implies that V(t) is bounded, since V(t)≥0 and Eq. (106) shows that V(t) is non-increasing. It follows from Eq. (101) that e₁ and e₂ are also bounded. Integrating both sides of Eq. (106) from t=0 to t=t₁, yields:V(t₁)≤V(0)−β∫₀^t1e₁^T(τ)e₁(τ)e₁(τ)dτ.  (107)From Eq. (107), e₁∈L₂ can be obtained. In addition, the observer error dynamics of Eq. (75) implies that ė₁∈L_∞. Based on this and using Barbalat's lemma, it follows that lim_t→∞e₁=0 and, as a result, lim_t→∞ė₁=0. Therefore, from Eq. (75):

$\begin{array}{cc}\underset{t\to \infty }{\lim }{D}_2(\tilde{H}\left(x\right)\Delta R_{\mathrm{sc}}-\tilde{H}\left(\hat{x}\mathrm{\_\Delta }{\hat{R}}_{\mathrm{sc}}\right)=0.& \left(108\right)\end{array}$ As lim_t→∞{circumflex over (x)}=x, Eq. (108) reduces to:

$\begin{array}{cc}\underset{t\to \infty }{\lim }{D}_2\left(\tilde{H}\left(x\right)\left({\mathrm{\Delta R}}_{\mathrm{sc}}-\Delta {\hat{R}}_{\mathrm{sc}}\right)\right)=0.& \left(109\right)\end{array}$ Moreover, the parameter estimation error converges to zero, i.e., Δ{circumflex over (R)}_sc→ΔR_sc if the persistent excitation condition of Eq. (100) holds, which is a typical requirement in system identification.

To demonstrate the performance of the disclosed adaptive observer, a simulation of gradual membrane fouling was developed. For the first 25 seconds, the fouling thermal resistance (ΔR_sc) was set to zero. Then, a gradual increase took place over the duration of 50 seconds, until the fouling thermal resistance reached 2:5×10⁻³ and stayed there for another 25 seconds. Referring to FIG. 30A, shown is an example of the actual and estimated profile of the thermal fouling resistance (ΔR_sc). It is clear that the adaptive observer is able to track the membrane fouling as it develops. The estimated and actual profile of the modules production rate (G) over time is shown in FIG. 30B, where the effect of membrane fouling on the production rate is clear. To confirm the response of the temperature distribution estimation, FIG. 30C shows the temperature spatial profiles (distribution) through the DCMD module for both the estimated and actual simulated temperature spatial profiles for both the feed and permeate sides of the module at the final time t=100 s.

Nonlinear Observer-Based Lyapunov Control

With no loss of generality, consider a DCMD module, where an example of a counter current DCMD module is depicted in FIG. 1. In this module, the vapor generated in the feed solution is forced to pass through the membrane dry pores to the permeate side, following thermodynamics rules. Accordingly, the heat transfer is driven by the trans-membrane heat flux and the latent heat of evaporation. In this model, there are two inputs to the DCMD process which are the feed and the permeate inlet temperatures, and two measured outputs for the DCMD process which are the feed and the permeate outlet temperatures.

Mechanical Model.

The mathematical mode of the DCMD model can be given by the following 2D ADE:

$\begin{array}{cc}\frac{\partial T\left(x,z,t\right)}{\partial t}+v_z^k\frac{\partial T\left(x,z,t\right)}{\partial z}={\alpha }_x^k\frac{{\partial }^{2}T\left(x,z,t\right)}{\partial x^2},& \left(110\right)\end{array}$ with k={f,p}, 0<t<τ_final, x∈[0,x_fm]∪[x_pm,X], 0<z<Z where k={f,p} refers to the feed and the permeate, respectively. The terms and, k={f,p} are the feed and permeate feed pumps flow rates along the x and z coordinates, respectively. The initial conditions of the model are:T(x,z,0)=T_finitial,T(x,z,0)=T_pinitial where the associated boundary conditions with the dynamic model (e.g., feed solution) are set to be a Dirichlet condition for the inlet boundary of each subsystem, and Neumann conditions for the rest of the boundaries like Eqs. (10)-(17),

$T\left(x,0,t\right)=T_{0_f}\mathrm{for}x\in \left[0,x_{\mathrm{fm}}\right],T\left(x,Z,t\right)=T_{0_p}\mathrm{for}x\in \left[x_{\mathrm{pm}},X\right],\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=0}=0,\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=X}=0,\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=x_{\mathrm{fm}}}=\left[J\left(T\right)H\left(T\right)+\frac{k_m}{{\delta }_m}\left(T\left(x_{\mathrm{pm}},z,t\right)-T\left(x_{\mathrm{pm}},z,t\right)\right)\right]\text{/}{k}_f,\frac{\partial T\left(x,z,t\right)}{\partial x}{|}_{x=x_{\mathrm{pm}}}=\left[J\left(T\right)H\left(T\right)+\frac{k_m}{{\delta }_m}\left(T\left(x_{\mathrm{fm}},z,t\right)-T\left(x_{\mathrm{pm}},z,t\right)\right)\right]\text{/}{k}_p,\frac{\partial T\left(x,z,t\right)}{\partial z}{|}_{z=Z}=\varphi ,\frac{\partial T\left(x,z,t\right)}{\partial z}{|}_{z=0}=\varphi ,$ where k_m is the average thermal conductivity of the membrane, δ_m is the membrane thickness, ϕ∈⁺ is a positive constant to represent the amount of flux that is released outside the process containers, and H(T) is the latent heat of evaporation that is responsible of the released sensible heat as given by Eq. (18). J(T) is the permeate mass flux to be transferred between the feed and the permeate solutions in the process based on Knudsen diffusion model for the mass transfer coefficient as given by Eq. (2). ΔP(T)=P₁(T)−P₀(T) is the vapor pressure gradient between the membrane sides. It is a function of the vapor temperature based on the Antoine equation.

Semi-Discretized Model.

The ADE model is semi-discretized (space only) using finite difference discretization scheme. Particularly, the convective and diffusive terms are solved with first-order upwind scheme, and central difference scheme, respectively. The upwind scheme warrants the stability and ensures the solution to satisfy the entropy condition, while the central difference scheme is suitable for second derivative diffusive term. The semi-discretization transforms the infinite dimensional partial differential equations model into a finite dimensional group of ordinary differential equations. It provides good model approximation by refining the mesh grid, and tuning the spatial step sizes. Furthermore, standard control strategies that are applicable for state space representations, can be extended to such systems. The relations between the mesh points in the process designated domain are significant to solve the model. These relations are shown in FIG. 2, where the discretization starts in a row-wise manner providing a solution for the feed then the permeate.

The semi-discrete version of the continuous 2D ADE model is a nonlinear high dimension system of ODEs because of the nonlinearity of the boundary conditions. The discrete 2D ADE model can be expressed as:

$\begin{array}{cc}\{\begin{array}{c}\stackrel{.}{T}\left(t\right)=\mathrm{AT}\left(t\right)+\mathrm{Bu}\left(T\right)+f\left(T\left(t\right)\right),\\ y\left(t\right)=\mathrm{CT}\left(t\right),\\ T\left(0\right)=T_0.\end{array}& \left(111\right)\end{array}$ where T∈^2n2×1 is the state vector, which represents the time evolution of the temperature distribution at each grid point inside the feed and the permeate designated domain T_f,T_p, respectively.

$\begin{array}{cc}T=\left[\begin{array}{c}{T}_f\\ T_p\end{array}\right],& \left(112\right)\end{array}$ where T_f,T_p∈^n2×1 and n is the number of discretization points in x and z axes, and T₀ is the initial condition of the state vector. Matrix A∈^2n2×2n2 represents the dynamics matrix which includes all dynamics of the feed and the permeate with their correspondent couplings on the boundaries, as in the following:

$\begin{array}{cc}A=\left[\begin{array}{cc}{A}_{f_{i,j}}& C_{\mathrm{jp}}\\ C_{\mathrm{pf}}& A_{p_{i,j}}\end{array}\right],& \left(113\right)\end{array}$ where i, j are the discretization grid indices.

$\begin{array}{cc}\{\begin{array}{c}{A}_{f_{i,j}\left(u_f\right)}=\left(-\frac{u_f}{\Delta z}-\frac{2{\alpha }_x^f}{\Delta x^2}\right)T_{i,j}+\left(\frac{u_f}{\Delta z}\right)T_{i,j-1}+\\ \left(\frac{{\alpha }_x^f}{\Delta x^2}\right)T_{i-1,j}+\left(\frac{{\alpha }_x^f}{\Delta x^2}\right)T_{i+1,j},\\ A_{p_{i,j}\left(u_p\right)}=\left(\frac{u_p}{\Delta z}-\frac{2{\alpha }_x^p}{\Delta x^2}\right)T_{i,j}-\left(\frac{u_p}{\Delta z}\right)T_{i,j+1}+\\ \left(\frac{{\alpha }_x^p}{\Delta x^2}\right)T_{i-1,j}+\left(\frac{{\alpha }_x^p}{\Delta x^2}\right)T_{i+1,j},\\ C_{\mathrm{fp}}=\frac{{\alpha }_z}{\Delta x^2}\frac{\Delta {\mathrm{xk}}_m}{{\delta }_m{k}_f},\\ C_{\mathrm{pf}}=\frac{{\alpha }_x}{\Delta x^2}\frac{\Delta {\mathrm{xk}}_m}{{\delta }_m{k}_p}.\end{array}\hspace{1em}& \left(114\right)\\ B=\left[\begin{array}{cc}\left(\frac{v_z^f}{\Delta z}\right)& 0\\ ⋮& ⋮\\ {\left(\frac{v_z^f}{\Delta z}\right)}_n& 0\\ 0& 0\\ ⋮& ⋮\\ 0& 0\\ 0& {\left(-\frac{v_z^p}{\Delta z}\right)}_1\\ ⋮& ⋮\\ 0& {\left(-\frac{v_z^p}{\Delta z}\right)}_n\end{array}\right],u\left(t\right)=\left[\begin{array}{c}{T}_{0_f}\\ T_{0_p}\end{array}\right],& \left(115\right)\end{array}$ where B is the input matrix and u(t) includes the inputs of the DCMD process. f(T)∈^2n2×1 contains the different boundary conditions of the feed and the permeate, especially the ones associated with the membrane boundaries. Some of the boundary conditions can be nonlinear and functions of the state vector T as shown in Eqs. (10)-(17). Vectors ç and η in the following equations represent the polynomial nonlinearity structure in f(T).

$\begin{array}{cc}f\left(T\right)={\left[{\zeta }_1\cdots {\zeta }_n{\eta }_1\cdots {\eta }_n\right]}^{\mathrm{tr}},& \left(116\right)\\ \zeta ={\left[0\cdots 0\frac{{\alpha }_x^k}{\Delta x}\frac{\left(T\right)H\left(T\right)}{k_f}\right]}^{\mathrm{tr}},& \left(117\right)\\ \eta ={\left[-\frac{{\alpha }_x^k}{\Delta x}\frac{\left(T\right)H\left(T\right)}{k_p}0\cdots 0\right]}^{\mathrm{tr}},& \left(118\right)\end{array}$ where tr refers to the transpose and (T) is the permeate mass flux to be transferred between the feed and the permeate solutions in the process. The two measured outlet temperatures in the output vector are:

$\begin{array}{cc}y=\left[\begin{array}{c}{T}_{\mathrm{fout}}\\ T_{\mathrm{pout}}\end{array}\right],& \left(119\right)\end{array}$ where y∈^2×1 and the output matrix is:

$\begin{array}{cc}C=\left[\begin{array}{ccc}\stackrel{\stackrel{n^2}{︷}}{0\cdots 0}& \stackrel{\stackrel{n-1}{︷}}{10\cdots 0}& \stackrel{\stackrel{n\left(n-1\right)+1}{︷}}{00\cdots 0}\\ 0\cdots \cdots & 0\cdots \cdots 0& 10\cdots 0\end{array}\right].& \left(120\right)\end{array}$ where C∈^2×2n2. The boundary condition vector f(T) has some properties, such as being nonlinear, continuous, smooth, bounded and Lipschitz bounded.

Nonlinearity Analysis.

The DCMD dynamic model is nonlinear, and the source of the nonlinearity results from the nature of the functions (T), H(T) and their multiplication in the boundary conditions. However, the present nonlinearity is Lipschitz continuous. This result can be proved by showing that each of H(T) and (T) is separately Lipschitz continuous, and then their multiplication is going to be directly Lipschitz continuous. For H(T), because it is a polynomial function that is smooth and infinitely differentiable, it is Lipschitz continuous. Then, all the constants can be collected in the expression of (T) and can be written in a simplified form as:

$\begin{array}{cc}\left(T\right)={c_1\left(\frac{1}{T}\right)}^{\left(\frac{1}{2}\right)}{10}^{\left(\frac{-c_2}{T}\right)},& \left(121\right)\end{array}$ where c₁,c₂ are constants, and T∈[a,b], b>a>0. Therefore, (r) is continuous and differentiable on this interval [a,b], and therefore it is bounded:∥(T)∥≤Υ,  (122)where Υ is a real number. As a result, there exists a Lipshitz constant λ that satisfies:∥f(T₁)−f(T₂)∥≤λ∥T₁−T₂∥.  (123)The boundedness of f(T) can be explained by the physical aspects and properties or the process.

Controller Design.

Nonlinear processes are challenging in stability and control, due to the high sensitivity to initial and forcing conditions. The operation of heat transfer, inside the MD module, is nonlinear and needs a sufficient temperature difference along the membrane boundaries to operate efficiently. Controlling this temperature difference to manage the water production inside the DCMD process will now be discussed.

State Transformation.

The structure of the DCMD process model in Eq. (111) is high dimensional order, and would result in many of its properties being mesh size dependent. For example, a 10×10 discretization mesh would generate more temperature states than a 5×5 mesh, and accordingly the model may not be controllable and observable for higher mesh sizes such as 18×18 mesh. In this disclosure, only the affecting states for control, such as the states of the membrane boundaries, are considered because of their significant role in driving the process.

In this context, the controllability and the observability of the states on the membrane boundaries can be examined. This would reduce the number of involved states hugely, so the following change of variables can be used:{tilde over (T)}=T,  (124)where {tilde over (T)}∈^n×1 is the temperature difference vector along the membrane boundaries, and can be defined as follows:{tilde over (T)}=[{tilde over (T)}₁ . . . {tilde over (T)}_n]^tr.  (125)∈^n×2n2 is a non-square transformation matrix that form the difference of each correspondent boundaries, such as in:

$\begin{array}{cc}=\left[\begin{array}{ccccccccccccc}0& \cdots & 1& 0& \cdots & 0& -1& \cdots & 0& \cdots & \cdots & \dots & \dots \\ 0& \cdots & \cdots & \cdots & 1& 0& \cdots & \cdots & -1& 0& \cdots & \cdots & \dots \\ 0& \cdots & \cdots & \cdots & \ddots & 0& 1& 0& \ddots & 0& -1& 0& \cdots \end{array}\right].& \left(126\right)\end{array}$ Consider the pseudo inverse ^pseudo=(^T·)^−1T, which is the inverse matrix of the non-square transformation matrix . The dynamics of {tilde over (T)} are:

$\begin{array}{cc}\{\begin{array}{c}\stackrel{.}{\stackrel{~}{T}}=\stackrel{.}{T},\\ =\tilde{T}+u+\left(T\right),\\ y=\tilde{T}.\end{array}& \left(127\right)\end{array}$ where =A^pseudo∈^n×n represents the dynamic matrix of the membrane boundaries, and has the following structure:

$\begin{array}{cc}=\left[\begin{array}{c}{1}_{}\\ 2_{}\\ ⋮\\ n_{}\end{array}\right],& \left(128\right)\end{array}$ where _i, i=1, 2, . . . , n are the rows of matrix , =B∈^n×2 is the new input matrix and has a structure such that:

$\begin{array}{cc}={\left[\begin{array}{cccc}{b}_1& 0& \cdots & 0\\ 0& \cdots & 0& b_2\end{array}\right]}^{\mathrm{tr}}& \left(129\right)\end{array}$ (T)=f(T)∈^n×1 is the nonlinear effect on the membrane boundaries with a structure of the following form:(T)=[₁(T) . . . _n(T)]^tr.  (130)=C^pseudo∈^1×n is the output matrix as defined earlier.

Accordingly, the structure of the DCMD process has been reduced, and therefore, it will be easier to examine and guarantee the controllability and the observability of the states on the membrane boundaries through the pairs () and (), respectively. Numerical simulations will show that the controllability and the observability conditions for the DCMD process are satisfied.

Nonlinear Lyapunov Control Design.

The efficiency of the DCMD process appears clearly when the production rate of the clean water is maintained almost constant with relatively low energy consumption. In this regard, Lyapunov control design can be used to maintain the temperature difference along the membrane boundaries to be almost constant, and equal to a desired level. The control will be nonlinear due to the present nonlinearity in the process. Hereby, nonlinear Lyapunov control guarantees the boundedness of the trajectories, and therefore, convergence for the states. Additionally, the control inputs should be constrained to be within an acceptable ranges for energy consumption considerations.

Satisfying such criteria and for practical purposes, we define the inlet temperatures for the feed and the permeate to be bounded with minimum and maximum temperature values equal to 20° C. and 60° C., respectively.0<U_min=20° C.<u(t)<U_max=60° C., ∀t>0.  (131)Following the reduced DCMD model in Eq. (127), Proposition 1 suggests a control law for the reference tracking error.

• • Proposition 1. Assume the temperature difference of the membrane boundaries obeys Eq. (127), then the following control law (u=[u₁ u₂]^tr):

$\begin{array}{cc}\{\begin{array}{c}{u}_1=\frac{1}{b_1}\left[-K_1\left({\tilde{T}}_1-y_r\right)-1_{}\tilde{T}-1_{}\left(T\right)\right]\\ u_2=\frac{1}{b_{21}}\left[-K_2\left({\tilde{T}}_n-y_r\right)-n_{}\tilde{T}-n_{}\left(T\right)\right]\end{array},& \left(132\right)\end{array}$

• • ensure the reference tracking, such that the temperature difference {tilde over (T)} follows the reference y_r, where K₁ and K₂ are positive control gains.

Proof.

The proof is divided into two parts. First, the derivation of the control laws is shown, and second, the convergence of the temperature difference to the reference y_r is proven. To begin, define the error function as follows:e={tilde over (T)}−y_r=[e₁e₂ . . . e_n]^tr,  (133)then restrict the subsystem to

$\left(\begin{array}{c}{e}_1\\ e_2\end{array}\right),$ since this is a boundary control problem; i.e. u₁ and u₂ act only on the first and the last dynamical equations. Next, define a Lyapunov function such as:

$\begin{array}{cc}\{\begin{array}{c}V\text{:}{n}^{}\to ,\\ V=\frac{1}{2}{e}_2^2+\frac{1}{2}{e}_n^2,\end{array}& \left(134\right)\end{array}$ for which the derivative can be expressed as:{dot over (V)}=e₁ė₁+e_nė_n,{dot over (V)}=e₁({dot over ({tilde over (T)})}₁−{dot over (y)}_r)+e_n({tilde over ({dot over (T)})}_n−{dot over (y)}_r),{dot over (V)}=e₁(₁{tilde over (T)}₁+b₁u₁+₁(T))+e_n(_n{tilde over (T)}_n+b₂u₂+_n(T)).   (135)To ensure that the derivative of the Lyapunov function in Eq. (135) is negative definite, we propose the control law in Eq. (132) which guarantees that:{dot over (V)}(t)=−K₁e₁²−K₂e_n²,  (136)where K₁ and K₂ are positive constants. Hence, the control law will be defined as in Eq. (132).

The second part of the proof shows the stability and convergence of the whole system using the control law defined in Eq. (132). It starts by substituting u₁ and u₂ into the dynamics of the error equation in Eq. (133), such as in the following form:ė=e+(T)+y_r,  (137)where (T) has a structure in the form of:(T)=[0(T) . . . _n−1(T)0]^tr,  (138)and and are such that:

$\begin{array}{cc}=\left[\begin{array}{c}0\\ 2_{}\\ ⋮\\ {n-1}_{}\\ 0\end{array}\right],=\left[\begin{array}{c}-{\tilde{K}}_1\\ 2_{}\\ ⋮\\ {n-1}_{}\\ -{\tilde{K}}_2\end{array}\right],& \left(139\right)\end{array}$ where K₁=[K₁ 0 . . . 0]∈^1×n, and K₂=[0 . . . 0 K₂]∈^1×n. The system can then be rewritten to be in the following form:{dot over (e)}=e+(T)+y_r,  (140)

Due to the structure of , which is derived from a finite difference scheme, and knowing that is Hurwitz, there exist M>0 and ω>0 such that:∥exp(t)∥<M exp(−ωt).  (141)Solving Eq. (139) yields the following:e(t)=exp(t)e(0)+∫₀^t exp((t−τ))(((T)+y_r)dτ,  (142)By taking the L₂-norm for the left and right sides of Eq. (142), the following can be determined:∥e(t)∥≤∥exp(t)e(0)∥+∫₀^t∥exp((t−τ))((T)+y_r)∥dτ.  (143)Letting 0≤∥e₀∥≤ϵ₀, then:∥e(t)∥≤M exp(−ωt)ϵ₀+M∫₀^t exp((τ−t)ω)∥((T)+y_r)∥dτ.  (144)With the boundness property of (T), it can be written:∥e(t)∥≤C₁ exp(−ωt)+C₂.  (145)where C₁ and C₂ depend on M, ω, ϵ₀ and the bound on (T), the convergence of the error to a small region, which can be made as small as possible by appropriate choice of the gains, is ensured.

Nonlinear Observer.

The nonlinear Lyapunov controller in Eq. (132) needs the availability of the temperature distribution along the membrane boundaries, which is not viable in practice. Nevertheless, the process model is observable and the present nonlinearity is continuous and Lipschitz continuous, therefore, designing an observer will help to estimate the necessary temperature states and feed them back to the controller.

Nonlinear Lyapunov-Based Observer Design.

A nonlinear observer that takes advantage of the well-modeled nonlinearity structure in the process can be used. It can be extended to include an online constructive algorithm for finding an adequate observer gain. In this context, consider the nonlinear Lyapunov-based observer in the following form:={circumflex over ({tilde over (T)})}+u+({circumflex over (T)})+L(y−{circumflex over ({tilde over (T)})}),  (146)where and {circumflex over (T)} are the estimates of {tilde over (T)} and T, and L is the observer gain matrix. Following this, the observer error dynamic equation can be formulated as:{dot over ({tilde over (T)})}(−L)T+(T)−({circumflex over (T)}),  (147)where T={acute over (T)}−{circumflex over ({tilde over (T)})} is the state estimation error. Proposition 2 gives sufficient conditions for exponential stability of the state observer error of Eq. (147).

• • Proposition 2. There exists an exponential stable observer of the form of Eq. (146), if there exists a symmetric positive definite matrix, P, such that he following standard algebraic Ricatti equation is satisfied:

$\begin{array}{cc}P+P{T}^{}+{\mathrm{PR}}_{1}P+Q_1=0,& \left(148\right)\\ \mathrm{where}{R}_1=\left({\lambda }^{2}I-\frac{C^{T}C}{ϵ}\right),Q_1=I\left(ϵ+1\right),ϵ>0,ϵ\in {+}^{}.& \left(149\right)\end{array}$

and the observer gain matrix is

$L=\frac{{\mathrm{PC}}^T}{2ϵ}.$

The performance of the nonlinear observer was tested with the convergence time and error. FIG. 31A shows the fast convergence of the estimated states using the nonlinear Lyapunov-based observer for different locations along the membrane length. Each curve represents the location on the membrane module. These locations are marked as a percentage of the total length of the membrane, starting from the top of the membrane to its bottom, therefore the curves in FIG. 31A illustrates the observer convergence on the top, the middle and the bottom of the whole membrane length. The fast convergence of the observer is guaranteed by using a smaller sampling time. In the same context, the error analysis shows the exponential convergence of the error to zero as depicted in FIG. 31B. Each curve represents the location on the membrane module, which is a percentage from the membrane total length.

Nonlinear Observer-Based Lyapunov Controller.

Following the design of the nonlinear Lyapunov-based observer, the temperature difference along the membrane boundaries is available at each instant of time. Proceeding further, the estimated states back can be fed to the nonlinear Lyapunov controller to operate, which is valid because the observer convergence is fast. The control law in Eq. (132) with the estimated temperature difference along the membrane boundaries in Eq. (146). The new control law will be:

$\begin{array}{cc}\{\begin{array}{c}{u}_1=\frac{1}{b_1}\left[-K_1\left({\hat{\tilde{T}}}_1-y_r\right)-1_{}\hat{\tilde{T}}-1_{}\left(\hat{T}\right)\right],\\ u_2=\frac{1}{b_2}\left[-K_2\left({\hat{\tilde{T}}}_n-y_r\right)-n_{}\hat{\tilde{T}}-n_{}\left(\hat{T}\right)\right].\end{array}& \left(150\right)\end{array}$ FIG. 32 is a schematic diagram illustrating an example of the structure of the nonlinear observer-based Lyapunov controller.

The convergence of the combined system can be examined by studying the stability of the augmented system of the dynamics of the tracking reference error e, and the observer error T. The augmented error system is shown in:

$\begin{array}{cc}\left(\begin{array}{c}\stackrel{.}{\stackrel{~}{T}}\\ \stackrel{.}{e}\end{array}\right)=\left(\begin{array}{cc}-L& 0\\ +K& \end{array}\right)\left(\begin{array}{c}\tilde{T}\\ e\end{array}\right)+\left(\begin{array}{c}\left(T\right)-\left(\hat{T}\right)\\ \Delta \left(T\right)+\left(\tilde{T}\right)\end{array}\right)+\left(\begin{array}{c}0\\ y_r\end{array}\right),& \left(151\right)\end{array}$ and can be further simplified to be:{dot over (ϑ)}=ϑ+Φ({tilde over (T)},{circumflex over ({tilde over (T)})})+ç,  (152)where

$\vartheta =\left(\begin{array}{c}\stackrel{\_}{T}\\ e\end{array}\right)\in {n\times 1}^{}$ represents the error state vector. The dynamics of the error are represented in

$=\left(\begin{array}{cc}-\mathrm{LC}& 0\\ +K& \end{array}\right)\in {n\times n}^{},$ where has a structure as follows:

$\begin{array}{cc}=\left[\begin{array}{c}{1}_{}\\ 0\\ ⋮\\ 0\\ n_{}\end{array}\right].& \left(153\right)\end{array}$ The nonlinear terms are in:

$\Phi \left(\tilde{T},\hat{\tilde{T}}\right)=\left(\begin{array}{c}\left(T\right)-\left(\hat{T}\right)\\ \Delta \left(T\right)+\left(\hat{T}\right)\end{array}\right)\in {n\times 1}^{},$ where Δ(T)∈^n×1 has a structure such as:Δ(T)=[(₁({circumflex over (T)})−₁(T))0 . . . 0(_n({circumflex over (T)})−_n(T))]^tr.In addition,

$\varsigma =\left(\begin{array}{c}0\\ \end{array}\right)\in {n\times 1}^{},$ and matrix

$K=\left[\begin{array}{c}{\tilde{K}}_1\\ 0\\ ⋮\\ 0\\ {\tilde{K}}_2\end{array}\right]\in {n\times n}^{}.$ Proposition 3 gives sufficient condition for stability of the augmented error matrix.

• • Proposition 3. The augmented error system in Eq. (152) is stable, if there exists two reals, M₁>0 and ω₁>0, such that:∥exp(t)∥<M₁ exp(−ω₁t),  (154)and if0<∥ϑ(0)∥≤ϑ₀,  (155)where ϑ>0 is a smaller scalar.

Proof.

Similar to the proof of Proposition 1, the stability of system in Eq. (152) is guaranteed if the matrix is Hurwitz. Solving Eq. (152) yields the following:ϑ(t)=exp(t)ϑ(0)+∫₀^t exp((t−τ))Φ)({tilde over (T)},{circumflex over ({tilde over (T)})})+çdτ.  (156)then the following can be determined taking the L₂-norm for the left and right sides:∥ϑ(t)∥≤∥exp(t)ϑ(0)∥+∫₀^t∥exp((t−τ))Φ({tilde over (T)}−{circumflex over ({tilde over (T)})})+ç∥dτ.  (157)If there exist M₁>0 and ω₁>0, then:∥ϑ(t)∥≤M₁ exp(−ω₁t)ϑ₀+M₁∫₀^t exp((τ−t)ω₁)∥Φ({tilde over (T)},{circumflex over ({tilde over (T)})})+ç∥dτ.   (158)With the boundness and the Lipschitz continuity of the nonlinearity, there exist constants such that:∥ϑ(t)∥≤C₃ exp(−ω₁t)+C₄,  (159)where C₁ and C₂ depend on M₁, ω₁, ϑ₀ and the bound on the nonlinearity, which ensure the convergence to a small region, that can be made as small as possible by appropriate choice of the gains.

Simulation and Results.

Extensive numerical simulations have been carried out on the dynamic model of the DCMD process. The simulations illustrated the behavior of the process, especially along the membrane boundaries to highlight the effect of the controller. All simulations were performed with real membrane parameter values, which are provided by the Water Desalination and Reuse Center (WDRC) at King Abdullah University of Science and Technology (KAUST). These values are listed in the table of FIG. 33A.

Model Simulation.

The production of the DCMD process with no imposed control is fluctuating, and therefore, the temperature difference along the membrane boundaries is not constant, and does not fall within practical values to provide fixed rate of permeate flux. The induced fluctuations in the temperature difference along the membrane boundaries lead to more energy consumption, and as a natural result, fluctuations in the water production rate.

Observer-Based Lyapunov Controller Simulations.

The overall structure of the whole DCMD process was simulated. The structure included the nonlinear observer-based controller, as well as the DCMD dynamic model. The role of the nonlinear controller was maintaining the difference temperature along the membrane boundaries to be almost fixed around a desired reference such as about 10° C. This difference was recommended from experimental perspective for high permeate flux production. The value of this reference was chosen to be within a range of 5° C. to 15° C., following practical advice for sufficient water production. Additionally, the controller gains K₁ and K₂ were set to be equal K₁=K₂=K. FIG. 33B shows the convergence of the temperature difference for different locations on the membrane boundaries to the desired reference. Each curve represents the location on the membrane module, which is a percentage of the membrane total length.

The convergence of the temperature difference was accomplished with an acceptable evolution of the control inputs within the predefined range. FIG. 33C depicts that time evolution of the control inputs. Reaching steady state phase, with no disturbance, the DCMD process experiences a stable behavior among all its designated domain. FIG. 33D illustrates that effect of the observer-based Lyapunov_controller on the process during the steady state phase. The time evolution of the error between the desired reference and the actual temperature difference along the membrane boundaries is shown in FIG. 33E. The results show the asymptotic stability of the error, which was achieved using the observer-based Lyapunov_controller.

Controller Gain Simulations.

The controller gains K₁ and K₂ have a significant role in smoothing out the performance of the DCMD process, where it reduces the potential oscillations during the transient phase. Therefore, the selection process of this constant is important, and needs a balance between the convergence time and the process behavior. The impact of the controller gain is shown in FIG. 33F, where it shows the effect on the control inputs time evolution. Consequently, low controller gain gives low oscillations in the control inputs transient phase, unlike high gain values, and reduces spikes that may appear in the response.

An observer-based Lyapunov control design was presented to manage the production of a DCMD process. The DCMD process was modeled with Advection Diffusion Equation model in two dimensions, where the model experienced a nonlinear behavior that originated from the associated boundary conditions. Due to the significance of the temperature difference along the membrane boundaries in driving the DCMD process, the dynamic model was reduced to include only the temperature distribution along the membrane boundaries. Afterwards, a nonlinear Lyaponuv control design was adopted to maintain the temperature difference around a sufficient desired reference. This allows the process be able to provide satisfying water production levels.

Following the needs of the controller, a nonlinear observer was designed to provide the nonlinear controller with the temperature distribution along the membrane surfaces at each time instance. The controller combined with the observer guarantees the convergence to the desired reference. Control inputs were constrained with maximum and minimum values to follow critical physical limitations as well as process requirements. Simulations were performed with real membrane parameter values.

With reference now to FIG. 34, shown is a schematic block diagram of an example of processing circuitry 1400 that may be used to implement various portions of the control of the membrane distillation (MD) process 100 of FIGS. 1 and 14 in accordance with various embodiments of the present disclosure. The processing circuitry 1400 includes at least one processor circuit, for example, having a processor 1403 and a memory 1406, both of which are coupled to a local interface 1409. To this end, the processing circuitry 1400 may be implemented using one or more circuits, one or more microprocessors, microcontrollers, application specific integrated circuits, dedicated hardware, digital signal processors, microcomputers, central processing units, field programmable gate arrays, programmable logic devices, state machines, or any combination thereof. The local interface 1409 may comprise, for example, a data bus with an accompanying address/control bus or other bus structure as can be appreciated. The processing circuitry 1400 can include a display for rendering of generated graphics such as, e.g., a user interface and an input interface such, e.g., a keypad or touch screen to allow for user input. In addition, the processing circuitry 1400 can include communication interfaces (not shown) that allow the processing circuitry 1400 to communicatively couple with other communication devices. The communication interfaces may include one or more wireless connection(s) such as, e.g., Bluetooth or other radio frequency (RF) connection and/or one or more wired connection(s).

Stored in the memory 1406 are both data and several components that are executable by the processor 1403. In particular, stored in the memory 1406 and executable by the processor 1403 are soft sensing application(s) 1415, an operating system 1418, and/or other applications 1421. Soft sensing applications 1415 can include applications that support, e.g., controllers for control of the operation of the MD process 100 and/or observers for estimation of states and/or characteristics of the MD process 100. It is understood that there may be other applications that are stored in the memory 1406 and are executable by the processor 1403 as can be appreciated. Where any component discussed herein is implemented in the form of software, any one of a number of programming languages may be employed such as, for example, C, C++, C #, Objective C, Java®, JavaScript®, Perl, PHP, Visual Basic®, Python®, Ruby, Delphi®, Flash®, LabVIEW® or other programming languages.

A number of software components are stored in the memory 1406 and are executable by the processor 1403. In this respect, the term “executable” means a program file that is in a form that can ultimately be run by the processor 1403. Examples of executable programs may be, for example, a compiled program that can be translated into machine code in a format that can be loaded into a random access portion of the memory 1406 and run by the processor 1403, source code that may be expressed in proper format such as object code that is capable of being loaded into a random access portion of the memory 1406 and executed by the processor 1403, or source code that may be interpreted by another executable program to generate instructions in a random access portion of the memory 1406 to be executed by the processor 1403, etc. An executable program may be stored in any portion or component of the memory 1406 including, for example, random access memory (RAM), read-only memory (ROM), hard drive, solid-state drive, USB flash drive, memory card, optical disc such as compact disc (CD) or digital versatile disc (DVD), floppy disk, magnetic tape, or other memory components.

The memory 1406 is defined herein as including both volatile and nonvolatile memory and data storage components. Volatile components are those that do not retain data values upon loss of power. Nonvolatile components are those that retain data upon a loss of power. Thus, the memory 1406 may comprise, for example, random access memory (RAM), read-only memory (ROM), hard disk drives, solid-state drives, USB flash drives, memory cards accessed via a memory card reader, floppy disks accessed via an associated floppy disk drive, optical discs accessed via an optical disc drive, magnetic tapes accessed via an appropriate tape drive, and/or other memory components, or a combination of any two or more of these memory components. In addition, the RAM may comprise, for example, static random access memory (SRAM), dynamic random access memory (DRAM), or magnetic random access memory (MRAM) and other such devices. The ROM may comprise, for example, a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or other like memory device.

Also, the processor 1403 may represent multiple processors 1403 and the memory 1406 may represent multiple memories 1406 that operate in parallel processing circuits, respectively. In such a case, the local interface 1409 may be an appropriate network that facilitates communication between any two of the multiple processors 1403, between any processor 1403 and any of the memories 1406, or between any two of the memories 1406, etc. The local interface 1409 may comprise additional systems designed to coordinate this communication, including, for example, performing load balancing. The processor 1403 may be of electrical or of some other available construction.

Although the soft sensing application(s) 1415, the operating system 1418, application(s) 1421, and other various systems described herein may be embodied in software or code executed by general purpose hardware as discussed above, as an alternative the same may also be embodied in dedicated hardware or a combination of software/general purpose hardware and dedicated hardware. If embodied in dedicated hardware, each can be implemented as a circuit or state machine that employs any one of or a combination of a number of technologies. These technologies may include, but are not limited to, discrete logic circuits having logic gates for implementing various logic functions upon an application of one or more data signals, application specific integrated circuits having appropriate logic gates, or other components, etc. Such technologies are generally well known by those skilled in the art and, consequently, are not described in detail herein.

Also, any logic or application described herein, including the soft sensing application(s) 1415 and/or application(s) 1421, that comprises software or code can be embodied in any non-transitory computer-readable medium for use by or in connection with an instruction execution system such as, for example, a processor 1403 in a computer system or other system. In this sense, the logic may comprise, for example, statements including instructions and declarations that can be fetched from the computer-readable medium and executed by the instruction execution system. In the context of the present disclosure, a “computer-readable medium” can be any medium that can contain, store, or maintain the logic or application described herein for use by or in connection with the instruction execution system. The computer-readable medium can comprise any one of many physical media such as, for example, magnetic, optical, or semiconductor media. More specific examples of a suitable computer-readable medium would include, but are not limited to, magnetic tapes, magnetic floppy diskettes, magnetic hard drives, memory cards, solid-state drives, USB flash drives, or optical discs. Also, the computer-readable medium may be a random access memory (RAM) including, for example, static random access memory (SRAM) and dynamic random access memory (DRAM), or magnetic random access memory (MRAM). In addition, the computer-readable medium may be a read-only memory (ROM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or other type of memory device.

It should be emphasized that the above-described embodiments of the present disclosure are merely possible examples of implementations set forth for a clear understanding of the principles of the disclosure. Many variations and modifications may be made to the above-described embodiment(s) without departing substantially from the spirit and principles of the disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

It should be noted that ratios, concentrations, amounts, and other numerical data may be expressed herein in a range format. It is to be understood that such a range format is used for convenience and brevity, and thus, should be interpreted in a flexible manner to include not only the numerical values explicitly recited as the limits of the range, but also to include all the individual numerical values or sub-ranges encompassed within that range as if each numerical value and sub-range is explicitly recited. To illustrate, a concentration range of “about 0.1% to about 5%” should be interpreted to include not only the explicitly recited concentration of about 0.1 wt % to about 5 wt %, but also include individual concentrations (e.g., 1%, 2%, 3%, and 4%) and the sub-ranges (e.g., 0.5%, 1.1%, 2.2%, 3.3%, and 4.4%) within the indicated range. The term “about” can include traditional rounding according to significant figures of numerical values. In addition, the phrase “about ‘x’ to ‘y’” includes “about to about ‘y’”.

Claims

1. A system for water desalination, comprising: a membrane distillation (MD) module including a feed side and a permeate side separated by a membrane having a membrane boundary layer, wherein the feed side further comprises a feed inlet and feed outlet and wherein the permeate side further comprises a permeate inlet and permeate outlet:temperature sensors configured to monitor feed solution temperatures at the feed outlet and permeate solution temperatures at the permeate outlet; andprocessing circuitry configured to estimate, during steady- and unsteady-state operation, feed solution temperatures along the membrane boundary layer, on the feed side, and permeate solution temperatures along the membrane boundary layer, on the permeate side of the MD module using an Advection-Diffusion Equation (ADE), wherein the processing circuitry utilizes, as inputs, the monitored feed solution and permeate solution temperatures at the feed outlet and permeate outlet, respectively,wherein the ADE includes a first part that describes the unsteady-state operation, a second part that describes a heat convection transfer between the feed side and the permeate side, and a third part that describes a heat conduction transfer between the feed side and the permeate side, andwherein the processing circuitry is further configured to control operation of the MD module comprising maintaining, at a desired reference value, a difference between the feed solution temperatures along the membrane boundary layer, on the feed side, and the permeate solution temperatures along the membrane boundary layer, on the permeate side, based on the estimation.

2. The system of claim 1, wherein the processing circuitry is further configured to estimate a membrane mass transfer coefficient of the MD module.

3. The system of claim 1, wherein the feed solution temperatures and permeate solution temperatures include boundary temperatures along the membrane boundary layer.

4. The system of claim 1, wherein the Advection-Diffusion Equation (ADE) is a two-dimensional ADE model including dimensions of the membrane.

5. The system of claim 4, wherein the ADE model includes a plurality of grid points extending across the feed side and the permeate side of the MD module.

6. The system of claim 5, wherein the feed solution temperatures and permeate solution temperatures include temperatures corresponding to at least a portion of the plurality of grid points.

7. The system of claim 1, wherein the MD module is a counter current direct-contact MD module.

8. A method for water desalination, the method comprising: monitoring feed solution temperatures at a feed outlet and permeate solution temperatures at a permeate outlet of a membrane distillation (MD) module;determining, during steady- and unsteady-state operation, a plurality of estimated temperature states along a membrane having a membrane boundary layer separating the feed side and the permeate side of the MD module using an Advection-Diffusion Equation (ADE),wherein the ADE utilizes, as inputs, the feed solution temperature monitored at the feed outlet and the permeate solution temperature monitored at the permeate outlet, wherein the ADE includes a first part that describes the unsteady-state, a second part that describes a heat convection transfer between the feed side and the permeate side, and a third part that describes a heat conduction transfer between the feed side and the permeate side; andcontrolling operation of the MD module comprising maintaining at a desired reference value, a difference between the feed solution temperatures along the membrane boundary layer, on the feed side, and the permeate solution temperatures along the membrane boundary layer, on the permeate side, based on the determining.

9. The method of claim 8, further comprising estimating a membrane mass transfer coefficient of the MD module.

10. The method of claim 8, further comprising adjusting characteristics of the feed solution or permeate solution based at least in part upon the feed solution temperature monitored at the feed outlet and the permeate solution temperature monitored at the permeate outlet.

11. The method of claim 8, wherein the estimated temperature states include temperatures corresponding to grid points extending across the feed side and the permeate side of the MD module.

12. The method of claim 8, wherein the ADE model is a two-dimensional ADE model including dimensions of the membrane.

13. A system for water desalination, comprising: a membrane distillation (MD) module including a feed side and a permeate side separated by a membrane having a membrane boundary layer;temperature sensors configured to monitor feed solution temperatures at a feed outlet and permeate solution temperatures at a permeate outlet; andprocessing circuitry configured to (1) apply an Advection-Diffusion Equation (ADE) model to estimate, during steady- and unsteady-state operation, feed solution temperatures and permeate solution temperatures along the membrane boundary layer, within the feed side and the permeate side, respectively, of the MD module using an Advection-Diffusion equation (ADE), wherein the feed solution temperatures monitored at the feed outlet and the permeate solution temperatures monitored at the permeate outlet are used as inputs to the ADE model, and (2) control operation of the MD module based at least in part upon the feed solution temperatures estimated within the feed side and permeate solution temperatures estimated within the permeate side,wherein the ADE model includes a first part that describes the unsteady-state, a second part that describes a heat convection transfer between the feed side and the permeate side, and a third part that describes a heat conduction transfer between the feed side and the permeate side, andwherein the processing circuitry is further configured to control operation of the MD module comprising maintaining, at a desired reference value, a difference between the feed solution temperatures along the membrane boundary layer, on the feed side, and the permeate solution temperatures along the membrane boundary layer, on the permeate side, based on the estimation.

14. The system of claim 13, wherein the processing circuitry is further configured to estimate a membrane mass transfer coefficient of the MD module.

15. The system of claim 13, wherein the feed solution temperatures and permeate solution temperatures include boundary temperatures along the membrane boundary layer.

16. The system of claim 13, wherein the estimated temperatures include temperatures corresponding to grid points extending across the feed side and the permeate side of the MD module.