Modeling of thermodynamic processes in gas systems

A review of numerical equation solvers for models of thermodynamic processes based on a high-pressure gas system.

Early system verification is key to avoiding unforeseen problems during the development process. For this reason, IMT uses mathematical methods not only to check the feasibility of complex systems but also to identify possible design problems at an early stage. The procedure usually involves describing the system mathematically in process models and then checking it using simulations.

Having a mathematical representation of the system also allows IMT to speed up the controller development process. The controller is emulated and connected to the process model. This offers the possibility to work on the controller without having to generate the controller for the device and also allows us to gain insight into what the controller is doing with the resolution of a single step.

In the first part, we will look at the modeling of a few typical components for high-pressure gas appliances and derive the equations that describe the system dynamics. In the second part, we will look at different solvers and show examples of different solutions depending on various parameters to show their advantages and disadvantages.

The example shown in this article is a typical interface to a high-pressure gas connector as it can appear in various applications such as control of actuators (e.g., air pistons), pneumatic tools (e.g., pneumatic drills), medical industry (oxygen therapy, ventilators), industrial process plants, etc.

The model (see Figure 1) consists of a pressure regulator and a control valve downstream. The pressure regulator is connected to an external pressure source (subscript 1) and provides a constant input pressure (subscript 2) to a control valve downstream.

The pressure regulator has input boundary conditions \( p_1 \) and \( T_1 \) and we assume an output boundary pressure \( p_3 \) for the control valve. The output of the pressure regulator (subscript 2) is only ideally constant, but depends on the flow through the system, therefore it is a state of the system.

We assume the following simplifications:

• Ideal gas
• Both the pressure regulator and the control valve are isenthalpic
• There is no gas mixing, i.e., there is only one single gas in the system

These simplifications are reasonable for many applications. If we deal with high-pressure differences though, a real gas shows a temperature change when the pressure drops significantly (for example because of an orifice). In this case, calculating with an ideal gas will not show that temperature change. If we were to simulate an industrial process where air piping needs to be filled with nitrogen for preservation, gas mixing would be present, thus also requiring more complex models.

The models presented in the sections below are simple but serve the purpose of simulating a thermodynamic network accurately enough to help with the control system development.

Boundary conditions

The boundary conditions are assumed to be constant. At the inlet side, we have \( p_1 \) and \( T_1 \). At the outlet side, we only need the pressure \( p_3 \), since we do not have reverse flow (\( p_1 \) > \( p_3 \)).

Step 1: Mathematical description of the model

First, we look at each single component and derive the mathematical representation. The component inputs, outputs, and states are identified, and the equations are defined.

Depending on the type of component, the equations are either:

• First order ordinary differential equations (ODE): For components with continuous states

\( \frac{dy}{dt}=\ f(t,y,u) \)

• or simple explicit equations

\( y\ =\ f(t,u) \)

where t is time, u is the input and y is the state/output.

The valve

Given an input and an output pressure, a valve returns a flow (see Figure 2). A simple valve model without a modeled valve actuator has no states.

• Inputs:
  • \( p_2 \), \( T_2 \), \( p_3 \), \( T_3 \), \( Hub \)
• Outputs
  • \( {\dot{m}}_{valve} \)
• States:
  • none

The mass flow equation can vary depending on the valve type, but typically takes the following form for gases:

\( {\dot{m}}_{valve}=\ 2\ast\ k\ast \sqrt{\frac{p_2 \ast \Delta p \ \ast \ p_3}{T_2 }} \ , when\ p_3>{\frac{p_2}{2}} \)

\( {\dot{m}}_{valve}=\ k\ast\ p_2\ast\sqrt{\frac{\rho_2}{T_2}}\ ,\ else \)

Mit = \( k(stroke) \) ~ \( K_V\left(stroke\right) \)

The flow coefficient \( K_V \) is the standardized way to relate the pressure drop with the flow across a valve. It depends on the valve stoke. If \( K_V \) is not available or the actual curve \( K_V\left(Hub\right) \) varies from the design (e.g., due to washing out), k can be derived from measurements. Since the valve is isenthalpic and the gas is ideal, there is no temperature change, thus \( T_3 \) = \( T_2 \).

The pressure regulator
The pressure regulator (see Figure 3) reduces a (varying) input pressure \( p_1 \) to an output pressure \( p_2 \). Ideally, the output pressure is controlled to the setpoint pressure. The pressure regulator consists of a spring opening the valve and the regulated pressure generating a counterforce closing the valve for pressures above the setpoint.

• Inputs:
  • \( p_1 \), \( T_1 \), \( p_2 \), \( T_2 \)
• Outputs:
  • \( {\dot{m}}_{valve} \)
• States:
  • none

A simple approach to model the pressure regulator is to use the same equations as for the valve, but with \( k\ =\ k\left(p2\right) \) (see Figure 4):

\( {\dot{m}}_{reg}=\ 2 * k(p_2) \ast \sqrt{\frac{p_1 * \Delta p * p_2}{T_1}} \ , \ when \ p_2 > \frac{p_1}{2} \)

\( {\dot{m}}_{reg}=\ k\left(p_2\right)\ast p_1\ast\sqrt{\frac{\rho_1}{T_1}}\ ,\ else \)

For an outlet pressure above the pressure regulator setpoint, the flow is 0. For pressure below a threshold, the pressure regulator is fully open. In between, the flow coefficient is interpolated. This lower threshold should not be too close to the set point pressure, otherwise, it can be difficult for an ODE solver to calculate a stable and converging solution.

The volume

Between the valve and the pressure regulator, there is some piping with a volume that needs to be modeled. It receives flow \( {\dot{m}}_{reg} \) from the pressure regulator and has outflow \( {\dot{m}}_{valve} \) through the control valve (assuming the valve is open and the pressure after the valve is lower). Additionally, energy transfer to/from ambient \( \dot{Q} \) can be included.

• Inputs:
  • \( {\dot{m}}_{reg} \), \( T_{reg} \), \( {\dot{m}}_{valve} \), \( \dot{Q} \)
• Outputs:
  • \( p_2 \), \( T_2 \)
• States:
  • \( p_2 \), \( T_2 \)

Instead of the temperature, the enthalpy can be used as a state (for example if the gas is steam and thus steam tables must be used). With the ideal gas equation, the conservation of mass, and the conservation of energy, we can derive the ODEs:

\( \dot{T}_2=\frac{1}{m\ast c_V}\left(\left[-{\dot{m}}_{reg}\ast c_V-{\dot{m}}_{valve}\ast R\right]\ast T_2+{\dot{m}}_{reg}\ast T_{reg}\ast c_p-\dot{Q}\right) \)

\( \dot{p}_2=\ \frac{R}{V}\left(\left[{\dot{m}}_{in}-{\dot{m}}_{out}\right]\ast T_2+\ m\ast{\dot{T}}_2\right) \)

The heat transfer \( \dot{Q} \) , if required, can be modeled for example with natural convection, forced convection or even radiation (for very high temperatures).

Step 2: Make the system solvable

The equation for the valve (and also for the pressure regulator) to calculate the mass flow

\( \dot{m}_{valve}=\ 2\ast\ k\ast \sqrt{\frac{p_2 \ast \Delta p \ \ast \ p_3}{T_2 }} \)

has a problem for \( p_3 \) approaching \( p_2 \). The derivative goes to infinity

\( \lim\underset{p_3\rightarrow p_2} \ {\left(\frac{d{\dot{m}}_{valve}}{dp_3}\right)\rightarrow\pm\infty} \)

which means that a small change in \( \Delta p \) around 0 leads to a proportionally large flow change. When calculating the Jacobian matrix for the more advanced ODE solvers, this can lead to stability problems.

To overcome this issue, the equation can be linearized for small \( \Delta p \) :

\( {\dot{m}}_{valve}=\ 2\ast k\ast\sqrt{\frac{p_2 \ast \Delta p \ast p_3}{T_2}} , \ for \ \Delta p ≥ \Delta p_{min} \)

\( {\dot{m}}_{valve}=\ 2\ast k\ast\sqrt{\frac{p_2\ast p_3}{T_2}} \ast c \ast \Delta p , \ for \ \Delta p < \Delta p_{min} \)

with \( c \) such that there is a continuous transition between the two equations.
Other approaches with higher order polynomials are also possible. For example, a second order polynomial not only allows a continuous transition, but also a smooth transition at \( \Delta p \) = \( \Delta p_{min} \).

Step 3: Bring the models together

We now have a set of equations that can be solved. The valve and the pressure regulator provide mass flows based on the states and the boundary conditions to the volume, the volume calculates the state derivatives, and the ODE solver can then calculate the new states based on the derivatives. Figure 6 shows the workflow for solving the system.

Steps:

1. The solver calculates the current states \( T_2 \) and \( p_2 \).
2. The mass flows \( {\dot{m}}_{reg}\left(p_1,T_1,p_2\right) \) and \( {\dot{m}}_{valve}\left(p_2,T_2,p_3,\ Hub\right) \) are calculated.
3. The derivatives \( T_2 \) and \( p_2 \) are calculated.

This can for example be implemented in Simulink („Simulink 9.2,“ The Mathworks Inc., Natick, Massachusetts, 2018), see Figure 7. It shows the system with an Euler forward solver.

Step 4: Solve the system

There are two main types of ODE solvers:

• Fixed step-size solvers
• Variable step-size solvers

The fixed step-size solvers assume a constant step size and calculate the next step based on the current and (depending on the method) past or intermediate steps. The variable step-size solvers on the other hand adjust the step-size to ensure the per-step error does not exceed relative and absolute limits. These solvers are more stable as the step size is reduced in transient conditions, but that also means that the calculation time is increased in such situations. The amount of time required for the calculations can therefore vary greatly depending on the current system state.

For simple applications, a fixed step-size solver can be sufficient. For more demanding systems though, a variable step-size solver is often the better choice. The next sections analyze a few different solvers to show their performance and stability given different system parameters.

Fixed step-size solvers

We will compare three methods:

1. Forward Euler (simplest Runge-Kutta method)
2. Heun’s method (2-stage Runge-Kutta method)
3. 3-stage Runge-Kutta method

For the initial conditions, we assume \( p_{2,0}=\ 1\ bar \) and \( T_{2,0} =\ 25\ ^{\circ}C \) and as boundary conditions, we have \( p_1=\ 5\ bar \), \( T_1=\ 25\ ^{\circ}C \) and \( p_3=\ 1\ bar \). For the valve, we assume a constant stroke.

We will analyze the following examples:

Example 1
For the first example, we set \( V\ =\ 3 \ {cm}^3 \) and \( dt\ =\ 5\ ms \). The result is shown in Figure 8. The methods Heun and 3-stage Runge-Kutta show an almost identical result for both the pressure and the temperature, whereas the temperature for the Euler forward method is a bit higher during the transient behavior. Euler forward is the simplest but also the most inaccurate of those three methods.

Example 2

For the second example, we set \( V\ =\ 2 \ {cm}^3 \) and \( dt\ =\ 5\ ms \). The result is shown in Figure 9. With the smaller volume, the Euler forward method starts oscillating at the beginning. All three methods are stable and converge.

Example 3

For the third example, we set \( V\ =\ 1,6 \ {cm}^3 \) and \( dt\ =\ 5\ ms \). The result is shown in Figure 10. With the even smaller volume, Euler forward and the 3-stage Runge-Kutta method are still stable, but they don’t converge anymore. Heun’s method seems to give good results as the solution is stable and converging, however when looking at the temperature at 0.1 seconds, it is obvious that the result cannot be correct as the temperature should converge towards the inlet boundary condition temperature of 25°C.

Example 4

For the fourth and last example, we set \( V\ =\ 1,6 \ {cm}^3 \) and \( dt\ =\ 1\ ms \). The result is shown in Figure 11. With the reduced step size, the solution is now converging again to a stable solution.

The number of calls to calculate the derivative depends on the method, but for all methods, it grows linearly with the amount of time that is simulated. This is a drawback of those simple methods. Even though the system states do not change a lot anymore, the step size is not increased.simuliert wird, linear an. Dies ist ein Nachteil dieser einfachen Methoden. Auch wenn die Systemzustände sich nicht mehr viel ändern, wird die Schrittweite nicht erhöht.

Variable step-size solvers

We compare two different solvers

1. ode45
2. lsode

We assume the same boundary conditions as for the fixed step-size solvers. The volume is set to \( V\ =\ 1,6 \ {cm}^3 \).

ode45

The ode45 solver is a standard Matlab („Matlab 9.5 (R2018b),“ The MathWorks Inc., Natick, Massachusetts, 2018.) linear multistep method ODE solver. It retains information from previous steps to gain efficiency in solving the ODE system, contrary to the fixed step-size solvers introduced in the previous section.

The simulation result of the thermodynamic system is shown in Figure 12. At the beginning, the step size is at around 1 ms, but then it is subsequently increased to around 5 ms.

The table below shows the number of calls to calculate the derivative.

lsode – Livermore Solver for ODE

The Livermore solver for ODE (A. Hindmarsh, "Ordinary Differential Equation Solvers," 13 02 2008. [Online]. Available: https://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html. [Accessed 04 08 2021].) is a solver for stiff (Backward Differentiation Formula BDF) and non-stiff (Adams methods) problems with a Jacobian matrix that can be user-supplied or internally generated. The solver provides solutions at exactly defined points in time and retains information from previous steps for the next steps.

When solving an ODE system as part of a larger simulation environment, it is usually required to solve the system stepwise and have solutions for a given point in time. In a HIL/SIL environment, there is a discrete controller that drives inputs to the thermodynamic system (e.g., valves). This controller has a fixed step size and reads the inputs and sets the outputs every x ms. Therefore, the thermodynamic process needs to provide the results every x ms as well. The lsode algorithm can provide that functionality. For this simulation, the valve was kept at a constant stroke, so that the result can be compared to the other simulations.

The simulation of the thermodynamic system with the lsode solver can be seen in Figure 13. The simulation was programmed and solved in MS Visual Studio 2015. The simulation was run with a step size of 5ms. However, the solver may perform multiple sub-steps that are not shown. In the beginning, a lot of sub-steps are required to get a result that fulfills the relative and absolute tolerances, however, towards the end, the next iteration can be calculated with a single step only. The result might look coarse because the sub-steps are not displayed in the figure.

The table below shows the number of calls to calculate the derivative.

If the step size is reduced to 1ms (see Figure 14), the result resembles more the result of the system solved with the ode45 solver as shown in Figure 12, where a lot of steps with a small step size were made in the beginning.

Comparison

A comparison of the different solutions is shown in Figure 15. The solution of the lsode solver looks very coarse compared to the other solutions, but this is simply because the solution is only shown for the required step size of 5ms, thus only showing the solution every 5ms, but not the solution for all sub-steps.

The table below shows the number of calls to calculate the derivative for each solver.

In this article, we first looked at the modeling steps of a simple thermodynamic system. The modeling leads to several equations that easily allow the calculation of the state derivatives.

In the second step, we looked at several solvers – fixed step-size and variable step-size solvers. The comparison shows the clear advantage of the variable step-size solvers as they can adapt the step size to the system behavior. Furthermore, the lsode variable step-size solver shows a clear advantage over the ode45 variable step-size solver, requiring significantly fewer steps, thus solving the system a lot faster. The ode45 solver is furthermore bound to Matlab, but the fixed step-size solvers can be programmed in any language, and the lsode solver can be integrated as a library in most common languages.

The selection of a suitable solver depends on the system requirements and the system setup. The fixed step-size solvers are very simple to implement and not bound to any specific language/tool whereas the variable step-size solvers are bound to a language/tool or need to be integrated as libraries but are more stable and faster.

These methods allow IMT to analyze the feasibility of complex thermodynamic systems and identify potential design problems. Furthermore, it allows IMT to work on controller design efficiently even when a prototype is not yet available.