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.