Substitution-Newton-Raphson method applied to the modeling of a vapour compression refrigeration system using different representations of the thermodynamic properties of R-134a

Substitution-Newton-Raphson method applied to the modeling of a vapour compression refrigeration system using different representations of the thermodynamic properties of R-134a

J. R. Figueiredo¹; R. G. Santos; C. Favaro; A. F. S. Silva; A. Sbravati

Energy Department Mechanical Engineering School State University at Campinas

ABSTRACT

This paper gives a detailed presentation of the Substitution-Newton-Raphson method, suitable for large sparse non-linear systems. It combines the Successive Substitution method and the Newton-Raphson method in such way as to take the best advantages of both, keeping the convergence features of the Newton-Raphson with the low requirements of memory and time of the Successive Substitution schemes. The large system is solved employing few effective variables, using the greatest possible part of the model equations in substitution fashion to fix the remaining variables, but maintaining the convergence characteristics of the Newton-Raphson. The methodology is exemplified through a simple algebraic system, and applied to a simple thermodynamic, mechanical and heat transfer modeling of a single-stage vapor compression refrigeration system. Three distinct approaches for reproducing the thermodynamic properties of the refrigerant R-134a are compared: the linear interpolation from tabulated data, the use of polynomial fitted curves and the use of functions derived from the Helmholtz free energy.

Keywords: Substitution-Newton-Raphson, Newton-Raphson, minimizing variables, sparse matrices, refrigeration system

Introduction

The literature on the solution of nonlinear systems of equations refers to two basic methods, the Successive Substitution method and the Newton-Raphson (NR) method. The Successive Substitution method is generally simple to program and demands low computer memory, but may lead to divergence unless the equations are appropriately ordered, whist the NR is more reliable, at least if the trial values are sufficiently close to the solution (Stoecker, 1989). When convergent, the NR possesses quadratic order, so accelerating the convergence rate as it approaches the solution, therefore providing very accurate solutions with small to moderate number of iterations (Ortega and Rheinboldt, 1970).

One difficulty of the NR method in its original form, the computation of the partial derivatives that constitute the jacobian matrix, can be removed using numerical estimates, so combining the accuracy of the NR approach with easy utilization. The user of these codes is only required to establish the functional relationships that constitute the physical model, and to provide trial values for the variables.

The present paper describes a method called Substitution-Newton-Raphson (SNR), which combines both methods, maintaining the convergence characteristics of the NR with almost as spare use of computer memory as the Successive Substitution.

The SNR method is suitable for a set of N unknowns and N equations forming a sparse non-linear problem. The set of unknowns is subdivided in two groups, one with n so-called "effective" variables, and the other with N-n "substitution" variables. Accordingly, the set of equations is subdivided in two groups, one with n "residual" equations and another with N-n "substitution" equations. The n effective variables must be chosen in such a way that the other variables are explicitly determined by means of the substitution equations. The remaining n equations are solved according to the NR method, enforcing their residual to vanish by manipulating the n effective variables.

Mathematically, the SNR method can be interpreted as a technique to reduce the mathematical dimension of a NR problem from N to n .

The SNR method was employed by Figueiredo (1980, 1985) while analyzing a designed solar driven ammonia-water absorption refrigeration system with three distinct physical models. The simplest model assumed steady state conditions with constant heat transfer coefficients; the second, also steady state, considered the variation of the heat transfer coefficients as the fluid transfer properties varied with temperature and concentration; the third model, transient, employed the SNR method at every instant forming an implicit method for time-wise integration. Despite the varying complexity of the physical models, involving at least 56 equations, only 8 to 9 effective variables were required in all cases.

Silverio (1999) employed the SNR procedure to simulate both steady state and transient operations of an existing ammonia-water absorption cycle for scale ice production, obtaining good agreement between theoretical and experimental results. Hydraulic pressure losses through the tubing, fittings and equipment were considered. In the simplest case there were about one hundred physically relevant variables, but the SNR approach allowed the problem to be solved with ten effective variables. In the transient case the number of effective variables was a dozen.

Sbravati (1999, 2000) modeled by similar means a designed solar-driven water-lithium bromide absorption cycle for air conditioning, using eight effective variables. Three types of expansion devices were considered: a float valve for flooded evaporators, a thermostatic expansion valve and a simple capillary tube for dry evaporators. Design programs were successfully constructed for the three cases, but simulation programs converged only for the two cases with valves; in the case of the capillary tube the simulation program failed to converge, in relation to the presence of shock waves.

Santos (1998) studied the problem of substituting R-12 by R-134a for a single stage compression cycle and a double stage compression cycle with flash vapor separation and intermediate refrigeration of the compressed vapor. The single stage cycle demanded three effective variables and the double stage, five. The thermodynamic properties of both working fluids were reproduced trough polynomial fitted curves. In a further work the pressure losses through tubing and fittings were introduced, and the constant adiabatic efficiency hypothesis was replaced by a real compressor curve, requiring no extra effective variable (Santos, 1999, Santos and Figueiredo, 1999).

Favaro (1998) simulated the use of a reduced number of cylinders of the compressor in a single stage compression cycle to accommodate low demands from the refrigerated chamber. The thermodynamic properties of the R-134a were represented through tables together with a linear interpolation for the intermediate values. The electric motor equation was also introduced into the modeling in order to access the effect of rotor slipping upon the motor and compressor rotation. Only two effective variables were needed in all cases.

Silva (1999) analyzed a double-effect heat pump computing the thermodynamic properties of the R-134a according to the Baehr and Tillner-Roth (1996) approach, whose complexity raised the number of effective variables to seven.

In the published works above referred, attention is paid to the physical systems, so that none of them intended to describe the numerical method. This paper gives a detailed presentation of the SNR method as a computational tool. After reviewing the Successive Substitution and the NR methods, the composed SNR method is presented formally, then exemplified by means of an algebraic problem and, finally, applied to the modeling of a single stage R-134a vapor compression refrigeration system, analyzing its steady-state response for varying ambient temperature. The physical model is the simplest that could be proposed without losing physical realism, but is comprehensive, embodying the relevant relations from thermodynamics, heat transfer and mechanics.

Attention is paid to the calculation of the thermodynamic properties, which is a strategic part in the computation of refrigeration systems. Three distinct approaches are employed: linear interpolation from tabulated data, polynomial fitted curves and functions derived from the Helmholtz free energy.

Nomenclature

Successive Substitution and Newton-Raphson Methods

Let us consider a system of non-linear equations:

or, in vector notation, f(x)=0, where and .

The Successive Substitution method requires each of the equations f_i (x₁, ...x_j, ...x_N) = 0 to be written in a form such as , where the asterisks indicate previous values. Starting with an approximate solution x^k, the hopefully improved solution x^k+1 is obtained by an iteration rule such as:

Adequate ordering of the equations is required for convergence. The scheme (2) was put in a form analogous to the Gauss-Seidel procedure for linear equations, with the variables updated immediately after being calculated. An alternative approach, analogous to the Jordan procedure, uses only values from the previous iteration, so that x^k+1 = g(x^k).

The NR method is obtained after expressing the function f in Taylor series around x^k and neglecting second and higher order terms, leading to the approximated linear form:

where J is the jacobian matrix, calculated at point x^k:

The presumably improved solution x^k+1 is obtained after solving the system (3) as:

Following a well-established practice, the computation of the derivatives in the jacobian matrix can be automatic by means of approximate numerical methods, such as the one-sided derivative:

For each problem the code demands the specification of the matrix dimension N and of the function f(x), which is settled in a specific subroutine. According to Eq.(6), the method requires calling such subroutine N+1 times per iteration, with arguments (x₁, ...x_j, ...x_N) and for j=1,...N . The increments are chosen sufficiently small for the one-sided derivative to produce good estimates. The central differencing could give better estimates of the actual derivatives for the same increments , but it would cost greater computing time, since the subroutine of function f(x) should be called 2N+1 times. To reduce problems of ill-conditioned matrices, the solution of the linear system by Gauss factorization uses partial pivoting with row interchange based on the relative weight of the elements.

Residual control can be used to minimize the risk of divergence when the initial estimates are far from the actual solution. The total residual of an approximate solution y^k is defined as the root mean square of its elementary residuals. If this residual increases from x^k to x^k+1, the solution x^k+1 is avoided and replaced by an intermediate estimate x^k+1 = (x^k + x^k+1)/2; a new residual test is done, and so on for five times.

The Substitution-Newton-Raphson Method

Assuming the non-linear system (1) with N variables to be sparse, the first step of the SNR method is to choose a reduced set of n variables X_j, which are recalled y_i:

where j(i) is a reordering function. The n variables x_j(i) must have been chosen so that the N-n remaining variables of the extended set can be obtained explicitly, in substitution fashion, by rearranging N-n equations of system (1):

The identities (7) and the substitution equations (8) define the function x=x(y) which, for prescribed values of the reduced set of effective variables y, determines the whole set of physically relevant variables x while satisfying the substitution equations (8). For the remaining n equations these prescribed values of y , and the corresponding values of the extended set x , will provide residuals, presented below with indexes reordered for programming convenience:

These residuals are forced to vanish through the NR procedure by manipulating the effective variables y. The SNR method consists in writing the subroutine for f(y) in the form f(x(y)), i.e., by first computing the function x(y) according to equations (7) and (8), then computing the residuals f(x) by means of the equations (9).

One could substitute the function x(y), given by equations (7) and (8), into the partitioned vector function f(x) presented by (9) to obtain f(y):

However, in the simulations of complex systems, the algebraic manipulations required for such substitutions may be cumbersome and risky. The strategy proposed by the SNR method is much simpler and more convenient. This method strongly relies on the numerical computation of the partial derivatives, since the analytic differentiation of f(x(y)) would involve a cumbersome sequential use of the chain rule of differentiation.

Mathematically, the SNR method can be seen as a pure NR method of order n. When the set of variables x, of dimension N, can be expressed in terms of a subset y, of dimension n<N, through an explicit function x=x(y), then the system f(x)=0 can be put as a problem of dimension n , in the form f(x(y))=0, which is mathematically identical to f(y)=0. In other words, the substitution process employing the N-n substitution equations is interpreted as a way to produce the functional relationship giving a n-dimensional residual vector f as a function of the n-dimensional set of effective variables y. In this sense the new method is only a strategy to reduce the mathematical dimension of the Newton-Raphson problem keeping all physical information and avoiding opportunities for algebraic mistakes. Important gains in computing performance are obtained by such reduction of the problem dimension: as the number of effective variables decreases, the computer memory requirements decrease, the CPU time decreases, the choice of the initial estimates becomes simpler and, more importantly, the risk of divergence due to the usual ill conditioning of big matrices is reduced.

Stoecker (1989) noticed that the NR method does not require the equations to be listed in any special order, but it may diverge if the trial values are too far from the final solution. Because of the substitution part of the SNR method, the ordering of the equations is not immaterial as in the NR, but there is no risk of divergence due to improper choice of order of the equations.

Also because of the substitution part, previous experience in dealing with sequential algorithms for simulation or design in some sort of system is very useful for the programmer of the SNR method for this system.

The SNR is adequate for using within any NR solver, the dimension of the problem being that of the reduced set of the effective variables. In Stoecker's (1989) Generalized System Simulation Program, for instance, all equations must be put inside subroutine EQNS, which is now divided in two parts: first x(y) is computed according to the Eqs. (7) and (8), second the residuals f(x) are computed by Eqs. (9).

Algebraic Example

A simple algebraic example may be useful to clarify the explanation above. Let us consider the non-linear system:

which has two solutions, one entirely positive, (1, 2, 3, 4, 5), the other being (-1,-2,3,-4,-5).

The problem can be put as a five-unknowns system f = [f₁,f₂,f₃,f₄,f₅]^T = 0 in multiple forms, two of them quite natural. The most immediate form is:

The first three Eqs. (12) cannot be written in simpler forms. But the last two Eqs. (11) can manipulated to produce simpler forms for f₄ and f₅, avoiding divisions by unknowns, so eliminating both the risk of division by zero during the iterative process and the strongly non-monotonic behavior of the function around x₂ =0, x₄ = 0 or x₅ = 0. The resulting systems will be given by:

The system can also be put as a two-equation system, either through some algebraic substitution or by employing the SNR method. Because the first three Eqs. (11) involve small number of unknowns, they are adequate to fix, for instance, variables x₁, x₃ and x₅ by means of variables x₂ and x₄, which form the set of effective variables y:

The other variables are found by substitution with the first three Eqs. (11):

For arbitrary values of the reduced set y, the wider set x is determined through Eqs. (14) and (15), and its substitution in the last two Eqs. (13) will generate a residual f:

The NR method is then employed to enforce the residuals (16) to vanish by manipulating the effective variables y and computing the residuals after obtaining x through Eqs. (14) and (15), i.e, in the form f(x(y)).

In this simple algebraic problem one could easily substitute the Eqs. (14) and (15) into (16) obtaining the problem directly as a two equation one, in the form f(y)=0, with:

However, as already emphasized, intensive algebraic substitutions may be difficult and risky in complex real problems. The SNR method employs Eqs. (14), (15) and (16), being mathematically equivalent to the form (17), but computationally different since there is no need of algebraic substitutions and because the round-off errors will possibly differ.

Alternative two-equation systems also exists. For instance, keeping other things unchanged, the fourth Eq. (11) can be chosen, instead of the third one, as substitution equation to compute x₃. The substitution stage of the method becomes:

Accordingly, the third and fifth equations in system (11) are left as residual equations, in the more convenient form:

The above algebraic example was solved with our code, employing the four alternative functional relationships presented, i.e.: 1- NR approach using the five residuals (12); 2- NR approach using (13); 3- SNR approach using the substitution equations (15) and the two residuals (16); and 4- SNR approach using (18) and (19). Computations are done without and with the residual control. The increments for the approximate derivatives were 10^-6 and the stopping criteria was that the maximum rms. residual should be below 10^-9.

Table 1 reproduces the behavior of the four functional relationships, presenting the number of iterations required for convergence, when it happened, departing from different trial values. Initially, no residual control is employed.

In the first row of the table, the trial values are close to the first exact solution, x⁽⁰⁾ = 0.9 x; the number of iterations is the same for all the functional relationships. In the second row the trial values are much smaller than the exact solution, x⁽⁰⁾ = 0.1x; the performances of the distinct functional relationships are very different, two of them not even allowing convergence. The third row uses trial values much bigger than the exact ones, x⁽⁰⁾ = 10x; all schemes converge, but with different number of iterations.

In the remaining tests the trial values are no longer a multiple of the solution. Perturbations with alternate directions are imposed to the variables X₂ and X₄, and no perturbation is imposed to the variables X₁ , X₃ and X₅ in the five-equation approaches. In the fourth row, trial values are close to the first exact solution, with = 0.9x₂ and =x₄/0.9; again all functional relations demanded five iterations to converge. In the fifth row the trial values are = 0.1x₂ and = x₄/0.1, and in the sixth row they are = 10x₂ and = x₄/1C; the number of iterations required to achieve convergence varies in each case, and no convergence is achieved with functional relationship (2).

Clearly, for trial values close to the solution not only the order of the equations is immaterial, but their form appears to be too. However, the specific form of the functional relation does affect the number of iterations for trial values far from the solution. An overall appreciation indicates the SNR approach 3 followed by the NR approach 2 as the best functional relations to represent the system (11). In this comparison the advantages of the SNR over the NR appears to be minor, but two aspects must be recalled. First, this table is not complete for an overall comparison of the computational effort since it does not show the memory and the time required per iteration, which are bigger for the NR than for the SNR. Second, the difference between both approaches may become dramatic for more complex systems, as will be shown later.

Comparing the two-equation approaches 3 and 4 one is led to the conclusion that the simpler and the more direct substitution equation is the best. Comparing the five-equation approaches 2 and 1, one concludes that it is wise to avoid divisions by unknowns. But if the two-equation approach 3 is compared with the five-equation approach 2, it seems that the introduction of divisions by unknowns in the sequential equations produced no harm. The computations are now repeated with the residual control. As it could be expected, the residual control has no effect in the cases where convergence is fast, but it proved effective in reducing the number of iterations in some cases and in enforcing convergence in other cases, although it does not provide absolute guarantee of convergence.

Modelling of a Compression Refrigeration System

The Model

A simple, but comprehensive, steady-state modeling, involving basic thermodynamic, mechanical and heat transfer relations, will be presented for the single-stage R-134a vapor compression refrigeration cycle sketched at Fig. 1, that is supposed to operate with a floating expansion valve.

Assuming negligible kinetic and potential energy variations between the inlet and outlet fluxes, the energy balances on each component yield:

The specific work and heat transfer rates are related to the respective powers via the mass flow rate:

Negligible pressure losses are assumed through evaporator, condenser, tubing, and fittings:

The working fluid is supposed to be a saturated liquid after the condenser:

Saturation conditions are assumed for the vapor leaving the flooded evaporator:

Superheated vapor is found after the compressor:

A mixture is necessarily found after the expansion valve:

The equal pressure and the saturation conditions at points 1 and 4 imply that:

A constant adiabatic compression efficiency h is assumed:

where:

For a reciprocating compressor, the mass flow rate is determined by its displacement rate and the vapor specific volume at its inlet:

The volumetric efficiency depends on the ratio R between the clearance volume and the cylinder volume and on the vapor states at the suction and discharge:

The heat transfer rates of the evaporator and condenser are modeled according to the effectiveness approach for uniform temperature in the internal side of each heat exchanger:

More rigorous treatment of the heat transfer in the condenser could be employed to cope with the superheated zone of this exchanger, but this simplistic approach is sufficient for establishing a realistic heat transfer relationship since the greatest part of the heat is latent, not sensible, and the lower heat transfer coefficient of the superheated zone of the condenser partially compensates for the higher temperature difference. On the other side, a more rigorous modeling of the volumetric efficiency can be obtained substituting the actual discharge specific volume v₂ in Eq. (43) for the specific volume of the isentropic discharge , both forms will be compared. (Stoecker, 1958).

Linear Interpolation from Tabulated Data

One approach for reproducing the thermodynamic relations (27) to (35) is storing the relevant parts of the thermodynamic tables and creating an algorithm for choosing the adequate interval and using a linear interpolation rule between the node values. This was done for the refrigerant fluid R-134a according to the data presented by Moran and Shapiro (1995). The linear interpolation from tabulated data have artificially constant derivatives on the intervals between nodes, so the jacobian is computed with some degree of inaccuracy. Such interpolation also fails to present continuous derivatives at the nodal points, what would pose a serious problem to using the NR with analytical derivation, but not to the numerical derivatives approach.

To increase the precision of the interpolation of the specific volume against pressure in the superheated zone in Eq.(36), the specific volume is linearly related to the inverse of the pressure, considering that the hyperbolic form of that function for the case of the ideal gas relation is approximately constant in the real vapor on the range of interpolation. Some tests cases confirmed that hypothesis, as described in detail by Favaro (1998).

In the design of refrigeration systems, one must specify the temperature of the refrigerated chamber T_cam, a reasonably high ambient temperature T_cam, and adequate temperature differences through the evaporator and the condenser, so that the evaporation and condensing temperatures T₁ and T₃ are established. The chamber and the ambient temperatures also determine the evaporator heat transfer rate .

The temperatures T₁ and T₃ are the basis for the thermodynamic design, allowing the various properties and the specific work or heat transfer rate in each component to be determined according to steps 1 to 19 of Table 2. The heat transfer rate defines the mass flow rate , allowing both and to be determined, as shown at steps 20 to 22. Then the mechanical and heat transfer designs are made, in the last steps. The mass flow rate and the vapor properties define the size of the compressor, and the heat fluxes on the evaporator and the condenser, together with the respective temperature differences, determine its size parameters.

So, design departs from specified temperatures T₁ and T₃ to obtain the equipment size parameters , and n; the problem of modeling off-design performance can be seen as inverse to the problem of design, departing from given equipment size parameters to obtain the operating conditions, represented by T₁ and T₃.

The modeling of the off-design performance requires solving the 26 Eqs. (20) to (45), on the extended set of variables x = [T₁, p₁, h₁, s₁, v₁, T₂, p₂, h₂, v₂, T₃, p₃, h₃, T₄, p₄, h₄, z₄, , h_vol, h₂^*, w^*, w, q_C, q_E, ]^T. The experience of designing a refrigeration system is useful for its inverse problem as well. The design has shown that the variables T₁ and T₃ have a key role because all the other intensive properties and variables can be obtained explicitly from them. Therefore, these temperatures are suitable for constituting the effective variables, i.e., y = [T₁, T₃]]^T.

Also, the sequence of operations required in the SNR for this case has much in common with the design calculations. All the intensive properties and specific fluxes are analogously determined after the effective variables T₁ and T₃, employing Eqs. (20) to (43) to obtain the various operation characteristics of the cycle . However, in the design, the mass flow rate was obtained through the evaporator demands whilst, in the simulation, it is determined by the compressor size. The values so found are substituted into Eqs. (44) and (45) to generate the residual f(T₁,T₃), in the form below, which must vanish through the NR method.

Polynomial Fitted Curves

The second approach for reproducing the fluid properties is provided by polynomial fitted curves, a detailed description of which is given by Santos (1998). For the saturation conditions all the relevant relations have been obtained with polynomials up to degree 5: P = p_sat (T), h = h_sat.liq.(T), h = h_sat.vap.(T), s = s_sat.vap.(T), v = v_sat.vap.(T). For the superheated region the following relations were put in polynomial form of degree 2 by parts: h = h_sup(p,s), h = h_sup(p,T), s = s_sup(p,h), v = v_sup(p^-1,T). No relations were obtained for T = T_sup (p,h) and v= v_sup(p,h), as given in Eqs. (35) and (36), which were substituted by two new equations in the form:

The design still followed the route previously outlined, except that the Eq. (36) is substituted by (49), and the explicit function T = T_sup(p, h) was substituted by a Newton-Raphson procedure for a single equation doing the same job using Eq. (48).

The modeling problem also bears much in common with the previous one, except for the following. Three effective variables were employed: the previously used T₁ and T₃ plus the new effective variable. Again the Eq. (36) is substituted by (49). The Eq. (43) is no longer satisfied during the substitution steps, and must be satisfied as a residual equation. So, a three-equation residual f(T₁, T₃), must vanish through the Newton-Raphson procedure; these are equations analogous to (46) and (47) together with a new residual equation that results from the equation (48):

Equations of State Based on the Helmholtz Free Energy

The third approach for the state equations was the methodology presented by Baehr and Tillner-Roth (1996), where the various properties are calculated by derivation of the Helmholtz free energy, expressed as function of the temperature and the density.

For the saturation properties the original method requires solving a three-equation non-linear system on the pressure, the saturated liquid density and the saturated vapor density, which proved strongly ill behaved, probably because the liquid density varies slowly with the pressure, and was abandoned. The authors offer alternative approaches called "working equations", which express the pressure and the saturated liquid properties as explicit functions of the temperature alone, similarly to the previous procedure. There are also "working equations" for the superheated region, which were not needed.

For the saturated vapor the functions h = h_sat.vap(T), s = s_sat.vap(T) and v = v_sat.vap(T), as expressed at equations (32) to (34), are no longer available. Neither the functions T= T_sup(p,h), v = v_sup(p,h) and h = h_sup(p,s) for superheated vapors, as in Eqs. (35), (36) and (41), which must be reproduced iteratively from the available functions p = p(T,v^-1), h = h(T,v^-1), and s = s(T,v^-1).

The design of the refrigeration system demands two iterative procedures to find the properties v₂ and , using two equations each; the NR with analytical computation of the derivatives was chosen.

Two different codes were written for system simulation off-design. In the first code, the design strategy was incorporated to the first NR code referred, so that two NR procedures with analytic derivation were run at every iteration of another NR code with numerical derivation. Further details of this approach are given by Silva (1999). The second simulation code, to be presented here, proceeds entirely according to the SNR, performing a single NR loop. Both approaches agreed at least up to the fourth significant digit, the difference being attributable to the stopping criteria.

Since, according to Baehr and Tillner-Roth methodology, all vapor properties are computed after temperature and density, the effective variables must include these two properties for all vapor states, even for the saturation states and the fictitious exit of the isentropic compressor. As a consequence, the problem was put as a seven effective variables problem, being the reduced set .

Instead of the six no longer applicable Eqs. (32) to (36) and (41), new equations must be settled, with the following forms:

Comparing to the simplest case, this method excludes six equations and introduces eleven new ones. Correspondingly, it introduces five new variables: p₁^*, p₂^*, s₂^*, T₂^* and v₂^*. There is now 31 equations, with equal number of physically relevant variables forming the extended set x = [T₁, p₁, p₁^*, h₁, s₁, v₁, T₂, p₂, h₂, v₂, T₂^*, v₂^*, p₂^*, s₂^*, T₃, p₃, h₃, T₄, p₄, h₄, z₄, , h_vol, h₂^*, w^*, w, q_C, q_E, ]^T. The new Eqs. (51) to (58) are included as substitution equations. The set of residual equations is formed by equations analogous to (46) and (47), plus the five residual equations below, which correspond to the no longer explicitly satisfied Eqs. (28) and (39) and to the new Eqs. (59) to (61):

The physical dimensions of both the effective variables and the residuals are not homogeneous. The effective variables are temperature or specific volume and the residuals are energy rate, pressure or entropy. The normalization of these variables is useful for adopting adequate increment of the numerical differentiation, for the stopping criteria and for the residual control, if used. Estimated maximum values at design conditions were chosen for absolute temperature, specific volume, energy rate, pressure and entropy.

Results

Initially, a refrigeration system was designed by the three approaches above considered, under the hypothesis: T_env = 35^oC, T_cam = -10^oC, DT_con = 15^oC, DT_eva = 15^oC, h = 0.70, R=0.045, =50 kW. The results were:

The differences between the three methods are moderate with respect to the compressor displacement rate, and very small with respect to the heat transfer areas.

Based on these results, the following mean parameters were chosen:

These parameters were taken as the designed refrigeration system parameters for the problem of modeling the steady-state operation when the ambient temperature varies from its projected value of 35^oC, maintaining the chamber temperature T_cam = -10^oC, the compressor clearance ratio R = 0.045, and the assumed compressor efficiency = 0.70. Two aspects will be taken into account: the effects of using the distinct representations of the thermodynamic properties, and these of using the two models of the volumetric efficiency.

The resulting evaporator heat flux is presented at Fig. 2, showing the steady, almost linear decrease of the evaporator capacity with the increasing ambient temperature. The coefficient of performance, as shown at Fig.3, follows the evaporator capacity decrease. The linear interpolation produces some mildly noisy results, and do not converge for ambient temperatures above 43ºC. The representations of the thermodynamic properties do not affect significantly the result. The modified model of the volumetric efficiency, using instead of V₂ in Eq. (43), leads to lower refrigeration charge, associated with lower mass flow, as will be seen.

The coefficient of performance is displayed at Fig. 3. Again the representations of the thermodynamic properties do not affect significantly the result, except for some mild oscillations of the linear interpolation. The second model of the volumetric efficiency increases the coefficient of performance by decreasing the compressor power at a rate that more than compensates for the decreased evaporator heat flux.

The maximum system temperature, T₂, that occurs at the compressor discharge, is presented at Fig. 4, showing a steady increase with the environmental temperature. This temperature appears to be the parameter more sensible to the different approaches for thermodynamic properties, presenting differences of at most 2.4ºC. A decrease from 1,0 to 1,4ºC is found for the model of the volumetric efficiency.

The mass flow rate does not vary significantly with the environmental temperature, allowing the magnified presentation of the differences between the approaches in Fig. 5. The curious oscillatory behavior of the linear interpolation with the first model of the volumetric efficiency reflects the approximation and distancing of the linear interpolation to the actual function as the state is close to or far from the nodal points in the thermodynamic table. The mass flow rate is the parameter most affected by the different models of the volumetric efficiency, decreasing with the second model, reflecting a lower volumetric efficiency, and this tendency strengthens with the temperature.

Table 3 informs the number of iterations required for convergence for several ambient temperatures, employing the Baehr and Tillner-Roth representation of the thermodynamic properties, using both SNR and NR approaches. The stopping criterion was that both the normalized rms residual and rms maximum variation between iterations should be bellow 10^-9 for the SNR, but was relaxed to 10^-8 for the NR. The increment for numerical differentiation was 10^-6 times the reference value for each variable. In all cases the trial value corresponds to the design condition, copied to six digits.

No residual control was needed for the SNR. For ambient temperature 35ºC, which is the design condition, the SNR still took two iterations, because of the very exigent stopping criteria; for ambient temperatures far from this initial condition the number of iterations increased to six at most, despite the very exigent stopping criteria, showing the excellent convergence performance of the SNR method. The second model of the volumetric efficiency maintained the same convergence characteristics.

The difference between the SNR and the traditional NR approach, with 31 equations, becomes dramatic for this case. The NR converged only for the case defined by ambient temperature 35ºC, using the residual control. Then it took 137 iterations with the criteria 10^-8, being incapable of achieving 10^-9. It could not converge for the case 35ºC without residual control, neither for the other conditions, even with the residual control.

Conclusion

The Substitution-Newton-Raphson method was shown as a very efficient computational tool for modeling most engineering systems described by large sparse nonlinear matrices. Constructively, it is composed by Successive Substitution strategies within the Newton-Raphson method. Mathematically, it is a technique to reduce the dimension of the Newton-Raphson problem.

Some features of the method were illustrated with an algebraic example. The scheme was subsequently applied to a simple but comprehensive steady-state model of a single-stage R-134a vapor compression refrigeration system, analyzing the performance of the system to a change in the ambient temperature, also comparing three distinct approaches for reproducing the R-134a thermodynamic properties: linear interpolation from tabulated data, polynomial interpolation techniques and Baehr and Tillner-Roth´s approach based on Helmholtz function. The results obtained with the three approaches are in good qualitative and quantitative agreement.

Recalling the compression refrigeration studies with the SNR method mentioned in the Introduction, significant improvements of the physical modeling, such as including pressure drops, super- or sub-heating, experimental compressor curves, modeling of electric motor variable rotation, etc., could be done requiring small or null increase in the order of the Jacobian matrix. The most significant variation of the number of effective variables and equations refers to the distinct methods of computing the thermodynamic properties, that lead to systems with 26, 26 or 31 equations and physical variables, that could be solved respectively with 2, 3 or 7 residual equations and effective variables.

The Baehr and Tillner-Roth´s approach was also employed in a test of convergence, where the advantages of the Substitution-Newton-Raphson over the Newton-Raphson method were shown as dramatic. The conditions of convergence of the traditional method were so narrow that the difference between both schemes very often means the difference between to converge or not to converge, that is the fundamental question.

• Baehr, H.D. and Tillner-Roth, R., 1996, "Thermodynamic properties of environmentally acceptable refrigerants", Springer-Verlag, Berlin.
• Favaro, C., 1998, "Analysis of off-design operation of compression refrigeration systems" (in Portuguese). Scientific Initiation Report, FEM/UNICAMP, Campinas.
• Figueiredo, J.R., 1980, "Design and modeling of a solar driven absorption refrigeration system" (in Portuguese) Master Degree Thesis, FEM, UNICAMP, Campinas.
• Figueiredo, J.R., 1985, "Static and dynamic modeling of a solar driven absorption refrigeration system" (in Portuguese) Proceedings of the VIII COBEM, Rio de Janeiro.
• Moran, M.J. and Shapiro, H.N., 1995, "Fundamentals of engineering thermodynamics" John Wiley & Sons, New York
• Ortega, J.M. and Rheinboldt, W.C., 1970, "Iterative Solution of Nonlinear Equations in Several Variables", Academic Press, New York.
• Santos, R.G.dos, 1998, "Energy effects of the substitution of R-12 for R-134a in compression refrigeration systems" (in Portuguese) Scientific Initiation Report, FEM/UNICAMP, Campinas.
• Santos, R.G.dos, 1999, "Development of the study of energy effects of the substitution of R-12 for R-134a in compression refrigeration systems" (in Portuguese) Scientific Initiation Report, FEM/UNICAMP, Campinas.
• Santos, R.G.dos and Figueiredo, J.R., 1999, "Development of the study of energy effects of the substitution of R-12 for R-134a in compression refrigeration systems" (in Portuguese) Proceedings of XV COBEM, Lindoia.
• Sbravati, A., 1999, "Design of a water-lithium bromide absorption air-conditioning system and analysis of its operation under off-design conditions" (in Portuguese) Scientific Initiation Report, FEM/UNICAMP, Campinas.
• Silva, A.F.S., 1999, "Design of a double effect heat pump and analysis of its off-design operation" (in Portuguese) Scientific Initiation Report, FEM/UNICAMP, Campinas.
• Silverio, R.J.R., 1999, "Simulation and analysis of an ammonia-water absorption system for scale ice production" (in Portuguese) Doctoral Thesis, FEM/UNICAMP, Campinas.
• Stoecker, W.F., 1958, "Refrigeration and air conditioning", Mc-Graw-Hill, New York.
• Stoecker, W.F., 1989, "Design of thermal systems", Mc-Graw-Hill, New York.

Publication Dates

• Publication in this collection
  27 Aug 2003
• Date of issue
  July 2002