Thermodynamic Process Models

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.

Schritte:

The solver calculates the current states $T_2$ and $p_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.
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:

Forward Euler (simplest Runge-Kutta method)
Heun’s method (2-stage Runge-Kutta method)
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 cm³	5 ms
2	2 cm³	5 ms
3	1.6 cm³	5 ms
4	1.6 cm³	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 cm³ 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 cm³ 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 cm³ 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 cm³ 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

ode45
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 [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 cm³.

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

	ode45
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 cm³, step-size = 5 ms.

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

	lsode
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 cm³, step-size = 1 ms.

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