Symbolic mathematics-based simulation of cylinder spaces for regenerative gas cycles

Introduction

The process configuration of a regenerative gas cycle always consists of a heat exchanger circuit and a combination of cylinder piston assemblies. In the sealed machines the working gas is warmed up and cooled down alternating by the heat exchangers to generate the pulsating process pressure. The pistons are mechanically connected to a gear; during a power stroke they push the gas through the heat exchangers, and they generate mechanical work by changing the whole gas volume. A configuration therefore consists of the elementary components which the working gas flows through in a certain order:

• cylinder space (displacement or working cylinder),

• heat exchanger (cooler or heater for external heat transfer),

• regenerator (internal heat transfer).

For the optimal design of these components different methods for modeling and simulation are used (cf. [ ⁴ ]). The differential equations for the conserved quantities: mass, energy, and momentum are usually not directly solvable, i.e., do not have closed form solutions. Instead, solutions must be approximated using numerical methods. For purposes of simplified calculations, one may assume that the process takes place from initiation to completion without a change in a special state quantity of the system. In this case also analytic solutions exist, such as the Schmidt model for Stirling engines [ ¹¹ ]. One more model that has occasionally been used is the model for the regenerator due to Hausen [ ² ].

For the symbolic mathematics-based simulation of the gas cycle the single components can be modeled separately, where every component has an independent solving method for itself. One advantage is that special boundary conditions for a component can simplify modeling to receive an efficient solution of the underlying differential equations. Another advantage is that it is possible to simulate a higher number of different process configurations (e.g. Stirling or Vuilleumier process [ ⁴ ]) with this modular construction. Thereby the components are arranged according to their number and location in the flow. The nK components are coupled by the whole energy balance of the process:

The process pressure sets again a boundary condition for the single components. Furthermore the components are coupled by the energy flow of the involved masses as well as without mass via the surrounding.

Simulation models for the component of the cylinder space exist for regenerative gas cycles and also for reciprocating compressors or internal-combustion engines. Analytic solutions for the cylinder space may be obtained by assumption of an adiabatic-, polytropic-, or isothermal process (cf. [ ¹ ]), as with the Schmidt model [ ¹¹ ]. In contrast, the studies of multidimensional effects like heat transfer [ ¹³ ] or the pressure loss are solely based on numerical methods. Other models are also using numerical methods simulating the cylinder spaces for piston engines based on ordinary differential equations, like [ ⁵ ] for Vuilleumier engines, [ ⁸ ] for diesel engines, or [ ¹⁰ ] for Stirling engines.

The analytical or semi-analytical methods are only partially suitable for aims of simulation-based optimization because gross simplifications for the real situation must be superinduced. The numerical approximations do not deal with this problem but with an analytical model the solution is held in form of an equation. It will then be easier to handle and one can clearly recognize the characteristics of the physical process. Additionally, for the simulation of regenerative gas cycles it is possible to calculate directly the periodical-steady state of the system. This paper aims at exploiting the benefits of an analytical calculation method without using gross simplifications for the modeling.

To determine the change of the state quantities one has to solve the differential equations for the conserved quantities under respect of the boundary conditions. In this connection, the change of the cylinder volume and the gas density gives the required information for the bordering components, like the emitting and receiving gas mass flow. For the multidimensional problem numerical methods may approximate the solution of the differential equations and perhaps yield useful information on analytic solution. In this case the investigations of the structural intrinsic properties of the system may be used, to provide the information for

• the determination of the influence parameters of the cylinder space,

• their effects and especially the influence of individual parameters change,

• the simplifications that result from fixing certain parameters for the simulation, and

• the conditions for the isothermal process, which is widely accepted in the field.

Finally, for simulation purposes analytic solutions are generally considered superior to numerical approximations.

Energy equation

A single cylinder is connected with the remaining cylinder spaces of a process configuration by heat exchangers or by regenerators. The influencing quantities of the energy balance of a semi-open system are shown in Fig.  ¹ . These are the transported heat and work via the surrounding and the enthalpy flow rate to connect the other modules by the interface x :

Fig. 1

It has to be distinguished whether it concerns a gas taking or gas delivering cylinder space according to the temperature at the interface. In the first case, the temperature is a function of the bordering component: the heat exchanger. In case of the gas delivering cylinder the temperature at the interface corresponds to the temperature inside the cylinder space. Therefore:

• Tx=Tin(φ) , a transient function for the gas taking up cylinder

• Tx=T , the temperature inside the cylinder space for the gas delivering cylinder

The working gas within the cylinder is considered as ideal mixed. The reason is that the gas streams alternately incoming and leaving the cylinder provide in each case an intense turbulence.

Furthermore, no basic differentiation of the energy balance is made between a displacement and working piston because that has only a constructive meaning. The appearing leaky mass flow about the piston seals (which can be predicted mathematically only inexactly) is extremely low at a displacement compressor due to the lower pressure differences, and the effects for the process can be neglected.

The energy equation for the semi-open system (cf. Figure  ¹ ) is expressed by:

With the thermodynamic work Pdt=-pdV to change the cylinder volume, the internal energy U=mcvT given by the total derivate and the enthalpy flow rate given by (2) Eq. ( ³ ) can be transformed into:

For regenerative gas cycles the working gas can be treated as an ideal gas, and with the differential form of the ideal gas law

Substituting (2) and (5) into (4) and applying the mass balance of the semi-open system dm=m˙dt , the differential equation for the gas taking up cylinder space becomes

The observed rate of heat flow via the heat transfer surface of the cylinder Aq is determined by the overall heat transfer coefficient k(φ) as well as by the temperature difference to the steady surrounding temperature TU . Inserting into (7) gives:

Correlations to the heat transfer of cylinder spaces are found in literature. In [ ¹² ] different functions to calculate the inside heat transfer coefficient are delivered for steady-state conditions. Beyond that also transient functions for the cylinder space are available in case of internal combustion engines, but not for regenerative gas cycles. Therefore, in [ ⁴ ] the function

To calculate the rate of heat transfer the temperature at the wall surface is accepted as steady, but in [ ⁵ ] a transient consideration was also implemented. Besides in the periodical-steady state, the gas temperature in the cylinder is characterized by a sine approximation that leads to the propagation of the penetration depth of the periodical temperature variation in the wall. The divergence of the heat transfer compared to a model with a steady inside temperature by metallic materials amounts to merely a few percent. Hence this approach is suggested, if necessary, in exclusively isolated cylinder rooms.

In the following the heat transfer surface of the cylinder consists only of the side area. The area of the top and bottom results in a thermally very well isolated piston for the avoidance of heat conduction loss. However, the crossing is already in the bordering heat exchanger, so these stationary heat loss terms can be regarded separately later in the whole energy balance of the process. In this case, the quotient of heat transfer surface and cylinder volume gives:

Taking the ideal gas law and Eq. ( ¹⁰ ) into account, an equivalent form of (8) can be found:

The differential equation depends upon three groups for heat transfer, pressure, and volume. They are:

Their meaning is as follows: μ is a heat transfer parameter between gas and surrounding, similar to a gas temperature specific number of transfer units (NTU); ε is the gas-dependent ratio of the pressure gradient to the pressure; finally, δ is the ratio of the volumetric flow rate to the present volume inside the cylinder space, which is usually not zero. Rewritten using the symbols μ , ε and δ , the Riccati equation for the gas taking up cylinder (11) becomes:

So far, the model focused upon the gas taking up cylinder. For the gas delivering cylinder we received from the energy balance (6) the Bernoulli differential equation:

In comparison to this equation an additional term which arises from the temperature difference of the incoming gas flow is considered with the Riccati equation. Therefore:

• δ=fφ in case of the gas taking up cylinder space, and

• δ=0 in case of the gas delivering cylinder space.

Solution of cylinder space temperature

A simplification of the solution to the underlying problem can be produced if special state changes are regarded. This is particularly the case if it acts around an adiabatic system ( μ=0 ) or acts around a stationary system pressure ( ε=0 ). With the first-named boundary condition especially the square term is omitted, and by integrating the function of the incoming gas temperature elementary solutions of the Riccati Eq. ( ¹² ) are received. Finally, the isentropic change in state (16) is a solution of (13) in the case of the adiabatic system.

The difficulty still exists in the consideration of the usually transient pressure progression, because the functions to be integrated always contain exponents of non-natural numbers. With this consideration an analytical solution for (14) is limited, and whether it gives a suitable solution generally depends naively on the function of the pressure. Therefore, simplifications deliver only solutions for special boundary conditions, but not for the general case. In particular, the pressure progression is given during the simulation as a rule not as a function, and will be received first by the whole energy balance (1). Therefore, another procedure must be applied to get the solution.

The non-linear Riccati equation may be put in linear-form to be more suitable for calculation. If we write (cf. [ ³ ] )

The square term for the temperature disappears, and the differential equation is compared with the original form out of higher order. Writing y1=z and y2=z′ (18) gives a first order system of linear differential equation:

Finally, the solution of the ordinary differential Eq. ( ¹⁷ ) gives ( z0≠0 )

The initial value for the temperature at φ=0 can be disregarded in case of the Riccati equation, because the cylinder space contains, on this occasion, still no volume. The intensive properties of the gas are not defined at that time. Only in the interval 0<φ≤π of fluid flow does the cylinder space take up the working gas. For the initial value it is finally possible to chose any (numerically sensible) temperature

In general, and in particular for (19), there is no closed form solution for a system of linear differential equations, and an approximation method such as Magnus expansion [ ⁶ ] may have to be used. This first equation of the Magnus expansion or Magnus series coincides exactly with the solution in the scalar case

The Matrix S̲̲ in (19) is not independent of φ and so (22) does not give the whole solution to solve the matrix initial value problem. If one insists on having an exponential representation, the exponent needs to be corrected with the rest of the Magnus series.

Another possibility to receive a correction is to regard the matrix S̲̲ in segments as steady. With (22) for each of these n time steps Δφ=φ/n we get

To compute the matrix exponential we have to diagonalise the matrix S by the eigenvalues and eigenvectors. The exponential can be obtained by just exponentiation every entry on the main diagonal. With the diagonal matrix D and eigenbasis V we get [ ⁶ ]:

The averaged components of S by (19) are

The symbolic solution for the Riccati equation is very efficient; particularly in comparison to a numerical approximation. Hence, it is applied in the following to solve the Bernoulli differential Eq. ( ¹³ ) for the gas delivering cylinder space. On this occasion, the matrix S is simplified first and all terms with the volume flow rate (the additional perturbation term of the non-linear Riccati equation δ ) are omitted. The components become:

Apart from that the solution of the matrix exponential (29) coincides with the gas taking up cylinder space.

In this case the initial temperature T0(φ=0) is not discretionary anymore; it is a function of the cylinder space temperature at maximum gas capacity. This is received again from the solution of the Riccati equation with φ=π for the gas taking up cylinder.

Plot of Solution

We now evaluate the magnitudes of influence during a half power stroke with the help of the simulation. For the cylinder space the Riccati Eq. ( ¹² ) based model is available. The equation depends on:

• the heat transfer parameter μ ,

• the pressure parameter ε ,

• the volume parameter δ ,

• the temperature of the surrounding area TU , and

• the temperature of the incoming gas Tin .

For the simulation of the gas delivering and the gas taking up cylinder, the proposed solution for the Riccati equation is used; the series expansion (25) and the analytic solution for the matrix exponential (29). The boundary conditions to calculate the progression of the temperature are presented in Table  ¹ .

Figure  ² shows the solution for the temperature profiles over the course of a half power stroke (interval Δφ=π) where the gas is delivered ( δ=0 ). It is evident that the pressure change, depending on the isentropic exponent, exerts the dominant influence on the cylinder space temperature. The contribution of the heat flow is lower even if the overall heat transfer coefficient is rising. If the regenerative gas cycle is well-designed, so that heat looses are small and the pressure changes are usually high ( ε≫μT ), the adiabatic model for the gas delivering cylinder is fairly appropriate. In this case, the equation for an isentropic change in state (16) may be used.

Fig. 2

By contrast, at the gas taking up cylinder corresponding to Fig.  ³ , a significant influence of the supplied enthalpy flow is observed. Thereby the contribution of the temperature of the incoming gas depends on the ratio of volumetric flow to the volume of the cylinder space δ . In the simplest case, the volumetric flow rate is constant; for regenerative gas cycles the volumetric flow usually oscillates with a sinusoidal function. For both cases, the influence of the flow rate results in:

Fig. 3

The Eqs. (30) and (31) are a formal statement of the well-known fact that at φ=0 , there is no gas inside the cylinder ( δ=∞ ). While the cylinder volume increases, the influence of the temperature of the incoming gas is reduced by decrease of the volume coefficient. Finally, the cylinder reaches the maximum gas capacity and the movement of the piston stops ( δ=0 ). In the following time interval the parameter for the volume stays at zero; the cylinder delivers the gas to the connected component.

Conclusion

The working gas inside a machine of a regenerative gas cycle is pushed through the heat exchangers and the regenerator by cylinder piston assemblies. For the optimal design of these components one requires the information of the gas mass flow emitting and receiving by the cylinder spaces. In that regard, account must be taken of the modeling of the mass change within any cylinder space, based on differential equations for the conserved quantities mass, energy, and momentum, which are theoretically as well as numerically difficult to treat.

The proposed model of the cylinder space describes the dynamic behavior already by an ordinary non-linear differential equation, in particular due to the assumption that the alternating incoming and outgoing gas mass flow ensures an intense turbulence. Depending on the selected boundary conditions analytical solutions can be obtained for some instances. In addition, by using a special mathematical transformation, the energy equation of the cylinder space could be obtained as a system of linear differential equations, where the general solution is based on the matrix exponential function. Of course, the energy equations could also be treated by any of the numerical algorithms to solve ordinary differential equations. The important point to be emphasized here is that the main reason for the usefulness of the exponential algorithms lies in the fact that it preserves structural intrinsic properties of the exact solution.

The approach is developed for simulating the dynamic behavior of the cylinder space, depending on:

• the pulsating process pressure,

• the heat transfer to the surrounding,

• the temperature of the incoming gas flow, and

• usually oscillatory volumetric flow from the gear in question.

For the temperature within the cylinder space the theoretical investigations showed an isentropic change in state for the influence of the process pressure depending on the isentropic exponent. To retain the isothermal state a heat transfer to the surrounding must be carried out, whereby the influence of the overall heat transfer coefficient would have to be very pronounced to compensate the pressure change. In addition, the temperature change by the pressure is pulsating. However, the heat flow has only one direction: cooling down or heating up. Finally, the temperature of the gas taking up cylinder depends on the temperature of the incoming volumetric flow, which is a function of the bordering component. Some people assume that the process takes place without an increase or decrease in a state quantity of the system like the temperature. The impact of the investigated parameters indicates this assumption is not always correct; the general solution appears indispensable for realistic simulation.

The information of the parameters can be provided to the featured cylinder space simulation continuously or discreetly, as this is present as a finite number of discrete time steps. The division of the time steps depending on required simulation accuracy is in this case virtually arbitrary since there is an analytical solution for the matrix exponential function. Thus an efficient symbolic solution, in comparison to the numerical approximation, of the underlying differential equation is obtained.

The presented model of the cylinder space can be used for design simulating software of regenerative gas cycles. By the modular construction of the simulation program the external boundary conditions of the cylinder space are provided via the coupling of the other modules. Number and arrangement of the individual modules of a regenerative gas cycle is determined by the process configuration and therefore the cylinder space module is suitable for simulating different machines. Furthermore, additional physical models such as for combustion or pressure loss may be integrated. Thus, the basic approach used for the symbolic mathematic based simulation of cylinder spaces opens up many fields of application like different regenerative gas cycles, reciprocating compressors or internal-combustion engines.