The basic MOC algorithm for a system

The Method of Characteristics - MOC

To get the full understanding, MOC should be explained in more detail. That's for a later post. For the time being I just refer to Wylie and Streeter and Wikipedia:

In mathematics, the method of characteristics is a technique for solving particular partial differential equations. Typically, it applies to first-order equations, though in general characteristic curves can also be found for hyperbolic and parabolic partial differential equation. The method is to reduce a partial differential equation (PDE) to a family of ordinary differential equations (ODEs) along which the solution can be integrated from some initial data given on a suitable hypersurface.

Explained below is the basic system algorithm for MOC in a system, not MOC itself. It's how it is programmed in LVTrans, and how it will be programmed in any eventual future versions.

System of hydraulic pipes

Each hydraulic pipe is modelled with MOC. A Pipe has a length and two ends. Each of these ends must be connected to a "non pipe" element, called NP. This NP element is typically a valve, a T-joint, a turbine, a reservoir but can also simply be an open or closed end of the Pipe. Thus for each Pipe, the NP represent the spatial boundary elements at each end.

Each NP can have one or more boundaries. These boundaries are typically Pipes, but can also be other NP elements. An example is a PID for a turbine. Many of these NP elements also have natural boundaries. For instance, a reservoir can be connected to one or several pipes but also have a surface level that is relatively constant seen in relation to waterhammer transients. This surface level is a natural boundary condition. The surface level can also vary, and can be measured and given as input for a digital twin.

The figure above shows a simple system of Pipes and NPs. Each Pipe has an arrow towards the NP elements. Why becomes apparent soon.

Basic element algorithm

Each pipe is divided into smaller sections. For each time step, the H and Q (head and flow) for each section is calculated based on the values for the previous time step. The boundaries are yet unknown, but the C+ and C- characteristics can be calculated. This is how far each Pipe can be calculated alone. All the internal values and the C+ and C+ characteristics at the boundaries. What is missing is the H and Q at each end.

Each NP element do almost the same. They start with the values from previous time step and the C+ and C- characteristics from the connecting pipes that just were calculated. All the inner variables and the H and Q at the boundaries can be fully resolved. These H and Q are then sent back to each pipe and make their calculations complete. Now the time step is fully resolved everywhere. The procedure is repeated for the next time step and so on.

Note that NP elements normally use Runge-Kutta or Newton-Raphson methods.

NP elements also can have manual inputs. For instance the opening of a valve. This is also updated for each time step.

MOC System algorithm

From the figure the full arrows are the C+ and C- characteristics at each end of the pipe. The back propagation of H and Q at the boundaries are the dotted lines. The PID does not send H, Q or C, but rather sends A and receives B, where A and B typically but not necessarily, are guide vane opening and turbine speed.

The MOC system algorithm therefore becomes:

1. For all Pipes and (PID).
1. Calculate the internal variables; H Q and predict the A in (PID) for T1 based on T0
2. Send C+, C- and A to each connected NP element
2. For all NP elements
1.  Calculate all internal variables and the boundary values for T1 based on T0 and values received from Pipe and (PID)
2. Iterate each NP-PID pair once more to correct the predicted A if needed*.
3. Send H, Q and B to each Pipe and (PID)
3. Update the boundary H and Q for each pipe and B for the (PID)
4. Start at 1 again for the next time step.

Programmatically there are many ways to to this depending on choice and the programming language. In "native" LVTrans these connections happens almost by itself. It's the way LabVIEW works by default. In a text based programming language this can be done by lists or arrays of Pipes and NPs. A major point is that the inner workings and transfer of values of each element is the same no matter how this is done. An independent library of MOC can thus be made. Here independent means independent from the actual implementation of the MOC system algorithm. 

*Note that the (PID) is not a Pipe with MOC inside. It does not have analytically correct C+ and C- characteristics for T1. The initial A it sends is a forward in time predicted value at T1 based on the time step T0 for all variables in (PID) including B. The B which is sent back from NP, is based on this predicted value of A. This B is used to update all variables inside (PID) and can be used to correct the predicted value of A at T1. This corrected value can be sent once more for a new iteration in the connected NP. Default is a predictor method. One extra iteration becomes a predictor-corrector method. An error analysis of this is found in the LVTrans manual, and is a part of choosing the correct time step.

It must be noted that the initial "error" comes from the fact that the B used when calculating the PID at T1 comes from the previous time step T0. This is also very much what happens in real life. The update frequency of a real PID is seldom better than 0.01 Hz for practical purposes. There's no reason it should be faster than this, because the physical frequencies for which it shall control is already several decades lower. What all this means is that with a time step in the simulation of roughly 0.05 s or lower, none of this matters - at all. This is also explained in detail in the LVTrans manual (I may also write a blog post at some time).  

In Wylie and Streeter the procedure is to gather all connected NP elements together and iterate them as one single block with regard to connected Pipe elements. In practice this is equivalent to an NP-NP pair iteration. It's a matter of choice, but the method of Wylie and Streeter may be more practical for blocks of NPs with more than two NPs. Such blocks never happens in LVTrans. With the exception of PID elements, it's probably more practical to create separate elements where larger collections of NP blocks are are pre-made.  

Conclusion

The MOC system algorithm is very simple. MOC is also superfast. Orders of magnitude faster than anything else while being analytically correct and with no numerical artifacts. With an object oriented approach for each elements, the same core library can be used in several different applications. This is very useful when different use cases are needed.

Ponderings about the future of open source transient simulation software

Intro

It's more than 25 years ago since LVTrans started. Lots have happened during that time. Neither AI or Python existed in some generally usable form 25 years ago, but today AI and Python are everywhere along with many other new and modern programming languages.

My first programming was done with interpreted BASIC on a Commodore 64. At computer courses at the university I learned Pascal and then FORTRAN. I liked them both. Then came C and C++ and ruled the world. I wrote a FEM software in C++ for my PhD, linking in FORTRAN matrix library. Compilers were expensive stuff. I used Watcom, which was the only one that really could combine C/C++ and FORTRAN. Today compilers are free and open source and can link in whatever other language you want.

C/C++ and FORTRAN

Even though FORTRAN has received more attention lately due to vector programming and AI/machine learning, it's still very specialized. Perfect for "hard numerical computing" but not particularly suited for anything else. C/C++ is still popular and can be used for everything. Compared with newer languages, C/C++ is complicated and not particularly suited for the casual programmer or engineers cheating in the field. I really like the cheer power, the brute force, of C/C++ but for the programming stuff I actually do, I always go in another direction, a simpler direction that gets the job done faster.

High level specialization

Another kind of programming has developed during the years, partly from academia, partly from the industry and goes way back. This is Matlab/Simulink, LabVIEW, Modelica, Scilab/Xcos and so on. Of these it's Labview first and foremost I have used but I have also used Scilab/Xcos a lot. Then there's Excel with Visual Basic/C or whatever "Visual" it has these days. It's a very powerful tool for a wide variety of tasks. All of them are perfect tools for the casual programmer and engineer who don't want to make an app, but rather want to get solutions to engineering problems.

Python

In comes Python. It has existed since long, but remained in the shadows for a long time. The largest effect in engineering and academia is that it has made Matlab etc. much less relevant. Matlab has still some considerable momentum due to its sheer size (and perhaps it's similarities with Python), but LabVIEW has gone backwards. It's no longer in the top 50 at the TIOBE index. Ladder Logic (PLC) and COBOL are much more popular. I never envisioned LabVIEW to become anything like C/C++ in popularity, but didn't really envision the ill fate it has received either. An open source program on a programming language that hardly exists anymore is a real tough hill to climb.

What to do?

The question is what to do. Does anything need to be done? LabVIEW is not likely to simply disappear from the face of the earth, but it will perhaps shrink in generality and focus more on the stuff it's actually designed to do. The choice of LabVIEW for LVTrans was more of a coincidence anyway. I was doing some LabVIEW work in the lab, and was asked to do a transient analysis of acoustics in gas risers. I said yes, and got myself a license of the commercial transient simulation software Flowmaster. It turned out that Flowmaster couldn't handle the particular problem, so I made a quick MOC solver for one pipe in LabVIEW that could. The Method of Characteristics (MOC) is very easy to program, for one pipe at least. The rest is history.

There have been several good points about this:

• LabVIEW is compiled to machine code, and is relatively fast. About 1/3-1/2 the speed of C in my experience. Sometimes even directly comparable to C.
• All graphics, graphs and charts are already there, professional looking and ready to use.
• The LabVIEW programming environment is used directly in setting up systems. It's very functional and easy to use, and additionally sets and organize the topology of the system.

This has enabled me to work as an engineer, doing lots of other stuff, while still maintaining and developing LVTrans, and using it in my work. There is no way this would have been possible in C/C++ for instance. Using Matlab (or Python for that matter) would be too slow execution vise. Last, but not least, LVTrans is also used as the numerical core in digital twins at Aker Solutions This is due to its speed and the native embeddedness of LabVIEW. It was just a matter of compiling it down to one of National Instrument's boxes, very simplistically speaking. 

It's safe to say that LabVIEW has served me well. It has enabled me to get the job done just perfectly. I have never been super happy with the solution though. The open source part has never taken off to put it mildly. Closed proprietary binary source code is no boost for open source (have to put in a "duh" here :-) ) 12-13 years ago I started making a C++ library. I never finished that, but I still have it for the core MOC solver. In the beginning, LVTrans actually had the core MOC solver programmed in (pure) C for 3-4 years. The long term plan was always to finish the C++ library and move on from there. Due to being busy with work, it never happened. Linking in the C++ library as DLL, LVTrans will be the same. LabVIEW would be used as the HMI only, while the numerical core could move on to other projects. It would be free and truly open source.

C++ or? 

Is C++ the right way forward today? In some ways yes, in other ways no. It's perfect for the job, but C++ is bloated and overly complicated for this simple job. There are other languages today that could be much better suited. Even professional programmers find C++ to be overkill for most of the work they do, or simply not suited and not efficient enough. I don't need efficiency as such, but there is a strong correlation between efficiency and having the right tool for the job.

This is where I'm right now. Is C++ the right tool for the job today, and if not, which other options are there. One of my sons is a lead programmer in IT now. He has opened my eyes to other stuff. Lots have happened in the last 25 years, lots more than just Python and AI. That's for another post.

Frequency Containment Reserve (FCR)

 Frequency Containment Reserve (FCR) is a new method and a service for containing the frequency on the grid during large and fast initial deviations.

FCR requires that the power-plant is qualified up front doing tests on site. This qualification is rather elaborate. It is also very costly because the plant will be off grid for roughly a week. A substantial amount of planning and preparations must be done to assure efficient tests, yet there is no guaranty up front that the plant is suitable. Statnett is responsible for this qualifications together with entso-e. A bit surprisingly, the qualification scheme describes no theoretical analysis, only on-site tests to be approved by Statnett.

Information about everything FCR can be found on Statnett web site here.

LVTrans is a transient simulation software. Since the plant is modelled in minute details, this enables transient simulations of all the FCR tests. Such a simulation will therefore be very useful. It will:

• Assess how likely the plant is able to run FCR, and if so:
• Determine the best possible PID parameter sets for FCR-N, FCR-D up and FCR-D down
• Determine the ranges of droop for FCR
• Determine the theoretical stability and performance of the plant when in FCR mode using the special FCP IT-Tool by entso-e/Svenska Kraftnät. The same tool as used in the tests on site.

Including FCR simulations using LVTrans will therefore reduce the economic risk associated with these tests dramatically. It will also enable the best possible starting point for the PID parameter sets. Less tuning on site. The low frequency sine-sweep of the tests can be dangerous for some plants since they may in a worst case scenario coincide with the resonant frequency of the waterway (surge shafts). Using LVTrans this will be analyzed and prevented.

All in all, using LVTrans up front as part of the planning and screening, will greatly enhance the possibilities for a successful, efficient and accurate FCR deployment, as well as prevent FCR deployment/tests on plants not suitable for it.

The FCR capability of LVTrans is not open source. It's much too special to have a broad interest, while at the same time requires a substantial amount of programming and testing to adhere to specifications. Interesting parties can get in contact for further information. As of today, four companies have purchased a license, enabling the development.

Links:

FCR in general

FCR and FCP IT-Tool

PID tuning - hybrid Skogestad

What could be called a hybrid, or combined, Skogestad method will consist of the following.

1. Find Ti according to Skogestad.

The plant should run at approximately 85% of nominal power on isolated grid. Set it in manual mode and adjust the guide vane until the speed is at or very close to 50 Hz (or 60 Hz as may be the case). The speed in pu should be close to 1.0.

Find the time it takes for the speed to reach 63% of the steady state value for a step in guide vane (injectors in this particular example, since it is a Pelton) of 0.1 pu from the value above. The steady state value after the step is noted as 1,059 pu. 63% of this change becomes a speed of 1.037 pu.

This gives a Ti of 7.1 s.

2. The PID is placed in auto mode on isolated grid and with this Ti = 7.1 s. Now adjust the SP n [pu] down and up again as in the previous trial and error method. Then increase Kp until the overshoot is somewhere between 5 and 20 % depending on how aggressive the PID should be.

In this example a Kp of 3.2 will do fine, but it certainly could be more aggressive.

That's it. The Ti is as good as it gets based on the physics of the plant, and the Kp is as good as one can hope for.

For tuning at idle no load, just set the P Grid [pu] to zero. This is usually done with only one injector.

A PID could be tuned by adjusting Kp until the overspeed is about 30% in step 2. Then increase Td until the overshoot is 5 to 20%. This is similar to the trial and error method. This last step, Td, has not been tested, but is in theory a straight forward extension.

PID tuning - modified Skogestad

As mentioned in the previous post, in some cases, particularly for servo positioning systems, it's a good idea to dive a bit deeper into the theories of tuning. One method I have found useful is the Skogestad method.

All the details can be found the linked paper. This is a modified form that I have found useful for hydro-turbines in combination with the trial and error method. The method consists of:

• A step in input, plotting the output and the input.
• Finding characteristic values from that plot.
• Using a single parameter called tau_c as tuning parameter.

The reason that this method is sometimes useful for hydro turbines, is that these systems can be very dynamic and oscillating. It can therefore be very difficult to identify the 15% and 30% overshoot from the trial and error method with any accuracy at all (although not entirely impossible). The Skogestad method (modified) will give usable values in these cases, at least for Ti.

The modifications to the method are:

• Disregard everything about time delay. Although one can often find something looking like a time delay, and it can be 2-3 seconds long, this is NOT a time delay in the normal sense, and should not be treated as one. Treating it as a time delay will in most cases produce unworkable PID parameters (instability or way too much overshoot).
• Use a tau_c of approximately 5 as a starting point (subject to tuning, see further down)
  

The disadvantage of the method, is it only works for PI. It does not work for PID unless approximating the system as a second order system (which a hydro power plant definitely is not, and a rather impossible thing to do in the field anyway). A good point of the method, is that Ti always comes out with a good value. It's really only tau_c (Kp in practice) that would eventually be adjusted to fit performance specs. A tau_c of 5 is found by me as a good starting point on idealized systems. It could however be way off on real systems.

Well behaved system

Looking first at a "well behaved" system. In hydro-power plant circumstances, this will mean a very simple system, an idealized plant.

The system is set in open loop, and a step in "kappa" (guide vanes) of 10% (0.1 pu) is done. This is seen as the green line, where we also can see the speed of the servo. The red line is the turbine speed in "per units", starting at 1.0. Here we can clearly see the speed going in the opposite direction during the first 2 seconds. This is fairly typical, but it's not really a time delay. The pink line shows the 63% rise in speed and the corresponding time for that rise (approximately 5.8 seconds).

This plot is all the information we need to tune according to Skogestad, which is rather neat. We find:

• Time constant tau_1 = 5.8 seconds
• Steady state gain k = delta(n)/delta(kappa) = 0.063/0.1 = 0.63

According to the rules of Skogestad, pretending this is some kind of first/second order system, disregarding time delay and setting tau_c = 5, we get:

• Kp = (1/k)(tau_1/tau_c) = (1/0.63)(5.8/5) = 1.84
• Ti = tau_1 = 5.8 s

That's it, super simple and super fast. Tuning this plant using the trial and error approach we get Kp = 1.9 and Ti = 6.3. From a practical point of view, this is more or less identical. Both will give fairly equal performance, and both well withing specs. It's not coincidental that a well behaved system has a generator inertia constant of approximately Ta = 6, and this corresponds well to the "optimal" Ti in this case.  

A realistic system

This is a real system. It's certainly on the "fuzzy" and oscillating side, but far from unusual. Doing the same as for the idealized system, we get the plot:

We find:

• Time constant tau_1 = 9.5 seconds
• Steady state gain k = delta(n)/delta(kappa) = 0.022/0.1 = 0.22

Then with tau_c = 5 as before:

• Kp = (1/k)(tau_1/tau_c) = (1/0.22)(9.5/5) = 8.63
• Ti = tau_1 = 9.5 s

This time Ti seems reasonable, and is in line with what the trial and error method gives. Kp on the other hand seems way off, and it is. It will produce a system behaving like this:

The system will oscillate forever, the PID is unusable. Note also the slow servo speed, needed for safety reasons on this plant to prevent water hammer. Since Ti is OK, we can simply use step 5 from the trial and error approach to get a working Kp. Doing this, I find that Kp = 1.4 is as "good as it gets", and still be within performance specifications. A step response in set-point of the closed system, certainly does not look good, but the physics of the system does not allow for anything better.

In the last plot below, I have introduced anti-windup. It doesn't improve the closed loop response, but makes a smoother servo run (less wear and tear) and will have a much better performance in terms of handling disturbances because Ti is reduced considerably. This is done by trial and error.

For real systems, some trial and error usually has to be done no matter what, but Skogestad method is very helpful in obtaining Ti for systems with difficult behavior. An open loop step plot will also reveal useful information of the system. In the next post, a combined approach is shown.

PID tuning the simple way

Tuning that works "every" time

Throughout the years I have tried every method I could find, from text book classical analytical s-plane methods to trial and error types. It's mostly turbine governors I have been working with, both in theory and in practice. The only method I have found that works almost every time, in theory and on site, is a trial and error method. Sadly I cannot find the copy of the original text, nor can I remember where I found it. I have written down the method though, ages ago.

The method requires that you are able to set the set-point arbitrarily and can read off the process value. A graphical representation is of much help rather than just numbers, although numbers works as well. The method is also described in the manual of LVTrans. The method requires no knowledge of the system, although such knowledge certainly helps.

The first thing to do is to stabilize the system. This can usually be done by setting Kp to 1 or less, and Ti to 10 or more and Td to 0. Then set the set-point (SP) to a desired value of which you are going to tune around, the operational point OP, and let the whole system settle. It's useful to have some integral value here to get the SP and the process value (PV) equal, but not strictly necessary. When this is done, and looks OK and steady state, the tuning can start.

1. Set Ti to a really large number, or deactivate further integration if possible. In parallel form, Ki can be set to zero. Keep Td = 0. Let the system stabilize if it isn't already.
2. Do steps in SP of approximately 10% up and down from the OP and observe the PV.
3. Increase Kp bit by bit until an overshoot of the PV of about 15% is observed. That is 15% of the size of the step in SP, not 15% of the total PV.
4. With the Kp from 3, decrease Ti bit by bit until the overshoot increases to around 30%.
5. For a PI governor: Decrease Kp bit by bit until the overshoot is 5-20% or looks OK - finished.
6. For a PID: Increase Td bit by bit until the overshoot is 5-20% or looks OK - finished.

That's it. An "optimal" PID will have some overshoot, but for some applications an overshoot may not be wanted. Then simply adjust Kp or Td in step 5/6 until little or no overshoot is seen. By coincidence these steps looks very much like the graphical animation in the Wikipedia article.

Note. This method makes it possible to tune most PID controllers in most systems, but doesn't necessarily create the "best possible parameters" whatever that may be. For turbines the specs for the governor are in frequency domain (a separate analysis) and the ability to follow changes in SP and changes in frequency. Any tuning that makes the system behave according to spec is good enough. If no formal specs exists, then this method results in no better or worse parameters than any other method. It does however create good and stable results almost every time, and with little knowledge of the system.

On a second note, the tuning of systems where servos saturate can be tricky. By saturation I mean the output is restricted in amplitude and speed. Some anti-windup algorithm should then be included for best behavior of the PID controller. This should be in effect for step 4 and 5/6 if it isn't already. Most commercial PID's have anti-windup implemented.

On a third note, a pure 1st order system may not give any noticeable overshoot in step 3. If that happens, then you know it's a first order system. A typical example is PID servo positioning. For such systems you can set Kp to almost "whatever you want", and adjust Ti to achieve a response with no, or desired overshoot. If performance is specified, you simply adjust to achieve the wanted performance. In this case, it is however a good idea to dive a bit further into the theories of PID tuning, starting for instance here.

Background

Tuning a PID governor is the process of determining the proportional gain, the integral gain and the derivative gain so the process behave the way you want it to. A key element is "behave the way you want it to".

If the stability is important, then there could be industry standards specifying the stability gains. Typically there are gain margins and phase margins that are standardized. These are obtained in frequency domain analysis, either directly or through FFT analysis. An industrial standard could also specify certain time parameters, how fast the PID should be, how accurate, the overshoot and so on. If no specifications are given, then the usual way is to get the PID so "optimal" as possible while keeping it stable.

There are therefore no universal descriptions of what a proper tuned PID actually is. It will be dependent on the application, the industry and the very process it is governing. The only common denominator is stability of the system as a whole. Theory and practice are usually different when working with PIDs, especially when it comes to derivative. This is sometimes referred to as realizable forms of the PID vs mathematical/theoretical forms. 

Two popular versions of the PID exists. The standard form and the parallel form. From Wikipedia the standard form is:

and the parallel form is:

where

As long as the relations between the constants (above) are remembered and taken into account, these two are equivalent mathematically. The tuning procedure shown above uses the standard form. The parallel form could also be used as long as the relations above are taken into account. Both these forms above happen to be mathematical forms rather than realizable forms, but that is of little relevance here.

Neglecting the convective terms

General

A thing that was done when deducing the governing equations was to neglect all convective terms. The details of this is rather elaborate, and involves lots of equations (different forms of Navier-Stokes equations). A full treatment is not done here, only the key points are included. The full analysis can be found in Appendix A of my thesis.

The basic principle is to make the Navier-Stokes equations from Part one non-dimensional in the same manner as was done in the Long wavelength approximation in the previous post.

Why would we neglect the convective terms in the first place? The reason is that they add considerable complexity in solving these equations. At the same time, if they happen to be unimportant for the solution nonetheless, it's only silly to include them.

We know that liquids are compressible. It's what makes transient analysis of liquid filled pipe interesting and important. However, we also know that liquids are only slightly compressible. Gas on the other hand is highly compressible. Can this compressibility be related in some way to some characteristic property, like for instance fluid velocity? Yes it can.

The key property describing the compressibility in relation to transient piping flow analysis is the Mach number. The Mach number M is defined as:

M = V/a

Where V is the fluid transport velocity and a is the wave propagation velocity (commonly known as the speed of sound). The speed of sound in water is typically around 1400 m/s, while the speed of sound in air is about 340 m/s at room temperature and atmospheric pressure. The (transport) velocity of the fluid in a pipe is for liquids around 1-10 m/s at max, while for gas it can be 100 m/s or more in some cases depending on application. For choked gas flow the Mach number is reached. It's safe to say that for liquids in pipes, the Mach number is very small, while this is certainly not the case in general for gas flow (although it could be the case in some specific circumstances as shown later).

Continuity

We make the NS equations non-dimensional carefully choosing appropriate characteristic properties. The continuity equation (in 2D cylindrical coordinates) pops out to become:

Here we see that the order of smallness for all convective terms is related to the Mach number squared. This is beneficial. If the Mach number already is small, a typical value is 0.004 for water in a pipe, then this squared becomes 0.000016, which is even smaller. A table for all the convective and space derivative terms is set up:

We see that the order of smallness for all the convective terms is very small, as long as M << 1, and they can be neglected. This means that when M << 1, then the convective terms can be neglected both for 1D and 2D equations in the continuity equation. Even if the Mach number is not super small, it has very little effect in the continuity equation. For typical flow in (large) gas lines, the Mach number is around 0.1 -ish, which means the order of smallness becomes 0.005 -ish (0.5 %) which is well inside the accuracy we can hope to achieve in a calculation.

Momentum

Doing the same for the momentum equation, the following non-dimensional equations in axial and radial direction pops out.

Now, this is unfortunate. The momentum equation is not related to Mach number at all. No matter what Mach number we use, the order of smallness of the convective terms becomes O(1) and cannot simply be neglected. There's nothing we can do about it, it's just the way things are. 

But, there are a few other things we can do. One thing that can be done, is in combination with the continuity equation to force the convective terms to be small. We say that the transport velocity is small and negligible. This is called the acoustic approximation, and is used for acoustic analysis (typically in 3D). In 1D this will simply produce the ordinary water hammer equations, but without the friction terms. A kind of side effect is that the ordinary water hammer equations can be used to calculate acoustics in gas pipes, and very accurately so, as long as the long wavelength approximation is taken into account. This was actually the first application of LVTrans, calculation of acoustics in gas risers.

What is done in my thesis is to apply the long wavelength approximation, L/D >> 1 to make the whole thing 1D. The continuity and momentum equations then become:

A.13 and A.14 is combined to obtain:

Now the Mach number is back, and we can compare the time derivative terms for smallness:

This shows that the in the combined set of equations for continuity and momentum, the residing convective term becomes insignificant when M << 1. This means that M << 1 does not in general remove convective terms in the momentum equation at all. It's only when applying the long wavelength approximation (1D flow) that the combined set combined will show the convective terms to be insignificant, and only when M << 1. the final non-dimensional momentum equation becomes:

Equations for (slightly compressible) gas flow

For the sake of completeness, it is worth considering the more general form of the non-dimensional equations where only the long wavelength approximation have been introduced (1D flow). We simply remove all the radial terms to obtain (removing all the subscripts for clarity):

These are two equations governed by two constants: The Mach number and the Reynolds number. It must be noted that the variation of density with respect to pressure is constant for these equations, the equation of state is over-simplified. The following non-dimensional variables are used in this post.

The long wavelength approximation

The long wavelength approximation is not as much a "mathematical trick" as it is a fundamental part of nature itself. The meaning of that is a the liquid in a pipe of some length compared to the diameter cannot only be treated as a 1D phenomenon, it essentially IS a 1D physical entity. Treating it in full 3D certainly can be done, but the 1D nature will completely overshadow any added wanted increase in accuracy from the 3D approach. The only valuable output from a 3D approach, at least for engineering purposes, will be the 3D approach showing it is essentially a 1D phenomenon.

This has all to do with inertia and the characteristic length and diameter of the pipe. A constraint is also that any externally forced frequency must be lower than first radial resonance frequency, which normally starts in the kHz range.

Why this is so is explained next. We look at equation 3.6 and 3.7 from part 1.

Equation 3.5 is not needed, because the constraint is frequencies below the first radial filled mode. In other words, what we are looking for is when we can use that constrain, and which parameters decides it. This is done using non-dimensional analysis, so we can see the relative magnitude of radial modes vs axial modes. The following non-dimensional variables are used:

These are inserted into B.1 and B.2 and we get:

These two equations are very similar. In fact, they are exactly equal, both mathematically and in magnitude (physically) when:

This means that if we have a pipe, and the L/D ratio is 1/2 (it's half as long as it's diameter), then radial and axial inertia will have equal effect on the dynamics of the system. Is a pipe where the length is only half as long as the diameter really a pipe? Hardly. So we introduce the constraint that L/D >> 1, or conversely that D/L << 1. We then see that the first term in Equation B.6 becomes very small (much less than 1), but nothing happens to the other terms in B.5 or B.6. The larger the ratio L/D becomes, the smaller the inertia term in B.6 becomes.

With L/D = 5, the effect or the radial inertia is already an order of magnitude less than the effect of axial inertia. Another way to see this is to derivate Equation B.5 with respect to r, and B.6 with respect to x, and combine them to get:

This equations shows the magnitude of the change of velocity in radial direction with respect to change in x and the magnitude of change of velocity in axial direction with respect to change in r. When inserting L/D  >> 1 (D/L << 1), the first term simply vanishes. Equation B.7 then becomes:

This equations say that the change of axial velocity with respect to r is zero, thus the axial velocity is constant across the cross section.

A key element of this analysis is of course non-dimensionality, which is a very powerful, and often necessary tool in all of engineering. Non-dimensionality enables us to extract the physical essence that we otherwise would not be able to see.

One question remains. What exactly is L/D in terms of a numerical value so that L/D >> 1? Is it 5? 100? 8.67? There is no simple answer. For instance, a complex piping system may have a few very long pipes, a few rather short ones and lots of pipes with lengths in between. We will then be more interested in the characteristics of the total system, than in each individual pipe. For the vast majority of systems, the longest pipes will decide the total characteristics. If the average L/D for all pipes is larger than say 10, then the 1D phenomenon will be the important one. If the average L/D is smaller than say 2, then this is hardly a piping system anymore, but a complex 3D shape and should be analyzed as such. The boundary conditions are also very important in this respect, as well as the reason for doing a transient analysis in the first place.

This is the case, not only for transient analysis, but also for steady state analysis. After a sharp 10 degree bend, the flow is not considered to be fully reorganized before after 10 diameters in length and so on. This doesn't mean that 1D analysis cannot be done. In most cases it only means that we have to be aware that the accuracy may not be as good as it normally would be. Another question is why would you do a transient analysis in the first place? As long as it is conservative and roughly accurate, it would be difficult to justify the added cost (and added other uncertainties) of a full 3D analysis, unless other considerations, typically FSI, also becomes part of the picture, and failure of the system is unacceptable.

Governing equation for transient pipe flow, part 4

From the last part we ended up with these two equations:

These equations are on a form that is suitable for FSI analysis when combining them with the equations of motion for piping systems (longitudinal, rotational and bending vibrations). They could also be used as is for pure hydraulic transients, but the normal starting point is to use pressure head, H, instead of the pressure, P, and velocity, V instead of flow rate Q (even though Q usually is eventually used in the end) and the diameter, D, instead of the radius, R.

Neglecting the Poisson coupling (sturdy piping structure), the last term in equation 3.19 and reorganizing we get:

Here the wave propagation velocity, a, becomes from the outset (for thin walled pipes):

For most cases this is a good approximation, but Wylie & Streeter have elaborated this for lots of different conditions in their books. Based on experience and practicalities however, the wave propagation velocity is a very fuzzy thing to determine. It can hardly be experimentally verified for any real piping system within +- 20%, mostly due to small, and changing, amounts of gas content, gas pockets and so on. Wylie & Streeter writes a lot about that too, and in the MOC (Method of Characteristics) it is used as a parameter that is changed within limits to discretize large systems. Even though it can be calculated theoretically at very high accuracy, that accuracy simply does not exist in real systems. Having a good starting point, is a good idea nonetheless though, because the difference between polymer based pipes and steel pipes can be more than 100% for instance.

What's missing now is the pipe friction. Diving into this, one discover it is super complex material. To this day, no good general transient friction term exists as far as I know. I have written a lot about it in my thesis. For steady state, Darcy-Weisbach and the Moody diagram is what is used. In lack of a better general non steady state model, Darcy-Weisback is also used for transient simulations, even though "everybody" knows this is far from correct. The best we can say about this, is that at least the calculations becomes conservative. Calculating water hammers that are smaller than in real life, does not normally happen. The exceptions are when FSI becomes important and for fluid cavity collapse. If this is not taken into account, then the calculations are usually not conservative anymore.

The Darcy-Weisbach equation:

Here f is the Darcy-Weisbach friction factor we normally obtain from the Moody diagram. The minus sign is there because of integration when positive x is in positive flow direction. This is simply the steady state result of the momentum equation above, and we can "reverse" it back into that equation to achieve exactly that. The final equations for liquid filled transient pipe flow analysis therefore becomes:

The wave propagation velocity, a, is as shown above.

Governing equation for transient pipe flow, part 3

 To get further from part 2, Hooke's law in cylindrical coordinates is introduced.

Then by combining Equation 3.13, 3.15 and 3.16 we get:

We have already discussed that the radial inertia term is unimportant (Long wavelength approximation). We can therefore substitute the tangential stress term for the quasi-stationary relation for a thin walled cylinder found in most text books, and assuming the pipe wall is much thinner than the pipe diameter.

A bit reorganizing, and Equation 3.17 becomes:

This is now fully 1D, and the stress term can be substituted using Hooke's law in 1D:

We now get the final fluid equations:

 

This is the same equations for transient pipe flow as can be found in any text book (when re-organized a bit), except the last term in Equation 3.19. This is the so called Poisson coupling of the fluid to axial motion of the pipe. The fluid wave propagation velocity is approximated by parenthesis before the second term (when multiplied by the fluid density), and is also essentially the same as found in most text books.

An interesting point is that when deducing the Water-hammer equations from Navier-Stokes equations, the FSI coupling pops out automatically, even for fully 1D equations. The only way to remove that coupling, is to force the axial movement of the pipe to be small compared with the pressure forces. We simply say that axial displacement of the pipe wall, u, is "small". In other words, it is assumed that the piping system is in practice 100% rigid. The other way is to set the Poisson ratio to be zero, which it isn't. The Poisson ratio of steel is a constant, and about 0.3.

For normal engineering purposes we simply say that the piping system is "rigid enough" for the transient fluid simulations to be correct. This is usually the case for the vast majority of systems. The only way the piping system is normally included, is in the calculation of the pressure wave velocity, see for instance Equation 3.18 and 3.19 where it can change the wave propagation velocity substantially. This will be correct also for thick walled pipes, except in cases where the pipe wall material becomes very "soft".

For the pipe itself, the axial motion ends up being described by:

The pressure is coupled to axial displacement by the Poisson ratio. When we say that the axial displacements are small compared with pressure forces, we can flip it around and say that the pressure forces are large compared with the axial displacements. Thus, vibrational analysis of the mechanical piping system without including the fluid transients, is therefore fundamentally wrong from the start. We see that the pressure forces are proportional to R/e, the radius of the pipe over the wall thickness. By increasing the wall thickness we will also reduce FSI effects, which fits well with intuition in this case.

For a full FSI analysis, junction coupling, bending and torsion also must be included, as well as all the clamping structures. It is a bit funny however, that transient simulations of the fluid inside the pipe can be done in most cases by disregarding the piping structure elasticity, but doing a vibrational analysis of the piping structure cannot normally be done correctly without including the fluid transients.

Governing equation for transient pipe flow, part 2

In the previous part, the so called extended water hammer equations was discussed (FSI), and a reduction of the Navier-Stokes equation was started.

Imagine a valve being slam shut at the left end. Water streaming from the right. The valve causes a pressure wave, a shoch wave (water hammer) propagating to the right at about 1200 m/s. The widening of the pipe wall causes axial stress, which will shrink the diameter through the poisson coupling. This causes the precursor wave propagating at about 5000 m/s.

The equations were 2D. We want them 1D. This is done by averaging the pressure and flow across the cross section. Mathematically this is done by multiplying with 2𝜋R, integrating with respect to r from 0 to R and dividing by 𝜋R^2.

V and P is the average velocity and pressure, and expressed as:

This is just a fancy way of saying V and P are the same all over the cross section of the pipe. This is of course not 100% correct physically, particularly not for the velocity V, but this is the simplification we must do to make them 1D.

But we had one more equation from the last part:

The assumption is that the characteristic fluid length scale, L is large compared to the characteristic diameter, D. That is L/D >> 1 under the constrain that the wavelength is "long". With this assumption, it can be shown that the radial inertia is negligible compared to the axial inertia. We can simply neglect the first term in equation 3.7. This is called "the long wavelength approximation", and is found to be correct for frequencies up to 63% of the first radial liquid filled pipe frequency (which normally is much higher than the first axial frequency). Proof of this is found in my thesis. Then we are left with the last term in equation 3.7 to be equal to zero, which simply states that the pressure is constant in the cross section. Thus "the long wavelength approximation" will also "one dimensionalize" the pressure, and remove one equation from the set.

On the pipe wall, the radial fluid velocity must be equal to the radial pipe velocity:

Inserting this, and we are left with only these two equations:

We have removed one equation. We have transformed them from 2D to 1D, but we have added a radial wall velocity term. This looks like a complication, but is a key point for the end result in part 3.