Collaboration:

• The Technical University of Denmark
• burgeon

ABSTRACT

The ideal Stirling working cycle has the maximum obtainable efficiency defined by Carnot efficiency, and highly efficient Stirling engines can therefore be built, if designed properly. To analyse the power output and the efficiency of a Stirling engine, numerical simulation programs (NSP) have been developed, which solve the thermodynamic equations.

In order to find optimum values of design variables, numerical optimization techniques can be used (Bartczak and Carlsen, 1991). To describe the engine realistically, it is necessary to consider several tens of design variables. As even a single call for NSP requires considerable computing time, it would be too time consuming to use conventional optimization techniques, which require a very large number of calls for NSP. Furthermore, objective and constraint functions of the optimization problem present some level of noise, i.e. can only be estimated with a finite accuracy. To cope with these problems, the multipoint explicit approximation technique is used.

STIRLING ENGINE WORKING PRINCIPLE

Stirling engines seem to have a future in small combined heat and power plants (co-generation plants). Two applications are currently under consideration at the Technical University of Denmark. The first one is an engine driven from biomass for co-generation plants in the range 150-1500 kW el. The other application is being developed for individual co-generation, when a small Stirling engine driven from natural gas could supply a single family house with heat and power with a high level of energy utilization.

The Stirling engine is based on a closed cycle, where the working gas is kept inside the cylinders and heat is added and removed from the working spaces through the cylinder wall. The engine has a working piston, which converts the gas pressure into mechanical power, and a displacer piston, which can move the working gas between the hot and the cold working spaces (Reader and Hooper, 1983).

The basic Stirling engine working cycle is illustrated in Figure 1. The cycle is divided into four parts:

From the p-V and T-s diagrams it can be seen that power is produced on the working piston.

If a regenerator is placed in the channel between the cold and the hot cylinder volumes, heat can be stored, when the gas flows from the hot to the cold cylinder volume and used for reheating the gas, when it flows back to the hot volume. The ideal Stirling working cycle will then have the maximum obtainable efficiency defined by the Carnot efficiency.

Exchange of heat through the cylinder walls is inadequate in a real Stirling engine, therefore, heat exchangers must be added. The Stirling engine differs from conventional internal combustion engines by having external combustion like in a boiler. All sorts of gas, liquid and solid fuels can be used depending on the external burner system only.

Stirling engines can be build to obtain high efficiencies, if designed properly. Heat exchangers, regenerators and cylinder volumes have to be optimized carefully for heat transfer, flow losses, dead volume, etc. Therefore, a numerical simulation program (NSP) is an important tool for the design of engines with maximum efficiency and power density.

STIRLING ENGINE SIMULATION PROGRAM

The computer program is developed for the design of Stirling engines. The thermodynamic program calculates power output, efficiency and heat flow from the equations. The Stirling machine is divided into five different volumes: two cylinder volumes, two heat exchanger volumes and one regenerator volume.

The heat exchanger volumes are subdivided into smaller control volumes with constant temperature in each volume as in the cylinder volumes, while the regenerator is treated as a single volume with a constant temperature gradient.

The first law of thermodynamics and the law of continuity are used together with conservation of mass and energy equations to establish a set of differential equations, which describe the change of temperature and pressure in the Stirling engine. The differential equations are solved by means of the fourth order Runge-Kutta method. The initial pressure and temperatures for steady state operation are found by iterative procedure.

The basic assumptions made in the analysis of the Stirling cycle are:

• perfect gas law applies in each control volume;
• the pressure is the same throughout the machine;
• the flow is one dimensional throughout the machine;
• the working gas is perfectly mixed in the volumes;
• the kinetic energy of the working gas is neglected;
• the heat conduction parallel to the flow is neglected in the gas and the regenerator matrix;
• power losses due to frictional pressure drop trough the heat exchangers and regenerator are added after the simulation is finished;
• parasitic losses caused by heat conduction, shuttle conduction, etc. are added to the results obtained by the simulation.

The temperature in the control volumes and the cycle pressure can then be calculated from the differential equations of conservation of energy and continuity. These equations are valid for all N control volumes in the Stirling machine. If the heat exchanger volumes are divided into eight subvolumes each, the regenerator into two subvolumes, we have a total of N = 22 control volumes resulting in 22 independent variables, as mass of gas in the machine is known.

The differential equations are solved as a initial value problem, where steady state initial values are found by iterations.

Calculation of the heat transfer coefficient is based on experimental results from Kays and London (1984). The pressure drop caused by the gas flow through the heat exchangers is calculated from the mass flow, and the power losses are found by integration. The friction factor for the calculation of the pressure drop is also taken from Kays and London (1984).

To verify NSP, the numerical results have been compared to the results of laboratory experiments with a 8 kW natural gas driven Stirling engine. The comparison showed that NSP allows to predict the Stirling engine performance accurately.

FORMULATION OF THE OPTIMIZATION PROBLEM

Consider a Stirling engine optimization problem:

Maximize F_o (x), x Î R_N subject to

F_j (x) £ C_j , (j = 1,..., M) and A_i £ x_i £ B_i , ( i = 1,..., N), where  

x = ( x₁ ,..., x_N ) is a vector of design variables which can be varied in the course of design procedure (bore, stroke, heat exchanger tube dimensions, etc.);

F_o (x) is an objective function, which provides a basis for choice between alternative acceptable designs (efficiency);

F_j (x) , (j = 1,..., M) are the constraint functions, which impose limitations on various characteristics of the Stirling engine (power output, specific power);

A_i and B_i are side constraints, describing physical upper and lower limitations on design variables.

Depending on the specific problem under consideration, functions F_j (x), (j = 0,..., M) can describe various Stirling engine response quantities, such as thermodynamic efficiency, power output, etc.

MULTIPOINT APPROXIMATION TECHNIQUE 

Introduction

The formulated optimization problem has the following characteristic features:

• an objective function and constraints are implicit functions of design parameters. To evaluate these functions means to use the numerical simulation program which involves a large amount of computer resources;
• function values present a certain level of noise, i.e. can only be estimated with a finite accuracy due to the iterative nature of NSP.

The prohibitive computational cost of the direct combination of numerical simulation programs with conventional methods of mathematical programming stimulated the idea of sequential approximation of the initial optimization problem (8) - (10) by explicit subproblems. Nowadays this concept has proven to be very efficient and it is widely used in structural optimization (see a recent review by Barthelemy and Haftka, 1993). It is based on the information obtained by numerical response analysis and first order design sensitivity analysis (i.e. values of functions and their derivatives) at a current point of the design variable space (single point approximations). At present, several first order approximation techniques are known, which are based upon the function value and its derivatives at the current and the last previous design points (two-point approximations) (Barthelemy and Haftka, 1993). The drawback of these approaches is that the convergence of optimization algorithms can be very slow if the objective/constraint function values present a considerable level of noise.

A different approach to complex engineering systems optimization (Schoofs, 1987) is to create approximate explicit expressions by analysing a chosen set of design points (design of experiments). This approach is based on the multiple regression analysis methods which can use information being more or less inaccurate. They are global in nature and allow to construct the explicit approximations valid in the entire design space, but restricted by relatively small optimization problems (up to ten design variables, Vanderplaats, 1989).

The remainder of the section presents a general iterative technique, which has been originated by complex engineering systems optimization problems. The simplified approximations of the original objective/constraint functions are obtained by the weighted least-squares method. It uses in each iteration the information about function values gained at several previous design points (multipoint approximations, Toropov, 1992) in order to improve the quality of simplified approximations and to reduce the total number of calls for complex engineering systems analyses. It combines the advantages of two above mentioned basic approaches.

Basic Approximation Concept

The approximation concept leads to the iterative approximation of the original functions F_j (x), (j = 1,..., M) by the simplified functions F_j^k (x). The initial optimization problem is then replaced with the succession of simpler subproblems as follows:

Find the vector x_*^k that maximizes the objective function F_o^k (x)

subject to F_j^k (x) £ C_j , (j = 1,..., M) and A_i^k £ x_i £ B_i^k , A_i^k ³ Ai , B_i^k £ Bi , ( i = 1,..., N) ,  

where k is the current iteration number. Note that the current move limits A_i^k and B_i^k define a subregion of the design space where the simplified functions F_j^k(x) can be considered as adequate approximations of the initial functions F_j (x). To estimate the order of their adequacy in comparison with the initial functions, the error parameters  r_j^k = | [ F_j (x_*^k) - F_j^k (x_*^k) ] / F_j (x_*^k) | , (j = 0, 1,..., M) can be evaluated. The next, (k+1)-th iteration is started from the obtained point x_*^k.

Multipoint Approximations

To construct the simplified expression F_j^k(x) in the above optimization subproblem, we shall implement the methods of regression analysis (Draper and Smith, 1981). These methods are intended for obtaining an expression that reflects the behaviour of an object considered as a function of its parameters, based on a set of experimental results. Here and in the remainder of this section, an experiment means a numerical experiment using NSP. Note that we do not intend to construct simplified expressions that are adequate in the whole of the search region determined by side constraints A_i and B_i because it takes too large number of numerical experiments in the case of a multiparameter problem. Therefore, we construct such expressions iteratively only for separate search subregions which are determined by move limits A_i^k and B_i^k at each k-th step of the optimization subproblem. Thus, the function F_j^k (x), (j = 0,..., M) give piece-wise approximation of the initial functions F_j(x). To simplify notation, we will suppress the indices k and j on the functions F_j^k (x). Assume that the functions of a current optimization subproblem are expressed in the following general form:

F = F (x, a).

The vector a = (a₀ , a₁ ,..., a_L)^T in the above expression consists of so-called tuning parameters, that is, free parameters the value of which is determined on the basis of numerical experiments at points located in the design variable space R^N in accordance with some design (plan) of experiments. Then the weighted least-squares method leads to the following problem:

Find the vector a that minimizes the function

G ( a ) = å w_p [ F ( x_p ) - F (x_p, a) ] ² (summation over plan points p = 1, ..., P ),

where p is the number of a current point in the plan of experiments, P is the total number of such points, w_p is the weight coefficient that characterizes the relative contribution of the p-th experiment information. The solution of this optimization problem is the vector a which makes up the simplified function F (x, a).

The quality of approximations depends strongly on values of weight coefficients w_p , they reflect the inequality of data obtained in different design points. Usually, the point that corresponds to the solution of a complex engineering systems optimization problem lies on a boundary of the feasible region. Therefore, the quality of approximations should not be the same within the current search subregion. On the contrary, it should be the best in a domain, which is located along the boundary of the feasible region. This is possible if the weights depend on the distance of a current point from the boundary of the feasible region, i.e. on the values of F (x) . Similarly, weights can reflect the difference in the contribution of the information given by different experiments depending on the objective function value at individual design points. Then the maximum value of the weight coefficient corresponds to the design point with the maximum value of the objective function.

Choice of the Simplified Expression Structure

To construct the simplified expressions F, it is necessary to define them as a function of tuning parameters a . Apparently, the efficiency of the optimization technique depends greatly on the accuracy of such expressions. The simplest case is a linear function of parameters a :

 F(a) = a₀ +   a₁ j₁ + ...+ a_L j_L .

If the results F (x_p) , (p = 1,...,P), of P numerical experiments are known, then the minimization problem of finding the values of tuning parameters a is equivalent to the solving of the linear system of L+1 normal equations with L+1 unknowns a:

[F ]^T [W] [F ] a = [F ]^T [W] f ,  

where [F ] is the rectangular matrix consisting of values of individual functions j_l (regressors) from the expression for F(a) above at every plan point; f is the vector containing information from experiments obtained at P plan points, i.e. values of implicit functions F ( x_p ); [W] is the diagonal matrix consisting of weight coefficients w_p.

Note that the structure of the linear simplified expression F(a) is rather general because the individual regressors j_l can be arbitrary functions of design variables.

The procedure described above can be further generalized by the application of intrinsically linear functions (Draper and Smith, 1981). Such functions are nonlinear, but they can be led to linear ones by simple transformations. For example, the multiplicative function

 F(a) = a₀ j₁^a1 ... j_L^aL

requires the logarithmic transformation

 F(a) = ln a₀ +   a₁ ln j₁ + ...+ a_L ln j_L .

Several particular forms of intrinsically linear functions were implemented for the solution of various complex engineering systems optimization problems (Toropov, 1989). The obtained results indicate that the multiplicative function is the most universal one. It should be noted, however, that the structure of the simplified function need not necessarily be linear or intrinsically linear. If there are arguments for the application of nonlinear expressions in the general form, then the problem of the least-square estimation would become a standard problem of nonlinear programming, which can be solved by any available technique. Moreover, the simplified functions F need not necessarily be explicit. There can be numerical procedures involved in their formulation, such as numerical integration. But, the basic requirements to such simplified models are:

• its description of the simplified process must depend on the same parameters x as the initial numerical model, presented by functions F(x);
• they have to contain some tuning parameters a, which are obtained by solving the least-squares fitting optimization problem;
• they have to be simple enough to be used in numerous repeating calculations;
• in order to achieve fast convergence of the algorithm, they have to be accurate enough in comparison with the original functions F(x);
• they have to be noiseless or, at least, the level of noise must not cause problems with convergence of an algorithm used to solve each optimization sub-problem.

Choice of a Search Subregion

After formulation of the simplified functions F (x, a), a current mathematical programming sun-problem is solved and the error parameters r_j^k in for the obtained point x_*^k are estimated. Then the task is to determine the move limits A_i^k and B_i^k for the next iteration. First, the condition

 r_j^k £ e_j^k , (j = 0,...,M)

 is checked, where e_j^kare small positive values which define the feasible accuracy of approximation of functions F_j (x) by the functions F_j^k (x) . If this condition is not satisfied even for one of the active constraints or objective function, (i.e. approximations are inaccurate), then the size of the search subregion of the (k+1)-th step must be reduced. When the accuracy conditions are satisfied for all active constraints and the objective function, we must decide upon the movement of the search subregion. If the point obtained, x_*^k, is located inside the k-th search subregion (none of the move limits are active), then this point can be considered as the current approximation of the solution x_* . In that case, the next search subregion should be reduced and the other conditions of the search termination should be checked. Otherwise, the search must be continued. This means that the search subregion must be moved in the direction x_*^k- x_*^{k - 1}. The search process is terminated when (i) the accuracy conditions are satisfied, (ii) none of move limits are active and (iii) the subregion has reached a required small size.

OPTIMIZATION RESULTS

The optimization technique has been used together with the Stirling engine simulation program in order to improve efficiency and power output. The efficiency has been maximized, while the power output has been controlled by two constraints. The first constraint forces the power output not to be less than 10 kW, and the second constraint keeps the specific power above a certain minimum.

The objective and constraint functions have been approximated by multiplicative expressions in the simplest form:

 F(a) = a₀ x₁^a1 ... x₁₈^a18.

Table 1 shows the results of the optimization of the Stirling engine, which has been used for tests in the laboratory. The optimization process allowed to improve the efficiency from 30.1% to 40.2% without changing the important design parameters like bore and stroke by more than 20%.

The number of calls for the numerical simulation program in all runs has not been greater than 200.

CONCLUSION

A numerical simulation program (NSP) for the calculation of power output and efficiency of Stirling engines has been combined with a numerical optimization program utilizing iterative multipoint explicit approximation technique. Previously the accuracy of the NSP had been verified by comparison of calculated results with results from the measurements on a 8 kW Stirling engine.

This numerical tool has been applied to the optimization of the Stirling engine used for laboratory tests. The optimization was carried out using 18 design parameters, the results showed that the efficiency could be increased by more than 10% with the same power output.

The combination of NSP and an advanced optimization technique has shown to be a very efficient tool for improving Stirling engine performance. Solutions were found, which had never occurred by the traditional "trial and error" approach. Also the optimization helped the designer to understand the fundamentals of Stirling engine performance in a completely new way, and the possibilities for improvements were found to be numerous.

Furthermore, the multipoint approximation technique, which has been used for the optimization, has showed its superior qualities, as the number of NSP calculations was reduced by the factor of ten compared to conventional optimization techniques. The capability of handling a large number of design variables without spending enormous amounts of computer time can result in various applications for the improvement of engineering designs.

REFERENCES

Bartczak, L. and Carlsen, H., 1991. An optimization study of Stirling engines based on advanced simulation. Proceedings, 5th International Stirling Engine Conference, Dubrovnik, pp.161-166

Barthelemy, J.-F.M. and Haftka, R.T., 1993. Approximation concepts for optimum structural design - A review. Structural Optimization, Vol. 5, pp.129-144

Draper, N.R. and Smith, H., 1981. Applied regression analysis, 2nd ed., Wiley, New York

Kays, W.M. and London, A.L., 1984. Compact Heat Exchangers, 3rd ed., McGraw-Hill, New York

Reader, G. and Hooper, C., 1983. Stirling Engines, E. & F.N. Spon, London

Schoofs, A.J.G., 1987. Experimental design and structural optimization. Ph.D. Thesis, Eindhoven University of Technology, Eindhoven, The Netherlands

Toropov, V.V., 1989. Simulation approach to structural optimization, Structural Optimization, Vol. 1, pp. 37-46

Toropov, V.V., 1992. Multipoint approximatiom method in optimization problems with expensive function values. Computational Systems Analysis 1992, A. Sydow, ed., Elsevier, Amsterdam, pp. 207-212

Vanderplaats, G.N., 1989. Effective use of numerical optimization in structural design, Finite Elements in Analysis and Design, Vol 6, pp. 97-112