A Numerically Validated Approach to Modeling Water Hammer Phenomena Using Partial Differential Equations and Switched Differential-Algebraic Equations

Hydraulic transients are an inherent occurrence in water distribution networks, arising from intentional or accidental alterations in the network configuration, which result in immediate structural changes. These sudden modifications can have profound consequences on flow regimes, encompassing issues such as pump malfunctions or even severe pipeline failures. Conventionally, the flow of water in pipes has been mathematically represented by a system of nonlinear hyperbolic balance laws, commonly referred to as partial differential equations (PDEs) [1], wherein sudden structural changes lead to large peaks and fast transients in the solution.

Switched differential algebraic equations (DAEs) have been proposed as a framework to model these rapid transients. The concept of switched DAEs was initially developed for modeling electrical circuits [2], enabling an accurate mathematical representation of peaks and fast transients using Dirac impulses and jumps. This study focuses on the phenomenon known as water hammer, occurring due to abrupt velocity changes in pipelines, leading to significant increases in pressure. Water hammer is typically triggered by the rapid closure of valves, pump shutdowns, or pump restarts. It has been suggested that pressure spikes observed in PDE simulations can be effectively approximated using an appropriate switched DAE model.

A DAE takes the following form:

where, $\left(E_{\mathbf{p}}, A_{\mathbf{p}}\right) \in \mathbb{R}^{n \times n}, \mathbf{p} \in\{1, \cdots, \mathbf{P}\}, \mathbf{P} \in \mathbb{N}$ and $\sigma: \mathbb{R} \rightarrow\{1, \cdots, \mathbf{P}\}$ is considered a piecewise constant, termed the switching signal. It is assumed to be right continuous and to have locally finitely many jumps.

For the existence and uniqueness of solutions in both switched and non-switched DAEs, reliance is placed on a notion of regularity for matrix pairs. This regularity condition ensures that the matrices involved in the DAE satisfy certain properties and exhibit well-behaved characteristics.

A commonly used regularity notion for matrix pairs is based on the concept of regularity indices. These indices provide information about the differential and algebraic aspects of the system, crucial for determining the solvability of the DAE. The differential index represents the highest order of differentiation appearing in the differential equations, while the algebraic index represents the highest order of algebraic dependence in the system.

A matrix pair $\left(E_{\mathbf{p}}, A_{\mathbf{p}}\right)$ is considered regular with index $(v, \mu)$ if the following conditions hold:

1. The matrix pair $\left(E_{\mathbf{p}}, A_{\mathbf{p}}\right)$ has a constant rank $V$ for all $p \in 1, \ldots, \mathbf{P}$. The matrix $\mathbf{A p}$ satisfies the rank condition rank $\left(\mathbf{A} \mathbf{p}-\lambda \mathbf{E}_{\mathbf{p}}\right) \leq \mu$ for all $p \in 1, \ldots, \mathbf{P}$ and all $\lambda \in \mathbb{C}$ [3].

2. The matrix $\mathbf{E p}$ satisfies the rank condition rank$\operatorname{rank}(\mathbf{E p})=\operatorname{rank}\left(\mathbf{E p}\left(\mathbf{A p}-\lambda \mathbf{E}_{\mathbf{p}}\right)\right)=v$ for all $p \in 1, \ldots, \mathbf{P}$ and all $\lambda \in \mathbb{C}$.

These regularity indices, $(v, \mu)$, provide important insights into the solvability and stability of the DAE. By analyzing these indices, conditions under which a unique solution exists can be determined, and effective numerical methods for solving the DAE can be developed.

Proposition 1 (QWF, c.f. [4]). A matrix pair $(E, A) \in \mathbb{R}^{n \times n} \times \mathbb{R}^{n \times n}$ is regular if and only if there exist invertible transformation matrices $S, T \in \mathbb{R}^{n \times n}$, which put $(E, A)$ into QWF

$(S E T, S A T)=\left(\left[\begin{array}{cc} I & 0 \\ 0 & N \end{array}\right],\left[\begin{array}{ll} J & 0 \\ 0 & I \end{array}\right]\right)$

where, $N \in \mathbb{R}^{n_2 \times n_2}$ with $0 \leq n_2 \leq n$ is a nilpotent matrix, $J \in \mathbb{R}^{n_1 \times n_1}$ with $n_1:=n-n_2$ is some matrix and $I$ is the identity matrix of the appropriate size.

The paper structure is as follows. Section 1 establishes the definition of the water network and its constituent elements as a graph, introducing the mathematical models for pipes, reservoirs, and valves. Section 2 provides a concise overview of the mathematical preliminaries concerning conservation laws and switched differential-algebraic equations (DAEs). In Section 3, a comprehensive analysis of a simple water network that demonstrates the occurrence of water hammer is conducted, deriving both the corresponding partial differential equation (PDE) model and the switched DAE model for this network. Finally, Section 4 presents a numerical comparison between the PDE and switched DAE models.

A finite undirected graph, denoted as $(V, E)$, provides an efficient means of representing the structure of water distribution networks. This graph-based approach effectively captures the organization of the network, where the set $V$ embodies the nodes and $E \subseteq V \times V$ corresponds to the edges that interconnect these nodes. In water distribution networks, each edge $e \in E$ signifies a distinct pipe, facilitating the transmission of water between nodes. Conversely, the nodes $v \in V$ comprise various network components, including junctions, pumps, valves, and reservoirs, each functioning as connection points or endpoints for the pipes. By employing this graph-based representation, the network's topology and connectivity can be accurately portrayed. The comprehensive framework formed by the nodes and edges enables modeling and examination of the interactions between diverse network components. This, in turn, permits the analysis and simulation of network dynamics, such as flow patterns and pressure fluctuations, while accounting for the interconnected pipes and associated components. Utilizing a graph-based representation yields valuable insights into the behavior of water distribution networks, facilitating a deeper understanding of the relationships between different components. This approach serves as a potent tool for analyzing network performance, optimizing its operation, and addressing challenges like water hammer and pressure variations [5].

In the following sections, a detailed overview of the modeling approach for the various components present in a water network is presented. The primary aim is to accurately capture the behavior and dynamics of these elements to enhance the understanding of the system as a whole. The models predominantly focus on pressure and flow variables, which are interconnected at the pipe endpoints. The specific type of node being represented determines the nature of this interconnection, as pressure and flow interact accordingly. Thorough analysis and modeling of these variables foster valuable insights into the functioning and performance of the water network.

To capture the intricate interactions within the water distribution system effectively, it is essential to consider the distinct characteristics of each network element. The nodes within the network encompass diverse entities such as junctions, pumps, valves, and reservoirs, collectively forming the set of nodes, denoted as $V$. Conversely, the pipes in the network are represented by the set of edges, denoted as $E$, forming a finite undirected graph $(V , E)$.

Each edge $e$ in $E$ corresponds to an individual pipe within the water network, while the nodes $v$ in $V$ function as connection points or endpoints for these pipes. The nodes play a pivotal role in regulating and controlling the flow of water throughout the network [6].

In the modeling approach, particular attention is devoted to the physical quantities of pressure and flow. The values of pressure and flow at the pipe endpoints are interconnected, with the specific relationship depending on the nature of the node being considered. Furthermore, the modeling of flow within the pipes accounts for water density. Although water is typically assumed to be an incompressible fluid with a constant density, the focus is on accurately capturing the water hammer effect. To accomplish this, the slight compressibility of water is acknowledged and incorporated into the modeling framework to capture the dynamic behavior of the system [7].

By integrating these essential elements and variables into the models, a robust framework for simulating and analyzing water distribution networks is provided. This approach enables the examination of the impact of various factors, such as abrupt velocity changes, valve closures, pump shutdowns, and restarts, on the occurrence of pressure peaks and water hammer phenomena. Through these simulations, the efficacy of the modeling approach in approximating and understanding the system's behavior under diverse operational conditions is demonstrated.

Water flow in pipes can be characterized by two distinct models, contingent on whether the compressibility of water is taken into account. These models cater to different purposes and offer varying levels of accuracy, making them suitable for diverse analysis scenarios. A brief overview of both models is provided in this section.

1. Compressible Flow Model: In scenarios where the compressibility of water is considered, the density and mass flow rate of water are treated as variables that vary spatially. This implies that as water flows through the pipe, its density can change due to pressure variations. Compressible flow models involve the utilization of partial differential equations, such as the continuity equation, Navier-Stokes equation (momentum equation), and energy equation, to describe the conservation of mass, momentum, and energy. Numerical methods, such as computational fluid dynamics (CFD), are employed to solve these equations and analyze water behavior under different conditions, including transient phenomena like water hammer [8].

2. Incompressible Flow Model: In practical engineering applications, it is often deemed acceptable to model water as an incompressible fluid when studying the qualitative behavior of flow, particularly in large networks like water distribution systems. In this model, water is assumed to have constant density regardless of pressure changes. The governing equations for incompressible flow are simplified by treating the density as a constant, eliminating it from the equations. The continuity equation, which ensures the conservation of mass flow rate along the flow direction, is the primary equation used in this model. These simplified equations are more straightforward to solve and provide a reasonable approximation for steady-state and low-speed flows [9].

It is noteworthy that although the incompressible flow model neglects density variation, it still captures many essential flow characteristics, including flow rates, pressure drops, and hydraulic behavior, in most practical scenarios. However, for accurate analysis of transient phenomena like water hammer or high-speed flows, the compressible flow model, which accounts for water's compressibility, becomes necessary. Both models possess their own advantages and limitations, and the choice between them hinges on the specific application and the desired level of accuracy. The selection of the appropriate model is determined by the nature of the problem and the level of detail required for analysis.

Following [2], the pressure law for compressible fluids is employed:

where, $K^\omega>0$ denotes the bulk modulus, $P_a>0$ represents the atmospheric pressure, and $\rho_a^\omega>0$ is the density. Furthermore, the bulk modulus is connected to the speed of sound $c>0$, as follows:

It is noted that $\beta^\omega:=1 / K$ is the compressibility coefficient. In this study, a pipe uniformly filled with water is examined, possessing a positive mass density function denoted as $(x,t)$ and a mass flux function denoted as $q^\omega(x, t)$, both defined on the interval $[0, L] \times+$. The objective is to model the compressible flow of water in the pipe using a balance law of the following nature:

with the pressure law $P: \mathbb{R}_{+} \rightarrow \mathbb{R}_{+}$ given by (2) and where $c_f>0$ is the friction against the pipe wall and $D>0$ is the diameter of the pipe. The initial condition for (4) is:

It is assumed that an initial flow function $q_0(x)$ and an initial pressure function $p_0(x)$ exist, defined on the interval $[0, L]$. The initial condition for water flow in a pipe is given implicitly in terms of pressure, reflecting its significance as a fundamental physical quantity in water networks. When pipes are connected to other elements in the network, additional boundary conditions are imposed to capture the behavior at the network boundaries. These boundary conditions depend on the specific configuration and characteristics of the interconnected elements [10], ensuring accurate representation of system dynamics. Including these boundary conditions is essential for modeling the flow dynamics and analyzing the behavior of water distribution networks. The inclusion of other network components introduces further constraints and equations that need to be satisfied for an accurate representation of the overall system dynamics. By expressing the initial condition in terms of the pressure, a more consistent and comprehensive representation of the interconnected pipes and their interactions with other network elements is ensured. This approach allows for capturing the complex dynamics and interactions within the water network, taking into account the pressure variations and their effects on the overall system behavior. Thus, the choice of an implicit initial condition in terms of the pressure provides a more suitable and physically relevant representation when considering the connections between individual pipes and other components in a water distribution network.

For an initial flow function $q_0:[0, L]$ and an initial pressure function $p_L:[0, L] \rightarrow \mathbb{R}^{+}$, the problem with an implicit initial condition given in terms of pressure rather than density is considered. This choice is motivated by the significance of pressure as a fundamental physical quantity, especially when the pipes are interconnected with other elements of the water distribution network. In such cases, additional boundary conditions are introduced as the individual pipes interface with the broader network [5]:

Within the water network graph, a reservoir functions as a node, contributing to the overall system dynamics. Characterized by unrestricted mass flow and a predetermined pressure value, the reservoir plays an essential role. When a pipe is connected to the left side of a reservoir, a boundary condition affecting the pressure behavior at the pipe's starting point is introduced, as demonstrated in the PDE model (4). This connection also imposes constraints on the variable $\mathrm{P}_0$ within the ODE model (6). These constraints govern the initial conditions and the reservoir's subsequent evolution, impacting the flow characteristics and pressure distribution throughout the water network. Understanding the interaction between a reservoir and its connected pipe is therefore crucial for comprehending the system's behavior.

To maintain simplicity, a modeling approach incorporating a valve with a reservoir, as depicted in Figure 1, is adopted. This integration facilitates accounting for the valve's behavior within the system while preserving a straightforward framework. The valve-reservoir combination imposes specific algebraic constraints that regulate their interaction. These constraints establish relationships between various variables and introduce additional conditions that affect the system's dynamics. By incorporating the valve and reservoir as a unified model, their interconnected effects can be effectively examined, and their implications on the overall behavior of the water network can be assessed.

case-I: open: $P_V^\omega=P_D$,

case-II: closed: $Q_V^\omega=0$,

where, $P_V^\omega$, $Q_V^\omega$ denotes the flow pressure through the valve and $P_D$ is the pressure of the reservoir located.

A simple water distribution system is investigated, consisting of an upstream reservoir with a fixed pressure $P_U$, a single conduit with a fixed length $L$, and a valve connected to a downstream reservoir with pressure $P$. $P_D$<$P_U$, as shown in Figure 2.

The described scenario, involving the abrupt closure of a valve, represents a typical situation in which water hammer can occur in a system. The objective is to evaluate the effectiveness of the PDE model and the switched DAE model in capturing the water hammer effect within this specific network architecture. To achieve this, solution concepts employed by both approaches to address this water network problem are discussed first. By examining and comparing various proposed solution methods, the capability of each model to accurately represent and explain the water hammer phenomenon can be assessed.

The balance law in the pipe is considered in conjunction with the pressure law (2).

$\Omega=[0, L], \quad \mathscr{U} \subseteq \mathbb{R}_{+} \times \mathbb{R}, \quad U(x, t)=\left(\begin{array}{c}\rho^\omega(x, t) \\ q(x, t)\end{array}\right)$,

and for $u=\left(\rho^\omega, q\right) \in \mathscr{U}$

$F(u)=\left(\begin{array}{c}q^\omega \\ \frac{\left(q^\omega\right)^\omega}{\rho^\omega}+P\left(\rho^\omega\right),\end{array}\right), \quad S(u)=\left(\begin{array}{c}0 \\ -c \frac{q^\omega\left|q^\omega\right|}{2 D \rho^\omega}\end{array}\right)$,

The initial condition is given by (5), where the initial density $\rho_0^\omega$ is induced by the typically considered initial pressure $p_0$ via the invertible pressure law. The boundary conditions in terms of the Figure 2 is given by

$\Psi_1(t,(\rho, q))=P(\rho)-P_U, \quad \forall t \in \mathbb{R}^{+}$,

It is crucial to emphasize that the discontinuity resulting from the switch affects only the boundary condition. Consequently, the well-posedness analysis can be performed independently for each mode. This is feasible because the partial differential equation (PDE) can effectively be “restarted” at $t=t_S$, using the final value of the solution on the interval $\left[0, t_S\right]$ as the initial condition. As a result, the conditions for well-posedness, such as strict hyperbolicity, genuine nonlinearity, and S-Lipschitz continuity [6], can be examined independently of the valve's state. This approach allows for assessing the well-posedness of each mode individually and investigating their behavior as self-contained entities, ultimately providing a more comprehensive understanding of the system's dynamics. With initial and boundary conditions imposed as constants, it is evident that Eq. (7) holds. The Jacobian of the flux function in the equation can be expressed as follows [11]:

$A(\rho, q)=\left(\begin{array}{cc}0 & 1 \\ -\frac{q^2}{\rho^2}+P^{\prime}(\rho) & 2 \frac{q}{\rho}\end{array}\right)$

and invoking the pressure law (2) we see that $P^{\prime}(\rho)=\frac{K}{\rho_a}>0$ independently of $\rho$. Hence the eigenvalues of $A(\rho, q)$ are

$\lambda_{1 / 2}(\rho, q)=\frac{q}{\rho} \pm \sqrt{\frac{K}{\rho_a}}$ .

Hence, $\lambda_1(u)>\lambda_2(u)$ for all $u=(\rho, q) \in \mathscr{U}$ and hence (4) with pressure law (2) is strictly hyper- bolic. The eigenvectors of $A(u)$ are

$r_{-}(u)=\left(\begin{array}{c}-1 \\ -\lambda_1(u)\end{array}\right), \quad r_{+}(u)=\left(\begin{array}{c}1 \\ \lambda_2(u)\end{array}\right)$.

Hence

$\begin{aligned} & \mathbf{D} \lambda_1(u) r_{-}(u)=\left[-\frac{q}{\rho^2}, \frac{1}{\rho}\right]\left(\begin{array}{c}-1 \\ -\frac{q}{\rho}-\sqrt{\frac{K}{\rho_a}}\end{array}\right) \\ & =-\frac{1}{\rho} \sqrt{\frac{K}{\rho_a}}<0 \quad \forall u=(\rho, q) \in \mathscr{U}\end{aligned}$

and analogously

$\mathbf{D} \lambda_2(u) r_{+}(u)=\frac{1}{\rho} \sqrt{\frac{K}{\rho_a}}>0 \quad \forall u=(\rho, q) \in \mathscr{U}$

Consequently, genuine nonlinearity is established. Clearly S is locally Lipschitz continuous. The wellposedness of the boundaries follows from [6] since for the open valve situation,

$D_U \Psi(u) r_{+}(u)=\left(\begin{array}{ll}P^{\prime}(\rho) & 0\end{array}\right)\left(\begin{array}{c}1 \\ \lambda_2(u)\end{array}\right)=P^{\prime}(\rho) \neq 0$

and for the closed valve

$D_U \Psi(U) r_{+}(u)=\left(\begin{array}{ll}0 & 1\end{array}\right)\left(\begin{array}{c}1 \\ \lambda_2(u)\end{array}\right)=\lambda_2(u) \neq 0$

it holds, as we are in a subsonic case i.e.,

$\frac{q}{\rho}<\sqrt{\frac{K}{\rho_a}}$.

Therefore, it is reasonable to expect that numerical simulations will yield accurate approximations of the solutions to the corresponding initial/boundary value problem for both the compressible and incompressible fluid models. Numerical methods, such as finite difference, finite element, or finite volume techniques, can be employed to discretize the governing equations and solve them numerically. In numerical simulations, the continuous domain of the pipe network is typically discretized into a finite grid or mesh. The governing equations are then approximated using discretization schemes, which transform the partial differential equations into a system of algebraic equations that can be solved numerically. The accuracy of the numerical approximation depends on factors such as the grid resolution, the choice of numerical scheme, and the time integration method [11].

It is important to emphasize that numerical simulations provide valuable insights and predictions, but they are approximations of the actual system behavior. The accuracy of the results relies on the accuracy of the mathematical model, the quality of the numerical discretization, and the appropriate selection of model parameters and boundary conditions. Validation of the numerical simulations is essential to ensure the reliability of the results. This can be achieved by comparing the numerical solutions with experimental data or analytical solutions, whenever available. Sensitivity analyses and convergence studies can also be conducted to assess the robustness and accuracy of the numerical methods.

In conclusion, numerical simulations play a crucial role in the study of water flow in pipe networks, enabling the investigation of complex system behavior, analysis of different scenarios, and optimization of network design and operation. Through numerical simulations, engineers and researchers can gain valuable insights into the dynamics of water distribution systems and make informed decisions for efficient and reliable water supply.

The quasi-stationary model (6), in conjunction with the valve-dependent boundary conditions (valve open on $\left[0, t_S\right]$ and valve closed on $\left[t_s, \infty\right)$), for a setup resembling the one illustrated in Figure 2, gives rise to a switched differential-algebraic equation (DAE) represented by:

where, $x=\left(Q, P_0, P_L\right)^{\top}$,

$\sigma(t)= \begin{cases}1, & t \in\left[0, t_S\right) \\ 2, & t \geq t_S,\end{cases}$

and

where, $c_1=\frac{A}{L}>0$ and $c_2=\frac{c_f}{2 D A \rho_a}>0$. It is worth noting that the switched DAE (8) contains a nonlinear term $g_\sigma(x)$, which makes it unsuitable for direct application of the distributional solution framework discussed in Section 1.1. Previous studies, such as LibeTren12, have examined nonlinear switched DAEs; however, these approaches inherently exclude the incorporation of Dirac impulses in $x$. Since capturing the water hammer effect requires the presence of Dirac impulses in the solution $x$ of (8), the evaluation of $g_\sigma(x)$ in general becomes unclear (e.g., the determination of the sine of a Dirac impulse).

Lemma 2. Consider the nonlinear initial-trajectory problem (ITP)

$\begin{aligned} x_{(-\infty, 0)} & =x_{(-\infty, 0)}^0 \\ (E \dot{x})_{[0, \infty)} & =(A \dot{x}+f+g(x))_{[0, \infty)}\end{aligned}$

where either $(E, A)=\left(E_1, A_1\right)$ or $(E, A)=\left(E_2, A_2\right)$ and $g(x)=g_1(x)=g_2(x)$ as in (9). Then for every initial trajectory $x^0 \in\left(\mathbb{D}_{\mathrm{pw}} \mathscr{C}^{\infty}\right)^3$ and every inhomogeneity $f$ induced by a piecewise-smooth function, there exists a unique solution of the ITP in the sense of Definition 1.

Case $(E, A)=\left(E_1, A_1\right)$. The DAE $E \dot{x}=A x+f+g(x)$ with $x=\left(Q, P_0, P_L\right)^{\top}$ and $f=\left(f_1, f_2, f_3\right)^{\top}$ reads as

$\begin{aligned} \dot{Q} & =c_1\left(P_0-P_L\right)+f_1-c_2 Q|Q|, \\ 0 & =P_0+f_2, \\ 0 & =P_L+f_3,\end{aligned}$

which is equivalent to $P_0=-f_2, P_L=-f_3$ and

The Ordinary Differential Equation (ODE) (6) ensures the existence and uniqueness of local solutions, even in scenarios where the right-hand side lacks global Lipschitz continuity. This guarantee hinges significantly on the presence of a negative sign in the quadratic term. This negative term, acting as a form of a stabilizing influence, curtails any tendency for the solution to grow unbounded. It provides a sort of damping mechanism, enabling the solution to stay well-regulated.

Indeed, the absence of global Lipschitz continuity could complicate the analysis. However, the unique structure of the ODE empowers us to affirm the existence and uniqueness of solutions across all times. The presence of the negative sign in the quadratic term precludes the solution from diverging or showing unbounded growth, thereby maintaining the solution's well-behaved nature.

In the context of the Initial Transition Problem (ITP), this theorem of existence and uniqueness carries over directly from the ODE. Even if there are discrepancies in the initial conditions, signified by $P_0$ and $P_L$, that result in a discontinuity in the solution, they only manifest as simple jumps in the solution without impacting its overall existence and uniqueness.

Therefore, despite the lack of global Lipschitz continuity, the unique characteristics of the ODE and the negative sign in the quadratic term ensure the global existence and uniqueness of solutions. The ITP manages any discrepancies in the initial values via simple jumps, thereby preserving the overall well-posedness of the problem.

Case $(E, A)=\left(E_2, A_2\right)$.

The DAE $E \dot{x}=A x+f+g(x)$ with $x=\left(Q, P_0, P_L\right)^{\top}$ and $f=\left(f_1, f_2, f_3\right)^{\top}$ reads as

$\begin{aligned} \dot{Q} & =c_1\left(P_0-P_L\right)+f_1-c_2 Q|Q|, \\ 0 & =P_0+f_2, \\ 0 & =P_L+f_3,\end{aligned}$

For the corresponding ITP with initial trajectory $x^0=\left(Q^0, P_0^0, P_L^0\right)$ we can directly derive the following unique solution:

$\begin{aligned} Q= & Q_{(-\infty, 0)}^0-f_{3[0, \infty)} \\ P_0= & P_{0(-\infty, 0)}^0-\left(f_2\right)_{[0, \infty)} \\ P_L= & P_{L(-\infty, 0)}^0+\frac{1}{c_1}\left(-\dot{Q}+c_1 P_0+f_1-c_2 Q|Q|\right)_{[0, \infty)} \\ = & P_{L(-\infty, 0)}^0+\frac{1}{c_1} Q^0\left(0^{-}\right) \delta_0 \\ & +\frac{1}{c_1} \frac{\mathrm{d}}{\mathrm{d} t}\left(f_{3[0, \infty)}\right)+\left(-f_2+\frac{1}{c_1} f_1+\frac{c_2}{c_1} f_3\left|f_3\right|\right)_{[0, \infty)},\end{aligned}$

where, $\delta_0$ denotes the Dirac impulse located at $t=0$.

Corollary 3. The switched nonlinear DAE (8), (9) has for every initial condition $Q(0)=Q^0 \in \mathbb{R}$ a unique solution in the sense of Definition given in the study [1]. In particular, the jump and the Dirac impulse in $P_L$ at $t_S$ are given by:

$\begin{aligned} & P_L\left(t_S^{+}\right)=P_U \\ & P_L\left[t_S\right]=\frac{1}{c_1} Q\left(t_S^{-}\right) \delta_{t_s}\end{aligned}$

The ultimate objective of this study is to establish a rigorous proof of convergence, in an appropriate sense, between the solution $(q, P(\rho)$) of the balance law given by Eqns. (2), (4), (5), (7), and the solutions $\left(Q, P_0, P_L\right)$ of the switched DAE represented by Eqns. (8), (9), as the compressibility coefficient $\beta^\omega$ approaches zero. In this context, the focus lies on analyzing the abrupt change and Dirac impulse in pressure resulting from the closure of the valve. Specifically, it is assumed that the initial condition specified in Eq. (5) ensures the solution of the PDE on the interval $\left[0, t_S\right)$ remains stationary. In other words, both the spatial and temporal components, denoted by $q(t, x)$, become constant once the valve is closed, and the dynamics within the pipe stabilize.

For the numerical simulations, a third-order HEOC (High-Order Essentially Non-Oscillatory Conservative) scheme, as described in the study [7], is employed. The upper portion of Figure 3 illustrates the results for the pressure values at the valve over the time interval $[0.4 s, 4 s]$ with initial values [12], [13], [14].

$q_0^\omega(x) \equiv 0, \quad \rho_0(x) \equiv 1.4115 \times 10^3$

and pipe parameters

$\begin{aligned} P_a^\omega & =1.01 \times 10^6, & \beta^\omega & =\frac{1}{K^\omega}=4 \times 10^{-9}, & & \rho_a=1000, \\ L & =5, & D & =0.5, & & c_f=0.02 .\end{aligned}$

A significant surge in pressure is observed immediately following the switching time, denoted as $t_S=0.5 \mathrm{~s}$. Subsequently, the pressure exhibits periodic settling to a new equilibrium value. The frequency of this periodic behavior is dependent on the length of the pipe, represented by $L$ (larger $L$ corresponds to a lower frequency), and the speed of sound, which is higher for pipes with smaller compressibility coefficients $\beta^\omega$.

To facilitate a comparison between the results, it is necessary to approximate the pressure value $p(t, L)$ as $t$ approaches infinity. Instead of running the simulation for an extensive duration, a settling time $\epsilon>0$ is employed, and the average of $p(t, L)$ is computed over the interval $\left(t_S+\epsilon, T\right]$, where $T>t_S+\epsilon$ denotes the total duration of the simulation.

$\bar{P}_L:=\frac{1}{T-\left(t_S+\varepsilon\right)} \int_{t_S+\varepsilon}^T p(t, x) \mathrm{d} t$.

With

$\varepsilon=1.5, \quad T=4$

we obtain

$\bar{P}_L \approx 8.23 \times 10^8$

The value predicted by the switched DAE is

$P_L\left(t_S^{+}\right)=P_U \approx 8.23 \times 10^8$

In Table 1 the relative error between $\bar{P}_L$ and $P_L\left(t_S^{+}\right)$ is presented for decreasing compressibility coefficients $\beta^\omega$.

For the largest value of $\beta^\omega$, the value $P_L\left(t_S^{+}\right)$ is found to be a very good approximation of $\bar{P}_L$, and the approximation improves for decreasing $\beta^\omega$.

To compare the peak in $p(t, L)$ immediately after the valve closure with the Dirac impulse $P_L\left[t_S\right]$ that occurs at the switching time, the approximation of a Dirac impulse can be employed. Recall that a Dirac impulse $\delta_{t_s}$ at $t_s>0$ can be approximated by a sequence of functions denoted as $\delta_{t_s}^\epsilon(t)$. These functions satisfy the conditions $\delta_{t_s}^\epsilon(t)=0$ for $t \neq\left[t_s, t_s+\epsilon\right]$ and $\int_{t_s}^{t_s+\epsilon} \delta_{t_s}^\epsilon(t) d t=1$.

$p(t, L)=\bar{P}_{t_S}^{\mathrm{imp}} \delta^{\varepsilon}(t)+\bar{P}_L, \quad t \in\left(t_S, T\right]$ ,

hence we can approximate the magnitude of the “smoothed-out” Dirac impulse occurring in the PDE model as follows:

$\bar{P}_{t_S}^{\mathrm{imp}}:=\int_{t_S}^{t_S+\varepsilon} p(t, L)-\bar{P}_L \mathrm{~d} t$.

The Dirac impulse induced by the switched DAE is (see Corollary 3):

$P_L\left[t_S\right]=\frac{1}{c_1} \boldsymbol{Q}\left(t_S^{-}\right) \delta_{t_s}=: P_{t_S}^{\mathrm{imp}} \delta_{t_s}$ .

In order to evaluate $P_{t_S}^{\mathrm{imp}}$ for comparison, we will proceed in the following manner: We will evaluate the steady state of (10) alongwith f given in (9) for $t \in\left(0, t_s\right)$:

$\dot{Q}=c_1\left(f_3-f_2\right)+f_1-c_2 Q|Q|$,

with $f_1=0, f_2=P_U, f_3=P_D$, we get:

$\dot{Q}=c_1\left(P_U-P_D\right)-c_2 Q|Q|$,

to get steady state put $(\dot{Q}=0)$ then substituting $c_1=0.0393, c_2=1.0188 e^{-04}$, and $P_U$, $P_D$ from numerical setup we get:

$Q\left(t_s^{-}\right)=1.2830 \cdot 10^6, \quad P_{t_S}^{\mathrm{imp}}=\frac{1}{c_1} Q\left(t_s^{-}\right)=3.2670 \cdot 10^7$ .

A comparison between $\bar{P}_{t_S}^{\mathrm{imp}}$ and $P_{t_S}^{\mathrm{imp}}$ for different values of the compressibility coefficient $\beta^\omega$ at tS is presented in Table 2. The accuracy of the approximation is not very high for large values of $\beta^\omega$. However, as the compressibility coefficient $\beta^\omega$ decreases, the accuracy of the approximation improves significantly.

It is crucial to note that the selection of $\varepsilon$ directly affects the accuracy of the approximation. This is demonstrated in Table 3, which showcases the relationship between the choice of $\varepsilon$ and the accuracy of the approximation.

However, it is important to emphasize that the qualitative trend of the error decreasing as the compressibility coefficient decreases remains consistent. While the precise values presented in Table 3 may differ, the general observation of a decreasing error with decreasing compressibility coefficient remains valid.

In this research, a novel method for water hammer modelling using switched differential-algebraic equations (DAEs) is proposed. This method was juxtaposed with a more complex compressible nonlinear system of balance laws. Validation of the model was achieved by performing numerical simulations of the partial differential equations (PDE) model, which substantiated the hypothesis that the switched DAE model could effectively approximate the PDE model under low compressibility coefficients.

Future endeavours will revolve around two crucial aspects. One concerns the establishment of a rigorous proof of convergence, which is expected to validate the switched DAE model's suitability for approximating the more intricate compressible nonlinear system of balance laws. It is postulated that such a formal proof will elevate the trustworthiness and applicability of the switched DAE model, thereby advancing its practical use in water hammer phenomena analysis.

Simultaneously, it is intended to broaden the investigation to encapsulate more complex scenarios, such as larger networks. It is well-documented that water hammer incidents often transpire within intricate piping systems, thus, necessitating the consideration of multiple pipes and nodes. By venturing into these larger networks, an enriched comprehension of water hammer dynamics in practical applications can be developed. This expansion will necessitate the incorporation of additional variables, for instance, junction conditions and network topology, to depict accurately the complex hydraulic behaviour observed in large-scale systems.

The ultimate aspiration is to enhance predictive capabilities and formulate effective strategies for mitigating the deleterious effects of water hammer in real-world scenarios. The amalgamation of the switched DAE model and numerical simulations of the PDE model proffers a promising blueprint for understanding and managing water hammer phenomena, irrespective of the system's complexity.