# Stability Analysis ofa Degenerate Hyperbolic System Modellinga Heat

Document Sample

```					         Stability Analysis of a Degenerate Hyperbolic
System Modelling a Heat Exchanger

Michael Hanke, K. Henrik A. Olsson and Magnus Strömgren
KTH – Numerical Analysis and Computer Science
SE-100 44 Stockholm
SWEDEN

August 11, 2005

Abstract
Mathematical modelling of a heat exchanger in a carbon dioxide heat pump, an evap-
orator, is considered. A reduced model, called the the zero Mach-number limit, is derived
from the Euler equations of compressible ﬂuid ﬂow through elimination of time scales
associated with sound waves. The well-posedness of the resulting partial differential-
algebraic equation (PDAE) is investigated by analysis of a frozen coefﬁcient linearisation
as well as by numerical experiments.
The linear stability analysis is done through transformation to a canonical form with
one hyperbolic component, one parabolic component and one algebraic component. Us-
ing this canonical form it is seen how to prescribe boundary and initial data and an
energy estimate is derived.
Numerical experiments on the nonlinear PDAE using a ﬁnite difference spatial dis-
cretisation support the linear stability analysis.

1 Introduction
A heat pump is basically an air conditioner that can switch between acting as cooler and
heater. It consists of several different components connected in a network: heat exchangers
(evaporator and condenser), compressor, valve and control units; see the schematic drawing
in Figure 1. Heat pumps have been around for a long time and traditionally the greenhouse
gas freon has been used as the refrigerant substance. Lately, concern for the environment has
led engineers to experiment with other refrigerant substances. In a recent physical exper-
iment with carbon dioxide as refrigerant in a heat pump system two different steady-states
with very different coefﬁcient of performance were observed [10]. To explain these obser-
vations one would like to derive a mathematical model of the system and make detailed
numerical simulations. However, such numerical simulations may be difﬁcult to realise for
a number of reasons:
• The refrigerant carbon dioxide has a complex state diagram.

• There are two-phase regions in the system, that is, the refrigerant is part ﬂuid and part
gas in some regions. This restricts the choice of primary variables in the mathematical
model.

1
condenser

valve                                        compressor

evaporator

Tair

Figure 1: Schematic drawing of part of a heat pump, illustrating how components are connected.

• The mathematical models of all components are based on the conservation of mass,
momentum and energy. In general, this would lead to the Euler equations of com-
pressible ﬂuid ﬂow. Special discretisation techniques are then necessary to handle the
possibility of discontinuities developing in the solution.

Especially the last property makes it hard to implement the model in system simulation
software such as, for example, gPROMS or Modelica.

In this report we consider a mathematical model of one of the heat exchangers: the evap-
orator. In Section 2 we derive from the Euler equations of compressible ﬂuid ﬂow a reduced
model called the the zero Mach-number limit. In Section 3 we investigate this degenerate
hyperbolic system by linearising the equations and freezing the coefﬁcients. We ﬁnd con-
ditions for how initial and boundary conditions should be imposed and prove an energy
estimate which shows that the linearised system is in fact weakly ill-posed. In Section 4 we
verify the relevance of the theoretical results for the complete nonlinear model by numerical
experiments. Finally, in Section 5 we present conclusions and directions for future work.

2 Derivation of a Reduced Heat Exchanger Model
Consider the mathematical modelling of one of the heat exchangers: the evaporator. The
evaporator is the device in the refrigerant cycle where energy is taken from the surrounding
air to evaporate the (ﬂuid) refrigerant. Here it is considered to be a long pipe with constant
cross section. In general it can therefore be modelled by the one-dimensional Euler equations
of compressible ﬂuid ﬂow:

∂          ∂
(Aρ) +    (Aρu) = 0,                              (1)
∂t         ∂x
∂            ∂
(Aρu) +      (A(ρu2 + p)) = R,                             (2)
∂t           ∂x
∂              ∂
(AEtot ) +     (A(Etot + p)u) = Q.                             (3)
∂t             ∂x

2
Here ρ is the density, u the velocity, p the pressure , Etot the total energy per unit volume
and A the constant cross section, while R models the momentum loss due to friction in the
pipe and Q models the energy exchange with the surrounding air. The length of the pipe is
denoted by L, and we let x = 0 at the inﬂow boundary and x = L at the outﬂow boundary.
The Euler equations are accompanied by constitutive relationships between pressure p, mass
speciﬁc enthalpy h, density ρ and temperature T :

ρ = f (p, h),                                      (4)
T   = g(p, h).                                      (5)

For one-component two-phase systems, the choice of (p, h) as thermodynamic variables is
advantageous since they uniquely determine the state. The more common choice (p, T ) is
not suitable since p is determined by T alone in the two-phase region of the phase space. As
a third variable it is convenient to use the mass ﬂow rate F = Aρu.
Assuming that the ﬂow in the heat pump system is characterised by very low Mach-numbers
everywhere except in the expansion valve, we can eliminate the time scales associated with
sound waves from equations (1)–(3). This is done by neglecting the kinetic energy contri-
1
bution to the total energy Etot = ρh − p + 2 u2 and also neglecting the inertial terms in the
momentum equation (2). The reduced model, called the zero Mach-number limit [6], reads:

∂          ∂
(Aρ) +      (Aρu) = 0,
∂t         ∂x
∂
(Ap) = R,
∂x
∂                 ∂
(A(ρh − p)) +     (Aρhu) = Q.
∂t               ∂x
Introducing F = Aρu and simplifying yields:

∂ρ ∂F
A +           = 0,
∂t    ∂x
∂p
A          = R,
∂x
∂p      ∂h     ∂h
−A    + Aρ    +F           = Q.
∂t      ∂t     ∂x
Finally, using the constitutive relation (4) we get a system of three equations in the unknowns
(F, p, h):

∂f ∂p     ∂f ∂h ∂F
A          +A       +       = 0,                              (6)
∂p ∂t     ∂h ∂t    ∂x
∂p
A      = R,                              (7)
∂x
∂p       ∂h     ∂h
−A     + Af    +F       = Q.                              (8)
∂t      ∂t     ∂x
The friction model used is
F2                  F2
R = −Lf F |u| = −Lf        sign(F ) = −Lf           sign(F ),
Aρ                Af (p, h)

3
while the energy exchange is given by

Q = Aexch α(Tair − T ) = Aexch α(Tair − g(p, h)).

Here Lf is a friction constant, Aexch is the area over which the energy exchange with the
surrounding air takes place, Tair is the temperature of the surrounding air, and α is the heat
conductivity coefﬁcient of the material [10].

3 Linear Stability Analysis
The mathematical properties of the reduced system (6)–(8) are different from those of the
original system (1)–(3). The lack of a time derivative in equation (7) makes the reduced sys-
tem a so-called partial differential-algebraic equation (PDAE); see, for example, [3], [7], [8].
We ask: Is this system well-posed? How should initial and boundary conditions be im-
posed? In this section we linearise the equations (6)–(8) around a constant state, freeze the
coefﬁcients, and undertake an analysis of the resulting system. We compute a canonical
form that provides insight into the properties of the system. Moreover, we prove an energy
estimate for the solution and in the process discover how initial and boundary conditions
should be imposed for the linearised system.

To make the following more readable we from now on denote partial derivatives by sub-
∂
scripts; for example, ∂t () = ()t .

3.1   Derivation of the Linearised Equations
We linearise equations (6)–(8) at a state (F c , pc , hc ) formally by introducing

F = F c + F ′,     p = pc + p′ ,      h = hc + h′ ,

Taylor expanding the relevant functions and neglecting quadratic terms in the, assumed
small, primed variables. The resulting system is linear and on the form

A(x, t)ut + B(x, t)ux + C(x, t)u = g(x, t).

Here u = u(x, t) = (F ′ , p′ , h′ )T . The structure of the system is important. We have
                                                                            
0 a12 a13                      1 0       0              0 c12 c13             g1
 0 0       0  ut +  0 b22 0  ux +  c21 c22 c23  u =  g2  .                              (9)
0 a32 a33                      0 0 b33                 c31 c32 c33            g3

All matrix coefﬁcients are functions of (x, t), as are the elements of the right hand side. For
completeness we state that
c
a12 = Afp ,           c
a13 = Afh ,     a32 = −A,     a33 = Af c ,     b22 = A,      b33 = F c ,
c
c12 = A(fpp pc + fhp hc ),
t
c
t
c
c13 = A(fph pc + fhh hc ),
t
c
t
c
c21 = −RF ,              c
c22 = −Rp ,
c
c23 = −Rh ,    c31 = hc ,    c32 = Afp hc − Qc ,
c              c33 = Afh hc − Qc ,
c
x                 t    p                 t    h
c
g1 = −Afp pc − Afh hc − Fx ,
t
c
t
c        g2 = −Apc + Rc ,
x                g3 = Qc + Apc − Af c hc − F c hc .
t         t        x

4
Here superscripts by c indicate that the function is evaluated at the state (F c , pc , hc ), so that
for instance f c = f (pc , hc ). Notice that the right hand side (g1 , g2 , g3 )T ≡ (0, 0, 0)T if and
only if (F c , pc , hc ) solves (6)–(8) exactly.
In order to further simplify the analysis we freeze the coefﬁcients of equation (9) at (x, t) =
(xc , tc ). The so obtained constant coefﬁcient linear system is on the form

Aut + Bux + Cu = g(x, t).                                   (10)

It is evident from the derivation that equation (10) is restricted to governing the behaviour
of a small perturbation of (F c , pc , hc ) in a neighborhood of (xc , tc ).

3.2   A Canonical Form
Computing the canonical, or characteristic, form of a system of differential equations can
lead to useful insights into the properties of the system. In the present case of equation (10)
the matrix A is not invertible. However, the matrix pencil (A, B) is regular and we can ﬁnd
matrices P, Q such that the transformation

(PAQ)(Q−1 u)t + (PBQ)(Q−1 u)x + (PCQ)(Q−1 u) = Pg,                           (11)

leads to the following canonical form:

I                J
vt +             vx + Dv = h.
N                I

Here I is the identity matrix, N is a nilpotent Jordan matrix, J is a regular Jordan matrix and
v = Q−1 u are the canonical variables. Several interesting conclusions can be drawn from
studying such a canonical form [8]. Therefore, we introduce matrices
                                       
0       a32 a−2 b−1 b33
33 22              a−1
33
P =  0 (a12 a33 − a13 a32 )a−1 b−1
33 22     0     ,
−2 −1               −1
1    −a13 a32 a33 b22 b33     −a13 a33

                                       
a13 a−1 b33
33                  0              1
Q=      0      a33 (a12 a33 − a13 a32 )−1 0  ,
1      −a32 (a12 a33 − a13 a32 )−1 0
and carry out the transformation (11). The canonical form that we ﬁnd reads
                  c          
1 0 0              u 0 0
 0 0 0  vt +  0 1 0  vx + Dv = h,                                     (12)
0 1 0               0 0 1

and we notice that the transport coefﬁcient in the ﬁrst equation is uc = F c /(Af c ). Hence,
we have one characteristic left from the Euler equations, as expected. We also note that the
block
0 0
N=
1 0

5
indicates that the system has differentiation index 2 with respect to time if the linear forcing
term Dv is neglected; see [8]. For the canonical variables Q−1 u = v = (v1 , v2 , v3 )T we have
the relations

v1 = h′ + a32 a−1 p′ ,
33                                                                 (13)
v2 =          a−1 (a12 a33 − a13 a32 )p′ ,
33                                                                 (14)
v3 =          F ′ − a13 a−1 b33 v1 .
33                                                       (15)

Equation (12) is a system of three coupled scalar differential equations:
3
c
v1,t + u v1,x +               d1j vj      = h1 (x, t),                                (16)
j=1
3
v2,x +         d2j vj      = h2 (x, t),                                (17)
j=1
3
v2,t + v3,x +               d3j vj      = h3 (x, t).                                (18)
j=1

Provided d23 = 0 the second equation (17) can be solved for v3 ,
2
1
v3 =            (h2 − v2,x −                  d2j vj ),
d23
j=1

and used to eliminate v3 from the other two equations. We ﬁnd
2
d13                                   d13                    d13
v1,t + uc v1,x −       v2,x +              (d1j −            d2j )vj   = h1 −       h2 ,
d23                                   d23                    d23
j=1
2
1          d21        (d22 + d33 )                                     d33                    d33       1
v2,t −       v2,xx −     v1,x −              v2,x +                (d3j −            d2j )vj   = h3 −       h2 −     h2,x ,
d23         d23            d23                                          d23                    d23      d23
j=1

or, more compactly written,
2
v1,t + uc v1,x − e12 v2,x +                  f1j vj     = h1 − e12 h2 ,                    (19)
j=1
2
v2,t − βv2,xx − e21 v1,x − e22 v2,x +                  f2j vj     = h3 − γh2 − βh2,x .               (20)
j=1

Here notably
1                      2F c                 c    c
f c fp + fh
β=        > 0, where d23 = Lf      sign(F c )             ,                                    (21)
d23                     Af c                 fc
which we have veriﬁed through numerical investigation of d23 = d23 (p, h) for (p, h) in the
phase space region relevant for the evaporator. Tracing our argumentation backwards we
see that equations (17) and (19)–(20) tell us something very interesting about the system (10),
namely that it consists of

6
• a hyperbolic equation for v1 ,

• a parabolic equation for v2 ,

• an algebraic relation for v3 .

This suggests that initial conditions should be imposed for v1 and v2 only. In terms of the
physical variables it is then clear that initial conditions should be imposed for h and p.

3.3     An Energy Estimate
We now show an energy estimate for the variables (v1 , v2 )T and then use (17) to ﬁnd an
estimate for v3 . In doing so, we will ﬁnd the combinations of boundary conditions that give
the desired estimate. The estimate depends on the fact (21) that β > 0.
Deﬁne a scalar product and norm by
L
(w, z) =              < w, z > dx,       < w, z >=               wj z j ,   ||w||2 = (w, w).
0                                             j

Then,

1d         v1           v1
,                   = (v1,t , v1 ) + (v2,t , v2 )
2 dt       v2           v2
= −uc (v1,x , v1 ) +e12 (v2,x , v1 ) − f11 (v1 , v1 ) − f12 (v2 , v1 ) + (h1 , v1 )
I
− e12 (h2 , v1 ) + β (v2,xx , v2 ) +e21 (v1,x , v2 ) +e22 (v2,x , v2 )
II                    III            IV
− f21 (v1 , v2 ) − f22 (v2 , v2 ) + (h3 , v2 ) − γ(h2 , v2 ) − β (h2,x , v2 ) .
V

Note that
1 2
I = (v1,x , v1 ) = [v1 ]x=L − (v1 , v1,x ) ⇒ (v1,x , v1 ) = [v1 ]x=L ,
2
x=0                                      x=0
2
II = (v2,xx , v2 ) = [v2,x v2 ]x=L − ||v2,x ||2 ,
x=0
III = (v1,x , v2 ) = [v1 v2 ]x=L − (v1 , v2,x ),
x=0
1 2
IV = (v2,x , v2 ) = [v2 ]x=L − (v2 , v2,x ) ⇒ (v2,x , v2 ) = [v2 ]x=L ,
2
x=0                                      x=0
2
V = (h2,x , v2 ) = [h2 v2 ]x=L − (h2 , v2,x ),
x=0

so that
1 d        v1            v1                   uc 2                          e22 2
,                   = [−      v1 + βv2,x v2 + e21 v1 v2 +    v − βh2 v2 ]x=L
x=0
2 dt       v2            v2                   2                              2 2
VI
2
+ e12 (v2,x , v1 ) − f11 ||v1 || − f12 (v2 , v1 ) + (h1 , v1 ) − e12 (h2 , v1 )
− β||v2,x ||2 − e21 (v1 , v2,x ) − f21 (v1 , v2 ) − f22 ||v2 ||2 + (h3 , v2 )
− γ(h2 , v2 ) + β(h2 , v2,x ).

7
Let us consider the boundary term VI. It holds
uc 2 x=L          1                                       1 (d33 − d22 ) 2 x=L
VI = −      [v1 ]x=0 − [v2     (h2 − v2,x − d21 v1 − d22 v2 )]x=L +
x=0                 [v2 ]x=0 ,
2                d23                                      2     d23
v3

and we would like to bound this expression in terms of boundary data. We proceed term
by term and ﬁnd that since v1 appears squared in the ﬁrst term, one Dirichlet condition for
this component should be imposed at inﬂow x = 0 because uc > 0. Similarly, the last term
indicates that one Dirichlet condition should be imposed for v2 at either inﬂow or outﬂow
x = L, depending on the sign of d33 − d22 . The middle term [v2 v3 ]x=L comes from the
x=0
parabolic equation (20) for v2 and presents a difﬁculty to us since we cannot, in general,
be sure of the sign of the product v2 · v3 at x = 0 or x = L. For the current purpose of
obtaining an energy estimate we therefore assume homogeneous Dirichlet conditions for v2
and v3 : a boundary condition for v2 is assumed to be imposed at either inﬂow or outﬂow
depending on the sign of d33 −d22 , and at the opposite end a boundary condition can then be
imposed for either v2 or v3 . This takes care of the second and third terms, bounding them by
zero. Not by necessity, but for ease of presentation we might now also assume homogeneous
boundary conditions for v1 so that the whole term VI is bounded by zero. The assumption of
only homogeneous data means no real restriction: the estimate obtained below also holds for
non-homogeneous Dirichlet boundary conditions, as long as they are constants. The more
realistic case of time dependent boundary conditions is however technically difﬁcult and not
considered here. Summing up, we can bound the expression VI in terms of boundary data if
constant Dirichlet conditions are imposed for
• v1 at inﬂow x = 0 and v2 at both ends x = 0 and x = L, or
• v1 and v2 at inﬂow and v3 at outﬂow (if d33 < d22 ), or
• v1 and v3 at inﬂow and v2 at outﬂow (if d33 > d22 ).
In terms of the physical variables (F, p, h)T we may interpret the relations (13)–(15) to sug-
gest that Dirichlet boundary conditions should be imposed for
• h at inﬂow x = 0 and p at both ends x = 0 and x = L, or
• h and p at inﬂow and F at outﬂow (if d33 < d22 ), or
• h and F at inﬂow and p at outﬂow (if d33 > d22 ).
Having a look at (9) where all three variables appear differentiated once with respect to
space, we conclude that for simulations the more natural choices are to impose boundary
conditions for the enthalpy h at inﬂow and mass ﬂow F and pressure p at opposing ends.
With homogenous Dirichlet boundary conditions as indicated above it holds
1d      v1         v1
,            ≤ |e12 | |(v2,x , v1 )| +|f11 | ||v1 ||2 + |f12 | |(v2 , v1 )| + |(h1 , v1 )|
2 dt    v2         v2
VII
+ |e12 | |(h2 , v1 )| − β||v2,x ||2 + |e21 | |(v1 , v2,x )| +|f21 | |(v1 , v2 )|
VIII
2
+ |f22 | ||v2 || + |(h3 , v2 )| + |γ| |(h2 , v2 )| + β |(h2 , v2,x )| .
IX

8
Moreover, since
σ2          1
||u||2 + 2 ||v||2
|(u, v)| ≤ ||u|| · ||v|| ≤
2          2σ
for all real σ = 0, we see that if we choose σ cleverly in VII, VIII and IX above (choose
σVII = β/(2|e12 |), σVIII = β/(2|e21 |) and σIX = 1) it holds:
1d        v1         v1              |e12 |2                              |f12 |                         1
,               ≤            ||v1 ||2 + |f11 | ||v1 ||2 +        (||v2 ||2 + ||v1 ||2 ) + (||h1 ||2 + ||v1 ||2 )
2 dt      v2         v2                β                                    2                            2
|(1 − |γ| − β)|                         |e12 |2               |e21 |2
+                     ||h2 ||2 +                     ||v1 ||2 +         ||v1 ||2
2                      2|(1 − |γ| − β)|               β
|f21 |                                           1
+         (||v1 ||2 + ||v2 ||2 ) + |f22 | ||v2 ||2 + (||h3 ||2 + ||v2 ||2 )
2                                              2
|γ|                             β
+      (||h2 ||2 + ||v2 ||2 ) + ||h2 ||2
2                              2
C1               C2              1             1         1
=       ||v1 ||2 +       ||v2 ||2 + ||h1 ||2 + ||h2 ||2 + ||h3 ||2 .
2                 2              2             2         2
Hence, we have shown that
d
(||v1 (·, t)||2 + ||v2 (·, t)||2 ) ≤ C3 (||v1 (·, t)||2 + ||v2 (·, t)||2 ) + ||h(·, t)||2
dt
and using Grönwall’s Lemma we get the estimate
t
||v1 (·, t)||2 + ||v2 (·, t)||2 ≤ eC3 t {||v1 (·, 0)||2 + ||v2 (·, 0)||2 +                         ||h(·, τ )||2 dτ }.   (22)
0

It remains to show an estimate for v3 . From (17) we have
2
1
v3 =     (h2 − v2,x −                d2j vj ),
d23
j=1

so that
||v3 ||2 ≤ C4 (||h2 ||2 + ||v2,x ||2 + ||v1 ||2 + ||v2 ||2 ).
Now the problem is to estimate ||v2,x ||2 . Upon differentiation of equations (19)–(20) with
respect to x and using artiﬁcial homogeneous Dirichlet boundary conditions for
• v1,x at inﬂow x = 0 and v2,x at both ends x = 0 and x = L, or
• v1,x and v2,x at inﬂow and v3,x at outﬂow (if d33 < d22 ), or
• v1,x and v3,x at inﬂow and v2,x at outﬂow (if d33 > d22 ),
we see from (22) that we can get the estimate
t
||v2,x (·, t)||2 ≤ eC3 t {||v1,x (·, 0)||2 + ||v2,x (·, 0)||2 +                    ||hx (·, τ )||2 dτ }.
0

Summing up, we have shown that
||v(·, t)||2 ≤ C(t){||v1 (·, 0)||2 + ||v2 (·, 0)||2 + ||v1,x (·, 0)||2 + ||v2,x (·, 0)||2
t
+           ||h(·, τ )||2 + ||hx (·, τ )||2 dτ }.                                                 (23)
0

9
The estimate (23) suggests that the solution is bounded by at most ﬁrst spatial derivatives
of the initial conditions and forcing functions, that is, the problem is weakly ill-posed. As a
consequence of this estimate we obtain immediately that the perturbation index of (10) is 1
with respect to time and 2 with respect to space. Having a look at (12) we see that the effect
of the linear forcing term (that is, essentially the friction) is to reduce the time index from 2
to 1, but at the same time increasing the space index from 0 to 2.
Finally, note that since the estimate (23) allows for exponential growth our considerations
do not answer the question of asymptotic long-time behaviour.

4 Numerical Stability Study
We now turn to a numerical investigation of how to simulate the nonlinear system. We want
to try different boundary conditions to see if the conclusions from the linear stability analysis
can be useful for a nonlinear simulation.
Thermodynamic properties of the carbon dioxide are computed using a piece of software
from Ruhr-Universität Bochum based on the equation of state by Span and Wagner [11].

4.1   Derivation of Numerical Model
The full nonlinear PDAE system consisting of equations (4)–(8) and the equation for the
internal energy can be written
∂ρ     1 ∂F
= −     ,
∂t     A ∂x
∂p        F |F |
0 = −A     − Lf         ,
∂x          Aρ
∂e     1 ∂F h Aexch α(Tair − T )
= −         +                 ,
∂t     A ∂x                A
0 = −ρ + f (p, h),
0 = ρh − p − e,
0 = −T + g(p, h),

with mass ﬂow rate F = Aρu, pressure p, enthalpy h, density ρ, internal energy e and tem-
perature T as dependent variables; and friction constant Lf , cross-section area A, surface
area Aexch , heat conductivity α and surrounding temperature Tair as constants. We get val-
ues for our model parameters from [10].
In order to arrive at a system suited for numerical simulation, we begin by scaling the
equations. The six dependent variables are scaled to become dimensionless and of order
one. Reasonable values for the dependent variables are a mass ﬂow rate of 0.095 kg/s, a
pressure of 4 MPa and an enthalpy of −250 kJ/kg. From the last two of these, corresponding
values for density and temperature are computed using the equation of state [11]. Finally,
an internal energy per unit volume is computed from the density, enthalpy and pressure.
The spatial variable is scaled by the length of the pipe, and the temporal variable is scaled
by the quotient of this length and the typical velocity (computed from mass ﬂow rate, pipe
area and density).
For the spatial discretisation, we use a uniform grid and derivatives are replaced by one-
sided ﬁrst order ﬁnite differences—toward the side where the boundary condition for the

10
dependent variable is given. We recall that this simple discretisation should work since the
time scales associated with the sound speed are eliminated. Moreover, for technical reasons
no shocks should develop. Note that we need to give numerical boundary conditions for F ,
p and h since they are differentiated with respect to x in the equation.
For time-stepping, we use a DAE solver by Hanke [4] with consistent initialisation as
described by Hanke and Lamour [5].
The functions for computing density and temperature from pressure and enthalpy are
continuous and monotone, but they are not differentiable at the saturated vapour line in
the pressure–enthalpy diagram. The jumps in the partial derivatives at this line are a big
problem for the DAE solver as the line comes close to a discretisation point during time-
stepping. Therefore, for each function we compute its values on a grid and approximate
it with an interpolating C 1 monotone quadratic spline surface, using a code by Beatson [1]
for the approximation by Beatson and Ziegler [2]. Note that an approximation that is not
monotone would be unphysical and might give spurious results in the time-stepping.
With the boundary conditions that the scaled mass ﬂow rate, pressure and enthalpy are
all equal to one at the inﬂow, we compute a steady state solution using the M ATLAB function
lsqnonlin (with the time derivatives set to zero). This is a function for solving nonlinear least
squares problems using a subspace trust region method based on a Newton method; it is
available in the M ATLAB Optimization Toolbox. Thus, we get an initial condition and also
values to use as boundary conditions at the outﬂow boundary when needed.
With the discretised variables F , p, h, ρ, e and T collected in a vector y, we can write our
system as
˙
M y = φ(y),
where M is a singular constant matrix and φ is a nonlinear function. We linearise at a steady-
state solution y0 by assuming a solution of the form y = y0 + z, which gives

˙     ˙
M z = M y = φ(y) ≈ φ(y0 ) + Jφ (y0 )z = Jφ (y0 )z

neglecting higher order terms in z. Starting from y0 we have for the deviation z:

˙
M z = Jφ (y0 )z ≡ Az,

where the constant matrix A is the Jacobian Jφ evaluated at the steady-state solution. We
can now study the stability of the discretised problem to small perturbations by computing
the eigenvalues of the matrix pencil A − λM , using the functions numjac, for evaluating the
Jacobian, and eig, for computing eigenvalues. The justiﬁcation of the eigenvalue analysis for
assessing the asymptotic stability of DAEs is given in [9].

4.2   Numerical Results
All computations are made in M ATLAB 6.1, with the usual double precision where the ﬂoat-
ing point relative accuracy is 2.22 × 10−16 . This was run under Linux on a PC with 1 GB
memory and a Pentium 4 3.00 GHz processor with 1024 kB cache.
We have used a spatial discretisation with a step 1/50, thus arriving at a system of size
294.
In the interpolation of the density and temperature functions, we have used a grid with
50 steps in the p-direction on the interval p ∈ [0.8, 1.2] and 100 steps in the h-direction on the

11
interval h ∈ [0.05, 1.1]. Note that these are scaled variables; we have (p, h) = (1.00, 1.00) at
the inﬂow and (p, h) = (1.00, 0.17) at the outﬂow in the stationary solutions.
For each of the variables F , p and h, we need to give a boundary condition at either the
inﬂow or the outﬂow. We try all possible combinations, and ﬁnd that for two of these eight
the pencil A − λM of the linearised system has all its ﬁnite eigenvalues in the left half of
the complex plane—meaning that the system is stable. It is clear from the structure of the
system that 4/6 of the eigenvalues are inﬁnite. We ﬁnd that 1/6 are located at the real axis
and the remaining 1/6 are spread in the complex plane; see Figure 2 for two examples of such
complex groups. Each of these 1/6 parts are either stable or unstable, as seen in Table 1. The

F      p      h   real stable           complex stable   real unstable      complex unstable   inﬁnite
in    in     in        1                      0                0                   1              4
in    in    out        0                      0                1                   1              4
in   out     in        1                      1                0                   0              4
in   out    out        1                      0                0                   1              4
out    in     in        1                      1                0                   0              4
out    in    out        1                      0                0                   1              4
out   out     in        1                      0                0                   1              4
out   out    out        0                      0                1                   1              4

Table 1: Number of groups of stable and unstable eigenvalues, located at the real axis or in a region
of the complex plane, for all the combinations of where to specify boundary conditions.

two combinations that give stable spectra are: F and h given at the inﬂow, p at the outﬂow;
and p and h given at the inﬂow, F at the outﬂow. For each of these two, we show the least
damped eigenvalues in Figure 2. The other ﬁnite eigenvalues lie on the negative real axis,
with magnitudes in the range from 104 to 108 .

200

150

100

50
Im{λ}

0

−50

−100

−150

−200
−400 −350 −300 −250 −200 −150 −100 −50    0
Re{λ}

Figure 2: Least damped eigenvalues for stable discretisations, with F , h at inﬂow, p at outﬂow (∗),
and p, h at inﬂow, F at outﬂow (o).

We also test the stability numerically by trying to do a time-stepping for each set of
boundary conditions. Both the relative and the absolute tolerance in the DAE solver are

12
set to 10−8 . Starting from a stationary solution, we run a simulation from t = 0 to t = 0.8
with a linear increase in the pressure boundary condition from 1 to 1.01. This works for the
two combinations that were seen to give stable spectra above; see statistics in Table 2 which
indicate that the one with F and h given at the inﬂow, and p at the outﬂow, is preferable.
For the other six combinations, the simulations crash. Longer time simulations, for 20 time

F      p     h   successful steps   failed attempts    function evaluations      solution time
in   out    in         164                 16                  1286                  2.47 s
out    in    in         357                 43                  7172                  9.88 s

Table 2: Statistics for simulations of the stable discretisations, with boundary conditions given at
inﬂow and outﬂow as indicated.

units with an increase in the pressure boundary condition by 10 per cent during the ﬁrst 8
time units followed by a constant pressure, lead to a new steady state and indicate numerical
stability. In all these simulations, the time-stepper identiﬁes the system as an index-1 DAE.

5 Conclusions
We have derived a mathematical model of an evaporator, called the the zero Mach-number
limit, from the Euler equations of one dimensional compressible ﬂuid ﬂow by eliminating
time scales associated with sound waves. The resulting degenerate hyperbolic system (6)–(8)
is a partial differential-algebraic equation (PDAE), and it is not obvious whether this system
is well-posed or not and how initial and boundary data should be imposed.
We analysed a frozen coefﬁcient linearisation (10) of the PDAE (6)–(8), transforming it to
a canonical form (12) which immediately suggested that one characteristic was left from the
Euler equations, as expected. Furthermore, we saw that under a reasonable assumption (21)
the canonical system (12) consisted of a hyperbolic equation, a parabolic equation and an
algebraic relation. This suggested that initial conditions should only be imposed for two (v1
and v2 ) of the three canonical variables (v1 , v2 , v3 )T . As for boundary conditions we found
three admissible combinations of Dirichlet conditions. Of these, the more natural choices
were to impose conditions for v1 at inﬂow and v2 , v3 at opposing ends. In terms of the
physical variables this could be interpreted to suggest that boundary conditions should be
imposed for the enthalpy h at inﬂow and mass ﬂow F and pressure p at opposing ends. With
these boundary conditions we were able to prove an energy estimate (23) for the solution.
The estimate suggested that the linear system was weakly ill-posed: the perturbation index
was found to be 1 with respect to time and 2 with respect to space.
In a numerical study of the nonlinear PDAE we found that numerical boundary condi-
tions should be imposed according to the results from the linear analysis, other combinations
lead to crashes in time-stepping. The time-stepper identiﬁed the system as index 1 with re-
spect to time.
We conclude that the present investigations has shown a good qualitative agreement
between the linearised frozen coefﬁcient system (10) and the full nonlinear system (6)–(8).
For future work it is important to know if the energy estimate (23) holds for the nonlinear
system, too. In the perspective of simulating the complete heat pump system, estimates
allowing for time-dependent boundary conditions are necessary.

13
References
[1] R. Beatson. Bivariate monotone interpolation by Beatson and Ziegler Mar 85. From
Netlib: http://www.netlib.org/a/beatson.

[2] R. Beatson and Z. Ziegler. Monotonicity preserving surface interpolation. SIAM Journal
on Numerical Analysis, 22(2):401–411, 1985.

[3] S.L. Campbell and W. Marszalek. The index of an inﬁnite dimensional implicit system.
Mathematical and Computer Modelling of Dynamical Systems, 5(1):18–42, 1999.

[4] M. Hanke. A new implementation of a BDF method within the method of lines. Journal
of Computational Methods in Sciences and Engineering, In press.

[5] M. Hanke and R. Lamour. Consistent initialization for nonlinear index-2 differential-
algebraic equation: large sparse systems in MATLAB. Numerical Algorithms, 32:67–85,
2003.

[6] S. Klainerman and A. Majda. Compressible and incompressible ﬂuids. Commun. Pure
Appl. Math., 35:629–651, 1982.

[7] W. Lucht, K. Strehmel, and C. Eichler-Liebenow. Indexes and special discretization
methods for linear partial differential algebraic equations. BIT, 39(3):448–512, 1999.

[8] W. S. Martinson and P. I. Barton. Index and characteristic analysis of linear PDAE sys-
tems. SIAM Journal on Scientiﬁc Computing, 24(3):905–923, 2002.

[9] R. März. Practical Lyapunov stability criteria for differential algebraic equations. In Numer-
ical Analysis and Mathematical Modelling, J. K. Kowalski and A. Wakulicz, editors,
volume 29 of Banach Center Publications, pages 245–266. 1994.

[10] P. G. Mehta, B. Eisenhower, and J. Oppelstrup. Bifurcation analysis of the CO2 steady
state model. Internal report, United Technologies Research Center, Hartford, CT, 2002.

[11] R. Span and W. Wagner. A new equation of state for carbon dioxide covering the ﬂuid
region from the triple-point temperature to 1100 K at pressures up to 800 MPa. Journal
of Physical and Chemical Reference Data, 25:1509–1596, 1996.

14

```
DOCUMENT INFO
Shared By:
Categories:
Stats:
 views: 6 posted: 6/15/2010 language: English pages: 14
How are you planning on using Docstoc?