# Timestep-dependent element interpolation functions for axial FSI analysis in pipelines

Igor Orynyak, Danylo Katerynych, Roman Mazuryk

► To cite this version:

Igor Orynyak, Danylo Katerynych, Roman Mazuryk. Timestep-dependent element interpolation functions for axial FSI analysis in pipelines. 2026. hal-05478034

HAL Id: hal-05478034

https://hal.science/hal-05478034v1

Preprint submitted on 26 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Timestep-dependent element interpolation functions for axial FSI analysis in pipelines

Igor Orynyak¹, Danylo Katerynych², Roman Mazuryk³

¹ Professor, Department of Applied Mathematics at National Technical University of Ukraine “Igor Sikorsky Kyiv Polytechnic Institute”, Peremohy str, 37, Kyiv 03056, Ukraine. ORCID: https://orcid.org/0000-0003-4529-0235. Email:igor_orinyak@yahoo.com.

² Master student, Department of Applied Mathematics at the National Technical University of Ukraine "Igor Sikorsky Kyiv Polytechnic Institute", Peremohy str, 37, Kyiv 03056, Ukraine. ORCID: https://orcid.org/0009-0002-1267-7425. Email: danilokaterinich@gmail.com.

³ Assistant Lecturer, Applied Mathematics Department at the National Technical University of Ukraine "Igor Sikorsky Kyiv Polytechnic Institute", Peremohy str, 37, Kyiv 03056, Ukraine. ORCID: https://orcid.org/0000-0003-4309-824X. Email: r.mazuryk.ua@gmail.com.

**Abstract.** This work is a continuation of a series of works by the authors, where dynamic problems (time derivatives) are first modified using a central difference scheme. Then the main function divided by the time step becomes part of a differential equation for an element, that is solved analytically (exactly in the case of one-dimensional problems, and approximately in the case of multidimensional problems). For Skalak's problem, this reduces to an ordinary fourth-order differential equation with a right-hand where the time step is a simple multiplier for the unknown function. First, the eigenfunctions of the homogeneous equation are found, which are represented as generalized Krylov functions with remarkable properties regarding their values and the values of their derivatives at the initial point of each element. All solutions (general homogeneous and partial inhomogeneous) are represented in an original way through these Krylov functions (introduced here for the first time) and are given in a form convenient for application of the transfer matrix method, TMM. The properties of the solution are presented in detail for the well-known Delft Hydraulics Benchmark Problem A. To verify the accuracy, we also implement the method of separation of variables (Fourier). A comparison of the results showed the high efficiency of the proposed solution scheme, where it is possible to understand in advance what accuracy can be achieved for a given time moment and the selected time step.

**Keywords:** Skalak problem, transfer matrix method, generalized Krylov functions, ordinary differential equation of the 4th order, separation of variables method, numerical amplitude dissipation

## Highlights

1. Implicit central difference scheme leads to ordinary 4th order differential equation.

2. Generalized Krylov functions are introduced as eigenfunctions of homogeneous equation.

3. Procedure of solution is organized according to transfer matrix method.

4. Method of separation of variables is implemented too for validation of the results.

5. Accuracy depends on the ratio of time step to square of smallest natural frequency.
# 1. Introduction

The problem of axial interaction between a straight thin-walled pipe and the fluid inside it was formulated back in 1955 and has since been named after Skalak [1]. This is a one-dimensional problem in the axial direction, which reduces to fourth-order differential equations both in the axial coordinate and in time, and in fact takes into account the influence of Poisson's ratio on the coupled dynamic behavior. If it were zero, we would have two separate problems for the pipe and the liquid, each of which would be a second-order problem. Methods for solving Skalak's problem can be roughly divided into two main categories: frequency domain analysis and time domain analysis. Frequency domain analysis primarily investigates the vibration characteristics, stability, and resonance, whereas time domain analysis focuses more on transient responses, with emphasis on transient pressures [2].

Consider the first category – concerning vibrations and resonance. It is assumed that the process is periodic in time with frequency $ω$. Then the general formulation reduces to a fourth-order differential equation in $x$. Solutions to Skalak's equation are widely represented in the literature, starting with the work of Wiggert et al. [2], where it is presented by a transfer matrix, i.e., a matrix that links the four parameters at the beginning of the segment (element) with the corresponding parameters at the end of the segment. Wiggert's solution involves as many as 4 matrices. The first is for fluid-solid interaction, FSI, problem, 4 parameters; the second is for vibration in the transverse direction, 4 parameters, with water considered as a distributed attached mass; the third is similar to previous one and relates for another transverse direction, 4 parameters; the fourth is for torsion around the axis, 2 parameters. The solution [2] is often referred to as the 14-equations model to emphasize that in addition to the 12 mechanical parameters, two additional liquid parameters are applied: pressure and displacement. Similar solutions were summarized in dissertations [3,4]. Tanterelli [3] considered springs, accumulators, and the extra mass, using linear lumped impedances, friction coupling, and discussed various boundary conditions at pipe ends; whereas in [4] the flexibility factor is employed to modify the bending stiffness of curved pipes.

By applying the separation of variables method, the obtained frequencies and vibration modes can be used to reproduce the entire dynamic process. However, in [5,6], it was concluded that the time-domain solutions derived from frequency-domain results are difficult and impractical, as they require a very large number of expansion terms.

Nevertheless, for large pipeline systems, the task often consists in determining the conditions of resonance or the response of the system to known dynamic actions. In these cases, it is appropriate to use harmonic analysis. For typical segments of the pipeline, analytical solutions are obtained that relate the main parameters of the task at any point of the segment with its beginning. A feature of such solutions is that the unknown frequency of oscillation, $ω$, is an unknown constant within the analytical solution, for example, $sin(ωax)$, where $α$ is a known constant and $x$ is the length coordinate. Such solutions are used in transfer matrix method, TMM, and are extremely effective for
analyzing complex branched pipelines with additional supports and complex constraints [7-10]. These approaches are often referred to as spectral element methods to emphasize the use of the combination of $ωx$ as an independent variable [11-13]. Other, less popular examples of frequency domain analysis are Component Synthesis Method [14] and Laplace transform method [15].

It should be noted that the TMM was initiated in Pestel's work [16], and the technique is discussed in detail there, in particular, - the procedures for breaking down the system into elements (which increases the number of unknowns) and eliminating (multiplying) matrices of neighboring elements (which reduces the number of unknowns). In piping systems, the TMM was introduced into the field of FSI by Chaudhry [17], who defined the transfer matrix as a matrix relating two state vectors. He presented three types of them. The field transfer matrix, FTM, relates state vectors at the two borders of a pipe and gives the solution of corresponding differential equations. The point transfer matrix, PTM, relates state vectors at the borders of two segments (elements), in particular, it is often used to take into account the curvature of the pipeline axis or the supports and lumped masses. The overall (global) transfer matrix links two ends of the longer segments of the piping system, which is obtained as a multiplication of the field and point transfer matrices. It is interesting to note that in [2], the FTM and PTM were already separated. The first is considered as a matrix of solutions to a system of equations for a segment, and the second is responsible for the rotation of the local coordinate system, which is important when considering a large pipeline system with curved segments. It should be noted that TMM is very popular because of its simplicity, numerical efficiency, and ease of implementation for solving problems of statics and vibrations of pipeline systems. Thus, it has been widely used for problems of statics of complex pipeline systems [18,19] or complex spatial beams [20,21]. Note that static solutions of the TMM do not consider the interaction of liquid and metal, since pressure is usually considered a constant value. However, the advantage of static methods [18,19] is that analytical solutions are very often obtained in curvilinear coordinates, where instead of four sets of equations solved by TMM, there are two – deformation in the plane – 6 equations and deformation out of the plane – 6 equations. This allows for an increase in the length of curved segments. Such analytical solutions, which are convenient for use in TMM, are rarely used for the analysis of vibrations of curved beams, which is explained by the complexity of their derivation and the justification of their convergence.

The second category of solutions concerns time domain analysis. According to reviews [5, 22], the most common discretization methods used to solve the governing differential equations are the Method of Characteristics (MOC), the Finite Element Method (FEM), the Finite Difference Method (FDM), or the Finite Volume Method (FVM). These methods can be used either individually or in combination (hybrid) with respect to the different (pipe or liquid) degrees of freedom [5].

The Finite Element Method is most commonly used in structural dynamics. However, it is rarely used for this particular problem. Examples include works [23, 24]. However, for steep wave phenomena, the FEM is less
appropriate than the MOC, because due to the partitioning of the pipe in FEM, only a limited number of fundamental frequencies can be described [24]. Therefore, in combined approaches, where it is necessary to reproduce the hydraulic wave front more accurately, the MOC is often used for the liquid and the FEM for the pipe [25, 26].

An FDM approach was suggested for treatment of wave phenomena in the curved pipe (so-called eight equation problem, where the axial and the bending deformation models are coupled) [27]. FDM (for the liquid) was formerly used mostly in combination with FEM (for the structure) [28,29]. FVM is mostly used for analysis of single water hammer phenomena [30], while there are only few examples of its application to the Scalak’s problem [31].

MOC is the most popular and efficient time domain approach for solving the water hammer related steep wave transients. As was noted in [32] “the method of characteristics (MOC) has undoubtedly been a wonderful bonus for the pressure surge community”. The background principle of the MOC is to transform partial differential equations into the ordinary ones along the characteristic lines related time with distance. For the first time, the Scalak problem was solved using the MOC approach in [33].

Early works suffer from numerical damping and phase error due to problems with meshing the distance-time plane and time march from the initial conditions. The main reason is that there are two different wave velocities: stress wave speed for the structure for the pipe and pressure wave speed for the liquid, the ratio of which is usually not an integer. So, it is hard to find the intersection of the characteristic line and the time grid, and three different techniques of grid processing are proposed [34]: space-line interpolation (SLI) [35, 36], time-line interpolation (TLI) [37,38], and wave speed adjustment (WSA) [39]. Besides the aforementioned single methods, the combined MOC-MOC method for both the liquid and the structure was shown to have high efficiency in [40]. Despite being the most accurate and effective, MOC is prone to multiple interpolation errors related to characteristic lines in managing complex boundary conditions [41].

Since the time of the review [5], the discrete time transfer matrix method (DTTMM) for solving the FSI fourteen-equations model in the time domain was presented in 2021 [42]. Although it looks a little cumbersome, it contains many attractive features. First, it is organized completely according to TMM, so the efficient procedure for solving any complex pipeline system can be constructed easily. Second, it is a direct time integration method [43], but contrary to the conventional FEM procedure, it does not contain the mass matrix. The time step is included in the stiffness matrix analogously as the frequency is used in spectral methods [11-13]. The idea is to consider the partial differential equations in matrix form and present the acceleration as the difference between the velocity at the current time and the previous time divided by the time step. Then the system of partial differential equations became the ordinary one with the right side dependent on the solution at the previous times. Then it was formally solved as the sum of homogeneous and non-homogeneous solutions. The solution contains no explicit expressions for eigenfunctions, nor was any analysis of them or the structure of the right-side terms performed. No hints about the
meshing and stability are given. But as it is known, the eigenfunctions, at least for flexural modes, might contain
exponential functions that are prone to instability, thus careful meshing is required. A similar problem in the
application of conventional TMM (for vibration analysis) was mentioned by many authors [4,7,8], where the
exponents related to different types (lateral, axial) of vibration may differ drastically. So, the problems with instability
were perhaps the reason that in subsequent analysis of liquid-structure interaction in pipelines the authors of DTTMM
[42] have returned to the well-known Laplace transform TMM [44, 45].

Nevertheless, the idea of [42] is very promising and very akin to our last investigations on direct time
integration of various problems, where the time step is treated as a factor in a field transfer matrix. In [46] we solve
the 1D heat distribution task; in [47] – the axial wave propagation in a rod (water hammer can be treated similarly);
in [48] – the wave propagation in the gas. Their main features are as follows: a) representing the second derivative of
time by the central difference method and obtaining ordinary differential equations with the right-hand side depending
on the previous solutions in time; b) getting the eigenfunctions for ODE and relating by TMM the state vector of any
point of each segment with the state vector at the first point; c) presenting the known state vector at two previous
moments of time as a polynomial of order k, and obtaining the zeroth particular solution (equal to zero at the first
point of the segment) of the nonhomogeneous equation; d) compiling the matrix and obtaining the continuous (with
respect to the length coordinate x) solution at the present time. The eigenfunctions have the form of $\exp(\alpha \frac{x}{\Delta t})$,
where $\Delta t$ is the time step and $\alpha$ is a known constant. Its accuracy and stability are governed by two factors. Stability
is guaranteed if $\alpha \frac{L}{\Delta t} < Z$, where L is the segment length and the number Z is related to the machine accuracy and can
be chosen as, for example, $Z \le 12$. The accuracy does not depend on the length of the segment L and is completely
determined by the time step if the following requirements for the number Z are fulfilled: $Z \le 2$, when $k = 1$; $Z \le 5$,
when $k = 2$; $Z \le 10$, when $k = 3$ [47].

Our goal is to implement this method to much more challenging and complicated Scalak problem. It is
interesting to note that the idea of inserting the time step as the parameter into space interpolation function was
expressed by well-known FEM scientist J.N. Reddy [49], and supported by K.J. Bathe et al. [50], where they appeal
to the analogy with the spectral method and therefore use terms in the element interpolation functions that change
very drastically to reflect all possible quick modes of vibration. A similar realization for the 1D rod problem is given
in [51].

However, our intention is much broader than the sole solution of the Scalak problem. Over the last several
years, we develop the multiphysical Method of Matched Section [52-56], which is similar to FEM. Like FEM, we
mesh the entire body on elements and give the solution within the element in a strong form sense as the approximate
solution of governing differential equations. With this goal, the particular (with respect to space coordinates)
differential equations are broken down into ordinary ones by introducing additional unknowns. The elements are matched along their common borders (sections), and such equations are treated as PTM, and the solutions within the element are presented as FTM. The advantage of the method is that the vibration problem [56] is performed on the same mesh, by the same algorithm, with the same unknowns as the static problem [53]; in fact, it realizes the spectral approach for a 2D plate there. The same can be said about the heat transient problems – its realization does not differ in principle from the steady-state problem solution [55]. The only requirement is to have 1D direct time integration solutions by TMM like the aforementioned solutions [47, 48].

Thus, we consider the presented solution of the Scalak problem as a prerequisite for treating more complicated 2D and 3D dynamic problems by MMS. For convenience, we name our method as timestep-dependent TMM, TDTMM. To be able to judge the accuracy of our method, we will compare our results with those obtained by Fourier analysis. The methodology of the Fourier solution was formulated by the first author and his former Ph.D. student Ya. Dubyk in [57,58]. Although this problem is relatively simple, to our knowledge, it has not yet been applied for description of the known experiments [59], which will be considered here to investigate the accuracy of our methods.

# 2. Setting equations. Homogeneous solution

## 2.1. Setting equations

The problem operates with four independent parameters: $p_{fl}$ is a water pressure, $\sigma_x$ is a stress in the pipe, $u$ is a liquid displacement, $v$ is a pipe metal displacement. Let us write down the basic equations of Scalak's problem:

$$ \frac{\partial p_{fl}(x)}{\partial x} = -\rho_{fl} \frac{\partial^2 u(x)}{\partial t^2} \qquad (1a) $$

$$ \frac{\partial \sigma_x(x)}{\partial x} = -\rho_p \frac{\partial^2 v(x)}{\partial t^2} \qquad (1b) $$

$$ \frac{du}{dx} = - \frac{p_{fl}a}{K} - 2\sigma_x \frac{\mu}{E} \qquad (1c) $$

$$ \frac{dv}{dx} = - \frac{\sigma_x}{E} - \frac{\varepsilon a p_{fl}}{K} \qquad (1d) $$

Where $\rho_{fl}, \rho_p$ - densities of liquid and metal, respectively, $\mu$ is the Poisson's ratio, $K$ is a modulus of compression of water, $E$ is a modulus of elasticity of metal, and constants $a, \varepsilon$ are given by the expressions:

$$ a = \left(1 + \frac{2KR}{tE}\right); \quad \varepsilon = \frac{\mu RK}{tE + K2R} < \mu/2 \ll 1 \qquad (1e) $$

Where $R$ is the radius and $t$ is the wall thickness of the pipe.

As in our previous works [47,48], we use the classical implicit central difference scheme for time derivatives, which we write in terms of equal time intervals $\Delta t$. For the brevity, we will write that $\Delta t = t$. Then, the time moment
i corresponds to the absolute time $T_i = i \cdot t$. Let's move on to capital letters to denote all parameters at time i. The classic central difference scheme is written as follows:

$$ \frac{\partial^2 u_{fl}(x)}{\partial t^2} = \frac{U_i(x) - 2U_{i-1}(x) - U_{i-2}(x)}{t^2} \qquad (2a) $$

$$ \frac{\partial^2 v_{fl}(x)}{\partial t^2} = \frac{v_i(x) - 2v_{i-1}(x) - v_{i-2}(x)}{t^2} \qquad (2b) $$

Then the governing equations take the form:

$$ \frac{\partial P_i(x)}{\partial x} + \rho_{fl} \frac{U_i(x)}{t^2} = \frac{U_{R,i-1}(x)\rho_{fl}}{t^2} \qquad (3a) $$

$$ \frac{\partial \Sigma_i(x)}{\partial x} + \rho_p \frac{V_i(x)}{t^2} = \frac{V_{R,i-1}(x)\rho_p}{t^2} \qquad (3b) $$

$$ \frac{dU_i(x)}{dx} = -\frac{P_i(x)a}{K} - 2\Sigma_i(x)\frac{\mu}{E} \qquad (3c) $$

$$ \frac{dV_i(x)}{dx} = -\frac{\Sigma_i(x)}{E} - \frac{\epsilon a P_i(x)}{K} \qquad (3d) $$

Where, for convenience, the right-hand sides of the equations (already known values) are denoted as:

$$ U_{R,i-1}(x) = 2U_{i-1}(x) - U_{i-2}(x); \quad V_{R,i-1}(x) = 2V_{i-1}(x) - V_{i-2}(x) \qquad (3e) $$

For convenience, differential equations (3c) and (3d) can be represented as follows:

$$ P_i(x) = -\frac{dU_i(x)}{dx} \frac{K}{a(1-2\mu\epsilon)} + 2 \frac{K\mu dV_i(x)}{a(1-2\mu\epsilon)dx} \qquad (4a) $$

$$ \Sigma_i(x) = -\frac{EdV_i(x)}{(1-2\mu\epsilon)dx} + \frac{\epsilon EdU_i(x)}{(1-2\mu\epsilon)dx} \qquad (4b) $$

We will need other combinations of differential equations. In particular, by differentiating (3c) and (3d) and substituting (3a) and (3b) into them, we obtain:

$$ \left( \frac{d^2 U_i(x)}{dx^2} - \frac{2\mu d^2 V_i(x)}{dx^2} \right) - m^2 U_i(x) = -m^2 U_{R,i-1}(x) \qquad (5a) $$

$$ \left( \frac{d^2 V_i(x)}{dx^2} - \frac{\epsilon d^2 U_i(x)}{dx^2} \right) - n^2 V_i(x) = -n^2 V_{R,i-1}(x) \qquad (5b) $$

Where we introduce the following notation:

$$ \frac{(1-2\mu\epsilon)}{c^2t^2} = m^2 \frac{(1-2\mu\epsilon)}{s^2t^2} = n^2 \qquad (5c) $$

And the classical wave propagation speeds separately for hydraulics and mechanics problems, i.e. c and s are defined as follows:

$$ \rho_{fl} \frac{a}{K} = \frac{1}{c^2}; \quad \rho_p \frac{1}{E} = \frac{1}{s^2} \qquad (5d) $$

From expressions (5a), we first find $\frac{d^2 V_i(x)}{dx^2}$, and then, differentiating (5a) twice and substituting $\frac{d^2 V_i(x)}{dx^2}$ into it, we obtain:

$$ \frac{d^4 U_i(x)}{dx^4} - \frac{n^2+m^2}{(1-2\mu\epsilon)} \frac{d^2 U_i(x)}{dx^2} + \frac{n^2 m^2}{(1-2\mu\epsilon)} U_i(x) = F(x) \qquad (6a) $$

$$ F(x) = - \frac{m^2}{(1-2\mu\epsilon)} \frac{d^2 U_{R,i-1}(x)}{dx^2} - \frac{2\mu n^2}{(1-2\mu\epsilon)} \frac{d^2 V_{R,i-1}(x)}{dx^2} + \frac{n^2 m^2}{(1-2\mu\epsilon)} U_{R,i-1}(x) \qquad (6b) $$

This is the basic differential equation, which we will solve as the sum of the partial and general solutions.

## 2.2. Fundamental solution

Consider the homogeneous linear equation (6), where the right-hand side $F(x)$ is set to zero, i.e.:
$$ \frac{d^4 U_i(x)}{dx^4} - \frac{n^2+m^2}{(1-2\mu\epsilon)} \frac{d^2 U_i(x)}{dx^2} + \frac{n^2 m^2}{(1-2\mu\epsilon)} U_i(x) = 0 \quad (7a) $$

Find its characteristic roots from the equation:

$$ \lambda^4 - \frac{n^2+m^2}{(1-2\mu\epsilon)}\lambda^2 + \frac{n^2m^2}{(1-2\mu\epsilon)} = 0 \quad (7b) $$

Find the discriminant of this equation:

$$ D = \frac{(n^2+m^2)^2-4n^2m^2(1-2\mu\epsilon)}{(1-2\mu\epsilon)^2} = \frac{(n^2-m^2)^2+4n^2m^2(1-2\mu\epsilon)}{(1-2\mu\epsilon)^2} \quad (7c) $$

Obviously, it is always positive as the sum of two positive numbers. Therefore, the squares of the characteristic roots are as follows:

$$ \lambda_{1,2}^2 = \frac{1}{2(1-2\mu\epsilon)} \left( (n^2 + m^2) \pm \sqrt{(n^2 + m^2)^2 - 4n^2m^2(1-2\mu\epsilon)} \right) \quad \lambda_1^2 > \lambda_2^2 \quad (7d) $$

Obviously, they are always positive. Therefore, the characteristic (eigen) functions are as follows:

$$ \Phi_1 = \exp(\lambda_1 x); \quad \Phi_2 = \exp(-\lambda_1 x); \quad \Phi_3 = \exp(\lambda_2 x); \quad \Phi_4 = \exp(-\lambda_2 x) \quad (7e) $$

In accordance with the general methodology for solving ordinary linear differential equations in mathematical physics problems, we introduce the concept of generalized Krylov functions, which have the properties listed in Table 1. These functions are given by the following expressions (8a):

<table><caption>Table 1 – Properties of Krylov functions</caption><thead><tr><th><i>i</i></th><th><i>K</i><sub><i>i</i></sub>(0)</th><th><i>K</i>'<sub><i>i</i></sub>(0)</th><th><i>K</i>''<sub><i>i</i></sub>(0)</th><th><i>K</i>'''<sub><i>i</i></sub>(0)</th></tr></thead><tbody><tr><td><i>K</i><sub>1</sub></td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td><i>K</i><sub>2</sub></td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td><i>K</i><sub>3</sub></td><td>0</td><td>0</td><td>1</td><td>0</td></tr><tr><td><i>K</i><sub>4</sub></td><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

$$ \begin{aligned} K_1(x) &= \frac{\lambda_1^2 \cosh \lambda_2 x - \lambda_2^2 \cosh \lambda_1 x}{\lambda_1^2 - \lambda_2^2} & K_2(x) &= \frac{\lambda_1^3 \sinh \lambda_2 x - \lambda_2^3 \sinh \lambda_1 x}{\lambda_1 \lambda_2 (\lambda_1^2 - \lambda_2^2)} \\ K_3(x) &= \frac{\cosh \lambda_1 x - \cosh \lambda_2 x}{\lambda_1^2 - \lambda_2^2} & K_4(x) &= \frac{\lambda_2 \sinh \lambda_1 x - \lambda_1 \sinh \lambda_2 x}{\lambda_1 \lambda_2 (\lambda_1^2 - \lambda_2^2)} \end{aligned} \quad (8a) $$

It is easy to verify the following differential dependencies between them:

$$ \begin{aligned} K'_{4}(x) &= K_{3}(x) & K'_{2}(x) &= K_{1}(x) \\ K'_{3}(x) &= K_{2}(x) + (\lambda_{1}^{2} + \lambda_{2}^{2})K_{4}(x) & K'_{1}(x) &= -\lambda_{1}^{2}\lambda_{2}^{2} \cdot K_{4}(x) \end{aligned} \quad (8b) $$

## 2.3. Solution for a homogeneous equation using TMM

For a numerical solution, we divide the entire pipe into a set of elements. For each element j, we use a local coordinate system starting at point x = 0, and designate the notions of all parameters at its beginning:

$$ U(x=0) = U_{j,0}; \quad P(x=0) = P_{j,0}; \quad V(x=0) = V_{j,0}; \quad \Sigma(x=0) = \Sigma_{j,0} \quad (9a) $$

Hereinafter, we omit the index j whenever possible, understanding that the calculation scheme applies to all elements. We will find expressions for all parameters in a form convenient for application to the TMM. First, we write down a homogeneous solution for the liquid displacement that satisfies the initial conditions (9a). Let us assume for it:

$$ U(x) = AK_1(x) + BK_2(x) + CK_3(x) + DK_4(x) \quad (9b) $$
First, we use condition (9a), from which we find, taking into account the properties of Krylov functions, given in Table 1, that $A = U_0$. Then, differentiating (9b) and taking into account the homogeneous equation (3c), conditions (9b) and (9d) and the properties of Krylov functions, we find $B = -\frac{P_0 a}{K} - 2\Sigma_0 \frac{\mu}{E}$. Differentiating (9b) again, we find, taking into account the homogeneous equations (3a) and (3b), the value of the coefficient C, and upon further differentiation, the value of the coefficient D. Collecting similar terms, we obtain an expression for $U(x)$. Similarly, we obtain expressions for the other three basic parameters. For homogeneous solutions, we use the upper dash and represent them in matrix form:

$$ \{\bar{U}(\mathbf{x})\} = \begin{pmatrix} \bar{U}(\mathbf{x}) \\ \bar{V}(\mathbf{x}) \\ \bar{P}(\mathbf{x}) \\ \bar{\Sigma}(\mathbf{x}) \end{pmatrix} = \begin{bmatrix} m_{11}(\mathbf{x}) & \dots & m_{14}(\mathbf{x}) \\ m_{21}(\mathbf{x}) & \dots & m_{24}(\mathbf{x}) \\ m_{31}(\mathbf{x}) & \dots & m_{34}(\mathbf{x}) \\ m_{41}(\mathbf{x}) & \dots & m_{44}(\mathbf{x}) \end{bmatrix} \begin{pmatrix} U_0 \\ V_0 \\ P_0 \\ \Sigma_0 \end{pmatrix} = [\mathbf{M}(\mathbf{x})]\{\mathbf{U}_0\} \quad (10a) $$

Where the coefficients of the matrix $[\mathbf{M}(\mathbf{x})]$ are given by the expressions in Table 2. According to the TMM methodology, the matrix $[\mathbf{M}(\mathbf{x})]$ has the following property:

$$ [\mathbf{M}(0)] = [I] \quad (10b) $$

Where $[I]$ is a unit diagonal matrix.

<table><caption>Table 2. Values of the coefficients of the homogeneous solution matrix.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>1</td><td>$K_1(x) + \frac{m^2}{(1-2\mu\varepsilon)} K_3(x)$</td><td>$\frac{2\mu n^2}{(1-2\mu\varepsilon)} K_3(x)$</td><td>$-\frac{a}{K} \left( K_2(x) + \frac{(m^2 + 2\mu\varepsilon n^2)K_4(x)}{(1-2\mu\varepsilon)} \right)$</td><td>$-\frac{2\mu}{E} \left( K_2(x) + \frac{(m^2 + n^2)K_4(x)}{(1-2\mu\varepsilon)} \right)$</td></tr><tr><td>2</td><td>$\frac{\varepsilon m^2}{(1-2\mu\varepsilon)} K_3(x)$</td><td>$\left( K_1(x) + \frac{n^2}{(1-2\mu\varepsilon)} K_3(x) \right)$</td><td>$-\frac{\varepsilon a}{K} \left( K_2(x) + \frac{m^2 + n^2}{(1-2\mu\varepsilon)} K_4(x) \right)$</td><td>$-\frac{1}{E} \left( K_2(x) + \frac{n^2 + 2\mu\varepsilon m^2}{(1-2\mu\varepsilon)} K_4(x) \right)$</td></tr><tr><td>3</td><td>$-\frac{K m^2}{a(1-2\mu\varepsilon)} \left( K_2(x) + \frac{m^2 K_4(x)}{(1-2\mu\varepsilon)} \right)$</td><td>$-\frac{2\mu n^2 m^2}{(1-2\mu\varepsilon)^2 a} \frac{K}{K_4(x)}$</td><td>$\left( K_1(x) + \frac{m^2}{(1-2\mu\varepsilon)} K_3(x) \right)$</td><td>$\frac{2\mu m^2 K}{E a (1-2\mu\varepsilon)} K_3(x)$</td></tr><tr><td>4</td><td>$-\frac{\varepsilon n^2 m^2 E}{(1-2\mu\varepsilon)^2} K_4(x)$</td><td>$-\frac{E n^2}{(1-2\mu\varepsilon)} \left( K_2(x) + \frac{n^2 K_4(x)}{(1-2\mu\varepsilon)} \right)$</td><td>$\frac{a \varepsilon E n^2}{K (1-2\mu\varepsilon)} K_3(x)$</td><td>$\left( K_1(x) + \frac{n^2}{(1-2\mu\varepsilon)} K_3(x) \right)$</td></tr></tbody></table>

# 3. Heterogeneous solution of the TMM

## 3.1. Methodology for finding a non-homogeneous solution

Since the general solution of a homogeneous equation satisfies the initial conditions and is represented by (10), the partial solution must satisfy zero initial conditions. This is achieved as follows.

1. Solutions at past moments of time (3d) are represented in the form of polynomials (for example, up to the third order). This will be explained in detail later:

$$ U_{R,i-1}(x) = A_0 + A_1 x + A_2 x^2 + \dots + A_k x^k \quad (11a) $$

$$ V_{R,i-1}(x) = B_0 + B_1 x + B_2 x^2 + \dots + B_k x^k \quad (11b) $$

2. A regular (classic) partial solution $\{\bar{U}_{cl,p}(\mathbf{x})\}$ that satisfies the differential equations (3a) - (3d) is sought. It is obtained by the well-known method of undetermined coefficients (in this case, in a polynomial expansion). For each $A_k$ or $B_k$, it is represented as
$$ \{\tilde{\mathbf{U}}_{cl,p}^k(\mathbf{x})\} = \{\tilde{\mathbf{U}}_{cl,p}^{A_k}(\mathbf{x})\} + \{\tilde{\mathbf{U}}_{cl,p}^{B_k}(\mathbf{x})\} = \begin{cases} \tilde{\mathbf{U}}_{cl,p}^k(\mathbf{x}) \\ \tilde{\mathbf{V}}_{cl,p}^k(\mathbf{x}) \\ \tilde{\mathbf{P}}_{cl,p}^k(\mathbf{x}) \\ \tilde{\mathbf{\Sigma}}_{cl,p}^k(\mathbf{x}) \end{cases} = A_k \begin{cases} \alpha_{U,0}^k + \alpha_{U,1}^k x + \cdots + \alpha_{U,k}^k x^k \\ \alpha_{V,0}^k + \alpha_{V,1}^k x + \cdots + \alpha_{V,k}^k x^k \\ \alpha_{P,0}^k + \alpha_{P,1}^k x + \cdots + \alpha_{P,k}^k x^k \\ \alpha_{\Sigma,0}^k + \alpha_{\Sigma,1}^k x + \cdots + \alpha_{\Sigma,k}^k x^k \end{cases} + B_k \begin{cases} \beta_{U,0}^k + \beta_{U,1}^k x + \cdots + \beta_{U,k}^k x^k \\ \beta_{V,0}^k + \beta_{V,1}^k x + \cdots + \beta_{V,k}^k x^k \\ \beta_{P,0}^k + \beta_{P,1}^k x + \cdots + \beta_{P,k}^k x^k \\ \beta_{\Sigma,0}^k + \beta_{\Sigma,1}^k x + \cdots + \beta_{\Sigma,k}^k x^k \end{cases} \tag{11c} $$

3. We are looking for a complete zero partial solution $\{\tilde{\mathbf{U}}_{z,p}(\mathbf{x})\}$, which should be equal to zero at $x = 0$. It is written as the sum of the ordinary partial solution and the supplemented solutions, $\{\tilde{\mathbf{U}}_{ad,p}(\mathbf{x})\}$. Then the complete zero partial solution, $\{\tilde{\mathbf{U}}_{z,p}^k(\mathbf{x})\}$ for each $A_k$ and $B_k$ is as follows:

$$ \{\tilde{\mathbf{U}}_{z,p}^k(\mathbf{x})\} = \{\tilde{\mathbf{U}}_{cl,p}^k(\mathbf{x})\} + \{\tilde{\mathbf{U}}_{ad,p}^k(\mathbf{x})\}, \tag{11d} $$

where the supplemented solution $\{\tilde{\mathbf{U}}_{ad,p}^k(\mathbf{x})\}$ is essentially the solution of the homogeneous equation (7a), and therefore can formally be represented as the general solution of the homogeneous equation (10), where instead of the specified initial conditions, some currently unknown constants are used, i.e.:

$$ \{\tilde{\mathbf{U}}_{ad,p}^k(\mathbf{x})\} = [\mathbf{M}(\mathbf{x})]\{\mathbf{X}^{A_k}\} + [\mathbf{M}(\mathbf{x})]\{\mathbf{X}^{B_k}\}, \tag{11e} $$

where $[\mathbf{M}(\mathbf{x})]$ is a known matrix given by (10), and the columns $\{\mathbf{X}^{A_k}\}$ and $\{\mathbf{X}^{B_k}\}$ contain unknown coefficients and are represented as follows:

$$ \{\mathbf{X}^{A_k}\} = \{X_1^{A_k}, X_2^{A_k}, X_3^{A_k}, X_4^{A_k}\}^T \qquad \{\mathbf{X}^{B_k}\} = \{X_1^{B_k}, X_2^{B_k}, X_3^{B_k}, X_4^{B_k}\}^T \tag{11f} $$

Equation (11d) is first used to find the unknown coefficients $\{\mathbf{X}^{A_k}\}$ and $\{\mathbf{X}^{B_k}\}$ in (11e) under the condition that at the zero point the $\{\tilde{\mathbf{U}}_{z,p}(\mathbf{x})\}$ solution is equal to zero, i.e.:

$$ \{\tilde{\mathbf{U}}_{cl,p}^k(\mathbf{0})\} + \{\tilde{\mathbf{U}}_{ad,p}^k(\mathbf{0})\} = 0 \tag{11g} $$

Then it is used in the general algorithm as a known value $\{\tilde{\mathbf{U}}_{ad,p}^k(\mathbf{x})\}$, which is a component of the complete zero partial solution $\{\tilde{\mathbf{U}}_{z,p}^k(\mathbf{x})\}$ (11d).

## 3.2. Examples of a heterogeneous solution

Let us formally demonstrate, using an example only for non-zero $A_k$, how formula (11g) works to find the unknowns $\{\mathbf{X}^{A_k}\}$. According to formulas (11c) - (11d), taking into account the known property (10b) of the matrix $[\mathbf{M}(\mathbf{0})]$, we obtain from equality (11g):

$$ X_1^{A_k} = -A_k \alpha_{U,0}^k, \quad X_2^{A_k} = -A_k \alpha_{V,0}^k, \quad X_3^{A_k} = -A_k \alpha_{P,0}^k, \quad X_4^{A_k} = -A_k \alpha_{\Sigma,0}^k \tag{12a} $$

Now let's explain this in more detail using a more complicated example $U_{R,i-1}(\mathbf{x}) = A_{R,2}x^2$. The usual partial solution for all four parameters of the problem is as follows:
$$ \{\tilde{U}_{cl,p}^{A_2}(\mathbf{x})\} = \left\{ \frac{2}{m^2} A_{R,2} + A_{R,2}x^2; -\frac{2\epsilon}{n^2} A_{R,2}; -\frac{K}{a(1-2\mu\epsilon)} A_{R,2}2x; \frac{E\epsilon}{1-2\epsilon\mu} A_{R,2}2x \right\}^T \quad (12b) $$

We are looking for the unknown coefficients {$\mathbf{X}^{A_2}$} from condition (11g), taking into account (10b):

$$ \{\tilde{U}_{cl,p}^{A_2}(\mathbf{0})\} + \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \{\mathbf{X}^{A_2}\} = \mathbf{0} \quad (12c) $$

From which we obtain the values of the unknown coefficients {$\mathbf{X}^{A_2}$} :

$$ \{\mathbf{X}^{A_2}\} = \left\{ -\frac{2}{m^2} A_{R,2} \quad \frac{2\epsilon}{n^2} A_{R,2} \quad 0 \quad 0 \right\}^T \quad (12d) $$

Then the complete zero partial solution $\tilde{U}_{z,p}^{A_2}(\mathbf{x})$ for each of the parameters according to (11d) will be as follows:

$$ \{\tilde{U}_{z,p}^{A_2}(\mathbf{x})\} = \begin{cases} \frac{2}{m^2} A_{R,2}(1 - m_{11}) + A_{R,2}(-x^2 + m_{12}) \\ 2A_{R,2}\left(-\frac{\epsilon}{n^2} - \frac{1}{m^2}m_{21} + \frac{\epsilon}{n^2}m_{22}\right) \\ 2A_{R,2}\left(-\frac{K}{a(1-2\mu\epsilon)}x - \frac{1}{m^2}m_{31} + \frac{\epsilon}{n^2}m_{32}\right) \\ 2A_{R,2}\left(\frac{E\epsilon}{1-2\epsilon\mu}x - \frac{1}{m^2}m_{41} + \frac{\epsilon}{n^2}m_{42}\right) \end{cases} \quad (12e) $$

Table 3 shows the complete zero partial solutions for $k = 0,1,2$, where the coefficients of matrix (10a) are used, i.e. $m_{i,j}(\mathbf{x})$.

<table><caption>Table 3 – Complete zero partial solution for all parameters at $k = 0,1,2$.</caption><thead><tr><th></th><th>k</th><th>{U<sup>A<sub>k</sub></sup><sub>z,p</sub>(x)}</th><th>{U<sup>B<sub>k</sub></sup><sub>z,p</sub>(x)}</th></tr></thead><tbody><tr><td rowspan="3">U&#773;</td><td>0</td><td>A<sub>R,0</sub>(1 - m<sub>11</sub>)</td><td>-B<sub>R,0</sub>m<sub>12</sub></td></tr><tr><td>1</td><td>A<sub>R,1</sub>x - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">K</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>A<sub>R,1</sub>(1 - m<sub>13</sub>) + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">εE</div><div style="padding-top: 1px;">1 - 2με</div></div>A<sub>R,1</sub>(1 + m<sub>14</sub>)</td><td>-<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μK</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>B<sub>R,1</sub>m<sub>13</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">E</div><div style="padding-top: 1px;">1 - 2με</div></div>B<sub>R,1</sub>m<sub>14</sub></td></tr><tr><td>2</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2</div><div style="padding-top: 1px;">m<sup>2</sup></div></div>A<sub>R,2</sub>(1 - m<sub>11</sub>) + A<sub>R,2</sub>(-x<sup>2</sup> + m<sub>12</sub>)</td><td>-<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μ</div><div style="padding-top: 1px;">m<sup>2</sup></div></div>2B<sub>R,2</sub>(1 - m<sub>11</sub>) - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2</div><div style="padding-top: 1px;">n<sup>2</sup></div></div>B<sub>R,2</sub>m<sub>12</sub></td></tr><tr><td rowspan="3">V&#773;</td><td>0</td><td>-A<sub>R,0</sub>m<sub>21</sub></td><td>B<sub>R,0</sub>(1 - m<sub>22</sub>)</td></tr><tr><td>1</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">K</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>A<sub>R,1</sub>m<sub>23</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">εE</div><div style="padding-top: 1px;">1 - 2με</div></div>A<sub>R,1</sub>m<sub>24</sub></td><td>B<sub>R,1</sub>(x - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μK</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>m<sub>23</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">E</div><div style="padding-top: 1px;">1 - 2με</div></div>m<sub>24</sub>)</td></tr><tr><td>2</td><td>2A<sub>R,2</sub>(-<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">ε</div><div style="padding-top: 1px;">n<sup>2</sup></div></div>-<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">1</div><div style="padding-top: 1px;">m<sup>2</sup></div></div>m<sub>21</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">ε</div><div style="padding-top: 1px;">n<sup>2</sup></div></div>m<sub>22</sub>)</td><td>2B<sub>R,2</sub>(<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">1</div><div style="padding-top: 1px;">n<sup>2</sup></div></div> + x<sup>2</sup> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μ</div><div style="padding-top: 1px;">m<sup>2</sup></div></div>m<sub>21</sub> - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">1</div><div style="padding-top: 1px;">n<sup>2</sup></div></div>m<sub>22</sub>)</td></tr><tr><td rowspan="3">P&#773;</td><td>0</td><td>-A<sub>R,0</sub>m<sub>31</sub></td><td>-B<sub>R,0</sub>m<sub>32</sub></td></tr><tr><td>1</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">K</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>A<sub>R,1</sub>(m<sub>33</sub> - 1) + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">εE</div><div style="padding-top: 1px;">1 - 2με</div></div>A<sub>R,1</sub>m<sub>34</sub></td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μK</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>B<sub>R,1</sub>(1 - m<sub>33</sub>) + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">E</div><div style="padding-top: 1px;">1 - 2με</div></div>B<sub>R,1</sub>m<sub>34</sub></td></tr><tr><td>2</td><td>2A<sub>R,2</sub>(-<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">K</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>x - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">1</div><div style="padding-top: 1px;">m<sup>2</sup></div></div>m<sub>31</sub> + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">ε</div><div style="padding-top: 1px;">n<sup>2</sup></div></div>m<sub>32</sub>)</td><td>2B<sub>R,2</sub>(<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μK</div><div style="padding-top: 1px;">a(1 - 2με)</div></div>x + <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2μ</div><div style="padding-top: 1px;">m<sup>2</sup></div></div>m<sub>31</sub> - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">1</div><div style="padding-top: 1px;">n<sup>2</sup></div></div>m<sub>32</sub>)</td></tr><tr><td>Σ</td><td>0</td><td>-A<sub>R,0</sub>m<sub>41</sub></td><td>-B<sub>R,0</sub>m<sub>42</sub></td></tr></tbody></table>
<table><tr><td>1</td><td>$\frac{\epsilon E}{1 - 2\epsilon\mu} A_{R,1}(1 + m_{44}) + \frac{K}{a(1 - 2\mu\epsilon)} A_{R,1}m_{43}$</td><td>$-\frac{E}{1 - 2\epsilon\mu} B_{R,1}(1 - m_{44}) - \frac{2\mu K}{a(1 - 2\mu\epsilon)} B_{R,1}m_{43}$</td></tr><tr><td>2</td><td>$2A_{R,2}\left(\frac{E\epsilon}{1 - 2\epsilon\mu}x - \frac{1}{m^2}m_{41} + \frac{\epsilon}{n^2}m_{42}\right)$</td><td>$2B_{R,2}\left(-\frac{E}{1 - 2\epsilon\mu}x + \frac{2\mu}{m^2}m_{41} - \frac{1}{n^2}m_{42}\right)$</td></tr></table>

At the end of this paragraph, we note that all four components of ordinary partial solutions are proportional to $A_k$ and $B_k$. Therefore, it is convenient to represent the complete zero partial solution as:

$$ \{\tilde{U}_{z,p}^k(\mathbf{x})\} = \sum_{k=0}^{n} A_{R,k} \begin{pmatrix} a_{Z,U}^k(\mathbf{x}) \\ a_{Z,V}^k(\mathbf{x}) \\ a_{Z,P}^k(\mathbf{x}) \\ a_{Z,\Sigma}^k(\mathbf{x}) \end{pmatrix} + \sum_{k=0}^{n} B_{R,k} \begin{pmatrix} b_{Z,U}^k(\mathbf{x}) \\ b_{Z,V}^k(\mathbf{x}) \\ b_{Z,P}^k(\mathbf{x}) \\ b_{Z,\Sigma}^k(\mathbf{x}) \end{pmatrix} \quad (12f) $$

# 4. General TMM algorithm

## 4.1. General solution

The entire pipeline is divided into elements $j$, each with a length of $l_j$. In total, we use $J$ elements. The complete solution for each section is found as the sum of the general homogeneous and complete zero partial solutions:

$$ \{\mathbf{U}^j(\mathbf{x})\} = \{\overline{\mathbf{U}}^j(\mathbf{x})\} + \{\tilde{\mathbf{U}}^j_{z,p}(\mathbf{x})\} \quad (13a) $$

Which, in turn, are given by the dependencies (10a) and (12f). For each element $j$, a state vector is introduced, which is a column vector of variables at its beginning and end $\mathbf{U}_0^j$ and $\mathbf{U}_e^j$:

$$ \{\mathbf{U}_0^j\} = \{\mathbf{U}^j(x = 0)\} \quad \{\mathbf{U}_e^j\} = \{\mathbf{U}^j(x = l_j)\} \quad (13b) $$

That is, for each element 8 unknowns are introduced, which means that the entire problem operates with 8J unknowns. For each element, we have 4 connection equations obtained from (13a) at $x = l_j$, that is:

$$ \{\mathbf{U}_e^j\} = [\mathbf{M}(l_j)]\{\mathbf{U}_0^j\} + \sum_{k=0}^{n} A_{R,k}^j \begin{pmatrix} a_{Z,U}^k(l_j) \\ a_{Z,V}^k(l_j) \\ a_{Z,P}^k(l_j) \\ a_{Z,\Sigma}^k(l_j) \end{pmatrix} + \sum_{k=0}^{n} B_{R,k}^j \begin{pmatrix} b_{Z,U}^k(l_j) \\ b_{Z,V}^k(l_j) \\ b_{Z,P}^k(l_j) \\ b_{Z,\Sigma}^k(l_j) \end{pmatrix} \quad (13c) $$

Thus, we have 4J connection (FTM) equations. Another 4(J - 1) conjugation equations (PTM) are set at the borders between elements. In this case, they are trivial:

$$ \{\mathbf{U}_0^j\} = [I]\{\mathbf{U}_e^{j-1}\} \quad (13d) $$

The problem is supplemented by 4 boundary conditions. With the correct formulation of the problem, we have 2 conditions on the left end of the system and 2 conditions on the right. Thus, the number of equations coincides with the number of unknowns.

Thus, a system of equations is formed and solved with respect to the values of the determining parameters (13b) at the edges of the elements. The values of the state vector functions at intermediate points are determined by formula (13a), taking into account (10a) and (12e).
## 4.2. Transition to the next time step

The result of each time step $i-1$ for element $j$ is solution (13a), which is formally represented as a combination of a polynomial of a pre-selected degree $k$ and a combination of functions $m_{g,h}(x)$. In turn, $m_{g,h}(x)$ are combinations of Krylov functions, which in turn are combinations of functions $ch\lambda_2x, ch\lambda_1x, sh\lambda_2x$ and $sh\lambda_1x$. Therefore, the solution for the liquid and the pipe displacements is formally written as follows:

$$U_i^j(x) = c_0^{j,i} + c_1^{j,i}x + \dots + c_k^{j,i}x^k + \gamma_{1,c}^{j,i}ch\lambda_1x + \gamma_{2,c}^{j,i}ch\lambda_2x + \gamma_{1,s}^{j,i}sh\lambda_1x + \gamma_{2,s}^{j,i}sh\lambda_2x \quad (14a)$$

$$V_i^j(x) = d_0^{j,i} + d_1^{j,i}x + \dots + d_k^{j,i}x^k + \epsilon_{1,c}^{j,i}ch\lambda_1x + \epsilon_{2,c}^{j,i}ch\lambda_2x + \epsilon_{1,s}^{j,i}sh\lambda_1x + \epsilon_{2,s}^{j,i}sh\lambda_2x \quad (14b)$$

The next step is to expand the functions $ch\lambda_2x, ch\lambda_1x, sh\lambda_2x$, and $sh\lambda_1x$ into polynomials of order $k$ [48]. Let us illustrate this with the example of the function $sh\lambda x$. Let us assume that we have expansion:

$$sh\lambda x = s_0^{\lambda} + s_1^{\lambda}x + \dots + s_k^{\lambda}x^k \quad (15a)$$

Multiply both parts of (15a) by $x^m$ and $0 \le m \le k$, respectively, and integrate over the entire length of the element. We obtain $k$ equations of the form:

$$\int_{0}^{l_j}(sh\lambda x \cdot x^m)dx = \int_{0}^{l_j}(s_0^{\lambda} + s_1^{\lambda}x + \dots + s_k^{\lambda}x^k)x^m dx \quad 0 \le m \le k \quad (15b)$$

Both parts are easily solved analytically, which allows to obtain a system of equations for the unknown coefficients $s_k^\lambda$. Similarly, expand the function $ch\lambda x$. Let us assume that we have:

$$ch\lambda x = c_0^{\lambda} + c_1^{\lambda}x + \dots + c_k^{\lambda}x^k \quad (15c)$$

Then, to determine the coefficients $c_k^\lambda$, we have the following system of equations, where the coefficients of the unknowns and free terms are easily found by analytical integration:

$$\int_{0}^{l_j}(ch\lambda x \cdot x^m)dx = \int_{0}^{l_j}(c_0^{\lambda} + c_1^{\lambda}x + \dots + c_k^{\lambda}x^k)x^m dx; 0 \le m \le k \quad (15d)$$

Replacing the exponential functions in expressions (14a) and (14b) with polynomial functions according to representations (15a) and (15b), we obtain a polynomial representation of the "past" solution at time $i$:

$$U_i^j(x) = \tilde{c}_0^{j,i} + \tilde{c}_1^{j,i}x + \dots + \tilde{c}_k^{j,i}x^k \quad (16a)$$

$$V_i^j(x) = \tilde{d}_0^{j,i} + \tilde{d}_1^{j,i}x + \dots + \tilde{d}_k^{j,i}x^k \quad (16b)$$

The last preparatory step for moving on to the next step in time is to apply formula (3e):

$$U_{R,i}^j(x) = 2U_i^j(x) - U_{i-1}^j(x); V_{R,i}^j(x) = 2V_i^j(x) - V_{i-1}^j(x) \quad (16c)$$

The algorithm for the proposed TDTMM is complete.

# 4. Fourier analysis for the reference problem [59].

This section is auxiliary. Its main purpose is to provide a reference analytical solution for the well-known Delft Hydraulics Benchmark Problem A [59], below Problem A, for simplicity. We will use this solution to compare the accuracy of the above-proposed TDTMM.

## 4.1. General solution of the TMM

In this section, we have to repeat some of the results of [2] for the sake of completeness. The Fourier method, like the full dynamic method, operates with the same four parameters and uses the same Skalak equations (1a-d). Let the sought functions be decomposed in time as periodic functions, i.e.:
$$u = U(x) \cos \omega t \quad p_{fl} = P(x) \cos \omega t \quad \sigma_p = \Sigma(x) \cos \omega t \quad v = V(x) \cos \omega t \quad (17)$$

Then, substituting (17) into Skalak's equations (1a-d), we obtain:

$$\frac{dP(x)}{dx} = \rho_{fl} \omega^2 U(x) \quad (18a)$$

$$\frac{d\Sigma(x)}{dx} = \rho_p \omega^2 V(x) \quad (18b)$$

$$\frac{dU(x)}{dx} = -\frac{P(x)a}{K} - 2\Sigma(x)\frac{\mu}{E} \quad (18c)$$

$$\frac{dV(x)}{dx} = -\frac{\Sigma(x)}{E} - \frac{\epsilon a P(x)}{K} \quad (18d)$$

Expressing all unknowns, for example, through $U(x)$, we obtain a 4th order differential equation with respect to water displacement:

$$\frac{d^4 U(x)}{dx^4} + \frac{(s^2+c^2)\omega^2}{s^2c^2} \frac{d^2 U(x)}{dx^2} + \frac{(1-2\epsilon\mu)\omega^4}{s^2c^2} U(x) = 0 \quad (19a)$$

The roots of the differential equation (19a) are:

$$\lambda_{1,2}^2 = -\frac{\omega^2}{2s^2c^2} \left( (s^2+c^2) \pm \sqrt{(s^2-c^2)^2 + 8\epsilon\mu s^2c^2} \right) \quad (19b)$$

Since for all materials $\mu < 0.5$ and $\epsilon < 0.25$, the squares of the roots (19b) are negative, and the roots themselves are imaginary. Thus, we obtain a set of eigenfunctions:

$$F_1(x) = \cos \lambda_1 x; \quad F_2(x) = \sin \lambda_1 x; \quad F_3(x) = \cos \lambda_2 x; \quad F_4(x) = \sin \lambda_2 x \quad (19c)$$

Let us form generalized Krylov functions $K_i$, $1 \le i \le 4$ from the four functions (19c), which have the properties listed in Table 1. Note that these functions are identical to the functions $C_i(x)$ [2]:

$$K_1(x) = \frac{\lambda_1^2 \cos \lambda_2 x - \lambda_2^2 \cos \lambda_1 x}{\lambda_1^2 - \lambda_2^2}; \quad K_2(x) = \frac{\lambda_1^2 \sin \lambda_2 x / \lambda_2 - \lambda_2^2 \sin \lambda_1 x / \lambda_1}{\lambda_1^2 - \lambda_2^2} \\ K_3(x) = \frac{\cos \lambda_2 x - \cos \lambda_1 x}{\lambda_1^2 - \lambda_2^2}; \quad K_4(x) = \frac{\sin \lambda_2 x / \lambda_2 - \sin \lambda_1 x / \lambda_1}{\lambda_1^2 - \lambda_2^2} \quad (19d)$$

Similar to (8b), we obtain differential dependencies between Krylov functions:

$$K'_1(x) = -\lambda_1^2 \lambda_2^2 K_4(x) K'_2(x) = K_1(x) \\ K'_3(x) = K_2(x) - (\lambda_1^2 + \lambda_2^2) K_4(x) \quad K'_4(x) = K_3(x) \quad (19e)$$

The general solution of the homogeneous equation (19a) can again be written in the form (9b). Write down the initial conditions at the zero point, i.e., the conditions:

$$U(x=0) = U_0; \quad P(x=0) = P_0; \quad V(x=0) = V_0; \quad \Sigma(x=0) = \Sigma_0 \quad (20a)$$

which, unlike the conditions (9a), refer to much longer segments. In the case of the Problem A [59] to be analyzed below, they refer to the entire pipe, i.e., only one segment is considered. Then taking into account the conditions (20a), dependencies (18) between all parameters, and the properties of Krylov functions (19d) and (19e) the solutions for $U(x)$, $V(x)$, $P(x)$, $\Sigma(x)$ can also be represented in matrix form convenient for application of TMM:
$$ \{U(\mathbf{x})\} = \begin{Bmatrix} \bar{U}(\mathbf{x}) \\ \bar{V}(\mathbf{x}) \\ \bar{P}(\mathbf{x}) \\ \bar{\Sigma}(\mathbf{x}) \end{Bmatrix} = \begin{bmatrix} n_{11}(\mathbf{x}) & \dots & n_{14}(\mathbf{x}) \\ n_{21}(\mathbf{x}) & \dots & n_{24}(\mathbf{x}) \\ n_{31}(\mathbf{x}) & \dots & n_{34}(\mathbf{x}) \\ n_{41}(\mathbf{x}) & \dots & n_{44}(\mathbf{x}) \end{bmatrix} \begin{Bmatrix} U_0 \\ V_0 \\ P_0 \\ \Sigma_0 \end{Bmatrix} = [N(\mathbf{x})]\{U_0\} \quad (20b) $$

Where the functions $n_{n,m}(\mathbf{x})$ are presented in Table 4.

<table><caption>Table 4. Values of the coefficients of the homogeneous solution matrix.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>1</th><td>$K_1(x) - \frac{\omega^2}{c^2} K_3(x)$</td><td>$-\omega^2 \frac{2\mu}{s^2} K_3(x)$</td><td>$\frac{a}{K} \left(-K_2(x) + \omega^2 \left(\frac{1}{c^2} + \frac{2\mu\epsilon}{s^2}\right) K_4(x)\right)$</td><td>$\frac{2\mu}{E} \left(-K_2(x) + \omega^2 \left(\frac{1}{c^2} + \frac{1}{s^2}\right) K_4(x)\right)$</td></tr><tr><th>2</th><td>$-\omega^2 \frac{\epsilon}{c^2} K_3(x)$</td><td>$K_1(x) - \frac{\omega^2}{s^2} K_3(x)$</td><td>$\frac{\epsilon a}{K} \left(-K_2(x) + \omega^2 \left(\frac{1}{s^2} + \frac{1}{c^2}\right) K_4(x)\right)$</td><td>$\frac{1}{E} \left(-K_2(x) + \omega^2 \left(\frac{1}{s^2} + \frac{2\mu\epsilon}{c^2}\right) K_4(x)\right)$</td></tr><tr><th>3</th><td>$\rho_{fl} \omega^2 \left(K_2(x) - \frac{\omega^2}{c^2} K_4(x)\right)$</td><td>$-\omega^4 \frac{\alpha \rho_p}{s^2} K_4(x)$</td><td>$K_1(x) - \frac{\omega^2}{c^2} K_3(x)$</td><td>$-\omega^2 \frac{\alpha}{s^2} K_3(x)$</td></tr><tr><th>4</th><td>$-\omega^4 \frac{\rho_{fl} \beta}{c^2} K_4(x)$</td><td>$\omega^2 \rho_p \left(K_2(x) - \frac{\omega^2}{s^2} K_4(x)\right)$</td><td>$-\omega^2 \frac{\beta}{c^2} K_3(x)$</td><td>$K_1(x) - \frac{\omega^2}{s^2} K_3(x)$</td></tr></tbody></table>

Where these symbols are used:

$$ \rho_{fl} \frac{2\mu}{E} = \frac{\alpha}{s^2} \rho_p \frac{\epsilon a}{K} = \frac{\beta}{c^2} \quad (20c) $$

## 4.2. Natural frequencies for Problem A

To verify the described method, we will compare it with the Problem A [59], which describes water hammer in a tank-pipe-valve system, Fig. 1. This system has the following parameters: pipe length $L = 20$ m, radius $R = 398.5$ mm, pipe wall thickness $t = 8$ mm, Young's modulus $E = 210$ GPa, pipe density $\rho_s = 7900$ kg/m³, Poisson's ratio $\mu = 0.3$, elastic modulus $K = 2.1$ GPa, water density $\rho_{fl} = 1000$ kg/m³, initial water flow velocity $V_0 = 1$ m/s. The pipe cross-sectional area $S_p$ and cross-sectional area $S_{fl}$ are equal to:

$$ S_p \approx 2\pi R t \approx 0.02 \, m^2, \quad S_{fl} = \pi R^2 = 0.4989 \, m^2 \quad (21a) $$

Calculate the other constants used in formulas (5g), (1d), and (20b):

$$ a = 2.00; \quad \epsilon = 0.07486; \quad s = 5155.8; \quad c = 1025.66, \quad \alpha = 0.07595; \quad \beta = 0.5914 \quad (21b) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Scheme of the model problem The Delft Hydraulics Benchmark Problem A</figcaption></figure>
Next, proceed to find the natural frequencies. To do this, apply zero conditions on the left edge for water pressure and pipe displacement to equations (20b):

$$P_0 = P(x = 0) = 0; V_0 = V(x = 0) = 0, \qquad (21c)$$

as well as the conditions of stress equality and force equilibrium on the right edge, i.e.:

$$U(L) = V(L); \quad P(L)S_{fl} + \Sigma(L)S_p = 0 \qquad (21d)$$

From the first two conditions (21c), we directly obtain that the state vector at the initial point is as follows:

$$\{U_0\} = \{U_0; 0; 0; \Sigma_0\}^T \qquad (22a)$$

Now let's apply the conditions on the right edge. Substituting into the matrix (20b) with (22a) and (21d) taken into account, we get two characteristic equations for obtaining roots (natural frequencies), as those for which the determinant of the system is zero:

$$U_0(n_{11}(L) - n_{41}(L)) = \Sigma_0(-n_{14}(L) + n_{44}(L)) \qquad (22b)$$

$$U_0(n_{21}(L)S_{fl} + n_{41}(L)S_p) = \Sigma_0(-n_{24}(L)S_{fl} - n_{44}(L)S_p) \qquad (22c)$$

Then search for the natural forms. Take $\Sigma_0 = 1$, and substituting this value into equation (22b), we find $U_0$ as the initial value of the form for the water displacement, $U_{\phi,0}$, i.e.:

$$U_{\phi,0} = \frac{(-n_{14}(L)+n_{44}(L))}{(n_{11}(L)-n_{41}(L))} \qquad (22d)$$

Then conditions (21c) are substituted into equation (20b), instead of $\Sigma_0$ we substitute 1, instead of $U_0$ - the value found in (22b), thus the eigenforms are given as:

$$\{\Phi_i(x)\} = \begin{Bmatrix} \Phi_{U,i}(x) \\ \Phi_{V,i}(x) \\ \Phi_{P,i}(x) \\ \Phi_{\Sigma,i}(x) \end{Bmatrix} = \begin{bmatrix} n_{11}(x) & \dots & n_{14}(x) \\ n_{21}(x) & \dots & n_{24}(x) \\ n_{31}(x) & \dots & n_{34}(x) \\ n_{41}(x) & \dots & n_{44}(x) \end{bmatrix} \begin{Bmatrix} U_{\phi,0} \\ 0 \\ 0 \\ 1 \end{Bmatrix} \qquad (22e)$$

Where each eigenform $\Phi_i(x)$ corresponds to the natural frequency $\omega_i$.

Let us list the first three natural frequencies and their corresponding periods $T = \frac{2\pi}{\omega}$. We have: $\omega_1 = 77.86$, $T_1 = 0.0807$; $\omega_2 = 199.7$, $T_2 = 0.0315$; $\omega_3 = 349.72$, $T_3 = 0.018$.

### 4.3. Orthogonality of eigenfunctions

In this section, we generally follow, with some modifications, the little-known work [58].

All eigenforms satisfy equations (18a-d). Write down equations (18a) and (18b), but replace all parameters with eigenforms that satisfy these equations. Multiply the first of them by the expression $S_{fl} \cdot \Phi_{U,m}(x)$, and the second equation by $S_p \cdot \Phi_{U,m}(x)$, where the lower index $m$ corresponds to another eigenfrequency. Integrate these expressions over the entire length of the pipe:

$$\int_0^L S_{fl} \frac{d\Phi_{P,n}(x)}{dx} \Phi_{U,m}(x) dx = \omega_n^2 \int_0^L S_{fl} \rho_{fl} \Phi_{U,n}(x) \cdot \Phi_{U,m}(x) dx \qquad (23a)$$
$$ \int_0^L S_p \frac{d\Phi_{\Sigma,n}(x)}{dx} \Phi_{V,m}(x) dx = \omega_n^2 \int_0^L S_p \rho_p \Phi_{V,n}(x) \Phi_{V,m}(x) dx \quad (23b) $$

Adding (23a) and (23b), we obtain the following combination:

$$ \int_0^L \left( S_{fl} \frac{d\Phi_{P,n}}{dx} \Phi_{U,m} + S_p \frac{d\Phi_{\Sigma,n}}{dx} \Phi_{V,m} \right) dx = \omega_n^2 \int_0^L \left( m_{fl} \Phi_{U,n} \Phi_{U,m} + m_p \Phi_{V,n} \Phi_{V,m} \right) dx \quad (23c) $$

Where $m_{fl}$ is the specific weight of the liquid, $m_p$ is the specific weight of the pipe material:

$$ m_{fl} = S_{fl} \rho_{fl}, \quad m_p = S_p \rho_p. \quad (23d) $$

Similarly, we obtain an expression with interchanged indices:

$$ \int_0^L \left( S_{fl} \frac{d\Phi_{P,m}}{dx} \Phi_{U,n} + S_p \frac{d\Phi_{\Sigma,m}}{dx} \Phi_{V,n} \right) dx = \omega_m^2 \int_0^L \left( m_{fl} \Phi_{U,m} \Phi_{U,n} + m_p \Phi_{V,m} \Phi_{V,n} \right) dx \quad (23e) $$

The right-hand sides of expressions (23c) and (23e) differ only in their multipliers—the squares of the natural frequencies. Our goal is to show that the difference between the right-hand sides of expressions (23c) and (23e) is zero. Then, in the case of $\omega_n^2 \neq \omega_m^2$, expressions of the form:

$$ \int_0^L (m_{fl} \Phi_{U,n} \Phi_{U,m} + m_p \Phi_{V,n} \Phi_{V,m}) dx = 0, \quad \text{if } n \neq m, \quad (23e) $$

should be equal to zero, which is the condition of orthogonality if $n \neq m$.

Let's prove this. Consider the left-hand side of expression (23c). Integrating it by parts, we obtain:

$$ E_{nm} = \left( S_{fl} \Phi_{P,n} \Phi_{U,m} + S_p \Phi_{\Sigma,n} \Phi_{V,m} \right) \Big|_0^L - \int_0^L \left( S_{fl} \Phi_{P,m} \dot{\Phi}_{U,n} + S_p \Phi_{\Sigma,m} \dot{\Phi}_{V,n} \right) dx \quad (24a) $$

The first expression on the right-hand side is equal to zero due to the boundary conditions (for the condition at the point $x = 0$ (21c) – this is obvious, because both the first and second terms contain the factor 0; and at the right edge, i.e., at the point $x = L$ (21d) – the sum of the two terms is equal to zero). Now let's consider the second term on the right-hand side (24a). Note that due to (18c) and (18d), it can be written as follows:

$$ E_{nm} = - \int_0^L \left( S_{fl} \Phi_{P,m} \left( -\frac{\Phi_{P,n}a}{K} - 2\Phi_{\Sigma,n} \frac{\mu}{E} \right) + S_p \Phi_{\Sigma,m} \left( -\frac{\Phi_{\Sigma,n}}{E} - \frac{\varepsilon a \Phi_{P,n}}{K} \right) \right) dx = E_{mn}, \quad (24b) $$

where it was taken into account the values of constants $\varepsilon$, $a$ (1e) and for areas (21a). Now it is obvious that the difference between the right-hand sides (23c) and (23e) is zero, so the expression (23f) is also zero. The statement is proven.

Thus, the generalized eigenfunctions:

$$ \Phi_i(x) = \begin{cases} \Phi_{U,i}(x) & \text{for fluid} \\ \Phi_{V,i}(x) & \text{for pipe} \end{cases} \quad (24c) $$

are orthogonal to each other with the following weight $\Phi_{w,n}(x)$:

$$ \Phi_{w,i}(x) = \begin{cases} m_{fl} \Phi_{U,i}(x) & \text{for fluid} \\ m_p \Phi_{V,i}(x) & \text{for pipe} \end{cases} \quad (24d) $$

When $n = m$, it is obvious that expression (23e) is nonzero. Let us call it the norm of the eigenform and denote it as $|M_{i,i}| = |M_i|$:

$$ |M_i| = \int_0^L \left( m_{fl} (\Phi_{U,i})^2 + m_p (\Phi_{V,i})^2 \right) dx \quad (24e) $$
### 4.4. Application of Fourier series to dynamics problems

Equations (1a-d) are considered with the generalization that the force equations (1a), (1b) take into account the action of an external force:

$$ \frac{\partial p_{fl}(x)}{\partial x} S_{fl} = -\rho_{fl} \frac{\partial^2 u_{fl}(x)}{\partial t^2} S_{fl} + S_{fl} q_{fl}(x) \quad (25a) $$

$$ \frac{\partial \sigma_p(x)}{\partial x} S_p = -\rho_p \frac{\partial^2 v(x)}{\partial t^2} S_p + S_p q_p(x), \quad (25b) $$

where $q_{fl}(x)$, $q_p(x)$ are some known external forces, including concentrated forces.

Let's show how to apply the method of decomposition into eigenforms to solve dynamics problems. Present the complete solution for displacements $\vec{\Pi}(x, t)$, as:

$$ \vec{\Pi}(x, t) = \sum_n A_i(t) \begin{cases} \Phi_{U,i}(x) & \text{for fluid} \\ \Phi_{V,i}(x) & \text{for pipe wall} \end{cases} \quad (25c) $$

where $A_i(t)$ are unknown functions of time for each individual form. Multiply equation (25a) by $\Phi_{U,i}(x)$ and equation (25b) by $\Phi_{V,i}(x)$, add them, and integrate over the entire length of the pipe. Due to the orthogonality conditions (23e) established above, we will obtain the following differential equation for each unknown function $A_i(t)$:

$$ \ddot{A}_i(t) + A_n(t) \omega_i^2 = \frac{q_i(t)}{|M_i|}, \quad (25d) $$

where the generalized force $q_i(t)$ for each mode $i$ is as follows:

$$ q_i(t) = \int_0^L (S_{fl} q_{fl}(x, t) \Phi_{U,i} + S_p q_p(x, t) \Phi_{V,i}) dx \quad (25e) $$

In the case of Problem A, the external forces are zero, and the generalized displacements $\vec{\Pi}(x, t)$ become as follows:

$$ \vec{\Pi}(x, t) = \sum_i \left( A_{i,c} \cos(\omega_i t) + \frac{A_{i,s}}{\omega_i} \sin(\omega_i t) \right) \begin{cases} \Phi_{U,i}(x), & \text{fluid} \\ \Phi_{V,i}(x), & \text{pipe} \end{cases} \quad (26a) $$

The unknown constants $A_{i,c}$ and $A_{i,s}$ are found from the initial conditions (at $t = 0$) for displacement and velocity. Let the initial position and initial velocity be known:

$$ \vec{\Pi}(x, 0) = \vec{\Pi}_0(x) = \begin{cases} U_{fl,o}(x) \\ V_{p,o}(x) \end{cases}; \quad \frac{d(\vec{\Pi}(x,t=0))}{dt} = \dot{\vec{\Pi}}_0(x) = \begin{cases} \dot{U}_{fl,o}(x) \\ \dot{V}_{p,o}(x) \end{cases} \quad (26b) $$

Then, equating (26a) with (26b) at time zero, we have that the coefficients must satisfy the following conditions:

$$ \sum_i A_{i,c} \begin{cases} \Phi_{U,i}(x), & \text{fluid} \\ \Phi_{V,i}(x), & \text{pipe} \end{cases} = \begin{cases} U_{fl,o}(x) \\ V_{p,o}(x) \end{cases} \quad (26c) $$

$$ \sum_i A_{i,s} \begin{cases} \Phi_{U,i}(x), & \text{fluid} \\ \Phi_{V,i}(x), & \text{pipe} \end{cases} = \begin{cases} \dot{U}_{fl,o}(x) \\ \dot{V}_{p,o}(x) \end{cases} \quad (26d) $$

To determine them, we multiply (26c) and (26d) by the system of orthogonal functions $\Phi_{w,i}(x)$ (24d).

Since both metal and water are in the zero position at time zero, it is obvious that:

$$ A_{i,c} \equiv 0 \quad (27a) $$
Now consider expression (26d), where, according to the conditions of the problem, $\dot{U}_{fl,o}(x) = 1$ and $\dot{V}_{p,o}(x) = 0$. Multiply both parts by the functions (24g). This gives the following solutions for $A_{i,s}$:

$$A_{i,s}|M_i| = |M_{fl,i}| \qquad (27b)$$

Where the partial mass of water for each it natural form $i$ is given by the formula:

$$|M_{fl,i}| = m_{fl} \int_{0}^{L} \Phi_{u,i}(x) dx \qquad (27c)$$

## 5. Results for The Problem A.

In [47], the problem of wave propagation only in a pipe (rod) was investigated using the proposed TDTMM, where the time step was included in the analytical solutions. The following features of the method were established there. It was found that the accuracy of the solution depends on two parameters. The first main parameter is the time step, and the second auxiliary parameter is the characteristic number, namely the maximum value of the argument of the exponential function included in all analytical solutions. Denote the value of the argument as $CH$:

$$CH = \lambda_1 \cdot l \qquad (28a)$$

Where we recall that $l$ is the length of the segment of pipe. Denote the maximum acceptable value of the argument as:

$$[CH_k] = f(k) \qquad (28b)$$

Where $k$ is the degree of polynomial approximation of the pipe and water displacement solutions.

Reducing the time step leads to a continuous improvement in accuracy (provided that $[CH]$ is chosen appropriately), i.e., to reduction of the solution dissipation [47]; while exceeding the acceptable value of $[CH]$ leads to the local distortions on the graphs.

These properties and limitations of solutions will be explored below.

### 5.1. Dependence of the acceptable characteristic number from the degree of the polynomial $k$

Let us write the characteristic roots $\lambda_1$ and $\lambda_2$ in a convenient form:

$$\lambda_1 = \frac{1}{tA}, \quad \lambda_2 = \frac{1}{tB}, \qquad (29a)$$

where $A$ and $B$ can be easily found from equation (7g) taking into account (5b). These $A$ and $B$ have the meaning of wave propagation velocity. For the values of the system parameters described in section 4.2, $A$ and $B$ take the following values: $A = 1024.71 \frac{m}{c}$, $B = 5280.51 \frac{m}{c}$.

Let's examine the maximum acceptable value of the characteristic number $CH$ in TDTMM for a second-degree polynomial $k = 2$. We choose a time moment such that the faster wave has traveled three-quarters of the pipe, i.e., approximately 15 meters. This moment of time is $T = \frac{15}{B} = \frac{15}{5280.51} = 0.0028 \approx 0.003 s$.

Consider the solution at this moment $T = 0.003 s$. Let's take a time step $t = 10^{-5}$. Then $\lambda_1 = 97.59$, $\lambda_2 = 18.94$. First, let's compare the results for the following three characteristic number values: $CN = 6$, $CN = 10$, $CN =$
14, which are realized by adjusting the length of each element at meshing (number of elements, $f$). After plotting the pressure graphs in Figure 2a, we see that the graphs of the solutions completely coincide on the macro scale. It may seem that the solutions obtained are completely identical. This may lead to the false conclusion that the criterion value $[CH_2]$ (i.e., at $k = 2$) can be taken as large as 14, but this is not the case. Let's take three characteristic number values: $CN = 6, CN = 7, CN = 8$ and construct a micro figure. In Figure 2a, we marked the micro figure area with a red circle. When constructing the micro pattern, Fig. 2b, we see that the graphs at $CN = 6$ and $CN = 7$ coincide between themselves, yet for the value $CN = 8$ there are small wiggles. Therefore, to be sure of the accuracy of the results obtained, we will conservatively fix the acceptable value $[CH]$ at $k = 2$ as $[CH_2] = 7$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Graph of the pressure at $T = 3 \cdot 10^{-3}$ s: a) for three characteristic numbers: $CN = 6, CN = 10, CN = 14$; b) micro figure area for three characteristic numbers: $CN = 6, CN = 7, CN = 8$</figcaption></figure>

Above, we examined the acceptable value of the characteristic number on the example of pressure. Now check whether the stress solution imposes additional restrictions on the value of the characteristic number. Construct stress graphs for the values of the characteristic number: $CN = 6, CN = 10, CN = 14$, Fig. 3a. We see that all graphs completely coincide. Now, for the values of the characteristic number: $CN = 6, CN = 7, CN = 8$, we construct micro stress graphs, Fig. 3b. Evidently, the graphs completely coincide. This indicates that stress does not impose additional restrictions, therefore the acceptable value of the characteristic number remains unchanged, i.e. $[CH_2] = 7$.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Stresses at $T = 3 \cdot 10^{-3}$ s: a) full graph for characteristic numbers $CN = 6, CN = 10, CN = 14$; b) micro figure area for three characteristic numbers: $CN = 6, CN = 7, CN = 8$</figcaption></figure>

Similarly, we will get the acceptable values of the characteristic number for $k = 1$ and $k = 0$. Summarizing these calculations, the characteristic values are:

$$[CH_2] = 7; [CH_1] = 2; [CH_0] = 0.000 \qquad (29b)$$

That is, it is impossible to approximate the solution with a zero-degree polynomial. This can be explained by the fact that Krylov functions are monotonic functions and cannot be represented as constants. To illustrate the choice of $[CH_1]$, Fig. 4 shows a part of the pressure graph for a first-degree polynomial for three characteristic numbers: 0.2, 2, 3. The graph shows that $k = 2$ matches well, while $k = 3$ gives slightly different results.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Local part of the pressure graph at $T = 3 \cdot 10^{-3}$ s, $k = 1$ and characteristic numbers $CN = 0.2, CN = 2$, and $CN = 8$</figcaption></figure>

## 5.2. Accuracy study depending on the time step, time windows

Herein and below TDTMM uses a second-order polynomial, $k = 2$, and a characteristic number of no more than 7.

As investigated in [47], in order for the wave to pass through the entire pipe and still keeps its amplitude, the time step must be at least 100 times smaller than the time it takes for the wave to pass. For the problem under
consideration, this time is equal to $T = \frac{L}{B} = \frac{20}{5280.51} = 0.0038$ seconds. Accordingly, the maximum reasonable time step is $3.8 \cdot 10^{-5}$. Below, we will show how the time step affects the result. We compare three steps: $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-6}$.

Consider the solution at the moment when the wave has traveled 15 meters, i.e., $T = 3 \cdot 10^{-3}$. Plot a graph of the pressure solution at this moment in time, which was obtained for three values of the time step, Fig. 5a. In Fig. 5, we see a clear pattern: the overall behavior of the graphs coincides, but with a larger time step value, the wave becomes more blurred. Compare the obtained solution with the solution using the Fourier method. Note that for the Fourier method, we use 1854 terms of the expansion. The resulting pressure graph is shown in Fig. 5b and compared with the TDTMM results for minimal time steps (the best TDTMM result). Evidently the graphs of the solutions in Fig. 5a and Fig. 5b coincide. Note that the Fourier method produces local artifacts (peaks) around jump discontinuities, which is a well-known Gibbs phenomenon [60].

<figure><img src="image_6.png" /><figcaption>Fig. 5. Pressure graph at $T = 3 \cdot 10^{-3}$ s calculated: a) by TDTMM for time steps $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 10^{-5}$, $\Delta t = 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>

Construct pressure graphs at a time moment equal to the sum of the previous time and the time it takes for the wave to travel through the entire pipe, i.e., $T = 3 \cdot 10^{-3} + 3.8 \cdot 10^{-3} = 6.8 \cdot 10^{-3}$ seconds, Fig. 6a. Also consider the solution for pressure using the Fourier method, comparing it with TDTMM at a time step equal to $\Delta t = 10^{-6}$, Fig. 6b. It is evident that the TDTMM results are in good agreement with the Fourier method, although by TDTMM the results are slightly smoothed at locations of small pressure jumps.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Pressure graph at time $T = 6.8 \cdot 10^{-3}$ s calculated: a) by TDTMM for time steps $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>

Let us consider another time point equal to $T = 6.8 \cdot 10^{-3} + 3.8 \cdot 10^{-3} = 1.06 \cdot 10^{-2}$ seconds. Construct pressure graphs by TDTMM for three-time step values, Fig. 7a, and a pressure graph using the Fourier method, Fig. 7b. As we can see in Fig. 7a, the dissipation becomes quite noticeable at a time step of $\Delta t = 3.8 \cdot 10^{-5}$ although the horizontality of graphs on certain length segments is still ensured. To demonstrate that reducing the time step can significantly improve the results of TDTMM, in Fig. 7b we show the results at a time step of $\Delta t = 4 \cdot 10^{-7}$.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Pressure graph at time $T = 1.06 \cdot 10^{-2}$ s calculated: a) by TDTMM for time steps $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>

Thus, it is evident that the smaller the time step, the steeper the wavefront remains on the graph, and the less the results are "smoothed out." However, as the time step increases, the difference in the pressure graphs compared to Fourier method becomes more noticeable.

Let's show the effect of the time step on the metal stress graphs obtained by TDTMM. For the same time step values, we will construct stress graphs at the same time points. Compare the results with the Fourier method. Fig. 8a shows the stress distribution with three different time steps at the time equal to $T = 3 \cdot 10^{-3}$. Fig. 8b compares the
TDTMM results with the Fourier method. Small deviations are observed only in the wavefront zone. As for the
"straight line" graph sections, they again coincide completely.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Stress graph at time $T = 0.003$ s calculated: a) by TDTMM for time steps $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-5}$, $\Delta t = 1 \cdot 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>

Now consider the time moment $T = 6.8 \cdot 10^{-3}$, which is constructed using the same time steps as in the
previous example. As we can see, there is more smoothing for stresses than for pressure, Fig. 9a. For example, local
small jumps in stresses around the point $x = 13.2$ are completely smoothed out at $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 10^{-5}$, and
are kept only at $\Delta t = 10^{-6}$, Fig. 9b.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Stress at time $T = 6.8 \cdot 10^{-3}$s calculated: a) by TDTMM for time steps $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 10^{-5}$, $\Delta t = 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>

Take another time point equal to $T = 1.06 \cdot 10^{-2}$ seconds. Construct pressure graphs using TDTMM for
three values of the time step, Fig. 10a, and those by the Fourier method, Fig. 10b. As we can see, Fig. 7a, at a time
step of $\Delta t = 3.8 \cdot 10^{-5}$, dissipation becomes quite noticeable, and even the horizontality of fairly long segments is
not ensured. The comparison of the Fourier results with TDTMM is performed at a time step equal to $\Delta t = 4 \cdot 10^{-7}$.
In general, the correspondence is good, but it is obvious that the "smearing" and attenuation for stresses is more pronounced than for pressure.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Stress at time $T = 1.06 \cdot 10^{-2}$ calculated: a) by TDTMM for time steps $\Delta t = 3.8 \cdot 10^{-5}$, $\Delta t = 10^{-5}$, $\Delta t = 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>

### 5.3. Accuracy study depending on the time step, point $X = 20$

Above, we constructed graphs along the entire pipe at a fixed point in time. Now fix the last point of the pipe $X = 20$ and construct graphs of pressure and stress in the time interval $T \in [0; 0.2]$. The time step is equal to $\Delta t = 1 \cdot 10^{-6}$, the number of terms in Fourier analysis is 1854. The correspondence between theoretical (Fourier) and numerical results are generally good for the pressure, Fig 11, and for the stress behavior, Fig 12. The general impression is that DTTMM better reproduces the pressure wave than the stress wave. Nevertheless, the DTTMM has a remarkable property to provide the better accuracy for the less time step.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Pressure graph at the last point $X = 20$ in the time interval $T \in [0; 0.2]$ calculated: a) by TDTMM with a time step of $\Delta t = 1 \cdot 10^{-6}$; b) by Fourier method (1854 terms)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Stress graph at the last point $X = 20$ in the time interval $T \in [0; 0.2]$ calculated: a) by TDTMM with a time step of $\Delta t = 1 \cdot 10^{-6}$ b) by Fourier method (1854 terms)</figcaption></figure>

## 5.4. Study of attenuation of TDTMM

A feature of all numerical methods for calculating dynamic processes is artificial amplitude decays (dissipation). Let's examine the behavior of dynamic processes (pressure and stress) at the last point $X = 20$. Take a time step $\Delta t = 4 \cdot 10^{-5}$ and construct graphs at the last point over a time interval $T \in [0; 120]$ seconds. Look how the wave attenuates over time. Fig. 13a shows the pressure values over time, and Fig. 13b shows the stresses over time. Both graphs resemble an exponential function. To verify this, take the logarithm of the local extrema and rebuild the graphs $\ln|P|$ and $\ln|\Sigma|$ for a time step of $\Delta t = 4 \cdot 10^{-5}$. In doing so, we show only the extreme values on the graphs.

They are shown in Fig. 14a (dependence of the pressure logarithm on time) and Fig. 14b (dependence of the stress logarithm on time). As expected, these graphs are linear.

<figure><img src="image_14.png" /><figcaption>Fig. 13. Last point $X = 20$ in the time interval $T \in [0; 120]$ with a time step of $\Delta t = 4 \cdot 10^{-5}$: a) graph of pressure; b) graph of stress</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14. Amplitude at the last point $X = 20$ in the time interval $T \in [0; 20]$ with a time step $\Delta t = 4 \cdot 10^{-5}$: a) graph of pressure; b) graph of stress</figcaption></figure>

Figures 13 and 14 show that the amplitudes of pressure and stress decrease according to the following law:

$$P = a_P \cdot e^{\frac{-t}{b_P}} \qquad \Sigma = a_\Sigma \cdot e^{\frac{-t}{b_\Sigma}} \qquad (30a)$$

Where $a$ and $b$ are certain constants, where $a$ corresponds to the extreme initial value of the parameter under consideration (pressure or stress), and the value of the constant $b$ determines how long it will take for the amplitude to decrease by a factor of $e$, i.e., the greater the value of $b$, the slower the attenuation. So far, the values of $a$ and $b$ in the two equations (30a) are considered to be different. By taking the logarithm of the left equation (30a), obtain the logarithm of the amplitude:

$$\ln|P| = \ln a_P - \frac{t}{b_P} \qquad (30b)$$

To find $b_P$ take two points from the graph, Fig. 13a, (15; 12.3021) and (20; 11.6967). Then $b$ is equal to:

$$b_P(\Delta t = 4 \cdot 10^{-5}) = \frac{20-15}{12.3021-11.6967} = 8.26 \qquad (31a)$$

Similarly, on the graph of the amplitude of the stress $\ln|\Sigma|$, Fig. 13b, take two points (15; 15.5172) and (20; 14.9119). Then $b_\Sigma$ is equal to:

$$b_\Sigma(\Delta t = 4 \cdot 10^{-5}) = \frac{20-15}{13.2105-12.9087} = 8.26 \qquad (31b)$$

As we can see, the values of $b$ for pressure and stress coincide.

Let's halve the time step $\Delta t = 2 \cdot 10^{-5}$ and plot the graph of the logarithm of the pressure at the last point $x = 20$, Fig. 15a. Plot a logarithmic graph for stress, Fig. 15b. Obviously, again for both dependencies we have straight lines, which can be represented as (30b). Again, substituting the two points of the graph, (15; 13.2105) and (20; 12.9087), into the equation of the straight line, it can be obtained for $b_P$:
<figure><img src="image_16.png" /><figcaption>Fig. 15. Amplitude at the last point $X = 20$ in the time interval $T \in [0; 20]$ with a time step $\Delta t = 2 \cdot 10^{-5}$ : a) graph of pressure; b) graph of stress</figcaption></figure>

$$b_P(\Delta t = 2 \cdot 10^{-5}) = \frac{20-15}{13.2105-12.9087} = 16.57 \quad (31c)$$

Similarly, on the stress amplitude graph $\ln|\Sigma|$, Fig. 15b, take two points (15; 16.4257) and (20; 16.1238). Then $b_\Sigma$ is equal to:

$$b_\Sigma(\Delta t = 2 \cdot 10^{-5}) = \frac{20-15}{16.4257-16.1238} = 16.56 \quad (31d)$$

Let's take a step in time that is half as small $\Delta t = 10^{-5}$ and construct a graph of the pressure and stress amplitude, Fig. 16.

<figure><img src="image_17.png" /><figcaption>Fig. 16. Amplitude at the last point $X = 20$ in the time interval $T \in [0; 20]$ with a time step $\Delta t = 10^{-5}$ : a) graph of pressure; b) graph of stress</figcaption></figure>

To find $b_P$, take two points from the graph(15; 13.662) and(20; 13.5111). Then $b$ is equal to:

$$b_P(\Delta t = 10^{-5}) = \frac{20-15}{13.662-13.5111} = 33.13 \quad (31e)$$

To find $b_\Sigma$, take two points from the graph,(15; 16.8702) and(20; 16.7262). We get:
$$b_{\Sigma}(\Delta t = 10^{-5}) = \frac{20-15}{16.8702-16.7262} = 34.72 \quad (31f)$$

The larger difference in the values of b is caused by errors in determining the corresponding parameters, taking into account the small attenuation.

Take a time step $\Delta t = 5 \cdot 10^{-6}$ and plot the logarithm of the pressure and stress amplitudes, Fig. 17.

<figure><img src="image_18.png" /><figcaption>Fig. 17. Amplitude at the last point $X = 20$ in the time interval $T \in [0; 20]$ with a time step $\Delta t = 5 \cdot 10^{-6}$: a) graph of pressure; b) graph of stress</figcaption></figure>

To find $b_P$ and $b_{\Sigma}$, take two points with abscises 10 and 5. Then we get:

$$b_P(\Delta t = 5 \cdot 10^{-6}) = \frac{10-5}{13.9946-13.9203} = 67.29 \quad (31g)$$

$$b_{\Sigma}(\Delta t = 5 \cdot 10^{-6}) = \frac{10-5}{17.2097-17.1354} = 67.29 \quad (31h)$$

Again, we see a good correspondence of the numbers b.

Present these numbers in tabular form, Table 5

<table><caption>Table 5. Dependence of the coefficient b on the time step and its inverse.</caption><thead><tr><th>Time step, Δt</th><th>Inverse step <sup>1</sup>&frasl;<sub>Δt</sub></th><th>b<sub>P</sub></th><th>b<sub>Σ</sub></th></tr></thead><tbody><tr><td>4 · 10<sup>-5</sup></td><td>2.5 · 10<sup>4</sup></td><td>8.26</td><td>8.26</td></tr><tr><td>2 · 10<sup>-5</sup></td><td>5 · 10<sup>4</sup></td><td>16.52</td><td>16.56</td></tr><tr><td>1 · 10<sup>-5</sup></td><td>10<sup>5</sup></td><td>33.13</td><td>34.72</td></tr><tr><td>5 · 10<sup>-6</sup></td><td>2 · 10<sup>5</sup></td><td>67.29</td><td>67.29</td></tr></tbody></table>

Based on the obtained data, construct a graph of the value b from the time step $b(\frac{1}{\Delta t})$, Fig. 18. Obviously, this is a straight line. We are looking for its equation in the form:

$$b = \alpha + \frac{1}{\Delta t} \beta \quad (32a)$$
<figure><img src="image_19.png" /><figcaption>Fig. 18. Graph of the dependence of b on $\frac{1}{\Delta t}$</figcaption></figure>

Connecting the first and last points, obtain $\beta = 3.37 \cdot 10^{-4}$, $\alpha = -0.11$. Obviously, consider the scale of the figure and the values of the coefficient b, the value of coefficient $\alpha$ can be taken as $\alpha = 0$. Therefore, the final coefficient of numerical attenuation, $\delta$, can be represented as:

$$ \delta(t) = -\ln \frac{P(t)}{P_0} = -\ln \frac{\Sigma(t)}{\Sigma_0} = t \cdot \frac{\Delta t}{3.37 \cdot 10^{-4}} \quad (32b) $$

Remind, that $\Delta t$ is the time step. Represent the value $\beta$ as:

$$ \beta = \frac{\beta_0}{(\omega_1)^2} \qquad \omega_1 = 77.86 \quad (32c) $$

Where the value $\omega_1$ is defined above and is equal to 77.86. Then $\beta_0 = 3.37 \cdot 10^{-4}(\omega_1)^2 = 2.042 \approx 2$, and the value $\delta(t)$ can be represented as follows:

$$ \delta(t) = t \cdot \Delta t \frac{(\omega_1)^2}{2} \quad (32d) $$

This formula allows us to control numerical damping, i.e., to understand the effect of time step on amplitude dissipation. For example, specifying the ultimate moment of time, $t = T_0$, and allowable numerical dissipation $[\delta]$, we can predict the required time step:

$$ \Delta t = (t) = \frac{2[\delta]}{(\omega_1)^2 T_0} \quad (32e) $$

Where $\omega_1$ is (approximate) first natural frequency, which could be assessed independently.

## Conclusion

The main novelty of the work lies in the application of the central difference method to time derivatives to get the integrable space related 4th order differential equation for the Skalak problem. However, the implementation of this seemingly simple idea lies in the original inclusion of the time step in the space differential equation, where the values of the functions at past moments of time become the right-hand side of this differential equation. It is solved
using specially constructed generalized Krylov functions, as the sum of the general solution of the homogeneous equation and the partial solution. They are obtained in a form convenient for application in the TMM. In fact, this work continues a series of our previous works on obtaining interpolation functions within an element that explicitly depend on the time step, but the complexity of this task greatly exceeds our previous works. Another important aspect of this work is the possibility of using its results in the method of matched sections for dynamic problems, which is an alternative to the finite element method. The achievements of this work are:

1. Obtaining specific TDTMM solutions and a detailed analysis of general and partial solutions using the generalized Krylov method, which is applied for the first time to such differential equations.

2. A detailed solution of the Skalak problem using the method of separation of variables (Fourier), which is used to compare and demonstrate the accuracy of the proposed TDTMM.

3. Demonstration of the accuracy of TDTMM, which, as shown in the work, depends proportionally on the time step, proportionally on the time moment under consideration, and proportionally on the square of the smallest natural frequency.

## Literature

1. Skalak, R., "An Extension of the Theory of Water Hammer," Trans. ASME, Vol. 78, No. 1, 1956.

2. Lesmez, M.W.; Wiggert, D.C.; Hatfield, F.J. Modal Analysis of Vibrations in Liquid-Filled Piping Systems. J. Fluids Eng.-T ASME. 1990, 112, 311-318.

https://doi.org/10.1115/1.2909406

3. Tentarelli, S. C. 1990. Propagation of noise and vibration in complex hydraulic tubing systems. Ph.D. Thesis, Lehigh University, Department of Mechanical Engineering, Bethlehem, U.S.A.

4. Jong, C. A. F. DE, 1994. Analysis of pulsations and vibrations in fluid-filled pipe systems. Ph.D. Thesis, Eindhoven University of Technology, Eindhoven, The Netherlands, ISBN 90-386-0074-7.

https://doi.org/10.6100/IR423649

5. Ferras, D., Manso, P. A., Schleiss, A. J., & Covas, D. I. (2018). One-dimensional fluid-structure interaction models in pressurized fluid-filled pipes: A review. *Applied Sciences, 8*(10), 1844.

https://doi.org/10.3390/app8101844

6. Hatfield, F.; Wiggert, D. Harmonic analysis of coupled fluid and piping. In Proceedings of the ASCE Engineering Mechanics Speciality Conference, West Lafayette, IN, USA, 23-25 May 1983; pp. 180-192.

7. Zhou, F., Liao, J., Chen, Z., Tan, X., &amp; He, L. (2024). The Mixed Variable Transfer Matrix Method and Its Application in Predicting the Frequency Domain Vibration Characteristics of Hydraulic Pipelines. *Applied Sciences, 14*(23), 10847

https://doi.org/10.3390/app142310847

8. Xu, Y., Johnston, D. N., Jiao, Z., & Plummer, A. R. (2014). Frequency modeling and solution of fluid-structure interaction in complex pipelines. *Journal of Sound and Vibration, 333*(10), 2800-2822.

https://doi.org/10.1016/j.jsv.2013.12.023

9. Guo, X., Ge, H., Xiao, C., Ma, H., Sun, W., & Li, H. (2022). Vibration transmission characteristics analysis of the parallel fluid-conveying pipes system: numerical and experimental studies. *Mechanical Systems and Signal Processing, 177*, 109180.

https://doi.org/10.1016/j.ymssp.2022.109180

10. Cao, Y. H., Liu, G. M., &amp; Hu, Z. (2023). Vibration calculation of pipeline systems with arbitrary branches by the hybrid energy transfer matrix method. *Thin-Walled Structures, 183*, 110442.

https://doi.org/10.1016/j.tws.2022.110442
11. Lee U. Spectral element method in structural dynamics. John Wiley & Sons; 2009. ISBN: 978-0-470-82374-3 (Hardcover)

12. Kim, T., Lee, B., & Lee, U. (2019). State-vector equation method for the frequency domain spectral element modeling of non-uniform one-dimensional structures. *International Journal of Mechanical Sciences*, **157**, 75-86.

https://doi.org/10.1016/j.ijmecsci.2019.05.004

13. Sun, Z., Yu, R., Wang, X., & Zhang, X. (2025). Dynamic analysis of liquid-filled clamp-pipe systems based on the spectral element method. *Archive of Applied Mechanics*, **95**(1), 1-28.

https://doi.org/10.1007/s00419-024-02094-2

14. Hatfield, F.J.; Wiggert, D.C.; Otwell, R.S. Fluid Structure Interaction in Piping by Component Synthesis. J. Fluid. Eng. 1982, 104, 318-325

https://doi.org/10.1115/1.3241840

15. Liu, G., Li, S., Li, Y., &amp; Chen, H. (2013). Vibration analysis of pipelines with arbitrary branches by absorbing transfer matrix method. *Journal of Sound and Vibration*, **332**(24), 6519-6536.

https://doi.org/10.1016/j.jsp.2013.06.020

16. F. Leckie, E. Pestel. Transfer-matrix fundamentals. International Journal of Mechanical Sciences, Volume 2, Issue 3, 1960, Pages 137-167.

https://doi.org/10.1016/0020-7403(60)90001-1

17. M.H. Chaudhry, Resonance in pressurized piping systems, PhD thesis, Department of Civil Engineering, University of British Columbia, Canada, 1970.

18. Arai, M., Kuroda, S., &amp; Ito, K. (2021). Static elastic-plastic analysis of a pipe structure by the transfer matrix method. *Journal of Pressure Vessel Technology*, **143**(1), 011801.

https://doi.org/10.1115/1.4047680

19. Orynyak, I. V., & Radchenko, S. A. (2004). A mixed-approach analysis of deformations in pipe bends. Part 3. Calculation of bend axis displacements by the method of initial parameters. *Strength of materials*, **36**, 463-472.

https://doi.org/10.1023/B:STOM.0000048394.98411.4f

20. Goto, Y., Morikawa, Y., & Matsuura, S. (1986). Direct Lagrangian nonlinear analysis of elastic space rods using transfer matrix technique. *Proc. of JSCE, Struct. Eng./Earthquake Eng*, **5**(1).

21. Tufekci, E., Eroglu, U., & Aya, S. A. (2017). A new two-noded curved beam finite element formulation based on exact solution. *Engineering with computers*, **33**, 261-273.

https://doi.org/10.1007/s00366-016-0470-1

22. Tijsseling, A. S. (1996). Fluid-structure interaction in liquid-filled pipe systems: a review. *Journal of Fluids and Structures*, **10**(2), 109-146.

https://doi.org/10.1006/jfls.1996.0009

23. Sreejith, B., Jayaraj, K., Ganesan, N., Padmanabhan, C., Chellapandi, P., & Selvaraj, P. (2004). Finite element analysis of fluid-structure interaction in pipeline systems. *Nuclear Engineering and Design*, **227**(3), 313-322.

https://doi.org/10.1016/j.nucengdes.2003.11.005

24. Mei, M., Cheng, H., Huang, H., Wu, Z., Yang, H., &amp; Wu, K. (2025). Dynamic analysis and shape sensing of pipes subjected to water hammer by the inverse finite element method. *Measurement*, **241**, 115701

https://doi.org/10.1016/j.measurement.2025.115701

25. Lavooij, C. S. W., & Tusseling, A. S. (1991). Fluid-structure interaction in liquid-filled piping systems. *Journal of fluids and structures*, **5**(5), 573-595.

https://doi.org/10.1016/S0889-9746(05)80006-4

26. Zanganeh, R., Ahmadi, A., & Keramat, A. (2015). Fluid-structure interaction with viscoelastic supports during waterhammer in a pipeline. *Journal of Fluids and Structures*, **54**, 215-234.

https://doi.org/10.1016/j.jfluidstructs.2014.10.016

27. Gale, J., & Tiselj, I. (2008). Godunov's Method for Simulations of Fluid-Structure Interaction in Piping Systems. *Journal of Pressure Vessel Technology*, **130**(3).
https://doi.org/10.1115/1.2937758

28. A-Moneim, M. T., & Chang, Y. W. (1978). Comparison of ICEPEL code predictions with straight flexible pipe experiments. *Nuclear Engineering and Design*, **49**(1-2), 187-196.

https://doi.org/10.1016/0029-5493(78)90061-4

29. Kulak, R. F. (1985). Three-dimensional fluid-structure coupling in transient analysis. *Computers & Structures*, **21**(3), 529-542.

https://doi.org/10.1016/0045-7949(85)90132-4

30. Sabbagh-Yazdi, S. R., & Mastorakis, N. E. (2007). Water hammer modeling by Godunov type finite volume method. *International Journal of Mathematics and Computers in Simulation*, **1**, 350-355.

31. Gale, J.; Tiselj, I. Applicability of the Godunov's method for fundamental four-equation FSI model. In *Proceedings of the International Conference on Nuclear Energy for New Europe, Bled, Slovenia, 5-8 September 2005*

32. Vardy, A. E., & Tijsseling, A. S. (2015, November). Method of characteristics:(Why) is it so good? In *Proceedings of the 12th International Conference on Pressure Surges, Dublin, Ireland* (pp. 18-20).

33. Wiggert, D.C.; Hatfield, F.J.; Stuckenbruck, S. (1987). Analysis of Liquid and Structural Transients in Piping by the Method of Characteristics. *Journal of Fluids Engineering, ASME*. 109(2), 161-165.

https://doi.org/10.1115/1.3242638

34. Huang, H., Guo, P., & Yan, J. (2025). Numerical discretization errors of fluid-structure interaction model in pressurized pipeline systems. *Applied Mathematical Modelling*, 11594.

https://doi.org/10.1016/j.apm.2025.115954

35. Ghidaoui, M. S., & Karney, B. W. (1994). Equivalent differential equations in fixed-grid characteristics method. *Journal of Hydraulic Engineering*, 120(10), 1159-1175.

https://doi.org/10.1061/(ASCE)0733-9429(1994)120:10(1159)

36. Zha, C., Pan, C., Sun, Z., &amp; Liu, Q. (2024). Reliability sensitivity analysis for water hammer-induced stress failure of fluid-conveying pipe. *Applied Mathematical Modelling*, 130, 51-65.

https://doi.org/10.1016/j.apm.2024.02.028

37. Goldberg, D. E., & Benjamin Wylie, E. (1983). Characteristics method using time-line interpolations. *Journal of Hydraulic Engineering*, 109(5), 670-683.

https://doi.org/10.1061/(ASCE)0733-9429(1983)109:5(670)

38. Xu, Y., & Jiao, Z. (2017). Exact solution of axial liquid-pipe vibration with time-line interpolation. *Journal of Fluids and Structures*, 70, 500-518.

39. Lavooij, C. S. W., & Tusseling, A. S. (1991). Fluid-structure interaction in liquid-filled piping systems. *Journal of Fluids and Structures*, 5(5), 573-595.

40. Vardy, A. E., & Alsarraj, A. T. (1989). Method of characteristics analysis of one-dimensional members. *Journal of Sound and Vibration*, 129(3), 477-487.

https://doi.org/10.1016/0022-460X(89)90436-7

41. Tang, L., Xiaoting, R., Jianshu, Z., &amp; Lina, Z. (2023). Riccati transfer equations for fluid structure interaction in liquid-filled piping systems. *Heliyon*, 9(5).

https://doi.org/10.1016/j.heliyon.2023.e15923

42. Gao, H., Guo, C., & Quan, L. (2021). Fluid-structure interaction analysis of aircraft hydraulic pipe with complex constraints based on discrete time transfer matrix method. *Applied Sciences*, 11(24), 11918.

https://doi.org/10.3390/app112411918

43. K.J. Bathe (2014). Finite Element Procedure. Prentice Hall, 1996. Second Edition. K.J Bathe, Watertown, Massachusetts.

44. Gao, H., Guo, C., Quan, L., & Wang, S. (2022). Frequency domain analysis of fluid-structure interaction in aircraft hydraulic pipe with complex constraints. *Processes*, 10(6), 1161.

https://doi.org/10.3390/pr10061161
45. Quan, L., Fu, C., Yao, R., & Guo, C. (2023). Dynamic Modeling and Parameter Identification of Double Casing Joints for Aircraft Fuel Pipelines. *Processes*, 11(9), 2767.

https://doi.org/10.3390/pr11092767

46. Batura, A., Bogdan, A., Orynyak, I., & Radchenko, S. (2014, July). Two axisymmetric models application for fast and conservative analysis of cold plumes in RPV at PTS condition. In *Pressure Vessels and Piping Conference* (Vol. 46001, p. V003T03A106). American Society of Mechanical Engineers.

47. Orynyak, I., Mazuryk, R., & Tsybulskyi, V. (2022). Semi-analytical implicit direct time integration scheme on example of 1-D wave propagation problem. *Mechanics and Advanced Technologies*, 6(2), 115-123.

https://doi.org/10.20535/2521-1943.2022.6.2.262110

48. Orynyak, I., Kostyushko, I., & Mazuryk, R. (2023). Semi-analytical implicit direct time integration method for 1-D gas dynamic problem. *Mechanics & Advanced Technologies*, 7(1).

https://doi.org/10.20535/2521-1943.2023.7.1.271273

49. J.N. Reddy, An Introduction to the Finite Element Method, 3rd edition, McGraw-Hill, New York, 2006.

50. Ham, S., & Bathe, K. J. (2012). A finite element method enriched for wave propagation problems. *Computers & Structures*, 94, 1-12.

https://doi.org/10.1016/j.compstruc.2012.01.001

51. Tang, T., Zhou, W., Luo, K., Tian, Q., & Hu, H. (2024). Dynamic modeling and analysis of discontinuous wave propagation in a rod. *Journal of Sound and Vibration*, 569, 117991.

https://doi.org/10.1016/j.jsv.2023.117991

52. Orynyak, I., & Danylenko, K., (2023). Method of matched sections in application to thin-walled and Mindlin rectangular plates. *Mechanics and Advanced Technologies*, 7(2).

https://doi.org/10.20535/2521-1943.2023.7.2.277341

53. Orynyak, I., Danylenko, K. (2024). Method of matched sections as a beam-like approach for plate analysis. *Finite Elements in Analysis and Design*. Volume 230, March 2024, 104103

https://doi.org/10.1016/j.finel.2023.104103

54. K. Danylenko and I. Orynyak, (2024). “Numeric analysis of elastic plane body static problem by the method of matched sections”, *Mech. Adv. Technol.*, vol. 8, no. 4(103), pp. 428–439,

https://doi.org/10.20535/2521-1943.2024.8.4(103).313412

55. Orynyak, I., Tsybulnyk, A., Danylenko, K., Oryniak, A., & Radchenko, S. (2025). Timestep-dependent element interpolation functions in the method of matched sections on the example of heat conduction problem. *Journal of Computational and Applied Mathematics*, 456, 116222

https://doi.org/10.1016/j.cam.2024.116222

56. Orynyak, I., Tsybulnyk, A., & Danylenko, K. (2025). Spectral realization of the method of matched sections for thin-plate vibration. *Archive of Applied Mechanics*, 95(2), 51.

https://doi.org/10.1007/s00419-024-02755-7

57. Dubyk, I. R., & Orynyak, I. V. (2016). Fluid-structure interaction in free vibration analysis of pipelines. *Scientific Journal of the Ternopil National Technical University*. 81(1).

58. Dubyk, I. R. (2016). Dynamic analysis of pipelines using the normal mode method taking into account interaction with the environment. *Naukoyemni Technologiyi*, (2), 242-248. (in Ukrainian)

https://doi.org/10.18372/2310-5461.30.10572

59. Tijsseling, A. S., & Lavooij, C. S. W. (1990). Waterhammer with fluid-structure interaction. *Applied Scientific Research*, 47, 273-285.

https://doi.org/10.1007/BF00418055

60. Gottlieb, D., & Shu, C. W. (1997). On the Gibbs phenomenon and its resolution. *SIAM review*, 39(4), 644-668.

https://doi.org/10.1137/S0036144596301390
**CRediT authorship contribution statement**

Igor Orynyak: Conceptualization, Methodology, Writing – original draft.

Danylo Katerynych: Software, Validation, Visualization

Roman Mazuryk: Algorithm, Writing– review &amp; editing

**Declaration of competing interest**

The authors declare that they have no known competing relationships that could have appeared to influence the work reported in this paper.