Correction Factors for the Use of 1D Solution Methods for Dynamic Laminar Liquid Flow through Curved Tubes

: The modeling of transient ﬂ ows of liquids through tubes is required for studies in water hammer, switched inertance hydraulic converters, and noise reduction in hydraulic equipment. While 3D gridded computational ﬂ uid dynamics (CFD) methods exist for the prediction of dynamic ﬂ ows and pressures in these applications, they are computationally costly, and it is more common to use 1D methods such as the method of characteristics (MOC), transmission line method (TLM), or frequency domain methods. These 1D methods give good approximations of results but require many orders of magnitude less computation time. While these tubes are typically curved or coiled in practical applications, existing 1D solution methods assume straight tubes, often with unknown deviation from the curved tube solution. This paper uses CFD simulations to determine the correction factors that can be used for existing 1D methods with curved tubes. The paper also presents information that can be used to help evaluate the expected errors resulting from this approximation.


Introduction
Modeling of fluid flows through curving pipes has a long history, both as examples of theoretical calculus solution methods as well as for practical engineering.The interested reader is referred to two review papers [1,2] for details, but a brief overview is presented here.Consider the curved tube shown in Figure 1, with the radius of curvature defined as  , the inner radius of the circular cross-section of the rigid tube as , and the curve angle is .Boundary conditions include the pressure,   , and volumetric flow rate,   , at the inlet and outlet (denoted with subscripts A and B, respectively).The pressure is assumed to be constant across the inlet or outlet.Given the initial pressure and velocity fields within the tube and two of these four transient boundary conditions, the problem is to determine the other two boundary conditions.Geometry and nomenclature for a curved tube, showing the curved tube in plan view and also the cross-section.The radius of the cross-section is defined as , the bend radius as  , and the bend angle as .
Although others had previously noted the effect of curved passageways on an open channel flow and air through pipes, steady-state solutions to this problem typically start with the 1927 works of Dean [3], who established much of the theoretical groundwork and nomenclature used in the subsequent research, including this paper.He defined what would become known as the Dean Number, as where  is the Reynolds number and / is the curvature ratio.This number characterizes the nature of the flow via the ratio of the square root of the product of the inertial and centrifugal forces to the viscous forces [2].
The steady-state flow solution results in a secondary vortex flow perpendicular to the main streamwise flow, with the strength and number of these vortices related to the Dean number.It has also been noted that these vortices can have a significant effect on the flow for some distance downstream of the pipe bend [2].For the sake of this paper, we will not concern ourselves with the internal details of the flow pattern but will concentrate on the inlet and outlet flows.
Although theoretical solutions exist for some simplified situations, the most reliable relations come from empirical correlations.The early part of the 20th century resulted in much of the careful experimental data still in use today.This includes references [4,5], which resulted in Hasson's empirical correlation for steady-state flow, formulated as where / is the ratio of the steady-state resistance to the flow of a curved pipe relative to the resistance of a straight pipe with the same parameters ([6] via [2]).Note that the solution for the laminar flow through a straight pipe is the classical Poiseuille flow solution: where  is the fluid dynamic viscosity and  is the tube length.Refer to [1] for a more comprehensive list of suitable correlations for various situations.The above correlation is for laminar flow, noting that transition to turbulent flow has been observed to occur at higher Reynold numbers than for straight pipes.While the steady-state solution for the flow through curved pipes is relatively well understood, the dynamic response due to wave propagation under the influence of viscous friction is much less studied.Starting with straight tubes, analytical solutions taking into account the frequency-dependent nature of viscous friction exist in the frequency domain [7][8][9].These solutions can be calculated very quickly at a very low computational cost but are only applicable for problems where the entire system is linear and time is invariant.In situations where there are other parts of the system model that require nonlinear dynamic models (common in hydraulic system design), it is common to use either the method of characteristics (MOC) or the transmission line method (TLM).The method of characteristics [10][11][12] divides the tube into an evenly spaced 1D grid along its length.The flow and pressure at any grid point can then be calculated, moving forward through time with a time step selected so that the propagating wave moves exactly one spatial grid distance in each time step.This method is somewhat more computationally expensive than the frequency-domain method, but allows for parameters that change in time, as well as along the length of the tube, and also may be connected to arbitrary dynamic models at its inlet and outlet.However, the required equations have not yet been developed for curved pipelines, and it appears that this may not be theoretically possible.
Another method of modeling dynamic responses is using the transmission line method.Developed by a number of researchers over recent years [13][14][15], this method uses a four-port model (flow and pressure at inlet and outlet) and applies a number of parameterized linear transfer functions and delays to connect the four ports, ignoring the internal pressures and flows.The parameters of these transfer functions are selected to give the desired overall input-output characteristics, typically to minimize the error with respect to an analytical solution.This is an approximate method, but the errors can be small, and it has the advantages of being relatively computationally efficient and integrating well with dynamic modeling software packages such as Simulink [16], Hopsan [17], etc.This method was originally designed for simple, straight, rigid-walled tubes, but the author of the present paper has collaborated with others to expand the concept to allow for straight tapered tubes, tubes with elastic pipe walls with varying thickness, as well as arbitrary cross-sections [18,19].
There has been relatively little work carried out, either experimentally or theoretically, on the dynamic response of general flow through curved tubes.One series of experimental results is available [20][21][22], but the apparatus had a very large bend radius and was designed to accentuate fluid-structure interactions, so it is difficult to separate this effect from the purely fluid response.Some work has been performed on specific applications, such as blood flow through the circulatory system [23], but these results are sufficiently specialized with respect to the pressure signal limitations and fluid properties to be inapplicable in most other general applications.
Some recent works have also used 3D computational fluid dynamics (CFD) to solve the transient flow of fluids through curved and other complex geometries, such as [24][25][26][27][28]. Unfortunately, these methods can be unstable and require many orders of magnitude more computation time and resources to complete than 1D methods.
The main objective of the present paper is to invest time in CFD results to determine the correlation for correction factors that will enable the use of high-speed 1D solution methods for curved tubes.This paper will also present information related to the expected errors when using these methods with or without these correction factors.

Materials and Methods
The objective of this paper is to provide correction factors so that existing 1D solution methods for transient pipe flows may be used to approximate the flows in curved pipes.The solution for a straight, rigid pipe with a laminar flow of Newtonian fluid and with a constant wave speed can be entirely defined by three parameters: the characteristic impedance, dissipation number, and time scale.The shape of the overall response is controlled by the dissipation number, while the scaling of the magnitude of the flow-pressure relationship is set by the characteristic impedance, and the scaling in time is controlled by the time scale.In this section, CFD methods are used to establish the factors required to correct these values for use with flow through curved tubes.
For the purpose of this paper, assume a laminar flow through a smooth, rigid tube.This ensures no structure-fluid interactions or compliance in the tube wall and ensures a linear, time-invariant dynamic system.Assume that the sonic speed of the fluid is large relative to the fluid velocities and assume small changes in pressure, which will allow for the wave propagation solutions to be superimposed over a variety of mean flows and to be calculated for arbitrary boundary conditions.The pressure distribution is constrained to be constant across the inlet and outlet, but the pressure distribution within the tube is not constrained, nor is the velocity profile constrained at any point.

Characteristic Impedance
The characteristic impedance of a tube is the ratio of the magnitude of a pressure wave traveling through a tube to the magnitude of the resultant flow wave.For a lossless straight tube, this is related to the ratio of the inertia per unit length to the compliance per unit length and is defined as where  is the fluid density,  is the wave speed, and   is the cross-sectional area of the tube (a subscript "s" is used to refer to quantities related to a reference straight tube in order to differentiate them from an arbitrarily curved tube under consideration).Note that this value is only related to the properties of the fluid and of the cross-section; there is no dependence on the length of the tube.This impedance would be measurable at the inlet of a tube of infinite length or at the inlet of a tube of finite length for the initial part of the response prior to the influence of returning reflected waves.This is shown in Figure 2 for a Method of Characteristics solution for the inlet flow,  , in a straight tube of finite length subjected to a step in inlet pressure of Δ .The ratio   /Δ is equal to 1 in the middle portion (around  /) of this plot for a lossless tube (small dissipation number  .A similar method is used to measure the characteristic impedance of curved tubes, using response calculated using a CFD simulation.The CFD responses are calculated using the OpenFOAM CFD package [29].The "sonicLiquidFoam" solver was used, which is a compressible flow solver intended for studies of wave propagation in liquids with laminar flows.It uses an equation of state with a density described by where subscripts 0 denote reference conditions and  1/ is the (constant) derivative of density with respect to pressure.The default finite volume discretization of Gaussian integration is used.
The tube was meshed with the grid, shown in Figure 3, for half of the pipe (symmetry along the central plane is exploited to reduce the solution size).Following a similar study by Fries [25][26][27], the pipe's cross-section is divided into a central core with equal cell spacing, while the outer 2/3 of the radius is graded, with smaller cells toward the pipe walls to allow for better resolution of near-wall effects.In each cross-section, there are a total of 290 cells.This cross-section is then extruded in either a straight line of total length  or revolved about an axis over an angle of  (with a centerline length of   ).The cells in the axial direction are spaced with a 1 mm spacing along the center line of the tube (approximately equal to the spacing in the radial and circumferential directions at the edge of the central core).OpenFOAM's data structures are only concerned with the connections between cells, which permits the situation of  2, which is physically impossible as the rotations would overlap in space, but allows for a longer length of tube and therefore more accurate solutions that minimize inlet effects.This can also be thought of as a helical coil in the limit as the helical pitch approaches zero.The baseline parameters for simulations are shown in Table 1.The above tube was initialized with zero velocity and an initial pressure of 100 kPa.Boundary conditions were applied with a constant pressure of 100 kPa at the outlet and 200 kPa at the inlet, corresponding to an inlet pressure step of Δ = 100 kPa.A no-slip boundary condition was applied to the pipe walls and a symmetry plane boundary condition along the central plane.It should be noted that one could instead apply a flow boundary condition to the inlet, which can improve solution stability but would require some other method to determine the non-parabolic velocity profile of the inlet velocity.
The solution was then solved with a time step of 5 10 s for a period of 5 ms.This was selected as a time that would allow for the sonic wave to traverse the length of the straight tube approximately five times.(The change in the effective length of the curved tube will change the number of times that the wave actually transits the tube).
The volumetric flow was calculated by integrating the normal velocity over the plane of the inlet and outlet (and doubled to take into account the other symmetric half of the tube).The characteristic impedance was then calculated as the inlet pressure step divided by the flow, measured at the first time that the outlet flow equals the inlet flow ( in Figure 4).The characteristic impedance correction factor is defined as where  is the geometrical characteristic impedance (Equation ( 4)) and  is the impedance measured from the CFD results.The above was then repeated for a range of curve radii.These calculations were performed on a desktop computer with two Xeon 4114 processors (a total 20 physical cores) with 128 GB of RAM.ESI OpenFOAM version v2312 was used, running on Ubuntu Linux 22.04.4.

Time Scale
The second parameter required for the use of 1D solution methods is the time scale /.In this case, it is assumed that the wave propagation speed, , is constant and assumed that changes in the response are due to changes in the effective length  .This change can be explained by the fact that some wave energy may preferentially travel from inlet to outlet following the shorter path length along the inside of the curve.
The same CFD results from Section 2.1 were used to determine this effect.In this case, we identify the time it takes for a wave to travel the length of the tube, reflect off the open boundary, and return.As shown in Figure 4, for a sample waveform,  is defined as the first time that the outlet flow,  , exceeds the inlet flow,  . is defined as the first time that the outlet flow exceeds the inlet flow after  2  /.The effective length can then be calculated using and the correction factor is defined as where   is the centerline length.

Dissipation Number
The dissipation number for a tube reflects the damping induced by the fluid viscosity.This is analogous to the damping ratio of a dynamic system: a dissipation number near unity will reach the steady-state value with little overshoot, while a system with a low dissipation number will take many oscillations before reaching steady-state, as shown in Figure 2, for MOC solutions for varying dissipation numbers.
For a straight tube, the dissipation number is defined as which is equal to where  is the steady-state laminar resistance of the tube; this second equation applies to curved tubes as well as straight tubes.
The estimate for the effective dissipation number of a curved tube is determined by using the same CFD model as described in Section 2.1, but with the end time extended to 10  /.As before, the dynamic flow response to a step in pressure is calculated.The dissipation for a method of characteristics (MOC) solution (available for download from [12]) was determined for a straight tube that best fits this response.The characteristic impedance and effective length relations developed in the previous sections were used to compensate for these effects.The iterative, non-derivative global optimization scheme fminsearch, provided with the Matlab 2023a software package, was used to minimize the error between the MOC straight tube response and the CFD curved tube response by varying the MOC dissipation number.The normalized root mean squared error between the straight tube MOC solution,  , and the curved tube CFD solution,  , was minimized, and defined as where  10  /.The effective dissipation number correction factor can then be defined as where  is the effective dissipation number and  is the dissipation number for a straight tube of the same centerline length (Equation ( 9)).

Validation
In order to have some confidence in the CFD solution, the CFD solution for a straight pipe was calculated with parameters in Table 1, which can be compared to the validated method of characteristics and frequency domain solutions.These results are shown in Figure 5, along with varied viscosity solutions to show the effect of dissipation number.The results have good general agreement.One area with some differences is during the rapid change in flow around / / 1, 2, 3, and 4 as the wave fronts are reaching the reflective boundaries.This may be due to the assumption in the 1D solutions that the tube is narrow with no waves moving in the radial direction.Since the 3D CFD solution allows some wave energy to move diagonally through the pipe, some energy will arrive later, leading to the wave front's arrival being less sharp.A mesh invariance test was also performed to ensure that the mesh had been selected appropriately.The mesh size was halved in each dimension, along with halving the time step to ensure a constant Courant number.The results are shown in Figure 6.The Courant number was also checked by maintaining the original grid spacing while reducing the time step, as shown in Figure 6.Again, there is good agreement, and one can have confidence in the time and spatial mesh.Finally, the linearity of the transient solution was verified by doubling the driving pressure step magnitude, Δ .As shown in Figure 7 for curved tubes, the resulting flow doubles, but the shape is nearly exactly the same, with deviations only noticeable for extreme bend ratios and low dissipation numbers.This validates the linearity assumption and demonstrates that the solution does not rely on the Dean number (as it would for static solutions).The nonlinear static solution would be important for very large pressure waves, but the results presented here can be viewed as accurate for small variations about an operating point.

Characteristic Impedance
Examples of inlet flow wave forms for varying bend radii are shown in Figure 8.Note that as / →1 (more extreme bends), the magnitude of the flow between the steps increases, meaning the characteristic impedance decreases.There is also a change in the timing of the steps, which will be addressed in the next section.The associated calculated characteristic impedance correction factors are shown in Figure 9. Data for double the tube radius is also included, which is nearly identical, demonstrating that the curve for the characteristic impedance correction factor is a function of the bend radius ratio, not the radius itself.This curve can be approximated by a fit of For better accuracy, a higher order polynomial could be fit to the data in A1.

Effective Length
As shown in Figure 8, the curve ratio clearly affects the arrival time of the steps in flow (after they have traveled the length of the tube).The calculated effective length correction factor is shown in Figure 10, which is very similar to the correction factor for characteristic impedance.The correlation was fit to a similar form: It is unclear at present whether the difference in values in this equation should be assumed to be due to simulation error and taken as identical to Equation (13) or whether this difference exists in reality.
It should be noted that the  correction factor is not equal to the change attributable to following the shortest path length through the inside of the curve, indicating that the wave traveling along this shortest path is "pulled back" by the neighboring waves traveling longer paths.

Dissipation Number
Examples of the inlet flow responses are shown in Figures 11-13 for a range of bend ratios and dissipation numbers, along with the associated 1D MOC solutions with the best-fit dissipation number.These demonstrate a relatively good fit, especially for the parts of the response between the steps with moderate rates of change and also for moderate bend ratios and dissipation numbers.The more extreme bend ratios exhibit an oscillation in the volumetric flow as the wave fronts reflect at the boundary.There is no fluid experimental data available to validate this effect for large bend ratios, but these oscillations are visible in some data from curved microwave antennas [30], which are governed by similar wave propagation equations.This may be due to the nonplanar wavefronts that are generated in the strongly curved tubes.An example of this nonplanar wave front is shown in Figure 14, with the white contour line showing that the wave has progressed faster in the middle of the tube than along the outer bend radius.Figure 15 shows a plot of the pressure distribution along a radial line.This shows that the assumption of constant pressure across the cross-section that is commonly made for straight tubes is not applicable here.It also shows the fact that constructive interference has occurred with the maximum pressure considerably above the inlet's driving pressure.
The calculated dissipation number correction factor is plotted in Figures 16 and 17, along with the RMS error in Figure 18.As shown in Figure 17, the correction factor can be significant for strongly curved tubes, especially if the dissipation number is low.In most cases, the corrected response provides a usefully accurate estimation of the flow, although it does not capture the oscillations immediately after each step for strongly curved tubes.
To provide an indication of the relative computation time required for the CFD and MOC solutions, the typical time required to calculate a single response in this section was 1600 s for CFD and 190 ms for the MOC, a speedup of over 8000 times.In addition, the CFD fully utilized the 20 cores of the computer, while the MOC code was not parallelized

Overall Effect
One additional piece of useful information is how much of an overall effect a curved tube has, and under what conditions a straight tube approximation can be used without correction factors.In this case, one can calculate the error between the CFD-calculated responses for curved and straight tubes with the same cross-section and centerline length.The resultant normalized root mean squared (RMS) error is shown in Figure 19.As one can see from this plot, it is possible to make the straight tube simplifying assumption for more extreme bends if the dissipation number is large.For example, for zero viscosity, the relative error reaches 1% at  / = 16.2, while for  0.1 the same error can be achieved with the more extreme bend ratio of  / = 3.5.Another feature that can noted from this graph is that the dissipation number makes little difference for  less than 1 10 , which is similar to the effect in straight tubes for most applications of interest.

Conclusions
The results in the previous can be used to provide guidance with regard to a number of important questions that may be relevant to those employing 1D models to speed up transient flow calculations relative to CFD solutions.The first is whether the curvature of the tube in a particular application can be ignored and safely modeled as a straight tube.This can be evaluated by referencing the error calculated in Figure 19, quantifying this effect for a range of curve ratios and dissipation numbers.
If it is determined that this level of error is unacceptable, one can then use the presented information to obtain better approximations using correction factors.In particular, the characteristic impedance can be corrected using Equation (13), the effective length corrected using Equation ( 14), and the dissipation number corrected using Figure 14  While these results are applicable for small pressure changes and laminar flows, it is left to future work to establish the exact bounds of applicability with regard to these factors and establish additional methods of operating outside of these bounds.It is also left to future work to experimentally verify the presented results.

Figure 1 .
Figure 1.Geometry and nomenclature for a curved tube, showing the curved tube in plan view and also the cross-section.The radius of the cross-section is defined as , the bend radius as  , and the bend angle as .

Figure 2 .
Figure 2. Inlet flow responses to a step in inlet pressure for a straight tube with an open (constant pressure) outlet for varying dissipation numbers.The (left) and (right) plots show the same data at different time-scale zoom levels.For straight tubes, this response is totally defined by the dissipation number, , the characteristic impedance,  , and the time scale /, where  is the length and  is the wave speed.

Figure 3 .
Figure 3. Example meshing for a portion of a tube with an extreme bend radius.Note the inner core with equal grid spacing and the outer section with cells graded toward the outer wall.A symmetry plane bisects the tube.

Figure 4 .
Figure 4. Sample CFD flow responses for inlet ( ) and outlet ( ), showing points used to measure  and  .

Figure 5 .
Figure 5.Comparison of flow responses (inlet flow at left and outlet at right) for straight tubes for CFD (solid lines) and MOC (dashed lines).

Figure 6 .
Figure 6.Mesh invariance test (left), showing little effect when the spatial mesh spacing is halved in each dimension while also halving the time step to maintain the same Courant number.(Right) the minimal effect of the Courant number is shown by reducing the time step (while maintaining a consistent spatial mesh).

Figure 7 .
Figure 7. Linearity test, showing the CFD flow response for a variety of bend ratios and dissipation numbers for the base case of Δ = 100 kPa (black), and for the case of Δ = 200 kPa (colored).The results are very similar and overlap for almost all cases.The only visible difference is for extreme bend ratios and small dissipation numbers (top right of the plot).The error between the two solutions, , is calculated using Equation (11).

Figure 8 .
Figure 8. Flow responses to a step in inlet pressure for  0, while varying bend ratio.

Figure 9 .
Figure 9. Characteristic impedance correction factor for a curved tube relative to a straight tube of the same dimensions and wave speed.The same data is plotted in each graph, with the exponent applied to the left plot's x-axis to show linearity.

Figure 10 .
Figure 10.Effective length correction factor for a curved tube relative to a straight tube of the same dimensions.This also shows the same  data as Figure 9 for comparison.
r/r b and ran on a single core, with associated speed improvements expected, it was to run on multiple cores.

Figure 14 .
Figure 14.Example of the pressure distribution on the central plane after the wave has traveled clockwise from the inlet for 15 μs, for an extreme bend radius of  / = 1.2 and  0.1.Note the nonplanar wave front, denoted by the white pressure contour at 150 kPa.Also, note the overpressure above the inlet pressure of 200 kPa, plotted along the green line in Figure 18.

Figure 15 .
Figure 15.The plot of the pressure distribution along the green line in Figure 17, showing the nonconstant pressure distribution across the cross-section, which exceeds the inlet pressure of 200 kPa.

Figure 16 .Figure 17 .
Figure 16.Measured dissipation number relative to the straight pipe dissipation number.Each "X" indicates a CFD simulation location.

Figure 18 .
Figure 18.Residual error between corrected MOC solution and CFD solution after optimizing the dissipation number.

Figure 19 .
Figure 19.Error between the straight and curved CFD flow solutions for a range of dissipation numbers and bend ratios.

Table A1 .
Results for the characteristic impedance correction factor, as plotted in Figure9.

Table A2 .
Results for the effective length correction factor, as plotted in Figure10.