Signal Processing for Transient Flow Rate Determination: An Analytical Soft Sensor Using Two Pressure Signals

Abstract

Accurate knowledge of the flow rate is essential for hydraulic systems, enabling the calculation of hydraulic power when combined with pressure measurements. These data are crucial for applications such as predictive maintenance. However, most flow rate sensors in fluid power systems operate invasively, disrupting the flow and producing inaccurate results, especially under transient conditions. Utilizing pressure transducers represents a non-invasive soft sensor approach since no physical flow rate sensor is used to determine the flow rate. Usually, this approach relies on the Hagen–Poiseuille (HP) law, which is limited to steady and incompressible flow. This paper introduces a novel soft sensor with an analytical model for transient, compressible pipe flow based on two pressure signals. The model is derived by solving fundamental fluid equations in the Laplace domain and converting them back to the time domain. Using the four-pole theorem, this model contains a relationship between the pressure difference and the flow rate. Several unsteady test cases are investigated and compared to a steady soft sensor based on the HP law, highlighting our soft sensor’s promising capability. It exhibits an overall error of less than $0.15$% for the investigated test cases in a distributed-parameter simulation, whereas the HP-based sensor shows errors in the double-digit range.

1. Introduction

Flow rate sensors are commonly employed in hydraulic engineering since knowing the flow rate is critical for detecting leaks and monitoring the wear of system components in both stationary and mobile hydraulic systems. Additionally, in the emerging field of predictive maintenance, knowing flow rates and pressure levels allows for calculating power, efficiency, and potential losses.

There are two main issues with current volumetric flow sensors: First, most must be installed directly into a pipe, often using turbines or rotors, which disrupt the flow and increase system complexity. Second, the physical components of these sensors possess inertia, making them ineffective for measuring flows above a particular frequency. Pump pulsation measurements, for instance, become problematic because they rely on flow rate sensors. The number of displacement units and the pump’s rotational speed influence the fundamental frequency of pump pulsations. Pumps also produce pulsations at harmonic frequencies, many of which exceed the capabilities of most state-of-the-art sensors. These limitations highlight the need for new flow measurement methods. An ideal flow sensor would be non-invasive and accurately measure transient flow conditions. Ideally, it would operate as a soft sensor, relying solely on signals from other sensors, particularly pressure signals, often available in fluid power systems, and are non-invasive.

The concept of calculating flow rates from pressure data is not new. Analytical methods for determining volumetric flow rates have been available since the 1900s. One of the most well-known formulas is the Hagen–Poiseuille law (HP) [1], which calculates flow rates based on the friction-induced pressure drop within a pipe. However, this law is limited to steady, laminar flows and neglects the effects that become significant at higher frequencies. The so-called Richardson effect [2] alters the velocity profile at higher frequencies. Another significant contribution is by Zhao et al., who created a method for measuring unsteady flow rates and velocities based on pressure differences, expressed through the transfer function between the mean flow velocity and the pressure gradient [3].

Hydraulic systems are often modeled analogously to electrical circuits, using resistances, inductors, and capacitances [4]. This approach is typically accurate and straightforward for many applications. However, at higher frequencies, the system behavior changes due to transient effects. For instance, the characteristic impedance, propagation operators, and wall shear stress within pipes have all been shown to depend on frequency [5].

Recent advances in transient flow rate determination include work by Brereton et al., who presented a method for arbitrary transients in laminar pipe flow, starting from an initially steady state. This method relates the flow rate to the pressure-gradient history without making assumptions about velocity profiles [6]. In 2008, Brereton et al. introduced an indirect method that links flow rate to the centerline velocity history [7]. Sundstrom et al. [8] enhanced frictional modeling, reducing errors in flow rate calculations for the pressure–time method, a technique commonly used to measure discharge in hydraulic machinery. Foucault et al. proposed a method in 2019 for time-resolved transient flow rate calculation based on differential pressure measurements [9]. Their approach utilized the kinetic energy to predict specific laminar flow conditions well in real-time, relying only on two coefficients. García García et al. focused on unsteady turbulent pipe flow, validating their method through the transient response of the velocity field due to external perturbations, using a pressure step function as a test case [10]. Further investigations into turbulence were conducted by García García et al. in 2022, studying the transition from turbulent to laminar flow without an increase in bulk velocity. The observed laminarization was explained using a derived mathematical model [11].

Urbanowicz et al. conducted a comprehensive review of analytical models for accelerated incompressible Newtonian fluid flow in pipes. They examined different approaches with imposed pressure gradients and flow rates, focusing on their mathematical complexity and applicability to laminar and turbulent flows. Although existing models predict specific laminar flow conditions well, they struggle with turbulent flow scenarios [12].

In 2023, Urbanowicz et al. developed an analytical solution for modeling laminar water hammer events in pipes, validated through numerical simulations and experimental measurements [13]. Later that year, Urbanowicz et al. published new analytical models for wall shear stress in water hammer scenarios. By employing quasi-steady and transient hydraulic resistance assumptions, these models broadened the range of validity and simplified the mathematical descriptions. Explicit analytical expressions for wall shear stress during water hammer events were derived and numerically validated [14]. In two separate studies, Bayle et al. explored wave propagation models for water hammer events in pipes. The first paper developed a rheology-based model for viscoelastic pipes and validated it using experimental data [15], while the second study provided a wave propagation model in the Laplace domain applicable to a variety of boundary conditions in pipe scenarios [16].

This paper aims to provide the mathematical foundation for a volumetric flow rate soft sensor based on two pressure signals by deriving an equation that accurately describes the flow rate in a pipe while considering the system’s transient behavior. The derived equation consists of the steady HP solution with an added term describing unsteady effects. This will be achieved by taking the Laplace transform (LT) of the fundamental fluid mechanics equations, solving them in the Laplace domain, and then deriving the inverse Laplace transform (ILT) of the equations in the Laplace domain to obtain a solution in the time domain. The Laplace transform has a significant advantage over the Fourier transform, as the Fourier domain is a subset of the Laplace domain; thus, the benefits of the Fourier transform are inherently included in the Laplace transform. Regarding the area of convergence, the Laplace transform exhibits a broader region of convergence, which is a half-plane compared to a unit circle for the Fourier transform. This extension results in the Laplace transform capturing both the system’s dissipation and oscillatory behavior, making it applicable to arbitrary pressure signals. A one-step integration method is utilized to decrease the computational costs of the soft sensor. The analytical model of the soft sensor based on two pressure transducers exhibits three main benefits compared to currently available sensors and numerical solvers: high accuracy of transient compressible flow, real-time capability, and, regarding application on a real test rig, no requirement of an actual volumetric flow rate sensor but two non-invasive pressure transducers. Analysis of a system’s high-frequency behavior allows the user to capture transients in the volumetric flow rate caused by fast changes in pressure. Valuable areas of application of this soft sensor include dynamic measurements of volumetric flow rates achieved by pumps and controlling hydraulic actuators through the correlation of the velocity of the piston to the volumetric flow rate times the cross-section. Additionally, the soft sensor allows for predictive maintenance of hydraulic systems by calculating the power from the pressure and the volumetric flow rate. Fast flow rate calculation allows a detailed simulative examination of a hydraulic system compared to the computationally demanding classical iterative solver.

Section 2 introduces the fundamentals of fluid mechanics and the Laplace techniques. Afterward, the analytical model for the soft sensor is derived and investigated in Section 2.3 to Section 2.8. The implementation of the soft sensor for the actual application is presented in Section 3. The soft sensor is validated for several test cases, presented in Section 3.2, and compared to a commonly used pressure-based soft sensor only relying on the law of HP to underline the high accuracy compared to the currently utilized methods and the relevance of incorporating unsteady phenomena in the analytical model. The results are shown in Section 4. Finally, the results are discussed in Section 5, and a conclusion is given in Section 6.

2. Materials and Methods

2.1. The Law of Hagen and Poiseuille and the Richardson Effect

Flow conditions can be characterized by the dimensionless Reynolds number, which gives a ratio between inertial and viscous forces. A flow is laminar if the Reynolds number $Re={\displaystyle \frac{\rho vD}{\eta }}$ is below approximately 2300. Here, $\rho$ is the density of the fluid, v is the axial fluid velocity, D is the pipe’s diameter, and $\eta$ is the dynamic viscosity. In this regime, the velocity profile across the pipe’s cross-section adopts a parabolic shape, with its vertex at the pipe center, and wave dissipation is disregarded. A critical aspect of this flow scenario is the no-slip condition at the pipe wall [17], which dictates that the fluid velocity at the wall is zero. Under these conditions, the Hagen–Poiseuille law [1] governs the flow:

$$Q={\displaystyle \frac{\pi R^4}{8\eta L}}\Delta p={\displaystyle \frac{1}{R_{\mathrm{H}}}}\Delta p$$

(1)

This formula provides the volumetric flow rate Q given the pressure difference $\Delta p$, with $R_{\mathrm{H}}$ representing the hydraulic resistance, dependent on the pipe radius R, fluid viscosity $\eta$, and pipe length L. The relevance of hydraulic resistance will be explored further in later sections.

In 1929, Richardson [2] conducted experiments on unsteady airflow through a pipe generated by sinusoidal pressure fluctuations. Using a hot wire anemometer, he measured the radial velocity profile. He found that, within the laminar flow range of Reynolds numbers, increasing pressure frequencies caused the velocity profile to deviate from a parabolic shape, as shown on the left side of Figure 1. In this figure, the Womersley number $Wo=R\sqrt{\omega /\nu }$ is the root of the normalized frequency of the harmonically oscillating pressure input. It is based on the assumption of an incompressible fluid [5], where R is the pipe radius, $\nu$ the kinematic viscosity, and $\omega$ the pressure variation frequency.

As the profile changes, it transitions to a flatter, more uniform shape in the center, with peaks near the wall and a sharp drop to zero at the boundary. This unsteady profile resembles those found in steady, turbulent flows. The right side of Figure 1 illustrates the variation in the radial velocity profile over the phase of the harmonically oscillating pressure input for a Womersley number where the Richardson effect becomes apparent ($Wo=30$).

Building on Richardson’s findings, Sexl [18] derived solutions for the instantaneous velocity profile as a function of the normalized radial position and the Womersley number. Stecki [5] performed the nondimensionalization of the axes in this context. The system equations developed in this paper allow for an analytical expression of the velocity profile transformation from parabolic to the one observed in Richardson’s experiments, thus extending the applicability of the presented equations beyond traditional laminar flow with a parabolic profile.

Since the early 20th century, various fluid behavior models have been developed, primarily from the fundamental Navier–Stokes equations. These models often incorporate simplifying assumptions to make them more computationally feasible [19,20,21,22]. Typical assumptions include liquid incompressibility, neglecting heat transfer effects, and the assumption of uniform pressure distribution across the pipe’s cross-section.

Stecki reviewed and categorized these models in his work [5,23], organizing them according to dimensionality and the effects considered. This paper’s general governing equations are those of the two-dimensional viscous compressible model. The volumetric flow rate equation is derived from this general framework.

2.2. Signal Processing Based on Laplace Techniques

Transforming functions from the time domain into the Laplace domain is a well-established technique in fields dealing with differential equations [24,25,26,27]. When this transformation is applied, all functions that depend on t are expressed in terms of s, where s is the Laplace variable. One key advantage of this approach is that time derivatives become simple multiplications by s, simplifying the differential equations. In the case of a differential equation that depends on one spatial operator and one temporal operator, the partial differential equation is reduced to an ordinary differential equation. Moreover, the Laplace transformation enables the analysis of arbitrary pressure signals, which is impossible with the Fourier method as it restricts solutions to harmonic, oscillating signals. This work aims to transform the system equations derived in the Laplace domain back into the time domain, ultimately obtaining an accurate equation that links the flow to two pressure signals in the time domain.

One needs to perform an inverse Laplace transformation to revert from the Laplace domain to the time domain. The function must meet the conditions for an ILT, notably that the function in the Laplace domain must converge to zero as $s\to \infty$. The ILT can be computed using a correspondence table after performing partial fraction decomposition for many Laplace domain functions. For example, Hsiao et al. calculated the hydrodynamic force and torque on a sphere using ILT techniques aided by integral tables [25]. However, in cases where partial fraction decomposition is not feasible, and no correspondence table exists, a general formula is needed. This is provided by the Bromwich contour integral (Equation (2)), which can be used to obtain the ILT of most functions [28]. In this transformation, $f\left(t\right)$ and $F^{*}\left(s\right)$ represent the same function in the time and Laplace domains, respectively, with * indicating the function in the Laplace domain. The poles of the function $F^{*}\left(s\right)$ are denoted as $s=s_j$.

The Bromwich contour integral does not need to be directly computed in many cases. It can be evaluated indirectly using Cauchy’s contour integral, the Residue theorem from complex analysis, or numerical methods. For instance, Young et al. used a numerical approach to determine the ILT and obtain a solution for three-dimensional time-dependent circulation in shallow lakes [24]. Manhartsgruber conducted further work for instantaneous volumetric flow rates in circular pipes [26,27]. He utilized the fast inverse Laplace transformation (FILT) to obtain a representation of the physical signal. The FILT inherently contains a truncation error and limits the validity range of the determined result to instantaneous volumetric flow rates.

The Residue theorem’s advantage over numerical methods is that it provides an exact formula, offering more profound insight into the investigated system. Additionally, it extends the range of validity compared to numerically derived solutions and avoids concerns about numerical instabilities.

The integral over the complex plane is replaced by the sum of the residues of the function at its poles. Residues are specific coefficients found in the Laurent series expansion of a function around a pole [29]. It is worth noting that the Residue theorem is especially useful for dealing with functions that have complex poles.

The following theorems and formulae (Equations (2)–(4)) regarding the fundamentals of the residue theorem and Laplace transformation can be found in the work by Krantz [29].

$$f\left(t\right)={\displaystyle \frac{1}{2\pi i}}{\int }_{\beta -\mathrm{i}\infty }^{\beta +\mathrm{i}\infty }{\mathrm{e}}^{st}{F}^{*}\left(s\right)\mathrm{d}s=\sum _{j=0}^N\mathrm{Res}\left({\mathrm{e}}^{st}{F}^{*}\left(s\right)\right){|}_{s=s_j}$$

(2)

The residues for poles of order $m>1$ and $m=1$ are calculated using the general formula (Equation (3)), which stems from the definition of the Laurent series.

$$\mathrm{Res}\left(F^{*}\left(s\right)\right){|}_{s=s_j}=\underset{s\to s_j}{\lim }{\displaystyle \frac{1}{(m-1)!}}{\displaystyle \frac{{\mathrm{d}}^{m-1}}{\mathrm{d}{s}^{m-1}}}\left[F^{*}\left(s\right){(s-s_j)}^m\right]$$

(3)

In the case of $m=1$, and when the function is expressed as a fraction with $X^{*}\left(s\right)$ as the numerator and $Y^{*}\left(s\right)$ as the denominator, the residues for simple poles are calculated using Equation (4). The denominator is differentiated, and the entire term is evaluated at the simple poles $s=s_k$:

$$\mathrm{Res}\left(F^{*}\left(s\right)\right){|}_{s=s_k}=\mathrm{Res}\left({\displaystyle \frac{X^{*}\left(s\right)}{Y^{*}\left(s\right)}}\right){|}_{s=s_k}={\displaystyle \frac{X^{*}\left(s_k\right)}{\frac{\partial Y^{*}\left(s\right)}{\partial s}{|}_{s=s_k}}}$$

(4)

For practical applications of these formulae, see Xie [30]. It is also crucial to remember that the product of two functions in the Laplace domain is equivalent to the convolution integral in the time domain [28]:

$$F^{*}\left(s\right)G^{*}\left(s\right)=f\left(t\right)\ast g\left(t\right)={\int }_0^{t}f(t-\tau )g\left(\tau \right)\mathrm{d}\tau$$

(5)

Furthermore, the inverse Laplace transform (denoted as ${\mathcal{L}}^{-1}$) of a function multiplied by s is given by:

$${\mathcal{L}}^{-1}\left\{s{F}^{*}\left(s\right)\right\}={\displaystyle \frac{\partial f\left(t\right)}{\partial t}}$$

(6)

2.3. Analytical Model for Signal Processing of Pressure Transducer Signals

The process of the derivation of an equation that can be used in a soft sensor to calculate the volumetric flow rate in dependency on the pressure difference over the length of a pipe is presented in Figure 2. First, the Navier–Stokes equations and the equation of conservation of mass are simplified by using assumptions, yielding a set of partial differential equations that describe the transient flow of fluid in a pipe. The assumptions are listed in Section 2.3. In the second step, the time dependencies in the differential equations are eliminated by taking the Laplace transform. Afterward, the differential equations are solved in the Laplace domain to obtain a Laplace domain solution for the volumetric flow rate. These two steps can be found in detail in a previously published work by Brumand et al. [31]. The resulting equation in the Laplace domain is shown in Section 2.3. Lastly, as shown in Section 2.4, Section 2.5, Section 2.6, Section 2.7 and Section 2.8, the inverse Laplace transform of the solution is derived to obtain a solution for the volumetric flow rate in the time domain. Section 3 treats the resulting time domain solution.

Generally, the Navier–Stokes equations govern fluid flow in three-dimensional space. However, solving the full Navier–Stokes equations is unnecessary for this case of liquid flow in a rigid pipe. A simplified version of the equations can be used to obtain a solution. The following assumptions are applied to simplify the Navier–Stokes equations:

• The flow is laminar (Reynolds number $Re\le 2300$), meaning the fluid layers do not mix. Consequently, the pressure remains constant across the pipe’s cross-section, and the pressure gradient in the radial direction is negligible.
• Both the flow and pipe geometry are axisymmetric, resulting in no gradient in the angular direction, i.e., $\left({\displaystyle \frac{\partial }{\partial \phi }}=0\right)$.
• The axial flow velocity $v_{\mathrm{x}}$ is much smaller than the speed of sound a, which governs the fluid’s pressure wave propagation speed. Therefore, $v_{\mathrm{x}}\ll a$, and the Mach number $Ma=v_{\mathrm{x}}/a\ll 1$, so supersonic effects can be ignored.
• The pipe radius is much smaller than the pipe length ($R\ll L$), meaning no pressure reflections occur at the pipe walls.
• The significant viscous effects in the motion equations are limited to those involving the radial distribution of axial velocity [5].
• Gravitational forces are negligible as the pipe is horizontal, keeping the forces constant over its length.
• Changes in fluid density due to vertical positioning are negligible because the pipe’s diameter is small.
• Heat transfer is ignored since the focus is on liquids, excluding gases [32].

Using these assumptions, the Navier–Stokes equations can be simplified, and the general equation for the volumetric flow rate in the Laplace domain was derived in prior works [5,33]. The solution for the volumetric flow rate at the outlet of the pipe ($Q_2^{*}$) can be expressed as follows:

$$Q_2^{*}={\displaystyle \frac{1}{Z_{\mathrm{c}}\sinh \left(\gamma L\right)}}{p}_1^{*}-{\displaystyle \frac{\coth \left(\gamma L\right)}{Z_{\mathrm{c}}}}{p}_2^{*}=W_1^{*}\left(\zeta \right)p_1^{*}-W_2^{*}\left(\zeta \right)p_2^{*}.$$

(7)

Inserting the definitions of the impedance $Z_{\mathrm{c}}$ and the propagation operator $\gamma$ into the equation and multiplying both sides by $R_{\mathrm{H}}$ yields Equation (8).

$$Q_2^{*}{R}_{\mathrm{H}}={\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sinh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}{p}_1^{*}-{\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\tanh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}{p}_2^{*}.$$

(8)

With:

$$\begin{array}{c}\hfill \gamma L={\displaystyle \frac{s}{a}}\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}=s{\displaystyle \frac{R^2}{\nu }}{\displaystyle \frac{\nu }{R^2}}{\displaystyle \frac{L}{a}}\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}=\zeta Dn\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}.\end{array}$$

(9)

$Dn$ is the dissipation number, which relates the time span of the axial propagation of pressure waves to the radial diffusion of axial momentum, and $\zeta$ is the normalized Laplace variable defined as the following:

$$Dn={\displaystyle \frac{\nu L}{R^{2}a}}\hspace{1em}\hspace{1em}\hspace{1em}\zeta =s{\displaystyle \frac{R^2}{\nu }}$$

(10)

Also, $R_{\mathrm{H}}$ is the hydraulic resistance introduced in Section 2.1. Note that $I_n\left(x\right)$ is the modified Bessel function of the first kind of order n. In the later analysis of the ILT criteria, it became necessary to introduce an expansion factor of $\zeta /\zeta$ to ensure all required criteria were satisfied. The denominator is incorporated into the weighting functions, while the numerator acts as a prefactor to the pressure. As a result, the time gradient of the pressures is considered in the time domain equation rather than the pressures themselves.

$$Q_2^{*}{R}_{\mathrm{H}}=W_1^{*}\left(\zeta \right){\displaystyle \frac{\zeta }{\zeta }}{p}_1^{*}-W_2^{*}\left(\zeta \right){\displaystyle \frac{\zeta }{\zeta }}{p}_2^{*}$$

(11)

Now, the $\zeta$ in the denominator will always be included in the weighting functions.

$$Q_2^{*}{R}_{\mathrm{H}}={\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sinh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}\zeta p_1^{*}-{\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\tanh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}\zeta p_2^{*}$$

(12)

2.4. Inverse Laplace Transformation of the Weighting Functions for Compressible Fluids

Unlike in the incompressible case [31], for compressible fluids, the weighting functions of the pressures at both ends of the pipe are not equal due to the finite wave propagation speed. Therefore, each of the functions must be analyzed separately. The prefactors of the pressure functions from Equation (8) will be written as weighting functions:

$$\begin{array}{c}\hfill W_1^{*}\left(\zeta \right)={\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sinh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}},\end{array}$$

(13)

$$\begin{array}{c}\hfill W_2^{*}\left(\zeta \right)={\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\tanh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}.\end{array}$$

(14)

Note the change from $\sinh \left(x\right)$ to $\tanh \left(x\right)$, which effectively is only a multiplication of $W_1^{*}$ with a $1/\cosh \left(x\right)$ term. The weighting function $W_1^{*}\left(\zeta \right)$ for a fixed dissipation number $Dn$ can be seen in Figure 3. Note the poles on the negative real axis starting from zero. For increasing values of $Dn$, the poles would move further away from the origin.

2.5. Investigation of the Necessary Criteria for an ILT

First, the existence of the inverse Laplace transforms for the weighting functions must be established. To do this, it is essential that the limits of the weighting functions approach zero as $s\to \infty$ and that the limits of the weighting functions multiplied by s remain finite as $s\to 0$. These conditions apply when using a normalized Laplace variable $\zeta$ instead of s. In what follows, the four required limits to confirm the existence of the inverse Laplace transforms for $W_1^{*}\left(\zeta \right)$ and $W_2^{*}\left(\zeta \right)$ will be examined. Firstly, the limits for $\zeta \to 0$.

$$\begin{array}{cc}\hfill & \underset{\zeta \to 0}{\lim }{W}_1^{*}\left(\zeta \right)\zeta \hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to 0}{\lim }8Dn{\displaystyle \frac{\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}{\displaystyle \frac{1}{\sinh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]}}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to 0}{\lim }8Dn1{\displaystyle \frac{1}{\infty }}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to 0}{\lim }{\displaystyle \frac{8Dn}{\infty }}=0\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & \underset{\zeta \to 0}{\lim }{W}_2^{*}\left(\zeta \right)\zeta \hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to 0}{\lim }8Dn{\displaystyle \frac{\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}{\displaystyle \frac{1}{\tanh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]}}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to 0}{\lim }8Dn1{\displaystyle \frac{1}{1}}=8Dn\hfill \end{array}$$

(21)

Since $Dn$ must be finite, the overall limit must also be finite. Next, the two limits as $\zeta \to 0$ are evaluated.

$$\begin{array}{cc}\hfill & \underset{\zeta \to \infty }{\lim }{W}_1^{*}\left(\zeta \right)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to \infty }{\lim }{\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sinh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to \infty }{\lim }{W}_1^{*}\left(\zeta \right)\zeta \underset{\zeta \to \infty }{\lim }{\displaystyle \frac{1}{\zeta }}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & =0{\displaystyle \frac{1}{\infty }}=0\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & \underset{\zeta \to \infty }{\lim }{W}_2^{*}\left(\zeta \right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to \infty }{\lim }{\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\tanh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & =\underset{\zeta \to \infty }{\lim }{W}_1^{*}\left(\zeta \right)\zeta \underset{\zeta \to \infty }{\lim }{\displaystyle \frac{1}{\zeta }}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & =8Dn{\displaystyle \frac{1}{\infty }}=0\hfill \end{array}$$

(29)

The outcomes of these calculations confirm that all the required conditions for an inverse Laplace transformation are satisfied.

2.6. Approximation of the Hyperbolic Sine Term

In the following, the weighting function $W_1^{*}$ will be examined, but it is important to note that the procedure for $W_2^{*}$ is analogous to that for $W_1^{*}$. This can be inferred from Equation (12), as only the poles of a function and its zeros, located at the same positions as the poles, are relevant for an inverse Laplace transform (ILT). Since cosh does not introduce additional poles or zeros, all formulas applied to $W_1^{*}$ can also be used for $W_2^{*}$.

$$W_2^{*}\left(\zeta \right)=W_1^{*}\left(\zeta \right)\cosh \left[\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}\zeta Dn\right]$$

(30)

With:

$$W_1^{*}\left(\zeta \right)={\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{\sinh \left[\sqrt{\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}\zeta Dn\right]\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}.$$

(31)

The order of the poles of a function is crucial when calculating the residues needed for an ILT. To determine the order of the poles, we use an approximation for sinh from Goodson [34], as shown in Equation (32). This approximation is ideal because it is expressed as a product, preserving the poles’ properties and the function’s zeros. In contrast, Taylor series expansions do not maintain the poles and are thus unsuitable [35]. Increasing n in Equation (32) corresponds to increasing the number of complex zeros of the sinh term.

$$\sinh \left(x\right)\approx {\sinh }_{\mathrm{approx}}(n,x)=x\left[\prod _{k=1}^n\left(1+{\displaystyle \frac{x^2}{k^2\ast {\pi }^2}}\right)\right]$$

(32)

By inserting $n=5$ and $x=\left[\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}\zeta Dn\right]$, Equation (32) becomes:

$$\begin{array}{c}\hfill {\sinh }_{\mathrm{approx}}\left(5,\left[\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}\zeta Dn\right]\right)=\sqrt{{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}}\zeta Dn\left(1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left({\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}}\right)\\ \hfill \left(1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(4{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}}\right)\left(1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(9{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}}\right)\left(1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(16{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}}\right)\\ \hfill \left(1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(25{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}}\right).\end{array}$$

(33)

It is important to note that the approximation ${\sinh }_{\mathrm{approx}}(5,x)$ is equivalent to the standard “small angle” approximation of $\sinh \left(x\right)=x$, which is derived from a truncated Taylor series and multiplied by the first five complex zeros of $\sinh \left(x\right)$.

$$W_{1,\mathrm{app}}^{*}\left(\zeta \right)={\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{{\sinh }_{\mathrm{approx}}(5,x)\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}$$

(34)

With this approximation, the weighting function can be separated into the incompressible part $W_{\mathrm{inc}}^{*}\left(\zeta \right)$, which is derived using the small angle approximation of sinh for $Dn\to 0$, and additional terms from the sinh approximation.

$$\begin{array}{cc}\hfill W_{1,\mathrm{app}}^{*}\left(\zeta \right)=& {\displaystyle \frac{8Dn\sqrt{{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}}{{\sinh }_{\mathrm{approx}}(5,x)\zeta \sqrt{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}}\hfill \\ \hfill =& W_{\mathrm{inc}}^{*}\left(\zeta \right){\displaystyle \frac{1}{\left(1+\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left({\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}\right)\left(1+\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(4{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}\right)\left(1+\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(9{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}\right)}}\hfill \\ \hfill & ·{\displaystyle \frac{1}{\left(1+\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(16{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}\right)\left(1+\frac{{\mathrm{I}}_0\left(\sqrt{\zeta }\right){\zeta }^{2}D{n}^2}{\left(25{\mathrm{I}}_2\left(\sqrt{\zeta }\right){\pi }^2\right)}\right)}},\hfill \end{array}$$

(35)

where $W_{\mathrm{inc}}^{*}\left(\zeta \right)={\displaystyle \frac{8{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}{{\zeta }^2{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}$. Note that $W_{\mathrm{inc}}^{*}\left(\zeta \right)$ can be obtained by taking the limit of ${\lim }_{Dn\to 0}{W}_1^{*}$ using ${\lim }_{x\to 0}\sinh \left(x\right)={\lim }_{x\to 0}\tanh \left(x\right)=x$, which are the Taylor series approximations of the hyperbolic terms. Although $W_{1,\mathrm{app}}^{*}$ will not be used to compute the residues, it demonstrates that the complex poles of sinh are simple. With this information, the ILT can be calculated from the residues of the function, which requires knowledge of the order and location of the poles. In conclusion, the weighting functions exhibit a double pole at $\zeta =0$, simple negative real poles from the incompressible part, and simple, complex conjugate poles that depend on $Dn$, which arise from the fluid’s compressibility effects.

2.7. Derivation of a Solution for a Fixed Dissipation Number for a Compressible Fluid

The previous chapter demonstrated that the poles of the weighting function depend on the dissipation number, making it impossible to derive a singular solution through an inverse Laplace transform (ILT) of the weighting function. This chapter derives the ILT for a fixed dissipation number of $Dn=0.0011$. This value is selected based on the pipe characteristics used in the test bench later in this project. Specifically, $Dn=0.0011$ corresponds to a pipe length of $L=1.5\mathrm{m}$, a diameter of $D=14\mathrm{mm}$, and standard fluid properties of HLP46, including a density of $\rho =860\mathrm{kg}/{\mathrm{m}}^3$, a kinematic viscosity of $\nu =46\times {10}^{-6}{\mathrm{m}}^2/\mathrm{s}$, and a bulk modulus of $K=14,000\mathrm{bar}$.

The general formula for calculating residues, provided in Equation (4), is now applied to determine the residue of the double pole at ${\zeta }_0=0$.

$$\begin{array}{cc}\hfill \mathrm{Res}\left(W_{\mathrm{inc}}^{*}\left(\zeta \right)\right){|}_{{\zeta }_0=0}=& \underset{\zeta \to {\zeta }_0}{\lim }{\displaystyle \frac{1}{(m-1)!}}{\displaystyle \frac{{\mathrm{d}}^{m-1}}{\mathrm{d}{\zeta }^{m-1}}}\left[W_{\mathrm{inc}}^{*}\left(\zeta \right){(\zeta -{\zeta }_0)}^m\right]\hfill \\ \hfill \stackrel{m=2}{=}& \underset{\zeta \to {\zeta }_0}{\lim }{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\zeta }}\left[W_{\mathrm{inc}}^{*}\left(\zeta \right){(\zeta -0)}^2\right]\hfill \\ \hfill =& \underset{\zeta \to {\zeta }_0}{\lim }{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\zeta }}\left[W_1^{*}\left(\zeta \right){\zeta }^2\right]\hfill \\ \hfill =& \underset{\zeta \to {\zeta }_0}{\lim }{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\zeta }}\left[\left({\displaystyle \frac{8{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}{{\zeta }^2{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}\right){\zeta }^2\right]\hfill \\ \hfill =& \underset{\zeta \to {\zeta }_0}{\lim }{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\zeta }}\left[{\displaystyle \frac{8{\mathrm{I}}_2\left(\sqrt{\zeta }\right)}{{\mathrm{I}}_0\left(\sqrt{\zeta }\right)}}\right]\hfill \\ \hfill =& 1\hfill \end{array}$$

(36)

The equation for the residues of simple, complex conjugated poles is given in Equation (5).

$$\mathrm{Res}\left(W_1^{*}\left(\zeta \right)\right){|}_{{\zeta }_k}={\displaystyle \frac{\mathrm{numerator}\left[W_1^{*}\left({\zeta }_k\right)\right]{\mathrm{e}}^{{\zeta }_k{t}_{\mathrm{n}}}}{\frac{\partial }{\partial \zeta }\left(\mathrm{denominator}\left[W_1^{*}\left({\zeta }_k\right)\right]\right)}}$$

(37)

The poles of the compressible component of the weighting function can be identified by solving $\sinh \left[x\right]=0$ along with Equation (32). This leads to the following expression:

$$1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{{\zeta }_k}\right){\zeta }_k^{2}D{n}^2}{k^2{\mathrm{I}}_2\left(\sqrt{{\zeta }_k}\right){\pi }^2}}=0,\mathrm{for}k\in \mathbb{N}.$$

(38)

Rearranging for ${\zeta }_k$, we obtain:

$${\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{{\zeta }_k}\right){\zeta }_k^2}{{\mathrm{I}}_2\left(\sqrt{{\zeta }_k}\right)}}={\left({\displaystyle \frac{k\pi }{Dn}}\right)}^2-1.$$

(39)

Assuming $Dn=0.0011$, the poles are located at:

• ${\zeta }_1=-38.30\pm 2818.21i$
• ${\zeta }_2=-53.95\pm 5658.55i$
• ${\zeta }_3=-65.96\pm 8502.53i$
• ${\zeta }_4=-76.08\pm 11348.4i$
• ${\zeta }_5=-85.00\pm 14195.4i$.

When the first five complex poles and their conjugates are substituted into the formula for the residues Res, the resulting values are as follows:

$$\begin{array}{cc}\hfill W_{1,\mathrm{comp}}& \left(t_{\mathrm{n}}\right)=\hfill \\ \hfill -& (0.0000180+0.00270i){\mathrm{e}}^{(-38.30-2818.21i)t_{\mathrm{n}}}+(-0.000018+0.00270i){\mathrm{e}}^{(-38.30+2818.21i)t_{\mathrm{n}}}\hfill \\ \hfill +& (0.0000066+0.00140i){\mathrm{e}}^{(-53.95-5658.55i)t_{\mathrm{n}}}+(0.0000066-0.00140i){\mathrm{e}}^{(-53.95+5658.55i)t_{\mathrm{n}}}\hfill \\ \hfill -& (0.0000036+0.00093i){\mathrm{e}}^{(-65.96-8502.53i)t_{\mathrm{n}}}+(-0.0000036+0.00093i){\mathrm{e}}^{(-65.96+8502.53i)t_{\mathrm{n}}}\hfill \\ \hfill +& (0.0000023+0.00070i){\mathrm{e}}^{(-76.08-11348.4i)t_{\mathrm{n}}}+(0.0000023-0.00070i){\mathrm{e}}^{(-76.08+11348.4i)t_{\mathrm{n}}}\hfill \\ \hfill -& (0.0000016+0.00056i){\mathrm{e}}^{(-85.00-14195.4i)t_{\mathrm{n}}}+(-0.0000016+0.00056i){\mathrm{e}}^{(-85.00+14195.4i)t_{\mathrm{n}}}.\hfill \end{array}$$

(40)

The method is based on a sum of poles with increasing imaginary parts. Thus, the equation can be expanded to include the system behavior corresponding to any desired frequency. The five poles considered were sufficient for the model to accurately capture the fluid flow behavior in all examined test cases. The curve-fitted ILT of the incompressible fluid is independent by $Dn$ and is presented in Equation (41):

$$\begin{array}{c}\hfill W_{\mathrm{inc}}^{*}\left(t_{\mathrm{n}}\right)=1-0.956788{\mathrm{e}}^{-5.783.t_{\mathrm{n}}}-0.03446{\mathrm{e}}^{-30.4713t_{\mathrm{n}}}-0.0006{\mathrm{e}}^{-{10}^{1.5}{t}_{\mathrm{n}}}\\ \hfill -0.0076{\mathrm{e}}^{-{10}^2{t}_{\mathrm{n}}}-0.0002{\mathrm{e}}^{-{10}^{2.5}{t}_{\mathrm{n}}}-0.0005{\mathrm{e}}^{-{10}^3{t}_{\mathrm{n}}}.\end{array}$$

(41)

The normalized time $t_{\mathrm{n}}$ corresponding to the normalized Laplace variable $\zeta$ is:

$$\begin{array}{cc}\hfill t_{\mathrm{n}}& ={\displaystyle \frac{\nu }{R^2}}t.\hfill \end{array}$$

(42)

A previously published paper provides further details about the derivation of the incompressible model describing the volumetric flow rate calculation based on true pressure signals [31].

Combining both equations, we obtain the final expression for the volumetric flow rate in the time domain for a specified dissipation number. The impact of incorporating the compressible residues into the incompressible equation is illustrated in Figure 4.

The complex conjugate poles introduce a step function behavior in the weighting function, which is attributed to the fluid’s compressibility. During the time it takes for pressure information to travel from one end of the pipe to the other, the weighting function for that pressure must remain zero. This is emphasized in the figure by the vertical line marking the moment in time that corresponds to the selected dissipation number $Dn$.

Notably, the dissipation number reflects the time required for the pressure to traverse the length of the pipe, expressed as $t_{\mathrm{n}}=L/a$. It is important to note that the vertical line at $t_{\mathrm{n}}=Dn$ coincides with the first step increase, indicating the moment when pressure information reaches the far end of the pipe. This occurs regardless of the other variables defined in $Dn$. Changes in the radius or kinematic viscosity of the fluid will also alter the normalized time in relation to $Dn$.

Furthermore, as illustrated in Figure 4, the weighting function oscillates around zero instead of remaining at zero. This oscillation results from the approximation of the hyperbolic sine term. Including additional terms in the sine approximation would progressively smooth out the step functions with each added term.

For completeness, the inverse Laplace transform (ILT) of the second weighting function, $W_2$, is shown in Figure 5 alongside the ILT of the weighting function $W_1$. These weighting functions are applied to the pipe’s corresponding pressures at both ends. Notably, $W_2$ does not exhibit the delay time characteristic of the first weighting function, as pressure information at the second end of the pipe instantaneously influences the volumetric flow rate at that end.

2.8. Derivation of a Solution for Arbitrary Dissipation Numbers

As previously mentioned, the inverse Laplace transform (ILT) of the weighting function for both ends of the pipe in the case of an incompressible fluid is independent of the dissipation number $Dn$ and was presented in Equation (41). The primary challenge in determining the general ILT of $W_1^{*}$ is that the poles depend on $Dn$. The kth pole is indirectly defined by the following formula, derived from Equation (32) for $k>0$:

$$1+{\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{{\zeta }_k}\right){\zeta }_k^{2}D{n}^2}{k^2{\mathrm{I}}_2\left(\sqrt{{\zeta }_k}\right){\pi }^2}}=0.$$

(43)

Rearranging for ${\zeta }_k$ yields:

$${\displaystyle \frac{{\mathrm{I}}_0\left(\sqrt{{\zeta }_k}\right){\zeta }_k^2}{{\mathrm{I}}_2\left(\sqrt{{\zeta }_k}\right)}}={\left({\displaystyle \frac{k\pi }{Dn}}\right)}^2-1.$$

(44)

Since it is not feasible to derive a solution applicable to all $Dn$, the calculations will focus on discrete values of $Dn$.

Using the software Maple [36], the first ten complex conjugate poles of the weighting functions were computed for dissipation numbers ranging from $1\times {10}^{-7}$ to $1\times {10}^{-1}$, with 100 increments per order of magnitude. This produced a matrix of dimensions 20 by 700. The positions of the poles in the complex plane for three different dissipation numbers are depicted in Figure 6. It can be observed that as the dissipation number decreases, the poles move further from the origin in the complex plane. Note that poles positioned further to the left on the negative half of the real axis correspond to faster-decaying residues, while those further away from zero on the imaginary axis are associated with faster-oscillating residues.

With the matrix of poles dependent on the dissipation number now in hand, we can calculate the residue for each pole across all considered dissipation number cases. The equation for the residues of simple poles is provided in Equation (45). Since all poles of the weighting function resulting from $\sinh \left(x\right)$ have been confirmed to be simple poles (as shown in Equation (35)), this equation can be utilized for the automated computation of the residues:

$$\mathrm{Res}\left(W_1^{*}\right){|}_{{\zeta }_k}={\displaystyle \frac{\mathrm{numerator}\left[W_1^{*}\left({\zeta }_k\right)\right]{\mathrm{e}}^{{\zeta }_k{t}_{\mathrm{n}}}}{\frac{\partial }{\partial \zeta }\left(\mathrm{denominator}\left[W_1^{*}\left({\zeta }_k\right)\right]\right)}}.$$

(45)

It is worth noting that in the exponent of $\mathrm{e}$:

$$s_{k}t=s_k{\displaystyle \frac{R^2}{\nu }}{\displaystyle \frac{\nu }{R^2}}t={\zeta }_k{\displaystyle \frac{\nu }{R^2}}t={\zeta }_k{t}_{\mathrm{n}}.$$

(46)

We can compute a matrix containing all residues using Maple and Equation (45). By combining this with the residues and poles from the incompressible case (as represented in Equation (41)), we can obtain the weighting function in the time domain for discrete values of $Dn$. The rounding of actual $Dn$ values to the nearest sampled values did not visibly impact the results.

3. Implementation of the Soft Sensor in the Time Domain

In this chapter, the solution of the equation for volumetric flow rate in the time domain is evaluated. This involves a detailed examination of the convolution integral obtained as a solution, which will be simplified in this section.

To begin, we revisit Equation (12). For clarity, this equation is restated as Equation (47). To determine the volumetric flow rate in the time domain, we need to revert this equation into that domain:

$$Q_2^{*}{R}_{\mathrm{H}}=W_1^{*}\zeta p_1^{*}-W_1^{*}\zeta p_2^{*}.$$

(47)

By applying the principle that the multiplication of two functions in the Laplace domain corresponds to the convolution of those same functions in the time domain, along with the fact that multiplying a function by the Laplace variable equates to taking the time derivative, we obtain:

$$Q_2\left(t_{\mathrm{n}}\right)R_{\mathrm{H}}=\left[{\int }_0^{t_{\mathrm{n}}}{W}_1(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial p_1(t_{\mathrm{n}}-\tau )}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}{W}_2\left(\tau \right){\displaystyle \frac{\partial p_2(t_{\mathrm{n}}-\tau )}{\partial t_{\mathrm{n}}}}d\tau \right].$$

(48)

Since the first term of both weighting functions is a constant, we can separate $W_1$ and $W_2$, placing the constant term in a separate integral:

$$\begin{array}{cc}\hfill Q_2\left(t_{\mathrm{n}}\right)R_{\mathrm{H}}=& \left[{\int }_0^{t_{\mathrm{n}}}1{\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau +{\int }_0^{t_{\mathrm{n}}}(W_1(t_{\mathrm{n}}-\tau )-1){\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right]\hfill \\ \hfill & -\left[{\int }_0^{t_{\mathrm{n}}}1{\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau +{\int }_0^{t_{\mathrm{n}}}(W_2(t_{\mathrm{n}}-\tau )-1){\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right].\hfill \end{array}$$

(49)

It is practical to reorganize the equations into constant and dynamic components. The constant integral can be evaluated as the pressure difference given by $p_1\left(t_{\mathrm{n}}\right)-p_2\left(t_{\mathrm{n}}\right)=\Delta p$.

$$\begin{array}{cc}\hfill Q_2\left(t_{\mathrm{n}}\right)R_{\mathrm{H}}=& \left[{\int }_0^{t_{\mathrm{n}}}1{\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}1{\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right]\hfill \\ \hfill & +\left[{\int }_0^{t_{\mathrm{n}}}(W_1(t_{\mathrm{n}}-\tau )-1){\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}(W_2\left(\tau \right)-1){\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right]\hfill \\ \hfill =& \left(p_1\left(t_{\mathrm{n}}\right)-p_2\left(t_{\mathrm{n}}\right)\right)\hfill \\ \hfill & +\left[{\int }_0^{t_{\mathrm{n}}}(W_1(t_{\mathrm{n}}-\tau )-1){\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}(W_2\left(\tau \right)-1){\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right]\hfill \\ \hfill =& \Delta p+\left[{\int }_0^{t_{\mathrm{n}}}(W_1(t_{\mathrm{n}}-\tau )-1){\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}(W_2\left(\tau \right)-1){\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right]\hfill \end{array}$$

(50)

The resulting equation now contains a stationary term that resembles the Hagen–Poiseuille law $(Q{R}_H=\Delta p)$ alongside dynamic terms for each pressure function.

It is convenient to define the dynamic components of $W_1$ and $W_2$ as distinct dynamic weighting functions:

$$W_{1,\mathrm{dyn}}^{\prime }=W_1-1,\hspace{1em}\hspace{1em}\hspace{1em}{W}_{2,\mathrm{dyn}}^{\prime }=W_2-1.$$

(51)

This leads to the expression:

$$Q_2\left(t_{\mathrm{n}}\right)R_{\mathrm{H}}=\Delta p+\left[{\int }_0^{t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{\prime }(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}{W}_{2,\mathrm{dyn}}^{\prime }(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right].$$

(52)

To clarify the impact of the dynamic components of the weighting functions, the negative signs from the summands of $W_{1,\mathrm{dyn}}^{\prime }$ and $W_{2,\mathrm{dyn}}^{\prime }$ are factored out and placed in front of the brackets:

$$Q_2\left(t_{\mathrm{n}}\right)R_{\mathrm{H}}=\Delta p-\left[{\int }_0^{t_{\mathrm{n}}}{W}_1,\mathrm{dyn}(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial p_1\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau -{\int }_0^{t_{\mathrm{n}}}{W}_2,\mathrm{dyn}(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial p_2\left(\tau \right)}{\partial t_{\mathrm{n}}}}d\tau \right].$$

(53)

With:

$$W_1,\mathrm{dyn}=-W_{1,\mathrm{dyn}}^{\prime },\hspace{1em}\hspace{1em}\hspace{1em}{W}_2,\mathrm{dyn}=-W_{2,\mathrm{dyn}}^{\prime }.$$

(54)

This outcome is reasonable since the Hagen–Poiseuille law describes an idealized flow condition that the equation aims to represent, specifically, incompressible and steady flow rather than compressible unsteady flow. As HP represents an idealization, the actual flow rate must increase over time to match the value dictated by HP due to the effects of viscous friction.

3.1. Application of the TVB Method to Obtain the Soft Sensor

From a computational perspective, performing a full convolution of functions over the entire time range is impractical for use in a soft sensor. Moreover, it is unnecessary, as the impact of pressure data diminishes over time, evidenced by the convergence of the weighting function to zero for large times. Vardy and Brown devised a method to compute convolutions efficiently between weighting functions and the time derivative of a discrete function [37]. This approach can similarly be applied to the convolution of a weighting function with a discrete function, as demonstrated in this paper. The weighting function is composed of exponential terms, and the discrete function is the difference between measured pressure signals.

The method’s efficiency lies in requiring only a few pressure data points while still delivering accurate results. The approach, termed the “Trikha–Vardy–Brown” (TVB) method, is an extension of a convolution simplification method initially introduced by Trikha [38]. In 1975, Trikha proposed a one-step scheme for the convolution integral, incorporating historical acceleration into parameters updated at each integration time step [38]. However, Trikha’s method had two limitations: first, the historical component lacked the generality needed for different types of flows, and second, it was only valid for very small time steps, which was impractical for real-world applications.

In 1983, Kagawa et al. improved Trikha’s technique by deriving a model with ten exponential terms, in contrast to the previously used three [39]. Although this approach required more memory, it enhanced Trikha’s method. By 1991, Suzuki et al. had developed a hybrid model, combining full convolution for short durations with Trikha’s algorithm for longer times, using a sum of five exponential terms [40]. Later, Schohl proposed an enhanced stepping formula to correct the first-time step limitations of Trikha’s method. Assuming uniform acceleration, Schohl’s approach also employed five exponential terms [41]. Similarly, Vítkovský et al. improved the first-step accuracy by introducing an empirical decay factor and developed advanced exponential sum weighting functions applicable to both laminar and turbulent flows [42].

In 2002, Ghidaoui and Mansoor presented an alternative to Trikha’s one-step method by using a weighting function derived from the work of Vardy and Brown [43], which they employed to study the evolution of specific integrals in unsteady flow [44]. More recently, Urbanowicz refined Schohl’s method by correcting the erratic recursive formula used to compute unsteady wall shear stresses [45]. His results were validated using the computationally expensive formula from Vardy and Brown in 2010 [46].

In this contribution, the TVB method was implemented in Matlab [47]. The derivation follows the authors’ approach but is applied to the current equations. The method is suitable for convolving a weighting function of the form $f\left(t\right)={\sum }_{i=1}^N{m}_i{e}^{-n_{i}t}$ with the time derivative of a second function expressed as $g\left(t\right)={\displaystyle \frac{\partial p\left(t\right)}{\partial t}}$.

Beginning with Equation (53), the flow rate can be divided into a stationary flow rate (according to the Law of Hagen and Poiseuille) and two dynamic flow rates corresponding to each pressure:

$$Q_2\left(t_{\mathrm{n}}\right)=Q_2,\mathrm{stat}\left(t_{\mathrm{n}}\right)-[Q_2,\mathrm{dyn},1\left(t_{\mathrm{n}}\right)-Q_2,\mathrm{dyn},2\left(t_{\mathrm{n}}\right)].$$

(55)

The convolution integral for the stationary volumetric flow rate can be readily solved and is already presented in Equation (53). Consequently, the TVB method must be applied solely to the dynamic volumetric flow rates. The following derivation will exemplify this method for $Q_{2,\mathrm{dyn},1}$, while the procedure for $Q_{2,\mathrm{dyn},2}$ would follow similarly.

As a starting point, the value of the integral evaluated up to the time $t_{\mathrm{n}}$ will be referred to as $Y_0\left(t_{\mathrm{n}}\right)$:

$$Q_{2,\mathrm{dyn},1}\left(t_{\mathrm{n}}\right)={\displaystyle \frac{1}{R_{\mathrm{H}}}}{\int }_0^{t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{*}(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau ={\displaystyle \frac{1}{R_{\mathrm{H}}}}{Y}_0\left(t_{\mathrm{n}}\right).$$

(56)

Next of interest is the time step coming after $t_{\mathrm{n}}$: $(t_{\mathrm{n}}+\Delta t_{\mathrm{n}})$. Inserting into the equation gives:

$$Q_{2,\mathrm{dyn},1}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}})={\displaystyle \frac{1}{R_{\mathrm{H}}}}{\int }_0^{t_{\mathrm{n}}+\Delta t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{*}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau .$$

(57)

The complete integral will be divided into two parts: the integral from 0 to $t_{\mathrm{n}}$, designated as $K_1$, and the integral from $t_{\mathrm{n}}$ to $t_{\mathrm{n}}+\Delta t_{\mathrm{n}}$, referred to as $K_2$.

$$\begin{array}{cc}\hfill Q_{2,\mathrm{dyn},1}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}})& ={\displaystyle \frac{1}{R_{\mathrm{H}}}}\left[K_1+K_2\right]\hfill \end{array}$$

(58)

With:

$$\begin{array}{cc}\hfill K_1& ={\int }_0^{t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{*}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau ,\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill K_2& ={\int }_{t_{\mathrm{n}}}^{t_{\mathrm{n}}+\Delta t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{*}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau .\hfill \end{array}$$

(60)

First, the second part of the integral, $K_2$, will be evaluated.

$K_2$ represents the change in volumetric flow rate during the interval from the current time step to the next one, which is $\Delta t_{\mathrm{n}}$ after the current step. Assuming a short time step, the change in the pressure gradient can be considered constant throughout the duration of $\Delta t_{\mathrm{n}}$. Thus, it follows that ${\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)=\mathrm{constant}$ for the entire period of $\Delta t_{\mathrm{n}}$. This allows the gradient to be factored out of the integral as a constant.

The constant value of the pressure difference gradient is determined by averaging over the time step, which involves taking the difference between the starting and ending values during the interval $\Delta t$. Therefore:

$$\Delta \dot{p}={\left({\displaystyle \frac{\partial \Delta p\left(\tau \right)}{\partial t_{\mathrm{n}}}}\right)}_{averagedover\Delta t_{\mathrm{n}}}.$$

(61)

Inserting this assumption into $K_2$ gives:

$$K_2=\Delta \dot{p}{\int }_{t_{\mathrm{n}}}^{t_{\mathrm{n}}+\Delta t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{*}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}}-\tau )d\tau .$$

(62)

The weighting functions are recognized as a sum of weighted, falling exponential terms, which can be expressed as:

$$W_{dyn}\left(t_{\mathrm{n}}\right)=\sum _{i=1}^N{m}_i{e}^{-n_i{t}_{\mathrm{n}}}.$$

(63)

Substituting this definition into the integral, the equation can be easily solved by integrating the exponential terms.

$$\begin{array}{cc}\hfill K_2& =\Delta \dot{p}{\int }_{t_{\mathrm{n}}}^{t_{\mathrm{n}}+\Delta t_{\mathrm{n}}}\left[\sum _{i=1}^N{m}_i{e}^{-n_i{t}_{\mathrm{n}}}\right]d\tau ,\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill K_2& =\Delta \dot{p}\sum _{i=1}^N{\displaystyle \frac{m_i}{n_i}}\left[1-e^{-n_i(\Delta t_{\mathrm{n}})}\right].\hfill \end{array}$$

(65)

Next, the first part of the integral, denoted as $K_1$, needs to be evaluated:

$$K_1={\int }_0^{t_{\mathrm{n}}}{W}_{1,\mathrm{dyn}}^{*}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau .$$

(66)

In this context, the authors of the TVB method employed a Taylor expansion of the weighting function around the point $(t_{\mathrm{n}}-\tau )$. As a result, $K_1$ becomes:

$$\begin{array}{c}\hfill K_1\approx K_{1,app}={\int }_0^{t_{\mathrm{n}}}{W}_{dyn}(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau +\Delta t_{\mathrm{n}}{\int }_0^{t_{\mathrm{n}}}{\left.{\displaystyle \frac{\partial W_{dyn}}{\partial t_{\mathrm{n}}}}\right|}_{(t_{\mathrm{n}}-\tau )}{\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau \\ \hfill +{\displaystyle \frac{{(\Delta t_{\mathrm{n}})}^2}{2!}}{\int }_0^{t_{\mathrm{n}}}{\left.{\displaystyle \frac{{\partial }^2{W}_{dyn}}{\partial t_{\mathrm{n}}^2}}\right|}_{(t_{\mathrm{n}}-\tau )}{\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau +\dots .\end{array}$$

(67)

It is important to note that the subsequent derivatives of the weighting function take the form of:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial W_{dyn}}{\partial t_{\mathrm{n}}}}& =-n_i{m}_i{e}^{-n_i{t}_{\mathrm{n}}},\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^2{W}_{dyn}}{\partial t_{\mathrm{n}}^2}}& ={n_i}^2{m}_i{e}^{-n_i{t}_{\mathrm{n}}},\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^3{W}_{dyn}}{\partial t_{\mathrm{n}}^3}}& =-{n_i}^3{m}_i{e}^{-n_i{t}_{\mathrm{n}}}.\hfill \end{array}$$

(70)

Due to this repetitive structure, the general derivative (of order k, where $k\in \mathbb{N}$) can be expressed as:

$${\displaystyle \frac{{\partial }^k{W}_{dyn}}{\partial t_{\mathrm{n}}^k}}=-n_i{\displaystyle \frac{{\partial }^{k-1}{W}_{dyn}}{\partial t_{\mathrm{n}}^{k-1}}}=C{\displaystyle \frac{{\partial }^{k-1}{W}_{dyn}}{\partial t_{\mathrm{n}}^{k-1}}}=\cdots =C^{k-1}{\displaystyle \frac{\partial W_{dyn}}{\partial t_{\mathrm{n}}}}.$$

(71)

Inserting these definitions of the derivatives of the weighting functions into the Taylor expansion yields:

$$K_{1,app}={\int }_0^{t_{\mathrm{n}}}{W}_{dyn}(t_{\mathrm{n}}-\tau ){\displaystyle \frac{\partial \Delta p}{\partial t_{\mathrm{n}}}}\left(\tau \right)d\tau \left[1+C\Delta t_{\mathrm{n}}+{\displaystyle \frac{{(C\Delta t_{\mathrm{n}})}^2}{2!}}+\dots \right].$$

(72)

Additionally, it is noteworthy that the first term in this equation is $Y_0\left(t_{\mathrm{n}}\right)$, as defined earlier.

$$K_{1,app}=Y_0\left(t_{\mathrm{n}}\right)\left[1+C\Delta t_{\mathrm{n}}+{\displaystyle \frac{{(C\Delta t_{\mathrm{n}})}^2}{2!}}+\dots \right].$$

(73)

Interestingly, the terms within the brackets following $Y_0\left(t_{\mathrm{n}}\right)$ resemble the series approximation of an exponential term with a factor C in the exponent. Consequently, the bracket can be rewritten as an exponential term instead (in this case, $C=n_i$):

$$K_{1,app}=Y_0\left(t_{\mathrm{n}}\right)e^{C\Delta t_{\mathrm{n}}}=Y_0\left(t_{\mathrm{n}}\right)e^{-n_i\Delta t_{\mathrm{n}}}.$$

(74)

Both integral parts are now simplified and can be substituted back into Equation (58):

$$\begin{array}{cc}\hfill Q_{2,\mathrm{dyn},1}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}})& ={\displaystyle \frac{1}{R_{\mathrm{H}}}}\left[K_1+K_2\right],\hfill \\ \hfill Q_{2,\mathrm{dyn},1}(t_{\mathrm{n}}+\Delta t_{\mathrm{n}})& \approx {\displaystyle \frac{1}{R_{\mathrm{H}}}}\sum _{i=1}^N\left[Y_{0,i}\left(t_{\mathrm{n}}\right)e^{-n_i\Delta t_{\mathrm{n}}}+\Delta \dot{p}{\displaystyle \frac{m_i}{n_i}}(1-e^{-n_i\Delta t_{\mathrm{n}}})\right].\hfill \end{array}$$

(75)

As a result, a one-step integration method for calculating the volumetric flow rate has been established. For each step, the flow rate from the previous time step is multiplied by a decaying exponential term, effectively reducing its value. The current time step is then evaluated independently of any prior flow rate or pressure values and added to the diminished value from the last time step.

The significant advantage of this approach is that only the values from the previous time step need to be stored. In contrast, without this simplification, the entire convolution integral would need to be evaluated for each time step. Since the bounds of the convolution integral range from 0 to $t_{\mathrm{n}}$, the length of the integral would increase indefinitely as time progresses, which would make the application of the soft sensor impracticable. With the implemented method, the analytical model can be used offline and online on a real fluid power system to constantly measure the volumetric flow rate. Especially for online use, the one-step integration method and the analytical model’s simple computational operations, like adding exponential functions, underline the soft sensor’s ability to perform real-time computation.

3.2. Test Cases

The following test cases are investigated to obtain a range of validity for the derived method. The volumetric flow rate computed by the soft sensor is validated against the results provided by the software DSHplus, which solves the partial differential equations for the conservation of momentum and mass with a finite difference scheme [48]. Furthermore, a soft sensor only utilizing the law of HP, representing a commonly used approach for pressure-based soft sensors, is investigated to underline the benefit of incorporating the unsteady flow rate behavior. The simulation model in DSHplus consists of a pipe with varying pressure inlet boundary conditions. Four different types of pressure signals are investigated: sine, sum of sines, sawtooth, and step signals. In

Table 1. Parameter set in the simulation for each test case.

Name	Variable	Unit	Value
Kinematic Viscosity	$\nu$	m2/s	$46\times {10}^{-6}$
Fluid Density	$\rho$	kg/m3	860
Bulk Modulus	K	Pa	$14\times {10}^9$
Diameter of the Pipe	D	m	$14\times {10}^{-3}$
Length of the Pipe	L	m	$1.5$
Hydraulic Resistance	$R_{\mathrm{H}}$	Pa/(m3/s)	$6.294\times {10}^7$

, the overall simulation parameters are displayed.

The comparison is performed for a pipe with standard hydraulic oil HLP46 and a set diameter and length. Based on these values, the hydraulic resistance can be computed. The change between the different presented simulations lies in the investigated pressure boundary conditions.

Table 2. Pressure conditions set for each test case.

Test Case	Initial Conditions ${\mathit{p}}_{\mathrm{IC},1,2}$ [bar]	Boundary Conditions $p_{\mathrm{BC},1}$ [bar]	Frequency f [Hz]
Sine	50	$50\pm 0.1$	1, 100, 1000
Sum of Sines	50	$50\pm [0.1,0.1,0.1,0.1]$	[$1,10,20,40$], [$10,100,200,400$]
Sawtooth	50	$50+0.2$	10
Step	50	$50+0.5$	N/A

displays the various types of investigated pressure boundary conditions. The pressure outlet is fixed; thus, only the pressure inlet varies.

The first investigated pressure signals are sine waves at different frequencies (1, 100, and 1000 Hz) to validate the analytical model for low and high frequencies. When a pump delivers oil into a system, it inevitably generates pressure pulsations. These pulsations occur at a fundamental frequency determined by the number of chambers in the pump and its rotational speed. In addition, integer multiples of this fundamental frequency are also excited. To replicate this behavior, the pressure boundary condition for the simulations is defined as a sum of sine waves oscillating at 1 Hz, 10 Hz, 20 Hz, and 40 Hz, as well as 10 Hz, 100 Hz, 200 Hz, and 400 Hz. Therefore, the second test case consists of two signals composed of sums of sines. Thirdly, a sawtooth function was examined as an example of a highly instationary boundary condition. Lastly, a step function was analyzed because this work focuses on the equation’s ability to calculate transient flow rates in a pipe. To validate its performance under such conditions, a step function increase in the pressure boundary condition will be examined. A step function represents the most unsteady case possible, theoretically involving a Fourier approximation as a sum of frequencies up to infinity.

4. Results

This section presents the agreement between the distributed-parameter simulation’s determined volumetric flow rate, the novel soft sensor, and a soft sensor based on the law of HP.

Table 3. Mean errors and standard deviation of the soft sensor for the investigated test cases.

Test Case	Frequency f [Hz]	Mean Error	Standard Deviation
Sine (Figure 7)	1	$0.11\%$	$0.076\%$
Sine (Figure 8)	100	$0.0084\%$	$0.38\%$
Sine (Figure 9)	1000	$0.14\%$	$1.3\%$
Sum of Sines (Figure 10)	[$1,10,20,40$]	$0.0049\%$	$0.067\%$
Sum of Sines (Figure 11)	[$10,100,200,400$]	$0.034\%$	$1.3\%$
Sawtooth (Figure 12)	10	$0.043\%$	$0.11\%$
Step (Figure 13)	N/A	$0.0051\%$	$0.059\%$

shows the performance of the soft sensor, characterized by the mean error and standard deviation of the volumetric flow rate calculated by the distributed-parameter simulation and the soft sensor.

Table 4. Mean errors and standard deviation of a sensor based on the law of HP for the investigated test cases.

Test Case	Frequency f [Hz]	Mean Error	Standard Deviation
Sine (Figure 7)	1	$29\%$	$52\%$
Sine (Figure 8)	100	$2700\%$	$1300\%$
Sine (Figure 9)	1000	$3800\%$	$1900\%$
Sum of Sines (Figure 10)	[$1,10,20,40$]	$57\%$	$91\%$
Sum of Sines (Figure 11)	[$10,100,200,400$]	$320\%$	$350\%$
Sawtooth (Figure 12)	10	$36\%$	$57\%$
Step (Figure 13)	N/A	$63\%$	$30\%$

displays the same errors for the HP-based soft sensor.

4.1. Sine Wave Pressure Signals

The first test cases consider sinusoidally varying pressure inlets. Frequencies of 1 Hz, 100 Hz, and 1000 Hz were selected for the sine waves in Figure 7, Figure 8 and Figure 9. HP is only shown in Figure 7 for the 1 Hz frequency, as the deviation increases with higher frequencies. As a result, scaling the plots becomes necessary, making it more challenging to visually compare the actual soft sensor output with the simulation model. The errors of the HP-based soft sensor are summarized in Table 4 and range from 29% up to 3800% for the lowest to the highest frequency.

The analytical model agrees well with the simulation data in all cases, with a mean error of less than 0.15%.

4.2. Sums of Sine Wave Pressure Signals

In the first following figure (Figure 10), the pressure boundary condition for the plots in this section is defined as the sum of sine waves oscillating at 1 Hz, 10 Hz, 20 Hz, and 40 Hz. In the following figure (Figure 11), the pressure boundary condition is set as the sum of sine waves at 10 Hz, 100 Hz, 200 Hz, and 400 Hz. Both figures demonstrate that the analytical model aligns well with the simulation data. The errors of the soft sensor are $0.0049$% and $0.034$%, respectively. The HP-based sensor exhibits errors of 57% and 320%, respectively.

4.3. Sawtooth Pressure Signals

This section presents the agreement between the soft sensor and the simulation data for the case of sawtooth pressure boundary conditions. The figure below (Figure 12) shows a 10 Hz sawtooth wave. The soft sensor demonstrates good agreement with the simulation data, with a mean error of $0.043$%, while the HP-based sensor exhibits an error of 36%. Additionally, the plot on the right in Figure 12 shows that the analytical model developed in this work also captures the compressible oscillations observed in the simulation data.

4.4. Step Function Pressure Signal

The last test case investigates a pressure step function. The two plots in Figure 13 display the same boundary condition at different zoom levels. The novel equation successfully captures the inductive behavior observed in the simulation data. Additionally, compressibility effects are evident in the step changes of the simulation data, and the equation reflects these steps. As shown in the right plot of Figure 13, the equation matches these steps through a Fourier approximation, oscillating around the constant values of the simulation data. The mean error of the soft sensor is $0.0051$%, while the HP-based sensor shows an error of 63%.

5. Discussion

The analytical soft sensor analysis based on two pressure signals demonstrated its promising ability to determine high-frequency flow rates. The model accurately computes the various flow rates using pressure signal information. Due to its analytical nature, the computation is fast and not delayed by iterative numerical calculations, unlike the distributed-parameter simulation model, DSHplus. Nevertheless, the results obtained from the soft sensor show good agreement with the numerically computed values from the simulation model. A soft sensor based on the HP law exhibits high percentile errors (up to 3800%) due to its neglect of transient effects.

The analytical model exhibits a small mean error of less than $0.15$% for all test cases, with most mean errors below $0.045$%. The standard deviation of the error further highlights the good performance of the soft sensor, with values ranging from $1.3$% to $0.59$%. The analytical model can accurately compute the correct flow rate for sine waves up to 1000Hz, with a mean error of $0.14$% and a standard deviation of $1.3$%. This high accuracy for transient pulsations represents a promising application of the soft sensor for the investigation of axial piston pumps, which exhibit pressure ripples in the range of 1000 Hz. Furthermore, the computation of volumetric flow rate based on the pressure difference allows for power monitoring of a pump, enabling the characterization of pulsation and a pump’s performance.

Regarding the sum of sine waves, the model provides precise values for several overlaying sine waves. For both cases, the mean error is below $0.04$%, with a mean error of $0.0049$% for the sum of sines with 1, 10, 20, and 40 Hz, and $0.043$% for the sum of sines with 10, 100, 200, and 400 Hz. The standard deviation values of $0.067$% and $1.3$% further underline the high accuracy of the soft sensor. These results demonstrate the soft sensor’s ability to capture different levels of transient effects in a system, enabling the detailed investigation of a hydraulic system’s dynamic behavior. Furthermore, the knowledge of pressure and volumetric flow rates for arbitrary systems allows for calculating hydraulic power for hydraulic systems with varying behavior. This is further highlighted by the good results of the sawtooth test case, with a mean error of $0.0043$% and a standard deviation of $0.11$%.

Regarding the step function, the soft sensor exhibits good agreement with the distributed-parameter simulation, with a mean error of $0.0051$% and a standard deviation of the error of $0.059$%. This accuracy allows the sensor to be applied to systems where sudden pressure changes, e.g., due to valve closing, resulting in a water hammer, are investigated. Notably, the step function contains a bandwidth of frequencies, which the analytical soft sensor accurately captures. Significantly, the investigation of the flow rate directly after the pressure step shows the model’s ability to capture compression phenomena, as demonstrated by the stair-like increase in the volumetric flow rate.

The novel soft sensor exhibits good agreement across all investigated test cases, considering various frequency ranges and transient pressure signals. In comparison, the HP-based soft sensor exhibits errors of up to 3800% and a mean error of 29% for the lowest frequency sine wave. This highlights the analytical model’s advancement over the typically applied HP law.

The demonstrated precision at high frequencies and compressible flow conditions underscores the soft sensor’s potential for deployment in various industrial applications, such as condition monitoring and control, where real-time and accurate computation of crucial system parameters is essential. Additionally, the sensor’s ability to track pressure and flow rate enables the calculation of hydraulic power, which can be utilized for predictive maintenance, offering more profound insights into system performance. In contrast, traditional numerical solvers often provide either real-time solutions or high accuracy, but not both simultaneously, due to the limitations inherent in numerical computation.

Future research will focus on extending the applicability of this soft sensor to turbulent flow conditions and validating its performance on a real-world test rig [49]. A key aspect of this research will be examining the sensor’s limitations through its application to practical setups. The soft sensor’s limitations arise from several factors. The underlying equation can be expanded to capture system behavior across any frequency range relevant to hydraulic systems. However, in real-world applications, the noise level of the measurement system sets the lower limit of detectable pressure differentials between the sensors, while the software’s numerical accuracy defines the minimum precision in simulations. By employing this soft sensor, significant advancements can be achieved in transient flow rate measurements.

6. Conclusions

This contribution demonstrates the capability of an analytical soft sensor to determine compressible volumetric flow rates without an actual flow rate sensor and only two pressure signals in a pipe. It begins with an introduction to current advances in the computation of transient flow rates, followed by a section on the fundamentals of fluid mechanics and the Laplace techniques to solve the governing equations to obtain an analytical model. Subsequently, the analytical model is presented, explaining the assumptions made and the inverse Laplace transformation. This transformation converts the solved equation from the Laplace domain to the time domain, enabling its use for flow rate computation. The developed model is investigated in the Laplace domain to obtain a deeper understanding, especially regarding its modeling of compressibility phenomena.

The findings of this research present a significant advancement in soft sensors for transient flow rate measurement. The soft sensor accurately computed the volumetric flow rate and was compared to numerically obtained results from the 1-D hydraulic simulation software DSHplus. Four test cases were investigated: sine waves, sums of sine waves, sawtooth waves, and a step function to validate the soft sensor. The soft sensor agreed well with each test case’s numerically computed values. The soft sensor was able to determine even high-frequency flow rates of up to 1000 Hz and capture the compression phenomena visible in the step case.