About Inverse Laplace Transform of a Dynamic Viscosity Function

A dynamic viscosity function plays an important role in water hammer modeling. It is responsible for dispersion and decay of pressure and velocity histories. In this paper, a novel method for inverse Laplace transform of this complicated function being the square root of the ratio of Bessel functions of zero and second order is presented. The obtained time domain solutions are dependent on infinite exponential series and Calogero–Ahmed summation formulas. Both of these functions are based on zeros of Bessel functions. An analytical inverse will help in the near future to derive a complete analytical solution of this unsolved mathematical problem concerning the water hammer phenomenon. One can next present a simplified approximate form of this solution. It will allow us to correctly simulate water hammer events in large ranges of water hammer number, e.g., in oil–hydraulic systems. A complete analytical solution is essential to prevent pipeline failures while still designing the pipe network, as well as to monitor sensitive sections of hydraulic systems on a continuous basis (e.g., against possible overpressures, cavitation, and leaks that may occur). The presented solution has a high mathematical value because the inverse Laplace transforms of square roots from the ratios of other Bessel functions can be found in a similar way.


Introduction
Water hammer is propagation of pressure waves in liquid-filled pipelines due to the change of flow velocity. The phenomenon of water hammer has fascinated generations of researchers in academia and industry. The high interest is reflected in an increasing number of scientific works on this subject. In most works [1][2][3], numerical solutions are used when analyzing this phenomenon. However, commonly used numerical methods have only been verified experimentally over the years. Hence, an analytical method would be very useful. It would enable the theoretical verification of these solutions in the full range of variability of the basic flow parameters. A detailed review of the literature carried out in Reference [4] indicated that a complete laminar flow analytical solution, which would correctly model relatively short-and long-lasting water hammer events in the time domain, has still not been presented (derived). Interestingly, such a solution is not known only for complex systems, but also for the simplest hydraulic system, such as a pressure reservoir-pipe-valve system ( Figure 1). As it has been found from the extensive review of the literature carried out for the purposes of this work, there is relatively little work strictly related to the analytical modeling of the water hammer phenomenon. Most of the theoretical work concerns solutions that have been determined by solving only one of the water hammer equations (equation of motion). Among such solutions, there are solutions for accelerated [5][6][7][8][9], decelerated [10], reverse [11], pulsating [12,13], and other types of flows [14][15][16][17][18][19]. Relatively recently, new solutions have been developed in which a variable viscosity has been introduced during the duration of unsteady pipe flow [20,21], and with the help of newly derived methods, the previously known laminar models can be extended to turbulent flows [22,23].
The pioneering research studies that concerned attempts to solve the system of equations of continuity and motion describing the problems of water hammer were the works of Joukowsky [24] and Allievi [25]. These authors did not take into account the friction forces; hence, they managed to obtain very simple analytical solutions. The first to become interested in friction was Wood [26], who derived a number of Laplace solutions for basic practical cases. His work was significantly expanded by Rich [27], who derived analytical formulas in Laplace and time domain for water hammer cases analyzed by Wood, i.e., with the assumption of no friction and friction varying in a quasi-steady manner. The unsteady (frequency-dependent) nature of hydraulic resistance was noticed firstly by Iberall [28]. Unfortunately, despite numerous attempts to account for the frequency-dependent friction into account in the late 1950s [29,30] and at the beginning of the 1960s [31,32], no one, except for Holmboe [33,34] and Zielke [35,36], succeeded in presenting a solution in the time domain. Unfortunately, Holmboe's solution was too simplified, considering only the first period of the water hammer event, and in addition, the author focused only on the cross-section at the valve. On the other hand, Zielke presented an inverse Laplace transform of the weighting function suitable for the numerical modeling of the wall shear stress, but he did not work on a complete analytical solution. In the 1970s, an interesting solution was presented by Karam [37,38]. This theoretical solution was limited to the step response at a downstream end in a semi-infinite fluid line, combined with a two-port representation of a finite line. The major feature introduced by Karam is two "filters" in a finite line model, representing a convolution of their arbitrary inputs with the unit impulse response at the equivalent location in a semi-infinite line. Fortunately, the promising Holmboe's model has received a decent correction by Muto and Takahashi [39].
The initial verification of this model presented in recent papers [4,40,41] showed its enormous usefulness, especially in modeling water hammer occurring in water pipe flows (relatively low-viscosity water supply systems). Unfortunately, in typical hydraulic systems, i.e., wherever the working fluid is hydraulic oil (water hammer number 0.05 > Wh > 0.5 [42]), computational compliance was not sufficient. Further attempts of the analytical solution were undertaken by Sobey [43], whose solution, unfortunately, does not take into account the unsteady friction, and Mei-Jing [44,45]. The Mei-Jing model is the other promising model. It has been recently corrected in Reference [4] so that the asymptotic pressure values coincide with the expected final pressure value, even in the range of high values of the water hammer number. In the same conference paper, the authors showed other possible mathematical forms of this solution that turned out to be mathematically close to Rich's quasi-steady model. The analysis of the Mei-Jing model showed its great potential, and, with a slight correction of the function describing the variability of the arguments, it seems that it will be (in the near future) the most accurate approximate model of the water hammer. It should be pointed out that recent research shows that the Atangana-Baleanu derivative and the Laguerre polynomial can be used to define a new computational technique for solving different types of differential equations [46]. The proposed scheme gives better convergence to the actual solution than the results available in the literature. It appears that use of the q-homotopy analysis algorithm in combination with the Laplace transform [47] can be a powerful tool to find new approximate solutions for this problem. In References [48][49][50][51], the authors proved that the finite element method is the most powerful scheme to compute the solutions of highly nonlinear ODEs.
Analytical solutions originally developed by Holmboe [33,34], which were later modified by Muto-Takahashi [39], are also the basis for the analysis of dynamic pressure and velocity waveforms, with the help of solutions commonly referred to as transmission line modelling. Interesting papers in this regard have been published by Oldenburger and Goodson [52], Brown [53], Viersma [54], and Ham [55]. These models are still very popular today, as exemplified by works of Hullhender [56], Johnston et al. [57], Krus [58], and Manhartsgruber [59]. The main role in these models plays the dynamic viscosity function. The inverse Laplace transform of it is the subject of this work. The presented method of the inverse Laplace transform of this physically important (in unsteady fluid pipe flow problems) dynamic viscosity function being the square root of the ratios of Bessel functions has not yet been presented in the literature (to the best knowledge of the authors). The novel method, after the expansion, will be helpful in finding similar inverse transforms from fractional ratios of Bessel functions.
Before we derive a new inverse Laplace solution, we briefly describe the problems of water hammer modeling. The necessary boundary conditions, which enabled Holmboe [33] to solve the problem in the Laplace domain, are here formulated. The frequency domain pressure solution is enhanced with the solution for the second basic parameter, i.e., the flow velocity. The complexity of a complete Laplace solution is the effect of nesting a square root of ratios of Bessel-based function in the arguments of hyperbolic functions. These nested functions are the main source accounting for the lack of a specific time domain analytical solution. The aim of this work is to show some new forms of the existing solutions from which it is possible to define a complete analytical solution in the Laplace domain. The presented solutions have not been discussed and published before. Our next step is finding an inverse Laplace transform of the dynamic viscosity function. The analytical time domain form of this function has not been investigated before. This dynamic viscosity function plays the most important role in the presented final solutions because it is responsible for dispersion and decay of pressure and velocity histories. A novel double infinite formulation for the dynamic viscosity function is derived. Inverse Laplace transform is performed with the help of the convolutional integrals' technique. In the new time domain solution, one can find a new type of unknown polynomial that is the result of the use of Calogero-Ahmed identities [60][61][62] for infinite series in the time domain, based on zeros of Bessel functions. In the near future, the novel solution will be used to find a complete analytical solution for any type of flow found in the analysis of unsteady pipe flows.

Water Hammer Equations
The equations that describe the water hammer phenomenon occurring in the simplest system presented in Figure 1 for the case of sudden valve closure are partial differential equations of the hyperbolic type. Holmboe [33], neglecting the dispersive term ∂ 2 v ∂x 2 , wrote the momentum equation in axial direction as follows: where v = v(x, r, t), and p = p(x, t).
The averaged continuity equation was written by Holmboe as follows: In the above continuity equation, the radial dependence was eliminated: and the mean quantities of pressure p and velocity v were defined. A typical phenomenon of water hammer takes place when the fluid flow is suddenly blocked (e.g., by the sudden closing of the valve (see Figure 1) at the cross-section x = 0), and then intense pressure pulsations occur with values significantly exceeding the values recorded during the steady state flow. There is a potential hazard: the possibility of unsealing, leakage, and, in the worst case, rupture of the pipe wall. When the valve is closed, the kinetic energy of the flowing liquid stream changes into pressure energy. At the upstream end of the closing valve, a wave of an increased pressure is generated which travels along the length of the pipe at the pressure wave speed, c. This wave, after reaching the pressure reservoir, is reflected and travels back toward the valve. After returning to the valve, it reflects again as a negative pressure wave, and then it travels to the reservoir, where the second rebound takes place, and the wave returns back to the closed valve. After the second arrival to the valve, the first period of the water hammer event is completed, and the phenomenon begins to proceed as at the initial phase, but under slightly different conditions, because during the first impact, the front of this wave was dispersed as a result of skin friction between the liquid and the pipe wall. In mathematical analysis, we need to know what happens at the beginning, i.e., before analyzed transient event. In our problem steady flow exists, a laminar flow take place, which velocity profile has a parabolic shape, as described by Equation (6). A linear pressure drop occurs at the same time along the length of the pipe.
Axial boundary conditions for the initial and transient laminar flow situations are as follows [33]: The initial condition of the velocity profile, at t = 0, may be given as follows: The radial boundary conditions (in pre-transient, as well as in unsteady flow) are as follows: (a) at r = 0; ∂v ∂r = 0; (b) at r = R; v(x, R, 0) = v(x, R, t) = 0 (no-slip condition).

Solution of Water Hammer Problem in Laplace Domain
Holmboe, after taking the Laplace transform of Equations (1) and (2), derived the pressure solution in the Laplace domain [33]: where Z = ρcv 0 is the Joukowsky pressure-rise formula. The equation for velocity solution was derived recently [4], with the help of Equation (7), and it reads as follows: In the above equations, the most important role is played by the dynamic viscosity function, M, which is a Laplace operator expressed in the following way [54,55]: where I 0 and I 2 are modified Bessel functions of the zero and second order, respectively. From Ham's dissertation [55], as well as from Brereton and Jing's paper [63], one can find the inverse Laplace transform of a function being an argument of this square root: where η i represents consecutive zeros of the first-kind Bessel functions of the second order, It is worth noting that the exponential infinite series contained in Equation (10) plays an important role in the numerical modeling of water hammer and, more precisely, in modeling the skin friction occurring between the flowing liquid and the walls of the conduit. This series was originally derived by Zielke [35,36] for the purposes of numerical calculations of the wall shear stress. To increase the convergence of this solution (limiting the number of terms), it can be written in the following exponential form [64]: where:t = t ν R 2 is a dimensionless time. In the literature on the modeling of unsteady pipe flows, it has been adopted as the weighting function. With the time and appearance of the effective solutions of convolutional integrals, the approximate solutions of this function have been shown by a number of authors in their works on numerical calculations of water hammer [64][65][66][67][68]. The use of the weighting function in the form of a finite sum presented by the formula in Equation (11) increases the efficiency of numerical calculations.
If the Laplace transform of the function described by the formula in Equation (10) is performed, and the following equivalent form of the sub-elemental function is obtained: The first two terms on the right-hand side of Equation (12) are quasi-steady solutions to the problem of fluid flow in a pressurized conduit (a linear friction model) that represent uniformly distributed resistance [58,69]. Thus one can write down the following: Let us graphically compare the courses of both complete solutions (Equations (9) and (13)) by taking into account the frequency-dependent friction and the quasi-steady one. For the purpose of the graphical representation of the obtained results, let us define the dimensionless Laplace transform operator asŝ = s R 2 ν . Then, by replacing the Laplace transform with the Fourier transform, i.e., substituting forŝ = jΩ, (where Ω represents the dimensionless frequency Ω = ω R 2 ν ), it is possible to graphically analyze the above solutions ( Figure 2). From Figure 2, it can be seen that the real part of the dynamic viscosity function, M, takes infinite values whenŝ → 0 ; otherwise, i.e., whenŝ → ∞ , it goes to 1. The imaginary part goes from minus infinity to the zero value whenŝ → ∞ . The quasi-steady solution behaves similarly asymptotically, although there are significant differences between the waveforms of the M and M q-s functions, especially in the dimensionless frequency range from about 1 to 1000. The Realis of the quasi-steady solution, as can be seen clearly from Figure 2c, takes the final value (≈1) already for Ω ≈ 30, while the unsteady solution has values greater than 1, still even for the frequency Ω ≈ 1000. The imaginary waveforms ( Figure 2d) confirm what is visible in the real-part waveforms, i.e., that both of these functions at different times tend to asymptotic (end) values. Hence, in Figure 2d, from the frequency Ω ≈ 10, it can be seen that these functions tend to the final values with a different dynamic.

Novel Method to Calculate the Inverse Laplace Transform of M Function
Let us use the following integral representation (special thanks for pointing this identity goes to Professor Alireza Ansari from Shahrekord University, Iran): where: and by working on the last term of the right-hand side of Equation (18), we achieve the following: where Γ(b, t) is the incomplete gamma function which for t = ∞ has a zero value, Γ(b, ∞) = 0, and one finally obtains the following: The above form based on the gamma function (Γ(k + 0.5) and Γ(k + 1)) of the solution is useful in our further derivations in this paper. The above Equation (20) can also be written in two other forms, respectively, dependent on binomial function and factorials: where n k is a binomial coefficient, which is defined as follows: We do not bother with the quotient of the gamma function, the square root of π, or the term 4ν presented in Equation (20) because they constitute a certain constant value depending on k. Thus, the problem of Laplace's inverse transformation is based on finding inverses from the following power function: The solution for k = 0 is δ(t) (Dirac delta function), while for k = 1, it is commonly known to be the following: For k > 1, the following approach based on convolutional integrals is proposed: To better understand what type of function we are dealing with after calculating the convolution, let us consider in detail a fairly simple case for k = 2. Let us simplify the problem by analyzing possible analytical solutions of this type of convolution: where a = µ 2 a ν The result of the above integral depends on the values of the coefficients a and b (always positive in our case). A simpler solution is obtained when a = b, as then we obtain the following: However, when a = b, one achieves the following: If, instead of a single product of exponential functions, we consider the product of infinite sums built on them, we get the following result for k = 2: The sum of two infinite series based on the differences of zeros of the Bessel function, i.e., the square bracket, can be written as follows: The infinite series of the above type was dealt with by Calogero [60,61] and then by Ahmed and Calogero [62] (the first author of this paper would like to thank Dr. Javier Garcia from the University of A Coruña, Spain, for pointing out this important relationship). In the first paper [60], Calogero proved that the infinite series based on Bessel functions satisfy the infinite set of equations: where h is the order of Bessel functions. The obtained results are related to the established connection between the motion of poles and zeros of special solutions of partial differential equations and solvable many-body problems. In the next paper [61], Calogero found that this infinite system (31), together with the Rayleigh sum rule, i.e., eliminates the scale invariance of (31) and characterizes the (squared) zeros of the Bessel function J h , namely, that the solution of (31) and (32) is unique. The author next proved that other relations of a similar type can be derived: Other solutions of this, as well as of a modified type (denominators multiplied additionally by µ a j,h ), were presented in a coauthored paper by Ahmed and Calogero [62]. For other interesting relations of similar types, please consult References [70][71][72][73].
Using the Calogero solution, Equation (32), the above sum for the zeros of the Bessel function of the first kind and zero order can be written as follows (in our case h = 0): The above identity significantly simplifies the obtained final solution for k = 2 in the time domain as follows: The obtained result can be easily verified numerically by going back to the Laplace domain with this solution: The solution for k = 2 multiplied by ν R 2 k has the following form: For the purposes of graphical comparative analysis, let us write it in a dimensionless form: In this way, it is possible to graphically analyze the obtained solution in the frequency domain ( Figure 3).
The above results of the simulation comparisons ( Figure 3) were obtained considering one million zeros of the Bessel function (l = 10 6 ). It was then noticed (Figure 4) that reducing the number of zeros taken into account resulted in a decrease in the convergence of the original function "Org" for large values of Ω. Such a decrease was not observed ( Figure 5) in the case of the new obtained solution "Ver"; here, in the analyzed range, i.e., up to the frequency Ω = 10 5 , the acceptable convergence was obtained (reVer graph, Figure 5a), taking into account only a thousand zeros (l = 10 3 ). This fact is undoubtedly influenced by the Calogero dependency.   The obtained solution in the time domain for k = 2 is the basis for finding a solution for k = 3. According to the formula in Equation (25), we should calculate the following convolution: The above notation means that we need to know the analytical solutions of the following integral: Consequently, it means determining the solutions of the following two convolution integrals: The solution to the first of the above integrals, B * I , when a = b is the function, is as follows: Meanwhile, when a = b, the solution is as follows: On the other hand, the solutions of the second integral, B * I I , are as follows when a = b: Meanwhile, when a = b, the solution is as follows: By considering their infinite series instead of single exponential expressions, the following solution is obtained for k = 3: The last term in the square bracket can be written as follows: Then we can obtain the following: Unfortunately, in the works of Calogero [60,61], as well as in the work of Ahmed-Calogero [62], one cannot find a solution for the function ∑ ∞ j = 1 appearing in the equation (Equation (49)). The following solution can be deduced from knowing that: In our case, more generally, we obtain the following: Moreover, we also notice the following: Thus, it is possible to use the formulas presented in the work by Ahmed-Calogero [62]: Then we obtain the following: Having the solutions of Equations (35) and (54), we can return to Equation (49) and derive its final time domain form: This equation reads in the Laplace domain as follows: By multiplying the formula in Equation (56)  , one can obtain a dimensionless solution: The comparison of the obtained solution for k = 3 (Equation (57)) is shown in Figure 6 below. Please note that, this time, the analyzed function is shown in semi-log graphs and not as in the case of k = 2 in a log-logarithmic scale. This change resulted from the nature of the variability of the analyzed solution; for k = 2, the real values did not exceed the abscissa axis, while for k = 3, such a drop occurs (see Figure 6a), hence the need for a different imaging and the use of a different ordinate axis scale. The above comparison ( Figure 6) was made by considering only one thousand zeros of the Bessel function. Further increasing the number of zeros of the Bessel functions taken into account does not increase the degree of agreement of the two functions ("Org" and "Ver").
For k = 4, it is easier to apply a slightly different method than for k = 3. This method differs from the proposed solution from the formula in Equation (25) in that our sought solution is determined as the convolutional integral from the solutions obtained for k = 2: Then we obtain the following: Remembering that the solution for k = 2 is described by the formula in Equation (36), we can search for solutions for k = 4 in the form of the following convolutional integral: The advantage of the above approach is that the final result in the time domain will be the result of three convolutional functions (not four as if one used approach compatible with Equation (25)). It also results from the analysis of the product of solutions for k = 2 in the Laplace domain, by means of which the following can be shown: The solution (60), as well as the return of it to the time domain with the help of the Equation (61), results, as can be seen, in the appearance of three convolution integrals. Assuming at the beginning, as for the cases of k = 2 and k = 3, the consideration of integrals for single exponential terms, one can write the following: where a = µ 2 i ν The solutions for these integrals are as follows: (c) C * I I for a = b (e) C * By considering their infinite series instead of single exponential expressions, the following solution is obtained for k = 4: Using the Calogero [61] formulas (Equations (A5) and (A8)) collected in Appendix A of this paper and the formula in Equation (54), one can simplify the obtained equation in its final form (detailed mathematical steps are summarized in Appendix C Equation (A20)): The Laplace transform of the above solution gives the following result: Multiplying the obtained formula by ν R 2 4 makes it possible to write this solution in a dimensionless form: The comparison of the obtained dimensionless solution for k = 4 is shown in Figure 7. The last solution that was derived in this paper concerns the case when k = 5. In order to obtain this solution, the convolutional integral from the solutions for k = 4 and k = 1 was used. This solution was obtained with the help of the formula presented in Equation (25) (it is possible to obtain the same final solution as a convolution from solutions for k = 2 and k = 3). The detailed procedure of obtaining the final result is similar to the one discussed in the previous cases (k = 2, 3, and 4); hence, it is not discussed here again. For the sake of clarity, however, it is worth writing down (see Equation (A21) in Appendix C) all the obtained infinite series of Calogero-Ahmed type from which the final coefficients of the last terms were determined.
It is possible that such knowledge will enable us to propose a general solution in the near future. In order to obtain the final result (and in particular the coefficient for the last term), it was necessary to apply a number of formulas for infinite series of the Calogero-type (Equations (2), (4), (5) and (14)- (18) from Reference [62], for which the values are summarized in Appendix A (in our case for h = 0)). Finally, the following result was obtained: Moreover, in the Laplace domain, we obtain the following: Multiplying the obtained formula by ν R 2 5 makes it possible to write this solution in a dimensionless form: The comparison of the obtained dimensionless solution for k = 5 is shown in Figure 8. The further determination of successive solutions for powers greater than 5 becomes problematic in terms of manual implementation, which is associated with the need to use more and more complex Calogero-Ahmad series solutions. These solutions are the result of successive convolutions. In addition, the literature lacks solutions to these series, so it becomes necessary to constantly determine new formulas or find a relationship that would enable the analytical determination of their mathematical forms.
The most important advantage of the presented novel technique is the possibility to fully analytically describe the dynamic viscosity function in an extended time range. The inverse Laplace transform of this important function was developed for the first time (to the best knowledge of the authors). With its help, it will be possible to find a complete solution to water hammer equations (investigations are in progress to achieve this goal) and to define accurate approximation solutions of these phenomena (traditionally preferred by practitioners due to their mathematical simplicity). Another advantage of the presented technique is that it can be easily extended (solution for larger k can be found) with the help of a machine-learning technique. Initial work is currently underway on the implementation of such a computational algorithm.

Comparison of M Function Transform in Time and Frequency Domain
The knowledge of the convolutional solutions for the first power expansions of the analyzed dynamic viscosity function allows us to present a collective solution. Looking at the constants that appear in the new notation of the analyzed function (double infinite summation form, Equation (20)), it can be seen that they are based on the values of the gamma function that can be expressed as a definite integral for (a) > 0 in Euler's integral form [74]: A time domain solution of such rational equations, i.e., was found in the previous section as an inverse Laplace transform for the first five powers of k (for collected successive solutions of this type, see Appendix B). Knowing the above solutions and the values of the constants Γ(k+0.5) Γ(k+1) · 4 k √ π calculated for the first five powers of k, we are able to write the final form of the dynamic viscosity function in time domain, i.e., the inverse Laplace transform of our searched function: Interestingly, the above sum can also be written in a different (shortened) form by ordering the coefficients with appropriate sums: Let us compare our obtained solution of the inverse transform of the dynamic viscosity function, M, with the inverse Laplace transform of the quasi-steady form of this function, M q-s , i.e., with the inverse transform of the function described by the formula in Equation (13): Now we are going to investigate asymptotic solutions, as well. With the help of the generalized Puiseux series, one can find a series expansion of the dynamic viscosity function, M (Equation (9)), for smallŝ → 0 : An inverse Laplace transform of the above function gives an approximation in the time domain that is valid for relatively large dimensionless times: where H(t) is a Heaviside step function, andt = t ν R 2 . Moreover, for largeŝ → ∞ , one can expand the M function by using the Hankel asymptotic expansion of the modified Bessel functions (p. 377 in Reference [64]): The use of the above expansion in our case results in the following: The series expansion for z = ∞ finally gives the following: The above expansion is similar to the one presented by Holmboe [33] and by Brown and Nelson [75]. Substituting for z = √ŝ , one obtains the following: The inverse Laplace transform of the above solution gives an approximate form (correct for relatively small dimensionless times) of the dynamic viscosity function in the time domain: In order to emphasize the influence of the number of power expansion terms on the fit of the obtained results, as presented in Figures 9 and 10, a comparison is made here for two cases of analytical solutions, i.e., for k = 0 . . . 3 and k = 0 . . . 5. Figure 9 shows a comparison of the histories obtained in the time domain (discussed briefly above). Meanwhile, in Figure 10, the solutions for k = 0 . . . 3 and k = 0 . . . 5 are compared with the output function (Equation (9)) and its quasi-steady counterpart (Equation (13)). Figure 9 compares the solutions in dimensionless form, i.e.,m = m R 2 ν .  The presented comparisons show that the use of the analytical solution based only on the first three expansions of power k = 0 . . . 3 is not sufficient to describe the analyzed dynamic function. The approximate form of the solution in the time domain, which is based on the first five power expansions k = 0 . . . 5, gives a much better fit. It seems that the next five exponential terms should guarantee a smooth transition in the time domain to the solution that is the inverse transformation of the dynamic viscosity function expansion into the Puiseux series (red dotted line). Knowledge of this transition will allow us to approximate a complex analytical solution and to propose a simplified analytical solution of water hammer. The presented comparison shows that the inverse Laplace transform with asymptotic expansion (blue line of line graph) is suitable for use in the dimensionless time domain in the range of 0 <t < 10 −1 . For dimensionless timest ≥ 1, this function can also be correctly modelled with the help of the quasi-steady solution (line of large green dots). The transitional range of this function (10 −1 ≤t ≤ 10 0 ), in which significant dynamics of its variability are observed, is the range that is responsible for the water hammer in oil-hydraulic systems (0.05 ≤ Wh ≤ 0.5). Hence, it is important to further extend the analytical solution proposed in this paper.
The obtained analytical solutions based on the initial power expansions k = 0 . . . 3 and k = 0 . . . 5 were also compared with the exact solution of the analyzed function in the frequency domain ( Figure 10). Both the comparison of the real part ( Figure 10a) and the imaginary part (Figure 10b) show that successive expansions extend the compatibility range of the new solution. Therefore, it is worth extending the analyzed analytical function (by successive powers of k) so that it represents a high accuracy in the range of dimensionless frequencies significant from the practical point of view, i.e., 10 0 ≤ Ω ≤ 10 1 .

Conclusions
The novel solutions presented in this paper are an important first step in the development of a complete analytical solution of water hammer equations; however, we consider them to be worthy of sharing with a wider group of researchers. A complete analytical solution can be obtained with the help of the approach presented in this paper through machine learning, as well as with the help of artificial intelligence (AI). To develop such solutions, it is necessary to generalize the series based on the zeros of Bessel functions of the Calogero-Ahmed type. To this day, we know too little about these infinite series. Their polynomial solutions for higher power expansion of analyzed dynamic viscosity function looks to be the key point to completely describe the theory of such Laplace inverses. The extended final solution that will be obtained in the future can be approximated by finite exponential series, exactly those that are currently used in weighting functions for effective numerical calculations of the wall shear stresses in unsteady pipe flows.
Further investigations aimed at finding the final simplified solution to the problem analyzed by the authors in this paper are underway. We should learn whether the solutions presented in recent works in the field of mathematical physics that were kindly proposed by the reviewers would help in the development of the simplified analytical solution of unsteady pipe-flow equations [76][77][78][79]. The initial attempt to find/formulate a generalized analytical formula for water hammer equations in the time domain by using the dynamical viscosity function shows that such a solution will probably be similar in form to the differential form of describing the fractional Bessel functions, as well as a derivative equation that describe the Laguerre polynomials (see Equation (25) in page 319 in Doetsch's book [80]). However, the missing important part will be the definition of new polynomials that appear during calculation (marked red terms in Appendix B of this paper).
The presented and discussed new method seems to be very promising, as it can be used not only in fluid mechanics but also in other basic sciences: mathematics, materials engineering, thermodynamics, physics, and medicine; it can also be used generally wherever there are solutions based on the similar functions discussed in this paper, i.e., exponential functions and series of the Calogero-Ahmed type.
We do hope that researchers will be able to further simplify the proposed time domain analytical solution. Perhaps they will also manage to simplify the final form, that is, to collapse our solution to the final form; in other words, they need to complete the opposite task to the one performed in Laplace domain, where we decomposed the original function into a doubly infinite series.

Appendix A
The Calogero-Ahmad series for zero-order first-kind Bessel functions (h = 0) used in this paper are as follows: (A4) (A7)

Appendix B
(a) Solutions in the time domain In the above solutions the values of final coefficients were found by using the Calogero-Ahmed identities [61][62][63] discussed in these papers. The coefficients highlighted in the red color form unknown polynomial series. If completely described, they should be termed as Calogero-Ahmed polynomials.

Appendix C
This appendix shows the intermediate forms of the selected solutions. The purpose of showing them is to give the reader an insight into how to find a mathematical formula for the higher powers of the presented analytical solution: