A Review of Numerical Solvers for Thermodynamic Process Models Based on a High-Pressure Gas System

Marcel Schoch, Senior Signal Processing Engineer
Reading time: 12 min

Insight in Brief

Thermodynamic process models are an important tool in the development of control loops and controllers for thermodynamic processes. They can be used to predict the system behavior regarding properties such as temperature/enthalpy, pressure or humidity. Incorporating such a model into a testing environment with the controller or the emulation of a controller allows to study the system dynamics and to design and tune controllers prior to working on the actual system.

This article deals with the following topics:

  • Mathematical description of basic components typically used in high-pressure gas systems
  • Challenges to solve such a system
  • Comparison of a few common differential equation solvers


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 speeding 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 to gain insight in 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.

Problem background

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.

Figure 1: Thermodynamic model with a pressure regulator, a control valve, and a volume in between.


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 accurate 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) \]


  • 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.


Figure 2: Pneumatic symbol of a valve with the model inputs (pressure and temperature) and the output mass flow.


  • Inputs:
    • p_2, T_2, p_3, T_3, stroke
  • 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 }} \ , if\ p_3>{\frac{p_2}{2}} \]

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

With = 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 be available or the actual curve K_V\left(stroke\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 counter force closing the valve for pressures above the setpoint.


Figure 3: Pneumatic symbol of a pressure regulator with the model inputs (pressure and temperature) and the output mass flow.


  • 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(p_2\right) (see Figure 4):

    \[ {\dot{m}}_{reg}=\ 2\ast k(p_2)\ast \ \sqrt{\frac{p_1 * \Delta p * p_2}{T_1}}} \ , \ if \ 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 \]


Figure 4: Pressure dependent flow coefficient for the valve model.


For an outlet pressure above the pressure regulator setpoint, the flow is 0. For a pressure below a threshold, the pressure regula-tor 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.


Figure 5: Volume with the model inputs mass flow, temperatures and heat flow and the states pressure and temperature.


  • 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 \ge  \Delta p_m_i_n \]

    \[ {\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_m_i_n \]

with c such that there is a continuous transition between the two equations.
Other approaches with higher order polynomial 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_m_i_n.



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.


Figure 6: Work flow to solve the system.



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


This can for example be implemented in Simulink [1], see Figure 7. It shows the system with a Euler forward solver.


Figure 7: Simulink implementation of the thermodynamic system with the 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 number Volume Step time
1 3 cm3 5 ms
2 2 cm3 5 ms
3 1.6 cm3 5 ms
4 1.6 cm3 1 ms


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.


Figure 8: System solved with fixed step-size solvers, V = 3 cm3 and dt = 5 ms



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.


Figure 9: System solved with fixed step-size solvers. V = 2 cm3 and dt = 5 ms


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. The 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.

Figure 10: System solved with fixed step-size solvers. V = 1.6 cm3 and dt = 5ms


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.


Figure 11: System solved with fixed step-size solvers. V = 1.6 cm3 and dt = 1 ms


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.

Euler forward Heun 3-stage Runge-Kutte
0.1 sec 100 200 300
1 sec 1’000 2’000 3’000
10 sec 10’000 20’000 30’000
100 sec 100’000 200’000 300’000


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.



The ode45 solver is a standard Matlab [2] 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.

Figure 12: System solved with Simulink ode45 solver. V = 1.6 cm3.


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

0.1 sec 289
1 sec 1’465
10 sec 13’111
100 sec 129’589


lsode – Livermore Solver for ODE

The Livermore solver for ODE [3] 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 is reading the inputs and setting 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 tol-erances, however towards the end, the next iteration can be calculated with a single step only. The result might look coarse be-cause the sub steps are not displayed in the figure.

Figure 13: System solved with the lsode solver. V = 1.6 cm3, step-size = 5 ms.


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

0.1 sec 221
1 sec 281
10 sec 302
100 sec 306


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.

Figure 14: System solved with the lsode solver. V = 1.6 cm3, step-size = 1 ms.



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 solu-tion every 5ms, but not the solution for all sub-steps.

Figure 15: Comparison of the system solved with different solvers.



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

Euler forward Heun 3-stage RK ode45 lsode
0.1 sec 100 200 300 289 221
1 sec 1’000 2’000 3’000 1’465 281
10 sec 10’000 20’000 30’000 13’111 302
100 sec 100’000 200’000 300’000 129’589 306


[1] “Simulink 9.2,” The Mathworks Inc., Natick, Massachusetts, 2018.
[2] “Matlab 9.5 (R2018b),” The MathWorks Inc., Natick, Massachusetts, 2018.
[3] 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].


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 a 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.

Are you interested in reading further articles?

Leave your e-mail address here and we will inform you as soon as we publish new contributions.

This might also interest you

Reading time: 3 min.

Ensuring successful outsourced engineering projects

How do you ensure success of your outsourced engineering project from the start?
Usability study for medical devices
Reading time: 3 min.

Usability study for medical devices

Preventing safety related use-errors or harm to the user using usability studies.
Reading time: 4 min.

Why event-based software architecture?

This article explains the event-based software architecture, its advantages, & possible disadvantages for the development...
Expert Blog Architekturprozess
Reading time: 4 min.

The systematic creation of a system and software architecture

This article deals with the development of a good system architecture...
Reading time: 5 min.

Additive manufacturing - reliable enough for medical technology?

For a customer's device, a dynamically rotating system was developed within a short period of time...
Reading time: 4 min.

Know what you don't know

Typical classification algorithms assign a class to every sample, even though that sample may be far off...