An implicit system of delay differential algebraic equations from hydrodynamics

. Direct spring operated pressure relief valves connected to a constantly charged vessel and a downstream pipe have a complex dynamics. The vessel-valve subsystem is described with an autonomous system of ordinary differential equations, while the presence of the pipe adds two partial differential equations to the mathematical model. The partial differential equations are transformed to a delay algebraic equation coupled to the ordinary differential equations. Due to a square root nonlinearity, the system is implicit. The linearized system can be transformed to a standard system of neutral delay differential equations (NDDEs) having more elaborated literature than the delay algebraic equations. First, the different forms of the mathematical model are presented, then the transformation of the linearized system is conducted. The paper aims at introducing this unusual mathematical model of an engineering system and inducing research focusing on the methodology to carry out bifurcation analysis in implicit NDDEs.


Introduction
The mathematical model in question describes the dynamics of a simplified mechanical model consisting of a vessel charged by a constant fluid flow rate, a direct spring operated pressure relief valve and a pipe delivering the fluid to the atmospheric pressure.Previous studies are presented in [7] and [8] mainly focus on the engineering aspects of the valve dynamics.The model construction and the derivation of the dimensionless mathematical model can be found in [7] in detail, but the various simplified forms of the mathematical model are summarized here with compressed parameters.
The dimensionless mathematical model of the vessel-valve subsystem is an autonomous system of ordinary differential equations (ODEs): ẏ1 (t) = y 2 (t) , (1.1) ẏ2 (t) = −2ζy 2 (t) − (y 1 (t) + δ) + y 3 (t) − y out , (1.2) ẏ3 (t) = β q − y 1 (t) y 3 (t) − y out , (1.3) with the state variables y 1,2,3 ∈ R, where dot refers to the derivative with respect to the dimensionless time t.The appearing parameters ζ, δ, β, q are positive according their physical meanings.The backpressure y out is usually the atmospheric backpressure that is zero for the vessel-valve subsystem.In the followings, we investigate the dynamics of the open valve for which the conditions y 1 > 0 and y 3 > 0 fulfil.The system can be linearized around its trivial solution y * 1,2,3 , which is given as the only physically relevant non-negative solution of the following algebraic equations depending on the parameters q and δ: For small perturbations η 1,2,3 ∈ R around the equilibrium, the linearised system is: The Navier-Stokes and continuity equations of the fluid flow through the pipe serve two dimensionless partial differential equations: where the new coordinate x denotes the dimensionless space coordinate along the pipe such that x = 1 corresponds to the end of the pipe, while y 4 is the gauge pressure distribution and y 5 is the fluid velocity distribution along the pipe.The parameters Γ 1 and Γ 2 are also positive and can be calculated from the physical parameters of the system and the fluid.The time delay τ connects the parameters Γ 1 and Γ 2 .The time delay is the time of the wave propagation along the pipe until the backward wave interacts with the valve dynamics.In this case, the valve backpressure is y out = y 4 (0, t) in the ODEs (1.1)-(1.3).The corresponding boundary conditions are: ) The first boundary condition prescribes that the fluid delivers to the zero overpressure.The second one expresses that the mass flow rate through the valve must be equal to the inlet mass flow rate to the pipe.The parameter Γ 3 > 0 can also be expressed from the physical parameters of the system and the flow.
The solutions of the partial differential equations in (1.7), (1.8) can be formulated in the so-called travelling wave solution form: where f : R → R, g : R → R are unknown functions.The travelling wave solutions can be substituted back into the boundary conditions.From the boundary condition (1.9), the connection is obtained, which can be further substituted into the boundary condition (1.10) leading to a nonlinear delay algebraic equation: where ϕ > 0 is called as pipe inlet parameter, and ϕ = Γ 1 Γ 3 τ/2.Note, that the differential algebraic equation is only one possible terminology for this kind of equation, see [3], while it can also be called as nonlinear difference equation, or shift map [4], or neutral renewal equation [1,2].Finally, the whole system of delay differential algebraic equations (DDAEs) can be derived by considering the coupling between the ODEs and the algebraic equation with This system is implicit, because of the square root nonlinearity in (1.18), thus, the algebraic equation cannot be arranged into standard explicit NDDE form [5,11]: where y : R → R 3 , u : R 3 → R 3 and v : R The derived implicit DDAE model is unusual in the field of engineering dynamics where the mathematical background is not yet elaborated.As it was referred to [1][2][3][4], even the terminology is not uniform and the numerical methods rarely handle neutral kind systems, especially implicit and/or algebraic ones.Paper [7] contains a thorough linear stability analysis, mode shape calculations and a numerical validation of the time-delay model of the vesselvalve-pipe system.The results presented in [8] mainly explain the nonlinear aspects of the safe operation of the pressure relief valves mounted on a vessel.Based on the results of [7], it is also shown how to use the attached downstream pipe consciously to avoid the harmful valve vibrations.At the end, it also presents nonlinear analysis of the time delay model of the vessel-valve-pipe system by means of DDE Biftool via creating an approximate system of retarded delay differential equations.Although the studies in [7,8] provided valuable information from engineering point of view, both analytical and numerical techniques have to be developed for the reliable handling of implicit nonlinear delay differential algebraic equations (DDAEs), which, to best of our knowledge, is missing from the literature.

Transformation to neutral delay differential equations
The presence of the pipe added one algebraic equation to the model of the vessel-valve subsystem.The trivial solution y * 1,2,3 for the vessel valve system is not changed by the presence of the pipe, but additionally the trivial solution f * appears for the DDAEs (1.15)-(1.18),extending the algebraic system of equations for the trivial solution as follows: The linearized system can be derived for the small perturbations η 1,2,3,4 ∈ R around the trivial solution: )

.4)
The vessel-valve dynamics around the equilibrium is described in the space of η 1,2,3 in (1.4)-(1.6).These state variables have actual physical meaning, while the pipe adds one new equation with the state variable η 4 originated in the wave propagation described by the function f , which cannot be directly associated to a physical measure.The most convenient interpretation from engineering point of view would be to present the vessel-valve-pipe system in a linear form in the space of the vessel-valve subsystem, thus in a delayed differential equation form with the variables η 1,2,3 .This way, the vessel gauge pressure and the disk vibrations remain in the model, and the effect of the pipe would appear through the delayed terms only.Also, the equilibrium of the pipe is detached from the equilibrium of the vessel-valve subsystem, while the τ → 0 case must be the transition between the vessel-valve-pipe and vessel-valve model.Indeed, the linearized system of DDAEs (2.1)-(2.4)can be transformed to a system of neutral delay differential equations (NDDEs) by means of some algebraic manipulation.The goal is to present a standard system of NNDEs in the space of the same state variables as in case of the vessel-valve subsystem presented in (1.4)-(1.6).
First, let us rearrange (2.4) to obtain Then take the linear combination of equation (2.2) and its τ-delayed form: (2.6) It can be done similarly for (2.3): (2.7) For both equations, (2.5) and its τ-delayed form appears.After their substitution, the following two equations are obtained in place of (2.2)-(2.3): (2.9) Finally, the whole system of NDDEs is presented in a matrix form: where η = [η 1 , η 2 , η 3 ] T , while A, B, C ∈ R 3x3 are given as follows:

.13)
As conclusion, a linear system of NDDEs is derived in the presence of the pipe with the state variables of the vessel-valve subsystem.It has a straight transition to the linearaized system of ODEs describing the dynamics of the vassel-valve subsystem in (1.4)-(1.6)with the limit of τ → 0. This limit case leads to the system with the identity matrix I ∈ R 3 , which is, indeed, an equivalent system to (1.4)-(1.6).

Linear analysis of the DDAEs
The characteristic equation corresponding to the linearized system of NDDEs (2.1)-(2.4)can be derived as leading to the form of with the coefficients (2.17) Since we temporarily extended the interval of the initial data to [−2τ, 0] in the steps of derivation in (2.6) and (2.7), the term a 3 + b 3 e −λτ appears as a multiplicator of ∑ 3 i=0 (a i λ i + b i λ i e −λτ ).The 2τ-delayed term η 4 (t − 2τ) was eliminated in (2.8) and (2.9), so the final set of initial data is defined again on [−τ, 0] and the factorized characteristic equation does not contain exponential terms e −2λτ explicitly.The linear stability analysis of the system of DDAEs can be carried out by direct substitution of an exponential trial solution into (1.15)-(1.18),see [7], leading to the characteristic equation ∑ 3 i=0 (a i λ i + b i λ i e −λτ ) = 0 with the exact same coefficients as in (2.17); moreover, the roots of a 3 + b 3 e λτ = 0 assign the spectral abscissa of the essential spectrum [9].Because of the neutral kind of the system, infinitely many characteristic roots may appear on the right-hand side of the complex plane causing the so-called essential loss of stability.The condition guarantees that the spectral abscissa of the essential spectrum is located on the left side of the complex plane avoiding the possibility of the essential loss of stability.For the coefficients in (2.17), this condition always fulfils.The associated difference equation of the NDDE [9]

Stability charts
The linear stability analysis detects the possible Hopf bifurcations in the system.The stability properties of the trivial solution can be determined for ∑ 3 i=0 (a i λ i + b i λ i e −λτ ) = 0 through the so-called D-separation method [6,10].The critical roots are the poor imaginary complex conjugate ±iω, where ω ∈ R + is the angular frequency of the possibly emerging periodic solution, which corresponds to the self-excited vibration of the valve.In this case, two solutions exist for ω depending on the parameters ζ , δ , q.The stability boundaries β cr (τ; ω) can be determined for both frequencies, depending on the time delay τ; for details and the closed form solutions, see [7].As a practically useful result, stability charts like the one in Fig. 2.1 can be presented to characterise the effects of the various sizes of the pipe and the vessel on valve stability.

Figure 2 . 1 :
Figure 2.1: Example stability chart for q = 6 , δ = 3 , ζ = 0.39 , ϕ = 48.2.The gray domains are stable, the white domains are unstable.The two types of lines represent the stability boundaries with the two distinct vibration frequencies of the self-excited vibrations.
+ b i e −λτ ) = 0. Consequently, in case of the transformations resulting (2.10), the characteristic roots of the associate difference equation are false roots of the characteristic equation itself.If the condition (2.18) fulfils, it leads to the same stability properties of the linear system of DDAEs (1.15)-(1.18)and NDDEs (2.10).