Gary Phetteplace August 1995

Optimal Design of Piping Systems

for District Heating

CRREL REPORT 95-17

80

RETURN TEMPERATURE (°C)

0

LOGARITHMIC MEAN

GEOMETRIC MEAN

ARITHMETIC MEAN

0.1 1.0

HEAT LOAD RATIO
Cover: Humpty Dympty falling off of the very high wall.

For conversion of SI units to non-SI units of measurement consult ASTM Standard

E380-93a, Standard Practice for Use of the International System of Units,

published by the American Society for Testing and Materials, 1916 Race St.,

Philadelphia, Pa. 19103.

This report is printed on paper that contains a minimum of 50% recycled

material.

Cover: Results from three return

temperature models.

Abstract

First, a method for determining the optimal size for a single pipe segment in a

district heating system is developed. The method is general enough to allow for

any set of economic or physical parameter values. In addition, any form of load

management, i.e., temperature or flow modulation, or both, can be accommodated by the integral form of the coefficients in the cost equation. An example

is presented that shows a 17% savings in life cycle costs over a design based

on a common rule of thumb. Next the heat consumer and his effects on the

piping system are studied. A new model is developed for the consumer’s heat

exchanger that uses the geometric mean temperature difference as an approximation for the logarithmic mean temperature difference. The new consumer

model is integrated into the previous single pipe model and, for a sample case,

its effect is determined. For systems having multiple pipes and consumers, the

constraints are first developed and then the general solution strategy. The

method makes use of the solution to the unconstrained problem as a starting

point for the constrained solution. Monotonicity analysis is then used to prove

activity of some of the constraints, and thus simplify the problem. Finally, the

branch-and-bound technique is shown to be suitable for finding a design with

discrete values for all the pipe diameters. A simple example is provided. In

addition, a method is also demonstrated for further refinement of the pipe

network to eliminate excessive throttling losses in the consumer’s control

valves. The method developed here should be feasible for designing the piping

networks for district heating systems of moderate size, and its major advantage

is its flexibility.
Gary Phetteplace August 1995

CRREL Report 95-17

U.S. Army Corps

of Engineers

Cold Regions Research &

Engineering Laboratory

Optimal Design of Piping Systems

for District Heating

Approved for public release; distribution is unlimited.
PREFACE

This report was prepared by Dr. Gary E. Phetteplace, Mechanical Engineer,

Applied Research Division, Research and Engineering Directorate, U.S. Army Cold

Regions Research and Engineering Laboratory. The work reported here was funded

in part by DA Project 4A762784AT42, Work Unit CO-M06, Low Temperature Heat

Distribution Systems.

Technical review of this report was provided by Dr. V.J. Lunardini (CRREL) and

Professor D.S. Wilde (Stanford University, Department of Mechanical Engineering).

The contents of this report are not to be used for advertising or promotional

purposes. Citation of brand names does not constitute an official endorsement or

approval of the use of such commercial products.

ii
iii

CONTENTS

Page

Preface ..................................................................................................................... ii

Nomenclature ......................................................................................................... v

Chapter 1: Introduction ........................................................................................ 1

Current design practice .................................................................................... 1

Optimization in district heating system design ............................................ 2

Chapter 2: Optimal pipe diameter for a single pipe segment ......................... 5

Cost of heat loss ................................................................................................. 5

Cost of pumping ................................................................................................ 6

Cost of pipes and pumps .................................................................................. 9

Cost of maintenance and repair....................................................................... 9

Total cost ............................................................................................................. 9

A simple example .............................................................................................. 12

Chapter 3: The consumer‘s heat load ................................................................. 17

Simple model for the consumer‘s space heating equipment ...................... 17

Design of a single pipe segment with a consumer model ........................... 27

Heat consumer constraints ............................................................................... 31

Chapter 4: Constraints on the design of systems with multiple consumers

and multiple pipes ......................................................................................... 33

System constraints ............................................................................................. 33

Differential pressure constraints ..................................................................... 33

Maximum absolute pressure constraints ....................................................... 35

Minimum absolute pressure constraints ........................................................ 38

Chapter 5: Optimal design of systems with multiple pipes and consumers 41

Some possible approaches ................................................................................ 41

Pipe and pipe junction labeling scheme ......................................................... 43

System objective function ................................................................................. 43

Solution strategy ................................................................................................ 45

Chapter 6: A simple multiple-pipe–multiple-consumer example .................. 63

Chapter 7: Conclusions and recommendations ................................................ 75

Summary ............................................................................................................. 75

Conclusions ........................................................................................................ 76

Recommendations ............................................................................................. 76

Literature cited ....................................................................................................... 79

Appendix A: Approximation of the friction factor ........................................... 81

Appendix B: Computer program listings ......................................................... 85

Abstract ................................................................................................................... 95

ILLUSTRATIONS

Figure

1. Assumed annual load curve...................................................................... 13

2. Load duration curve ................................................................................... 14

3. Schematic of a hydronic heating system radiator .................................. 18

4. Errors from approximating the logarithmic mean temperature

difference with arithmetic and geometric mean temperature

differences .................................................................................................. 22
iv

Figure (cont‘d) Page

5. Results for the return temperature models ............................................. 25

6. Mass flow rate function with and without consumer models

Tr = 60 oC .................................................................................................... 28

7. Mass flow rate function with and without consumer models

Tr = 55 oC .................................................................................................... 29

8. Hypothetical pressure distribution under high and low flow

conditions and absolute pressure constraints....................................... 52

9. Method A ..................................................................................................... 56

10. Method B ...................................................................................................... 57

11. Slope of the heat loss cost term as a function of pipe diameter ........... 58

12. Slope of the pressure loss as a function of pipe diameter ..................... 58

13. Hypothetical branch-and-bound problem .............................................. 61

14. Simple multiple-pipe–multiple-consumer system................................. 63

TABLES

Table

1. Pressure drops and costs for discrete pipe sizes under maximum flow

conditions ................................................................................................... 16

2. Errors from approximating the logarithmic mean temperature

difference.................................................................................................... 22

3. Return temperatures and flow rates calculated with the consumer

models for n1 = n2 = 1.3 ............................................................................ 26

4. Pressure drops and costs for discrete pipe sizes under maximum

flow conditions with and without the consumer model..................... 30

5. Assigned values for nodes of example in Figure 14 .............................. 64

6. Assigned values for the pipe segments in the example of Figure 14 .. 64

7. Parameter values and optimal independent diameters for example

of Figure 14 ................................................................................................ 65

8. Standard metric steel pipe sizes ............................................................... 66

9. Discrete bounding diameters and variable costs for the example ....... 66

10. Pressure losses for the pipe segments in the example........................... 67

11. Heating plant pressure increase required by eq 4-2 for each

consumer .................................................................................................... 67

12. Pipe size combinations for the example .................................................. 69

13. Constraint satisfaction and costs for the remaining combinations ..... 70

14. Pressure levels in the piping network ..................................................... 73
NOMENCLATURE

a coefficient in friction factor equation determined by curve fitting

(dimensionless)

A1 empirical parameter related to pump costs ($/pump)

A2 empirical parameter related to pump costs ($/W)

A3 parameter related to pipe costs ($/m)

A4 parameter related to pipe costs ($/m2)

A5 defined parameter related to pumping energy at the design conditions

($ m5+b+c)

A6 defined parameter related to the density and dynamic viscosity of the

supply and return fluids at the design conditions (m3+c sc/kg1+c)

A7 defined parameter related to the density and dynamic viscosity of the

supply and return fluids (m6+c sc/kg2+c)

A7,d A7 evaluated at the design condition (m6+c sc/kg2+c)

A8 defined parameter related to the fluid density (m6/kg2)

A9 defined parameter related to maintenance and repair costs ($/m)

A10 defined parameter related to heat losses (m)

A11 defined parameter related to pipe costs (m)

A12 defined parameter related to pumping energy (mb+1)

A13 over-design factor for the consumer’s radiators (dimensionless)

A14 empirical coefficient related to the annual load curve (dimensionless)

A15 empirical coefficient related to the annual load curve (dimensionless)

A16 defined parameter related to heat losses ($/[°C hr])

Aη empirical coefficient related to pumping (dimensionless)

Ahe expression that relates the fluid properties and physical properties of

the heat exchanger to the pressure drop at a given flow rate (kg1–β/

m–s2–β])

Am&r annual maintenance and repair rate as a fraction of initial capital costs

(dimensionless)

At number of hours per year (8760)

AF approach factor for the heat exchanger (dimensionless)

b exponent in friction factor equation determined by curve fitting

(dimensionless)

c exponent in friction factor equation determined by curve fitting

(dimensionless)

Ccpe cost of pumping energy dissipated in the consumer’s heat exchanger

and control valve ($)

Ce cost of electricity ($/Wh)

Ch cost of heat ($/Wh)

Chl cost of heat losses ($)

Cm&r cost of maintenance and repair ($)

Cp cost of the pipes ($)

cp specific heat of water at constant pressure (kJ/kg °C)

Cpe cost of pumping energy ($)

Cpp capital costs of pipes and pumps ($)

v
Cfixed fixed cost of pipes and pumps, and the maintenance and repair on this

portion of their costs ($)

Cpv diameter variable cost of pipes and the maintenance and repair on that

portion of pipe cost ($)

Cpev diameter variable cost of pumps and pumping energy attributable to

piping pressure losses, and the maintenance and repair on that portion

of pump costs ($)

Cpvc variable cost of pumps attributable to the pressure losses at the

consumer ($)

Cpumps cost of the pumps ($)

 

Cpt ′ total diameter variable pipe costs for the system ($)

Cst total system cost ($)

Ct total system owning and operating cost ($)

 

Ct′ portion of total system cost that is dependent on the pipe diameter ($)

d inside diameter of pipe (m)

Do outer diameter of insulation (m)

f friction factor (dimensionless)

g acceleration of gravity (9.8 m/s2)

gl equality constant numbers l (dimensions vary)

hl inequality constant numbers l (dimensions vary)

Hp burial depth to pipe centerline (m)

I1 defined integral parameter related to heat loss costs ($)

I2 defined integral parameter related to pumping energy ($ m5+b+c)

I3 defined integral parameter related to pumping energy and maintenance and repair costs ($ m5+b+c)

ki insulation thermal conductivity (W/m °C)

ks soil thermal conductivity (W/m °C)

L pipe length (m)

m ̇ mass flow rate (kg/s)

 

m ̇ d maximum (design) mass flow rate (kg/s)

n total number of heat consumers

n1 empirically determined exponent in radiator equation (dimensionless)

n2 empirically determined exponent in the radiator equation (dimensionless)

np number of pumps

npi total number of pipe segments (measured in supply and return pipe

pairs)

Pa atmospheric pressure (≈ 105 N/m2)

Pasa minimum safety margin above atmospheric pressure (N/m2)

Php,s absolute pressure in supply pipe at heating plant (N/m2)

Php,r pressure in the return line at the inlet to pump (N/m2)

PI pressure at the inlet to the pipe segment (x=0) (N/m2)

Pmax maximum absolute pressure for the piping system being used (N/m2)

PNPSH minimum allowable pressure at the pump inlet due to NPSH requirements (N/m2)

vi
PPa actual pumping power required, including pump and pump driver

inefficiencies (W)

PPf frictional pumping power, exclusive of pump and pump driver inefficiencies (W)

Ps absolute pressure in supply pipe at point in question (N/m2)

Psaf minimum allowable safety margin on saturation pressure requirements (N/m2)

PVFe present value factor for electrical energy (dimensionless)

PVFh present value factor for heat (dimensionless)

PVFm&r present value factor for maintenance and repair costs (dimensionless)

Px,sat saturation pressure of the liquid at point x within the pipe segment

(N/m2)

Px pressure at point x (N/m2)

q heat output from the radiator (W)

Qhl rate of heat loss (W)

Re Reynolds number for the pipe flow (dimensionless)

Ro overall resistance to heat transfer (W/m °C)

RR relative roughness of pipe (dimensionless)

t time of year (hr)

tu equivalent full load utilization time (hr)

T water temperature (°C)

Ta indoor air temperature (°C)

Tao air temperature at radiator outlet (°C)

Tavg average of supply and return temperature (°C)

Tg soil temperature (°C)

Tm mean soil temperature (°C)

Tma arithmetic mean temperature difference (°C)

Tmg geometric mean temperature difference (°C)

Tml logarithmic mean temperature difference (°C)

Tp pipe outer surface temperature (°C)

Tr return temperature (°C)

Ts supply temperature (°C)

Tsmin minimum supply temperature required by next consumer (°C)

v flow velocity (m/s)

w1 weight of term 1 (dimensionless)

w2 weight of term 2 (dimensionless)

x position along the pipe in the direction of flow with x = 0 being defined

as the inlet end to the pipe segment in question (m)

z elevation at point in question relative to heating plant (m)

Greek

β exponent yielding the appropriate mass flow rate dependency for the

heat exchanger (dimensionless)

ε absolute roughness of the piping (m)

εa relative approximation error for the arithmetic mean temperature

difference (dimensionless)

vii
εg relative approximation error for the geometric mean temperature

difference (dimensionless)

γ ki

/ks (dimensionless)

ηp pump efficiency (dimensionless)

ηpm efficiency of the motor driving the pump (dimensionless)

μ dynamic viscosity (Pa s)

ρ fluid density (kg/m3)

ρd fluid density at design conditions (kg/m3)

∆Pcv pressure drop in the control valves (N/m2)

∆Pcvm minimum pressure drop in the control valve (N/m2)

∆Pcvs,i slack variable for consumer control valve pressure losses (N/m2)

∆Pd total pressure drop (supply and return) at design flow rate (N/m2)

∆Phe pressure drop in the heat exchangers (N/m2)

∆Php pressure increase across the pump (N/m2)

∆Pncv,i non-control-valve pressure losses (N/m2)

∆Pr,j pressure loss in the servicing return line j (N/m2)

∆Pr pressure drop in the return piping (N/m)

∆Ps,j pressure loss in supply pipe j (N/m2)

∆Ps pressure drop in the supply piping (N/m2)

∆Ps&r combined pressure loss of supply and return (N/m2)

∆Tra smallest temperature difference between fluids in the consumer’s heat

exchanger (°C)

∆Tsa greatest temperature difference between fluids in the consumer’s heat

exchanger (°C)

∆xi insulation thickness (m)

(∂z/∂x) partial derivative of the elevation of the pipe with respect to its

position (dimensionless)

(dP/dx)d hydrodynamic pressure gradient (N/m3)

(dP/dx)h hydrostatic pressure gradient (N/m3)

Subscripts

d design maximum load condition

i consumer index

j pipe segment index

sp straight pipe heat exchanger

he conditions within the heat exchanger or physical parameters of the

heat exchanger

0 “design” condition for the consumer’s heat exchanger, usually the

maximum load condition at maximum supply temperature

1 condition for the consumer’s heat exchanger of actual supply

temperature with the flow rate as determined under the design

condition

2 any actual operating condition in the consumer’s heat exchanger

viii v
CHAPTER 1: INTRODUCTION

District heating is the practice of heating multiple buildings from a single heating

plant. Heat is conveyed to the buildings by means of hot water or steam. District

heating systems offer enormous potential for energy conservation, in addition to the

advantages of fuel flexibility and reduced environmental impact. For these reasons

district heating has been used extensively in Europe with favorable results. For

example, in Denmark district heating serves 42% of the demand for space and hot

tap water heating (NRC 1985).

In the United States district heating is much less widespread, accounting for

about 4% of the space and hot tap water heating (NRC 1985). A few cities have

systems, as well as a number of college campuses and other large institutions. With

approximately 6000 miles (10,000 km) of district heating piping in place (Segan and

Chen 1984), the Department of Defense is the single largest user of the technology

within the United States. A major barrier to more widespread use of district heating

in the United States is the high capital cost of the piping required to convey the heat

to the buildings. The piping system is most often the major cost of district heating.

However, the lack of development of district heating in the United States is often

attributed to “institutional barriers.” Such barriers, where they truly exist, would be

significantly weakened if not removed should the economics become more favorable.

CURRENT DESIGN PRACTICE

Because the hot water or steam piping networks represent such a major portion

of the capital costs, they also represent an opportunity for significant cost savings by

optimal design. Despite this, in practice little effort is expended either here or in

Europe to ensure that proposed designs reduce costs as much as possible. Currently,

most designs are based on previous experience and often may be far from optimal.

Rules of thumb are commonly used, as are design guides developed for other

purposes, such as for plumbing within buildings.

To achieve an optimal design with minimum life cycle costs, all major costs

associated with constructing and operating the system must be considered. Capital

costs for piping and installation vary widely and must be determined for each case.

Operating costs strongly depend on the nature of the load and the load management

strategy adopted. For these reasons, it is impossible to develop guidance that can be

applied universally to obtain designs that are sufficiently close to the lowest life

cycle cost.

As with most areas in the practice of engineering, computer-aided design

methods are becoming more widespread in district heating system design. The use

of such methods allows the rapid evaluation of many alternate designs, a formidable

task if carried out without such methods. A number of computer-aided design

methods are available for thermal piping networks (Reisman 1985, Rasmussen and

Optimal Design of Piping Systems for

District Heating

GARY PHETTEPLACE
Lund, undated, Cowiconsult 1985, Hart and Ponsness 1992). Most of the computeraided design methods that have been developed are proprietary and thus any

optimization that they profess to make is not open to inspection or modification. An

optimal design model that is open to inspection and modification is therefore

needed. The objective of this work is the development of such a model.

OPTIMIZATION IN DISTRICT HEATING SYSTEM DESIGN

Determination of pipe sizes is one of the major decisions that the designer of a

district heating system faces. Other critical decisions that must be made are heat

source, distribution media, distribution temperatures and load management strategy. Many of these factors have been addressed by previous studies. The emphasis

of this work will be on determining pipe size.

Currently, pipe sizes are usually determined on the basis of simple criteria, such

as maximum pressure loss or maximum flow velocity. A number of investigators

have addressed the issue of pipe size determination, trying to improve on these

simple criteria. Aamot and Phetteplace (1976) presented a method that relies on

establishing the ratio between the heat losses and the pumping cost and then finds

the lowest cost pipe diameter by minimizing the sum of capital, heat loss and

pumping costs. Their work only addressed a single pipe segment and did not

include the effect of varying load over the yearly cycle. Szepe and Calm (1979)

presented a model for single pipe segments that neglected heat losses and time

varying loads, but used geometric programming theory to achieve additional

insight into their simplified problem. In later work, Phetteplace (1981) included the

effect of annual load variations, but only single pipe segments were addressed.

Frederiksen (1982) provided a detailed analysis of the heat generating station and

the consumer’s systems, but simplified the transmission network to a single supply

and return pipe.

A number of investigators have addressed multiple pipe networks. Of course, a

great deal of work has been done for water distribution systems where the problem

is much simpler owing to the lack of heat losses and load variation with temperature

as well as mass flow rate. Marconcini and Neri (1979) described a model that

calculates the flow rate, pressure and temperature in networks of steam pipes. They

discussed the effect that pipe diameter has on operation, but did not offer any

methodology for selecting diameter values.

Stoner (1974) discussed models that are capable of modeling either steam or water

networks. Although the models do not determine optimum diameters, he gave a

procedure for achieving an optimal design by sensitivity analysis, but did not discuss how this process would be accomplished for networks of more than one pipe.

Zinger et al. (1976) described a computer program for calculating flows and

pressure levels in branched networks of hot water pipes. Their program accounts for

pressure drops in the consumers’ equipment and throttling devices placed in the

network. Diameters are assumed to be known and they did not discuss how to

determine them.

Morofsky and Verma (1979) developed a feasibility analysis and costing tool for

district energy systems, not intended for detailed design. They found the appropriate pipe sizes by finding those that absorbed all of the available pressure difference.

They started the search for pipe size at the smallest available discrete pipe diameter

and then calculated pressure losses. If the pressure losses were more than the

available pressure difference, they increased the pipe size to the next discrete size

and repeated the calculation. They proceeded in this fashion until they reached a

discrete pipe diameter that did not result in pressure losses greater than the available

pressure difference.

2
McDonald and Bloomster (1977) discussed a model for laying out and sizing the

piping network for a city heated with geothermal water. Pipe diameter is determined using a “simple search” of feasible pipe sizes by minimizing the sum of the

annual capital cost, heat loss cost and pumping cost. They provided no information

on how to handle network constraints or consider annual load variations.

Bøhm (1986) noted that, in the case of consumers directly connected to the

network, the “classical” approach of determining the optimal diameter by finding

the minimum of the sum of the capital, heat loss and pumping costs results in

pressures that are too high at the heating plant. He suggested the use of Munser’s

(1980) method, which proportions the total available pressure loss in a network

using the equation

 

1

0

1

1 3

1 1

∆

∆

P

P

L

L m

i m

n =

∑ 

 

=  i i /  ̇  ̇

(1-1)

where∆P1 = pressure loss in pipe number 1 (N/m2)

∆P0 = total pressure loss in the pipe network (N/m2)

L1 = length of pipe number 1 (m)

Li = length of pipe i (m)

 1 m ̇ = mass flow rate in pipe number 1(kg/s)

 i m ̇ = mass flow rate in pipe i (kg/s)

n = total number of pipes

i = pipe index.

Equation 1-1 is intended for use on “linear networks” that do not have branches.

Koskelainen (1980) developed a method that is able to solve for optimal diameters

in a branched network. His method consists of successively assuming that the

objective function and constraints locally are linear and repeatedly solving the

problem with a linear programming algorithm. He gives an example where his

“optimal” network has a cost that is 16.4% less than one sized using a head loss

design rule.

In this work we develop a rational design method that yields the optimal pipe

sizes for an application based on case-specific parameter values. This method allows

for the inclusion of all major costs and can account for such factors as escalation of

energy prices, seasonal energy costs, increases in heat losses over system life,

variation in seasonal heat demand, load management strategy, the effect of the heat

consumer, etc. Each of the major constraints on the design of a realistic district

heating network is derived and considered. This method is felt to be practical for

sizing much of the piping of a district heating system.

We begin our study in Chapter 2 by first finding a suitable method for determining the optimal size for a single pipe, independent of any others. In developing

this method, we endeavor to keep the formulation as simple as possible, yet

complete and accurate enough for design calculations. We make use of geometric

programming theory to identify a lower bounding problem that can be used to guide

us to our solution. At the end of Chapter 2 is an example that shows a 17% saving

in life cycle cost.

In Chapter 3 we study the heat consumer and the effect he has on the piping

system. We develop a new model for the consumer’s heat exchanger, which uses the

geometric mean temperature difference as an approximation for the logarithmic

mean temperature difference, thus allowing for an explicit expression for return

temperature. We integrate this consumer model with our single pipe model of

Chapter 2 and show what effect the consumer has on the system.

In Chapter 4 we develop the constraints for systems with multiple pipes and

3
consumers. Both absolute and differential pressure constraints are derived and

where possible strategies are given to allow for constraint satisfaction at all points

implicitly without considering every point in the system.

After a brief review of general methods for constrained nonlinear optimization

techniques at the beginning of Chapter 5, our general solution strategy is developed

for systems with multiple pipes and consumers. The method makes use of the

solution to the problem, unconfined by the network constraint requirements.

Monotonicity analysis is used to prove activity of some of the constraints and thus

simplify the problem somewhat. The result is used as a starting point for two

methods proposed to find a solution to the constrained problem with continuous

values for some of the pipe diameters. Finally, the branch-and-bound technique is

used to find a design with discrete values for all the pipe diameters.

In Chapter 6 we work a simple example with only four consumers and seven pipe

segments. The example illustrates the use of our method and also shows how the

branch-and-bound technique can be used to quickly eliminate candidate designs.

In Chapter 7 is a summary of our results and offers some conclusions and

suggestions for further study.

Because of the inordinate number of variables and parameters involved in the

analysis that follows, in choosing symbols for them, an attempt has been made to

make their meaning as intuitive as possible. Where accepted symbols exist they have

been used to the maximum extent possible. Where it has been mathematically

convenient to represent quantities that may have no particular physical significance

by a symbol, subscripted A’s have been used for sums, products and quotients and

I’s have been used for integrals.

4
CHAPTER 2: OPTIMAL PIPE DIAMETER FOR A

SINGLE PIPE SEGMENT

To find the optimal diameter for a single pair of supply and return pipes, we need

to consider the costs involved and minimize their sum with respect to the pipe

diameter. The cost minimization is done for the life cycle of the system using a net

present value approach. Some types of heat distribution systems may have a salvage

value, while others will, in fact, have a disposal cost associated with the end of their

useful lifetime. Since these will in general be mild functions of the pipe diameter,

they will not significantly affect the optimal pipe diameter and thus will not be

treated here. With these limitations in mind, our objective function, the total life

cycle cost, becomes

min. Ct = Chl + Cpe + Cpp + Cm&r (2-1)

where Ct = total system owning and operating cost ($)

Chl = cost of heat losses ($)

Cpe = cost of pumping energy ($)

Cpp = capital costs of pipes and pumps ($)

Cm&r = cost of maintenance and repair ($).

Now let’s look at each of the costs in eq 2-1 in detail, starting with the cost of heat

losses.

COST OF HEAT LOSS

The basic form of the heat loss cost is

 

C PVF C Q t hl h h hl yr = ∫ d (2-2)

where Ch = cost of heat ($/Wh)

PVFh = present value factor for heat (dimensionless)

Qhl = rate of heat loss (W)

t = time of year (hr [0 ≤ t ≤ 8760]).

In the most general case, the cost of heat Ch can be a function of time because of

seasonal usage rates. The rate of heat loss Qhl will also be a function of time over the

yearly cycle. In fact, deterioration of the thermal insulation will result in increasing

heat losses as the system ages. This can not be incorporated directly into the

formulation as given above, but could be allowed for by using an appropriate

escalation factor in the present value factor for heat costs PVFh.

The only variable defined above that is dependent on our decision variable, the

pipe diameter d, is the heat loss rate itself Qhl. For a single buried pipe, the

relationship is

Qhl = L(Tp – Tg)/Ro (2-3)

where Tp = pipe outer surface temperature (°C)

Tg = soil temperature (°C)

Ro = overall resistance to heat transfer (W/m °C)

L = pipe length (m).

The dependence on d is from the overall thermal resistance Ro. This resistance is

5
found by adding the resistance attributable to the soil to that resulting from the pipe

insulation (Phetteplace and Meyer 1990). After simplification the following result is

obtained

Ro = ln[ (Do/d) (4Hp/Do)

γ

 ]/2ki (2-4)

where γ = ki/ks (dimensionless)

ki = insulation thermal conductivity (W/m °C)

ks = soil thermal conductivity (W/m °C)

Do = outer diameter of insulation (m)

Hp = burial depth to pipe centerline (m).

In this form it becomes easy to see how each factor affects this parameter. The

(4Hp/Do)

γ

 factor represents the contribution of the soil to the overall thermal

resistance. If γ << 1, that is, if the soil conductivity is much greater than the insulation

thermal conductivity, then this factor will be close to unity and the overall thermal

resistance reduces to the thermal resistance of the insulation alone.

To obtain a simpler form for the cost of heat loss, we make the following

assumptions:

1. That the soil temperature at the pipe depth varies sinusoidally over the yearly

cycle about a mean temperature.

2. That the cost of heat is constant over the yearly cycle. This does not limit us to

fixed heat cost over the life of the system, since escalation factors may be used.

3. That the outer surface temperature of the carrier pipe is equal to the temperature of the carrier medium.

The result of these assumptions is the following form for the cost of heat loss

Chl = I1/ln(A10/d) (2-5)

where I1 = PVFh L 4πki (∫ Ch Tavg dt – At Ch Tm) ($)

A10 = Do(4Hp/Do)

γ

 (m)

Tavg = (Ts+ Tr)/2 (°C)

Tm = mean soil temperature (°C)

Tr = return temperature (°C)

Ts = supply temperature (°C)

At = number of hours per year (8760).

COST OF PUMPING

Now let’s consider the pumping costs. A cost is associated with the electrical

energy input to drive the pumps. The portion of this energy that results in frictional

heating of the fluid in the pipes is recovered as heat. In general the value of the heat

recovered will, of course, be less than the value of the electrical energy input to drive

the pumps. It can be significant, however, and therefore it has been included here.

Thus, we have the following for the pumping cost

Cpe = PVFe ∫ CePPa dt – PVFh ∫ ChPPf dt (2-6)

where PVFe = present value factor for electrical energy (dimensionless)

Ce = cost of electricity ($/Wh)

PPa = actual pumping power required, including pump and pump driver

inefficiencies (W)

yr yr

6
PPf = frictional pumping power, exclusive of pump and pump driver

inefficiencies (W).

The first integral term represents the total cost of electrical energy input to drive

the pumps. The second integral term is the value of heat recovered in frictional

heating of the fluid. The actual pumping power and the fluid frictional portion are

related as follows

PPa = PPf /ηpηpm (2-7)

where ηp = pump efficiency (dimensionless)

ηpm = efficiency of the motor driving the pump (dimensionless).

The pumping power of a closed system with return lines will not be affected by

elevation differences within a network and therefore they need not be considered

here. Elevation differences will, however, become a factor in determining the

absolute pressure level within a network. A constraint will arise owing to absolute

pressure limitations of the piping. This will be addressed later.

Now we assume that the product of the pump and motor efficiency can be

expressed as a function of the fraction of maximum volumetric flow. A similar

approach was used by Phetteplace (1981) based on data from Gartman (1970). This

gives an expression of the form

ηpηpm = Aη ( 

m ̇ /ρ) (ρd/ 

m ̇ d ) (2-8)

where Aη = empirical coefficient (dimensionless)

 

m ̇ = mass flow rate (kg/s)

 

m ̇ d = maximum (design) mass flow rate (kg/s)

ρ = fluid density (kg/m3)

ρd = fluid density at design conditions (kg/m3).

The frictional pumping in the supply or return line is given by

PPf = 2 (2/π)2 f L ρ–2 

m ̇ 3 d–5 (2-9)

where f is a friction factor (dimensionless).

Using the above expression for both the supply and return pipes, we substitute

the results, along with our earlier result for PPa, and our expression for PPf, into our

original expression for the pumping energy cost and simplify to obtain

 

C dA C m

A m

PVF

PVF C Am f t

yr

pe e d

d

h

e

= − h 8 d 









 −

∫ 5

11

ρ 3

η ρ

 ̇

 ̇  ̇ (2-10)

where A11 = (4/π)

2 PVFe L (m)

A8 = 

ρ ρ s r

− − ( ) + 2 2 /2 (s and r subscripts denote supply and return conditions

respectively) (m6/kg2).

Now we would like to find a simple function to approximate the friction factor f

over a range of interest. A simple power function relationship would be desirable to

keep the number of terms to a minimum and thus not complicate the above

expression further. The form of such a function is suggested by the dimensionless

groups of the Moody diagram (Jeppson 1976). A method of finding an approximat7
ing function by converting to logarithmic variables and using a least-squares curve

fit was developed (Appendix A) to fit an equation of the form

f = a (ε/d)b Rec (2-11)

where a, b and c = coefficients determined by curve fitting (dimensionless)

ε = absolute roughness of the piping (m)

Re = Reynolds number for the pipe flow (dimensionless).

As an example, the following coefficients are obtained over the parameter range

given

a = 0.119

b = 0.152

c = –0.0568

for

50 ≤ T ≤ 130

0.5 ≤ v ≤ 4.5

0.050 ≤ d ≤ 0.770

where T = water temperature (°C)

v = flow velocity (m/s).

When compared to the Colebrook and White equation (Jeppson 1976), the

maximum error of this approximation is 6.9%, with the average error over the range

given being only 1.1%. If more accuracy is required, a much better result could be

obtained by narrowing the parameter ranges. Some examples of results for other

parameter sets are given in Appendix A. The coefficients will be carried for the

general case in the derivations following to allow for values obtained with other

parameter sets.

By expressing the Reynolds number as a function of the quantities previously

used in the formulation, our equation for the friction factor becomes

f = a (4/μπ)

c εb d–(b+c)

 

m ̇ c (2-12)

where μ is dynamic viscosity (Pa s).

Now if we substitute this result into our expression for the cost of pumping

energy and simplify, we obtain

Cpe = I2d –(5+b+c) (2-13)

where

 

I A C m

A m

PVF

PVF

C Am t

yr

c bc

2 12

3 5 = − 







 ∫  ( ) e d + ++

d

h

e h 7 d  ̇

 ̇  ̇ $ ρ

η ρ m

A7 = [(ρ –2μ –c)s + (ρ–2μ–c)

r

]/2 (m6+c sc/kg2+c)

A12 = a (4/π)2+c εb PVFe L (mb+1).

8
COST OF PIPES AND PUMPS

Now we need to find expressions for the capital cost of the pipes and pumps. In

general, for the entire system the pump capital costs will be assumed to be of the

following form

Cpumps = A1np + A2(

 

m ̇

d /ρd) ∆Pd (2-14)

where A1 = empirical constant ($/pump)

A2 = empirical constant ($/W)

np = number of pumps

∆Pd = total pressure drop (supply and return) at design flow rate (N/m2).

The total pressure drop at maximum flow conditions is given by

∆Pd = a εb (4/π)

2+c A6 

m ̇ d

2+c L d–(5+b+c) (2-15)

where A6 = [(ρ–1 μ–c

)

d,s+ (ρ–1 μ–c

)d,r]/2 (m3+c sc/kg1+c).

So, our pump cost becomes

Cpumps = A1 np + A5 d –(5+b+c) (2-16)

where A5 = A2a εb(4/π)2+c A7,d 

m ̇ d 3+c L ($ m5+b+c)

A7,d = A7 evaluated at the design condition (m6+c sc/kg2+c).

For the capital cost of the supply and return piping, including installation, we

assume the following form

Cpipes = (A3 + A4d)L (2-17)

where A3 = empirical constant ($/m)

A4 = empirical constant ($/m2).

COST OF MAINTENANCE AND REPAIR

The cost of maintenance and repair is assumed to be of the following form

Cm&r = PVFm&r Am&r Cpp (2-18)

where Am&r = annual maintenance and repair rate as a fraction of initial capital cost

(dimensionless)

PVFm&r = present value factor for maintenance and repair costs (dimensionless).

TOTAL COST

With each of the component costs defined, our expression for the objective

function, the total cost Ct

, becomes

9
 

Ct′ = I1/[ln (A10/d)] + I3 d –(5+b+c)

 + A9 d (2-19)

where 

Ct′ = Ct

 – [(1+ PVFm&rAm&r)(A1np+ A3L)]

I

3 = I2 + (1+PVFm&rAm&r) A5

A9 = (1+PVFm&rAm&r) A4L.

Minimizing 

Ct′ with respect to d is, of course, equivalent to minimizing the

original total cost function Ct. Therefore, we can work with 

Ct′ for convenience. If

we neglect the first term, which represents the cost of heat losses, we have a

geometric programming problem (Papalambros and Wilde 1988) with zero degrees

of difficulty. Without specifying the parameter values, we see from inspection that

the weights of the two remaining terms will be

w1 = 1/(6 + b + c) and w2 = (5 + b + c)/(6 + b + c).

With heat losses neglected, at the optimum pipe diameter the variable costs

associated with pumping are 1/(6 + b + c) of the total variable costs. The variable

costs attributable to pipe capital and maintenance costs are the remaining portion.

Here, the variable costs represent that portion which is a function of our decision

variable, the pipe diameter. Also note that the pumping costs include the variable

portion of the capital cost of the pumps and the maintenance associated with that

portion, as well as the pumping energy costs.

Considering a more specific case, if the values of parameters b and c found in the

example given for eq 2-11 (b = 0.152, c = –0.0568) are used, we find the following

values for the weights

w1 = 16.4%

w2 = 83.6%.

These results vary very little over the range of values found for b and c in

Appendix A. Thus, when heat losses are neglected, we find this very simple solution

is applicable in most cases.

Once values for the remaining parameters are known, the pipe diameter is found

by using the equations given above and the two terms of the objective function

remaining. The resulting expression is

d = [(5 + b + c)(I

3/A9)][1/(6+b+c)]. (2-20)

It should be noted that this solution obtained using geometric programming theory

also could have been easily obtained using classical differential methods, as used

later. The advantage of the geometric programming method is that it ensures that a

global rather than a local minimum has been found. Differential methods only ensure

a local extremum and require the evaluation of second order terms to determine the

nature of the extremum, i.e., maximum or minimum.

To arrive at this simple expression for the pipe diameter, we have neglected the

heat losses. Because the cost of heat losses will always be greater than zero, we have

constructed a lower bounding problem for our original problem by neglecting this

cost, i.e.

I

3 d–(5+b+c)

 + A9 d ≤ I

1/[ln(A10/d)] + I3 d–(5+b+c)

 + A9 d.

The cost of any design that includes heat losses can never be less than the same

10
design excluding the cost of heat losses. And, since we have found the optimum

design (lowest cost) neglecting heat losses, we now know that no design can achieve

a lower cost when heat losses are included. This simple result can be very useful. It

may be possible to find a design, not necessarily known to be optimal, whose cost

including heat losses is acceptably close to that of the optimal design for the lower

bounding problem.

The solution to the complete problem including heat losses is slightly more

complicated, but is easily obtained. To find the extremum of the total variable cost

function 

C′

t , we simply take its partial derivative with respect to d and set the result

to zero. Before proceeding to do so, however, we must take note of the value of A10

in the heat loss term being a function of the pipe diameter. This results from the outer

diameter of the pipe being a function of pipe diameter and insulation thickness. The

appropriate insulation thickness is determined by a separate optimization procedure that would consider the insulation and jacket material costs and the cost of heat

loss. As a result of this separate “sub-optimization,” the insulation thickness

becomes a function of the pipe diameter.

For a given set of operating conditions and economic data, the optimal insulation

thickness can be found as a function of the pipe diameter. Here, for the sake of

simplicity, we will assume that the insulation thickness is fixed. We then find an

expression for A10 as follows

A10 = (d + 2∆xi)1–γ (4Hp)γ (2-21)

where ∆xi is the insulation thickness (m). In turn we approximate this expression by

one of the following form

A10 = (d1–γ + (2∆xi)1–γ) (4Hp)γ. (2-22)

For a typical set of parameter values

ki = 0.030 W/m °C

ks = 1.3 W/m °C

∆xi = 0.050 m.

This approximation is within 2% for values of d from 0.025 to 1.0 m. Using this

approximation for A10, we obtain the following equation for 

C′

t

 

C′

t = [I1/ln((4Hp)

γ

 (d–γ + (2∆xi)1–γ d–1))] + I3d–(5+b+c) + A9d. (2-23)

If we take the partial derivative of 

C′

t and set the result to zero, we have

0 = {I1/d [ln((4Hp)

γ

(d–γ+ (2∆xi)

1–γ d–1))]2}{1–((1–γ)/(1+(2∆xi

)

1–γ d γ–1))}

– (5 + b + c)I3d–(6+b+c) + A9. (2-24)

This equation can not be solved explicitly for d. A solution can be obtained by

using a root-finder technique. An initial estimate needed for the solution can be

found by using the value of d obtained from the solution to the lower bounding

problem, which neglects heat losses. The cost associated with this optimal design

that neglects heat losses also provides us with a global lower bound on the actual

cost. In the next section, we consider the solution to a representative case.

11
A SIMPLE EXAMPLE

To illustrate the application of the method developed in the previous section, an

example using realistic parameter values is presented in this section. In addition, the

results obtained are compared to those obtained by using a common rule of thumb.

Before we can proceed with the calculation of the parameters, we need to define the

economic and technical assumptions upon which our solution is based.

Economic assumptions

The assumptions concerning economic conditions of an application are the most

controversial in an analysis of this type. Here, we will endeavor to select conditions

that are felt to be representative of the majority of applications rather than any

specific application. The analysis is quite sensitive to these economic assumptions,

so the reader is cautioned that specific information must be obtained before

applying these results. This can not be overemphasized.

The most significant cost to be considered is the capital cost of the piping system.

This cost is also highly variable, depending on the piping system used and, above

all, the site conditions. Our pipe cost equation is from Phetteplace (1981) and has

been adjusted for inflation from 1980 to 1988 using factors from R.S. Means Co.

(1987). The following expression results

Cp = 218 L + 2180 Ld. (2-25)

The capital costs of pumps is also taken from Phetteplace (1981), with cost adjustments made as indicated above to arrive at

Cpumps = 1060 np + 0.242 ∆Pd 

m ̇ d/ρd. (2-26)

The present value factors (PVF’s) are assumed to be equal for electrical cost, heat

cost, and maintenance and repair costs. As noted earlier, these PVF’s could be

modified to allow for escalation of energy, labor or material costs and could even be

modified to allow for increasing heat loss over the life of the system. For simplicity,

however, such modifications are not made here. We assume an interest rate of 10%

per annum and a system lifetime of 25 years in calculating the PVF’s using the

following expression (CRC Press 1987)

PVF = [1 – (1+i)

–n]/i = [1 – (1.10)–25]/0.10 = 9.08. (2-27)

We assume that both heat and electricity costs are constant over the year,

although the formulation allows for varying rates over the yearly cycle. Again, these

costs are highly variable, dependent mainly on the sources of the energy and the

values given are not for a specific application, although it is felt that they are

representative. The assumed costs are

Ce = 7.0 × 10–5 $/Wh ($0.07/kWh)

Ch = 3.4 × 10–5 $/Wh ($10/106Btu).

The rate of maintenance and repair on the system is also taken from Phetteplace

(1981) as 2% of the capital cost per year. This factor is assumed to apply to both the

piping system and the pumps. Again, note that maintenance and repair increasing

with component age could be easily accounted for using escalation factors to adjust

PVFm&r , although this is not done here.

12
Technical assumptions

The major technical assumptions that we make are related to the heat load

characteristics and the method by which the heat load is met. In district heating

systems, the amount of heat supplied can be varied to accommodate varying

demand by adjusting either the supply temperature or flow rate. However, certain

constraints imposed by consumer equipment and minimum temperature requirements must be observed. In larger systems, both the supply temperature and

flow rate are varied over the course of the year. In small systems, which must adopt

simpler control strategies, often only the flow rate is varied. For the sake of

simplicity, we assume the latter here and assume that the supply and return

temperatures remain fixed over the yearly cycle. This is never actually the case, but

for system design, it is felt that this is an appropriate simplifying assumption for a

first analysis. Ideally, for example, the supply temperature, return temperature and

flow rate at the consumer would be determined by the heat transfer characteristics

of his heat exchanger under the prevailing load. A model that simulates the

consumer’s heat exchanger is developed in Chapter 3, but here it is not considered.

The actual heat load in district heating systems has several major components. A

detailed treatment of the heat load would be difficult and is not warranted for design

purposes. For an excellent treatment of the actual heat loads in operating district

heating systems, see Werner (1984). The assumption we make here for design

purposes, that the heat load can be approximated as sinusoidal, is supported by the

data presented by Werner (1984) as well as by the data of Phetteplace et al. (1981).

We assume here that the heat load varies sinusoidally from a minimum of 15% of its

maximum value to its maximum value. The assumed minimum load of 15% would

result primarily from hot tap water use and heat losses from the pipelines. Thus,

with our assumption of constant supply and return temperatures, the mass flow rate

as a function of time is

 

m ̇ / 

m ̇ d = 0.575 + [0.425 cos(2πt/8760)]. (2-28)

The ratio of the mass flow rate to its maximum value as determined by eq 2-28 is

shown in Figure 1. With this simple function for the load curve, it is easy to determine

the shape of what is normally referred to as the “load duration curve.” Since eq 2-

Figure 1. Assumed annual load curve.

1.0

0.8

0.6

0.4

0.2

0

d m/m· ·

2190 4380 6570 8760

t (hr)

13
Figure 2. Load duration curve.

1.0

0.8

0.6

0.4

0.2

0

d m/m· ·

2190 4380 6570 8760

t (hr)

28 is an even function about its midpoint of 8760/2 = 4380 hours, the number of hours

for which the load will exceed any given level is simply twice the number of hours

from the time of maximum load (t = 0) until the time that load would have occurred.

This yields an equation of the same form as eq 2-28, except that the period of the

function is now twice as long, i.e.,

 

m ̇ / 

m ̇ d = [0.425 cos(2πt/17,520) + 0.575]. (2-29)

The form of the resulting load duration curve is shown as Figure 2. This result is

similar to the shape of empirical load curves determined by other investigators, such

as Werner (1984), except that it over-predicts the number of hours at which high

loads occur. Since this will result in a conservative design, this is a suitable

approximation for design purposes. The equivalent full load utilization time is

another factor used to evaluate the load in district heating systems. The equivalent

full load utilization time is the amount of time at which the load would have to be

at its maximum value to result in the total heat supplied being equal to that supplied

over the actual yearly cycle. This time can be easily found by integrating eq 2-28 over

the yearly cycle

 

t tt u = (0.425 cos(2 /8760) + 0.575) d ∫

0

8760

π (2-30)

where t

u is the equivalent full load utilization time (hr). Carrying out this integration

yields a value of tu = 5037 hours. This value is near the upper end of the range of

measured values reported by Bøhm (1988) for a number of Danish district heating

systems.

A number of other technical parameters need to have values assigned to them for

us to proceed with this sample calculation. The following values selected are felt to

be well within the range of what could reasonably be expected in an actual design:

Aη = 0.90 (dimensionless)

Tm = 6.4°C

ki = 0.030 W/m °C

14
ks = 1.3 W/m °C

Hp = 1.0 m

∆xi = 0.050 m

ε = 5 × 10–5 m

a = 0.119 (dimensionless)

b = 0.152 (dimensionless)

c = –0.0568 (dimensionless)

Ts = 120°C

Tr = 60°C.

Application

For the application, we consider a main portion of a distribution system that

would serve a large number of consumers. We assume a maximum heating load,

including pipeline heat losses, of 25 MW. At the temperature difference specified

above, this would require a “design” or maximum flow of 

m ̇ d = 100 kg/s. We also

assume that the length of the pipeline is 1000 m. Note, however, that the length used

will not affect the diameter determined, as the length could be factored out of each

of the variable terms in the objective function eq 2-23, and the calculation would then

be done on a unit length basis. To arrive at a realistic total cost, which includes the

cost fixed with respect to d, the calculations here are for the system length specified

above. We also assume that only one pump is associated with the system.

Solution

For the problem described above, we arrive at the following values for the

parameters in the objective function:

γ = 0.0231 (dimensionless)

A9 = $2.58 × 106/m

I

1 = $8.56 × 104

I

3 = 44.1 $ m5.095

A1 = $1060/pump

A3 = $2.18 × 105.

The calculation of the above parameters is straightforward with the exception of

I3. The integral in the I

3 parameter was evaluated numerically by a FORTRAN

program adapted from Ferziger (1981), which uses Romberg integration. The

program is included in Appendix B.

Before solving eq 2-24 to determine the optimum diameter, we first find an

approximate solution using eq 2-20, which neglects the heat losses. From eq 2-20 we

solve for the diameter directly, obtaining d = 0.216 m. Using this value of d as an

initial estimate, we can proceed to solve eq 2-24. We know that the solution to eq 2-

24, which includes heat losses, will be a smaller diameter than the solution to eq 2-

20, which does not include heat losses, since heat losses are an increasing function

of the diameter. Various “root finder” methods can be used to find the solution to

eq 2-24. Guided by the value obtained above, a simple trial-and-error method was

used here, which yielded a solution to three significant digits with several function

evaluations. The optimal diameter d was found to be 0.208 m. The total cost for this

design is Ct

 = $1.11 × 106. In the following section, this result will be compared to one

obtained using a common design rule of thumb.

Comparison with a design based on a rule of thumb

Ideally, an analysis similar to the one above would be used to size all major district

heating pipes. In reality, however, most systems are designed on the basis of rules

of thumb that have evolved from practice. Although such rules of thumb may prove

15
Table 1. Pressure drops and costs for discrete pipe sizes

under maximum flow conditions (pipe data from Marks

1978).

Nominal Inside diameter

pipe size schedule 40 ∆Pd Ct

(in.) (in.) (m) (Pa/m) ($ × 106)

— 8.187 0.208 340 1.111

8 7.981 0.203 384 1.112

10 10.020 0.255 120 1.178

12 11.938 0.303 50 1.305

adequate in some cases, they lack the flexibility to account for varying conditions,

most notably economic. Because these rules of thumb are based on designs proven

only to be functional, they cannot profess to yield least life cycle cost designs. To see

how the results of the above example would compare with a rule of thumb based

design, we consider a very common design rule of thumb used in Europe for systems

in this temperature range: that the pressure loss in the piping not exceed 100 Pa/m.

For this example, standard schedule 40 pipe sizes are used.

To apply the above rule of thumb, we simply calculate the pressure loss that

would result at maximum flow conditions using increasing pipe size until we find

a size that satisfies the rule. This calculation is done using eq 2-15 given earlier. The

results are shown in Table 1. We see from Table 1 that a 12-in. (300-mm) pipe would

be necessary to satisfy the rule of thumb. The pressure loss for the 10-in. (250-mm)

pipe exceeds the 100-Pa/m level by over 20% and therefore would probably be

considered unacceptable.

Now we need to determine what discrete pipe diameter would be recommended

by the procedure outlined in the previous section. The optimal nondiscrete diameter

was found to be 0.208 m or 8.187 in. We see from Table 1 that this lies between the

inside diameters of the 8- and 10-in. (200- and 250-mm) nominal pipe sizes.

To determine which to use, we simply calculate the cost of each alternative using

eq 2-19. These results are also included in Table 1. We see from these figures that the

total life cycle cost of the 8-in. pipe is about 6% less than the 10-in. pipe and thus the

8-in. pipe should be selected. We also note that the life cycle cost of the 8-in. pipe is

only 0.1% greater than that of an optimal 8.187-in. inside diameter pipe, if such a pipe

were available.

If we compare the cost of the 8-in. pipe, which our method recommends, to the

12-in. pipe required by the maximum pressure drop rule of thumb, we find that the

life cycle cost of the rule based design is 17% greater. This great saving in life cycle

cost is also accompanied by an even greater 30% reduction in capital costs (sum of

eq 2-16 and 2-17). As the financing of a new district heating system is often a barrier

to implementation, such large reductions in capital costs could make a system

feasible where it might not be otherwise.

We have arrived at an optimal pipe size that promises to save 17% in life cycle cost

over a rule of thumb based design. This result is consistent with the results of others

(Bøhm 1986, Koskelainen 1980) who have compared optimized designs with rule of

thumb based designs. In determining this pipe size, we have not considered any

constraints on the selection, other than it be a commercially available size. Of course,

in reality, other constraints exist. Before this method could be used to design an

entire system, the constraints that arise from interconnection of the pipes need to be

considered. Constraints also arise because of the consumer’s equipment and minimum temperature requirements. Other constraints are associated with the limitations of the piping system and the plant that supplies the heat. These constraints will

be considered in the following chapters.

16
CHAPTER 3: THE CONSUMER’S HEAT LOAD

In the design and subsequent operation of a district heating system, the characteristics of the load can be very significant. The load will not only dictate the

combination of supply temperature and mass flow rate necessary for its satisfaction,

but the heat exchanger equipment used at the consumer will also determine the

return temperature of the water. Lowering return temperature is desirable because

it results in larger temperature differences and thus lower mass flow rates, pumping

energy expenditure, and possibly smaller pipes. The importance of this issue is

evidenced by many district heating utilities in Europe having taken significant

actions to achieve large temperature differences. Thus, it is essential that our design

methodology for the distribution piping system account for the characteristics of the

consumer’s load and the constraints that result.

The primary type of heat load for most district heating systems is space heating.

In some cases industrial process loads can also be significant. In most cases where

buildings rely on a district heating system for space heat, they also use the system

to heat hot water. Here we will develop simple models for space heating loads only.

SIMPLE MODEL FOR THE CONSUMER’S SPACE HEATING EQUIPMENT

As noted above, in addition to the maximum magnitude of the load placed on the

district heating system by the consumer, several other characteristics of the load are

important. The way in which the load varies is of primary importance. This was

discussed in Chapter 1 and will be addressed in more detail later. The other major

way in which the load affects heat distribution systems is through the response of

the consumer’s heat exchanger to changes in supply temperature. To address this

issue, we need a model for the consumer’s heat exchanger equipment. We will

develop such a model in this section.

In district heating systems using hot water, the water-to-air heat exchangers of the

consumers can either be directly connected to the network or indirectly coupled by

a heat exchanger. Each type of connection has its advantages and limitations. For the

sake of simplicity, we will assume that the buildings are directly connected in this

work. To address indirect systems, it would be necessary to either develop alternate

models or attempt to modify the model for a direct system developed below.

The normal radiator common on many residential and light commercial hydronic

heating systems can be classified as a cross flow heat exchanger with one of the fluids

mixed (water) and the other fluid (air) unmixed, as described by Kays and London

(1964). Although the term “radiator” is commonly used for these heat exchangers,

they function via both convective and radiative heat transfer within the temperature

ranges normally encountered in practice. A schematic representation of this cross

flow heat exchanger is shown in Figure 3.

Because the water is considered to be ideally mixed, its temperature is assumed

to be uniform in the direction of air flow at any point along the heat exchanger. As

the water moves through the heat exchanger, it varies from the supply temperature

Ts

 at the water inlet to the return temperature Tr at the outlet. The incoming air

temperature Ta is assumed to be constant along the length of the heat exchanger.

However, the outgoing air temperature Tao will vary along the length of the radiator

owing to the decline in water temperature.

Although it would be possible to describe the performance of a radiator using

traditional approaches, such as those described by Kays and London (1964), simpler

equations have been proposed. These are based on experimental results for such

heat exchangers, an example being the equation given by Bøhm (1988)

17
q2/q0 = [(Tml)

1/(Tml)0]

n1 [(Tml)

2/ (Tml)1]

n2 (3-1)

where q = heat output from the radiator (W)

Tml = logarithmic mean temperature difference (°C)

n1, n2 = empirically determined coefficients (dimensionless).

and the subscripts denote the following operating conditions

0 = “design” condition for the radiators, usually the maximum load condition

at maximum supply temperature

1 = condition of actual supply temperature with the flow rate as determined

under the design condition

2 = any actual operating condition.

Equation 3-1 uses the logarithmic mean temperature difference Tml and experimentally determined constants to predict heat exchanger performance. The logarithmic mean temperature difference is defined as

 

T

TT TT

TT TT

T T

TT TT m1 sa ra

sa ra

s r

sa ra

= ( ) − − − ( )

( ) − − − ( ) = −

ln ln ln / [ ] ( ) − ( ) − . (3-2)

One problem that results from using the logarithmic mean temperature difference is that an explicit expression for either the supply temperature Ts, or the return

temperature Tr, cannot be obtained from the expression for the logarithmic mean

temperature. This limits the extent of closed form analysis and ultimately, when

calculations are required, it forces solution by iterative numerical methods. As a

solution to these problems, the use of the arithmetic mean as an approximation for

the logarithmic mean was proposed by Soumerai (1987). The arithmetic mean

temperature difference for this case is defined as

Tma= [(Ts – Ta) + (Tr – Ta)]/2 = (Ts + Tr – 2Ta)/2. (3-3)

The arithmetic mean temperature difference has the advantage that it can be used

to find a very simple explicit expression for either the supply temperature Ts or the

return temperature Tr given the value of the arithmetic mean temperature difference. The disadvantage of using the arithmetic mean temperature difference as an

approximation for the logarithmic mean temperature difference is the error induced

by this approximation. As Soumerai (1987) points out, within certain ranges of the

temperatures involved, the resultant errors are usually acceptable, given the other

uncertainties in heat transfer engineering. Soumerai (1987) recommends the use of

the arithmetic mean as an approximation for the logarithmic mean in cases where

the approach factor AF is equal to or greater than 0.5. The approach factor for this

type of heat exchanger is given by

Air Temperature

Supply

Temperature

Return

Temperature

Figure 3. Schematic of a hydronic heating system radiator.

18
AF = (Tr – Ta)/(Ts

– Ta) . (3-4)

In the case where the above criterion for the approach factor is met, the error of

approximation is always less than 4%. The arithmetic mean always overestimates

the logarithmic mean and thus any estimates of the heat transfer based on the

arithmetic mean will overestimate the actual heat transfer that will be achieved. This

could result in undersized heat exchangers, assuming that no other margin of safety

is included, which is of course seldom the case.

As an alternative approximation to the logarithmic mean temperature difference,

the use of the geometric mean temperature difference is proposed. The geometric

mean temperature difference for this type of heat exchanger is defined as

Tmg = (Ts – Ta)

1/2 (Tr

 – Ta)1/2 = ( 

Ta

2 + Ts

Tr – TaTs – TaTr)

1/2. (3-5)

The geometric mean temperature difference, like the arithmetic mean temperature difference, has the advantage that an explicit expression for the supply or return

temperature can be obtained from it. The geometric mean temperature difference,

however, is a much better approximation of the logarithmic mean temperature

difference than is the arithmetic mean temperature difference, as will be shown

below.

To simplify the analysis, we introduce the following expressions

∆Tsa = Ts

 – Ta (3-6)

∆Tra = Tr – Ta (3-7)

AF = ∆Tra/∆Tsa (3-8)

where ∆Tsa = greatest temperature difference between fluids (°C)

∆Tra = smallest temperature difference between fluids (°C)

AF = approach factor for the heat exchanger (dimensionless).

Two limiting cases of heat transfer set the range of values possible for the

approach factor AF. The first case is the case where no heat transfer takes place in the

heat exchanger. In this case the temperature of the water flowing through the

radiator will not decrease, and thus the supply and return water temperatures will

be equal and the approach factor becomes unity. The other limiting case occurs when

the maximum amount of heat transfer occurs in the heat exchanger, in which case

the return temperature equals the air temperature and approach factor becomes

zero. Thus, we have the following range of values for the approach factor AF

0 ≤ AF ≤ 1 . (3-9)

Now we can examine the errors that can result from each of the approximations

presented above over the entire range of possible approach factors. First, we define

the relative error of each of the approximations

εa = (Tma/Tml) – 1 (3-10)

εg = (Tmg/Tml) – 1 (3-11)

where εa is a relative approximation error for the arithmetic mean temperature

difference (dimensionless) and εg is a relative approximation error for the geometric

mean temperature difference (dimensionless). Then, by combining eq 3-2, 3-3, 3-6,

19
3-7, 3-8 and 3-10, we can arrive at the following expression for the relative error of

the arithmetic mean temperature difference εa in terms of the approach factor AF

εa = { [(AF + 1) ln(AF)]/[2(AF – 1)] } – 1. (3-12)

Similarly, we combine eq 3-2, 3-5, 3-6, 3-7, 3-8 and 3-11 to arrive at an expression for

the relative error of the geometric mean temperature difference εg in terms of the

approach factor AF

εg = {[(AF)

1/2 ln(AF)]/[AF – 1]} – 1. (3-13)

Now we can study the approximation errors for the arithmetic and geometric

mean temperature differences over the range of possible approach factors by

examining eq 3-12 and 3-13 respectively. It is immediately obvious that the error

from the arithmetic mean temperature difference approximation εa becomes infinite

as AF approaches zero. However, it is not clear what the error from the geometric

mean temperature difference approximation becomes as AF approaches zero. To

determine what value εg approaches as AF approaches zero, we use l’Hôpital’s rule.

It states that

 

limit ( )/ ( ) limit ( )/ ( ) fx gx f x g x

x x

( ) = ( ) ′ ′ → → λ λ (3-14)

where f(x) and g(x) are some functions of x that both approach either zero or infinity

when x approaches the value λ. To apply this to the error expression for the

geometric mean temperature difference, we let

x = AF

f(x) = f(AF) = ln(AF)

g(x) = g(AF) = (AF–1)/AF1/2 = AF1/2 – AF –1/2

εg = f(AF)/g(AF) – 1 .

Taking the first derivatives of f(AF) and g(AF), we have

f ′(AF) = 1/AF

g′ (AF) = 0.5 AF –1/2 + 0.5 AF –3/2 .

Now we can determine the value that εg approaches as AF → 0 from

limit limit limit g

AF AF

f AF g AF f AF g AF → → → ( ) + = ( ) = ( ) ′ ′ =

0 0

ε 1 ( )/ ( ) ( )/ ( )

AF 0

limit 

AF

AF AF →+ =

0

05 1 0 1 2 ( /[ . ( )]) /

Thus, we find that the error from approximating the logarithmic mean temperature

difference with the geometric mean temperature difference εg reaches –100% as the

approach factor AF goes to zero. Although this is a very high relative error, it is still

much better than that of the arithmetic mean temperature difference approximation,

which becomes infinite at the same condition. Of course, in reality this limiting case,

where heat transfer is at its maximum value and the approach factor becomes zero,

is never achieved. As we will now show, the errors attributable to using the

geometric and arithmetic mean temperature difference approximations for the

20
logarithmic mean temperature difference are always less than the values for an

approach factor of zero.

The other limiting value for the approach factor is unity. In this case no heat

transfer occurs. First, let’s examine what happens to the error from the arithmetic

mean approximation εa. We must again use l’Hôpital’s rule, proceeding as before

f(AF) = (AF – 1) ln(AF)

g(AF) = 2(AF –1)

f′ (AF) = 1 + ln(AF) + 1/AF

g′ (AF) = 2 .

In the limit as AF approaches unity, we have

 

limit limit limit g

AF AF AF

f AF g AF f AF g AF →→ → ( ) + = ( ) = ( ) ′ ′ =

11 1

ε 1 ( )/ ( ) ( )/ ( )

 

limit

AF

AF AF → ( ) + + =

1

( ln( ) / )/ 1 1 21 .

So, we see that the error induced by using the arithmetic mean temperature

difference as an approximation for the logarithmic mean temperature difference

approaches zero as AF approaches unity. Now let’s look at what happens to the error

for the geometric mean temperature difference as AF approaches unity. Again we

see that l’Hôpital’s rule is needed and we proceed as follows

f (AF) = AF1/2 ln(AF)

g (AF) = AF – 1

f′ (AF) = AF –1/2 [(ln(AF)/2) + 1]

g′ (AF) = 1 .

In the limit as AF approaches unity we have

 

limit limit limit g

AF AF AF

f AF g AF f AF g AF →→ → ( ) + = ( ) = ( ) ′ ′ =

11 1

ε 1 ( )/ ( ) ( )/ ( )

 

limit

AF

AF AF →

− ( ) [ ] + =

1

1 2 21 1 / (ln( ))/ ) .

Thus, we find that the error for the geometric mean temperature difference

approximation to the logarithmic mean temperature difference also approaches

zero as AF approaches unity. The errors resulting from using the arithmetic and

geometric mean temperature differences as approximations for the logarithmic

mean temperature difference are shown in Figure 4. Some numerical values for these

errors are also given in Table 2.

Several important observations can be made by studying Table 2. First, we note

that the error from approximating the logarithmic mean temperature difference

with the arithmetic mean temperature difference is always positive. Since the heat

transfer is proportional to the logarithmic mean temperature raised to some positive

power, using the arithmetic mean temperature difference as an approximation will

always over-predict the actual heat transfer. Also note that, as we have shown

analytically, the arithmetic mean temperature difference approaches infinity as the

AF goes to zero and approaches zero as AF goes to unity. For the geometric mean

temperature difference, the error resulting from using it as an approximation for the

21
300

200

100

0

Approach Factor

Approximation Error (%)

Arithmetic

Error

Geometric 

Error (abs.)

10–3 10– 2 10 – 1 100

Figure 4. Errors from approximating the logarithmic mean temperature difference with arithmetic and geometric mean temperature differences.

Table 2. Errors from approximating the logarithmic mean

temperature difference.

Approach Arithmetic Geometric Ratio of

factor mean error mean error errors

AF εa(%) εg(%) εg/εa

0.0 + ∞ –100.0 – 0.00

0.0001 361.0 –90.8 – 0.25

0.001 246.0 –78.1 – 0.32

0.01 135.0 –53.4 – 0.40

0.1 40.7 –19.1 – 0.47

0.2 20.7 –10.0 – 0.48

0.3 11.8 –5.79 – 0.49

0.4 6.90 –3.41 – 0.495

0.5 3.97 –1.97 – 0.497

0.6 2.17 –1.08 – 0.498

0.7 1.058 – 0.528 – 0.4992

0.8 0.415 – 0.207 – 0.4997

0.9 0.0925 – 0.0462 – 0.4999

1.0 0.0 0.0 → – 0.500...

logarithmic mean temperature difference is always negative. Thus, the predicted

heat transfer using this approximation would always be conservative; i.e., it would

under-predict the actual heat transfer. We also note from Table 2 that, as we had

shown analytically, the error resulting from the use of the geometric mean approximation approaches 100% in magnitude as AF goes to zero and approaches zero as

AF approaches unity.

The ratio of the error from using the arithmetic mean and geometric mean

approximations is also given in Table 2. Because the error from the arithmetic mean

approximation becomes infinite and the error from the geometric mean approximation approaches –100% as the approach factor AF goes to zero, their ratio approaches

zero at that point. Thus, the geometric mean approximation is infinitely better than

the arithmetic mean approximation at that point. Since neither approximation is

acceptable near that point, this observation is of little use. However, it is of interest

to note that the ratio of errors approaches 1/2 as AF approaches unity. Although this

22
is seemingly apparent from Table 2, we can prove this analytically by using

l’Hôpital’s rule. In this case it becomes necessary to take successive derivatives up

to the third derivatives in order to arrive at an expression that is not indeterminate.

This is an acceptable application of l’Hôpital’s rule, however. The analysis is

presented briefly below.

 

ε

ε

g

a

= − +

[ ] + − +

AF AF AF

AF AF AF

1 2 1

1 21

/ ln( )

( )ln( )/

Let

εg = f(AF) = AF1/2ln(AF) – AF + 1

εa = g(AF) = [(AF + 1) ln(AF)/2] – AF + 1

f ′(AF) = [AF1/2ln(AF)/2] + AF –1/2 – 1

g′ (AF) = ln(AF)/2 + 1/2AF – 1/2 .

In the limit as AF→ 1, we see that both f ′(AF) and g′(AF) approach zero; thus, we still

have an indeterminate expression. Applying l’Hôpital’s rule to that expression

f ′′(AF) = – [AF –3/2ln(AF)]/4

g′′(AF) = (AF –1– AF –2)/2 .

Again, we see that both f ′′(AF) and g′′(AF) approach zero as AF approaches unity

and we are left with another indeterminate expression. Once again we take derivatives so that we can apply l’Hôpital’s rule

f ′′′(AF) = [(3AF –5/2ln(AF))/8] – AF –5/2/4

g′′′(AF) = AF –3 – AF –2/2

 

limit limit ( )/ ( ) ( / )/( / ) / ε ε ( ) g a + = ( ) ′′′ ′′′ = − =−

→ →

f AF g AF

AF AF

14 12 12

1 1

And we now have our desired result.

The significance of this result is that we now know that the error from using the

geometric mean temperature difference as an approximation to the logarithmic

mean temperature difference is always 50% or less of the error that would result

from using the arithmetic mean temperature difference. For applications where the

use of the logarithmic mean temperature difference is undesirable, use of the geometric mean temperature difference will result in errors of less than 5% for values

of AF greater than 0.33. As Soumerai (1987) points out, given the other uncertainties

in heat exchanger design calculations, errors of this magnitude are certainly acceptable. Most heat exchanger designs will have approach factors greater than 0.33 and

thus our findings here should be applicable in the majority of cases.

Now that we have this approximation for the logarithmic mean temperature

difference, we can construct a simple model for the consumer’s heat exchanger using

it. Equation 3-1 will be used to construct our model. For our model two possible cases

exist, dependent on the values of the empirical parameters n1 and n2. Here, we will

only address the simpler case where n1 = n2. This is the result that occurs for “high

radiators” according to Bøhm (1988), in which case n1 = n2 = 1.3. In that case eq 3-1

becomes

 

qq T T

n

2 0 2 0

1

/ = 

( )( ) /   m1 m1  . (3-15)

23
Our model of the consumer’s heat exchanger should accept as input the supply

temperature to the radiators and the heat load q2. The values of the operating

parameters at the design “0” condition will also be needed. The output from the

model will be the return temperature. We would like this model to be as accurate as

possible while being in a simple form, thus allowing ease in both analytical and

numerical procedures involving it. As stated earlier, it is not possible to obtain an

explicit model if the log mean temperature difference is used. Thus, we will proceed

using our geometric mean approximation for the log mean temperature difference.

Making this substitution, we have

 

qq T T

n

2 0 2 0

1

/ / = ( )( ) 



 

 mg mg  . (3-16)

Note that the log mean temperature difference at the design 0 condition is also being

approximated with the geometric mean temperature difference. Since all of the

temperatures are known at the design condition, we could have evaluated the log

mean temperature difference and used the result here and still achieved an explicit

result. However, to ensure that no error occurs in the resultant model at the design

condition, we have used the geometric mean approximation. This will also reduce

the errors at the “off-design” (2) condition. The same procedure has been adopted

for the model using the arithmetic mean temperature difference.

To obtain our model for the return temperature as a function of the load and

supply temperature, we simply solve eq 3-16 for the return temperature at the 2

condition (actual load). The result is

 

( ) ( ) ( )( / ) /

T T TT T qq n

r a s a mg 2 2

1

0

2

2 0

2 1

=+ − 



 



 − . (3-17)

We can also obtain a model for the return temperature using the arithmetic mean

approximation to the log mean temperature difference. It is

 

( ) ( )( / ) ( ) / T T T qq T n r a ma s 2 02 0

1

2 = + 2 1 { } [ ] − . (3-18)

To evaluate the performance of our models that use approximations to the log mean

temperature difference, we need a model that uses the log mean temperature

difference. As noted earlier, this model will be implicit and thus will require solution

by an iterative numerical method of some type. The model can be arranged in several

forms for numerical solution, one being

 

( ) ( ) /exp ( / ) ( ) /( ) /

T T TT qq TT T n

r a s a s r m1 2 2 20

1

2 0

1

=+ − − 



 

















− . (3-19)

A number of iterative methods can be used to solve this implicit equation for the

return temperature (Tr)2. Most iterative methods are very sensitive to the quality of

the initial estimate. Here, we are fortunate to have the geometric mean temperature

difference approximation that can be used to obtain the initial estimate. Figure 5

below shows some of the results for the three models developed. In addition more

detail as well as numerical values are given in Table 3. It is clear from Figure 5 that

the model using the arithmetic mean temperature difference is unacceptable for

most values of the load ratio q/q0, while the model using the geometric mean

temperature difference is acceptable over the entire range of values given for q/q0.

The results for the model using the log mean temperature difference were obtained

24
by using the geometric mean approximation as an initial estimate, then proceeding

with a simple “iterative improvement” method where the previous estimate of (Tr)2

was substituted into the right-hand side of eq 3-19, which yielded the next estimate.

This procedure is repeated until successive estimates of (Tr)

2 varied by no more than

the prescribed tolerance.

For each of the consumer models, we can easily develop an expression to calculate

the flow rate for any given load condition by starting with the heat balance for the

radiator. We assume that the mass of the radiator is negligible so that with the

gradual temperature changes typical of these systems, conditions are very close to

steady state. Treating the radiator as a control volume, we have

q = 

m ̇ cp(Ts – Tr

) (3-20)

where cp is the specific heat of water at constant pressure (kJ/kg °C).

Thus, our mass flow rate relative to the mass flow rate at the design condition is

given by

 

m ̇ / m ̇ d = (q/qd)(Ts – Tr

)d/(Ts – Tr)2. (3-21)

Notice that we have used the d subscript to denote the design condition for the

piping system rather than the 0 subscript used to denote the design condition for the

consumer’s radiators. If both the network piping and the consumer’s radiators are

designed for the same maximum load condition, then it would not be necessary to

distinguish between these two conditions. However, in most cases this will not be

the case. For the piping network, little or no over-design is desirable in order to keep

costs at a minimum. In fact, diversity of demand between consumers will allow the

network to be designed for a total maximum demand of less than the sum of the

individual demands, as will be discussed later. The consumer’s radiators, on the

other hand, will always be somewhat oversized. In addition to the normal conservatism in design, quick recovery from night setback and other off periods also favors

significant over-design. Relative mass flow rates calculated using eq 3-21 for each of

the consumer models are given in Table 3. The results in Table 3 assume the same

design condition for the piping network and the consumer’s radiators. Some

examples with differing design conditions will be given later.

80

60

40

20

0.2 0.4 0.6 0.8 1.0

Heat Load Ratio

Return Temperature (°C)

Logarithmic Mean

Geometric Mean

Arithmetic Mean

0

Figure 5. Results for the return temperature models. Top family is for Ts =

80°C; middle family is for Ts = 90°C; bottom family is for Ts

 = 100°C.

25
Table 3. Return temperatures and flow rates calculated with the consumer models for n1 = n2 = 1.3.

Heat GMTD AMTD LMTD

Supply loss rtn. rtn. rtn. GMTD AMTD

temp. ratio temp. temp. temp. GMTD AMTD LMTD GMTD AMTD LMTD GMTD AMTD m/md m/md

(°C) (q/qo) (°C) (°C) (°C) AF AF AF error error m/md m/md m/md error error

100 1 63.75 60.00 62.75 0.55 0.50 0.53 –0.0160 0.0438 0.5369 0.5517 0.5000 –0.0277 0.0687

100 0.9 57.20 50.66 55.56 0.47 0.38 0.44 –0.0295 0.0883 0.4051 0.4206 0.3648 –0.0383 0.0994

100 0.8 51.04 41.07 48.73 0.39 0.26 0.36 –0.0473 0.1572 0.3121 0.3268 0.2715 –0.0471 0.1300

100 0.7 45.27 31.21 42.31 0.32 0.14 0.28 –0.0700 0.2625 0.2427 0.2558 0.2035 –0.0541 0.1614

100 0.6 39.94 21.01 36.39 0.25 0.01 0.20 –0.0974 0.4227 0.1887 0.1998 0.1519 –0.0590 0.1948

100 0.5 35.06 10.41 31.09 0.19 –0.12 0.14 –0.1277 0.6653 0.1451 0.1540 0.1116 –0.0611 0.2309

100 0.4 30.68 –0.70 26.57 0.13 –0.26 0.08 –0.1550 1.0263 0.1089 0.1154 0.0794 –0.0594 0.2708

100 0.3 26.86 –12.47 23.04 0.09 –0.41 0.04 –0.1657 1.5411 0.0780 0.0820 0.0533 –0.0522 0.3158

100 0.2 23.68 –25.21 20.81 0.05 –0.57 0.01 –0.1379 2.2113 0.0505 0.0524 0.0319 –0.0376 0.3675

100 0.1 21.27 –39.58 20.03 0.02 –0.74 0.00 –0.0618 2.9763 0.0250 0.0254 0.0143 –0.0157 0.4271

95 1 66.67 65.00 66.20 0.62 0.60 0.62 –0.0071 0.0181 0.6944 0.7059 0.6667 –0.0165 0.0400

95 0.9 59.68 55.66 58.62 0.53 0.48 0.51 –0.0181 0.0505 0.4948 0.5097 0.4575 –0.0301 0.0753

95 0.8 53.11 46.07 51.38 0.44 0.35 0.42 –0.0335 0.1034 0.3668 0.3819 0.3270 –0.0411 0.1086

95 0.7 46.96 36.21 44.55 0.36 0.22 0.33 –0.0540 0.1874 0.2775 0.2914 0.2381 –0.0500 0.1420

95 0.6 41.27 26.01 38.22 0.28 0.08 0.24 –0.0798 0.3195 0.2113 0.2233 0.1739 –0.0567 0.1770

95 0.5 36.07 15.41 32.49 0.21 –0.06 0.17 –0.1099 0.5258 0.1600 0.1697 0.1256 –0.0606 0.2147

95 0.4 31.40 4.30 27.55 0.15 –0.21 0.10 –0.1397 0.8438 0.1186 0.1258 0.0882 –0.0605 0.2563

95 0.3 27.32 –7.47 23.62 0.10 –0.37 0.05 –0.1569 1.3163 0.0841 0.0887 0.0586 –0.0547 0.3034

95 0.2 23.92 –20.21 21.03 0.05 –0.54 0.01 –0.1378 1.9610 0.0541 0.0563 0.0347 –0.0408 0.3579

95 0.1 21.35 –34.58 20.05 0.02 –0.73 0.00 –0.0651 2.7253 0.0267 0.0272 0.0154 –0.0177 0.4216

90 1 70.00 70.00 70.00 0.71 0.71 0.71 0.0000 0.0000 1.0000 1.0000 1.0000 0.0000 0.0000

90 0.9 62.52 60.66 62.00 0.61 0.58 0.60 –0.0084 0.0216 0.6429 0.6550 0.6135 –0.0188 0.0457

90 0.8 55.47 51.07 54.33 0.51 0.44 0.49 –0.0210 0.0600 0.4486 0.4634 0.4110 –0.0330 0.0837

90 0.7 48.88 41.21 47.06 0.41 0.30 0.39 –0.0388 0.1243 0.3260 0.3405 0.2869 –0.0444 0.1199

90 0.6 42.79 31.01 40.27 0.33 0.16 0.29 –0.0625 0.2300 0.2413 0.2542 0.2034 –0.0533 0.1570

90 0.5 37.21 20.41 34.09 0.25 0.01 0.20 –0.0917 0.4013 0.1789 0.1894 0.1437 –0.0592 0.1966

90 0.4 32.21 9.30 28.68 0.17 –0.15 0.12 –0.1231 0.6756 0.1305 0.1384 0.0991 –0.0611 0.2401

90 0.3 27.84 –2.47 24.30 0.11 –0.32 0.06 –0.1460 1.1017 0.0913 0.0965 0.0649 –0.0571 0.2895

90 0.2 24.20 –15.21 21.30 0.06 –0.50 0.02 –0.1363 1.7139 0.0582 0.0608 0.0380 –0.0441 0.3470

90 0.1 21.45 –29.58 20.07 0.02 –0.71 0.00 –0.0686 2.4741 0.0286 0.0292 0.0167 –0.0201 0.4152

85 1 73.85 75.00 74.17 0.83 0.85 0.83 0.0044 –0.0112 1.8471 1.7931 2.0000 0.0292 –0.0828

85 0.9 65.79 65.66 65.75 0.70 0.70 0.70 –0.0006 0.0014 0.9351 0.9370 0.9306 –0.0020 0.0048

85 0.8 58.20 56.07 57.61 0.59 0.55 0.58 –0.0102 0.0268 0.5843 0.5970 0.5531 –0.0218 0.0533

85 0.7 51.11 46.21 49.86 0.48 0.40 0.46 –0.0249 0.0734 0.3985 0.4131 0.3609 –0.0366 0.0943

85 0.6 44.54 36.01 42.59 0.38 0.25 0.35 –0.0458 0.1545 0.2829 0.2966 0.2449 –0.0482 0.1343

85 0.5 38.54 25.41 35.91 0.29 0.08 0.24 –0.0733 0.2924 0.2037 0.2152 0.1678 –0.0566 0.1762

85 0.4 33.15 14.30 29.99 0.20 –0.09 0.15 –0.1053 0.5231 0.1454 0.1543 0.1132 –0.0609 0.2219

85 0.3 28.45 2.53 25.11 0.13 –0.27 0.08 –0.1331 0.8992 0.1002 0.1061 0.0728 –0.0591 0.2737

85 0.2 24.53 –10.21 21.65 0.07 –0.46 0.03 –0.1331 1.4715 0.0631 0.0661 0.0420 –0.0476 0.3346

85 0.1 21.56 –24.58 20.11 0.02 –0.69 0.00 –0.0722 2.2228 0.0308 0.0315 0.0183 –0.0229 0.4078

80 1 78.33 80.00 78.49 0.97 1.00 0.97 0.0020 –0.0192 13.2420 12.0000 ∞ 0.0938 ∞

80 0.9 69.60 70.66 69.90 0.83 0.84 0.83 0.0043 –0.0108 1.7827 1.7315 1.9268 0.0287 –0.0809

80 0.8 61.38 61.07 61.29 0.69 0.68 0.69 –0.0015 0.0036 0.8553 0.8594 0.8454 –0.0048 0.0116

80 0.7 53.70 51.21 53.02 0.56 0.52 0.55 –0.0127 0.0343 0.5190 0.5323 0.4862 –0.0256 0.0631

80 0.6 46.58 41.01 45.22 0.44 0.35 0.42 –0.0302 0.0931 0.3450 0.3591 0.3078 –0.0409 0.1079

80 0.5 40.08 30.41 37.99 0.33 0.17 0.30 –0.0551 0.1996 0.2380 0.2505 0.2016 –0.0524 0.1529

80 0.4 34.25 19.30 31.52 0.24 –0.01 0.19 –0.0866 0.3876 0.1650 0.1748 0.1318 –0.0596 0.2012

80 0.3 29.15 7.53 26.07 0.15 –0.21 0.10 –0.1181 0.7112 0.1113 0.1180 0.0828 –0.0606 0.2559

80 0.2 24.90 –5.21 22.08 0.08 –0.42 0.03 –0.1277 1.2357 0.0691 0.0726 0.0469 –0.0512 0.3203

80 0.1 21.69 –19.58 20.16 0.03 –0.66 0.00 –0.0757 1.9714 0.0334 0.0343 0.0201 –0.0262 0.3991

Averages –0.0700 0.6902 0.0360 0.1899

Table Nomenclature: LMTD = log mean temperature difference (°C); GMTD = geometric mean temperature difference (°C);

AMTD = arithmetic mean temperature difference (°C); m/md = ratio of mass flow rate to design condition mass flow rate.

Assumed conditions: at “0” condition Ts = 90°C, Tr = 70°C; for all calculations Ta= 20°C, n1= n2 = 1.3.

Tolerance for iterative calculation of Tr with LMTD model was < 0.01°C.

26
As expected, the results for the normalized mass flow rates in Table 3, which use

the geometric mean approximation for the log mean temperature difference, are

much better than those obtained using the arithmetic mean approximation. Clearly,

the errors induced by the arithmetic mean approximation are unacceptable for any

load–supply temperature condition that deviated significantly from the design

condition for the radiators. Over the range of supply temperatures and heat loads

given in the Table 3, the average error in the return temperature obtained using the

geometric mean approximation is only 7% compared to 69% for the results obtained

using the arithmetic mean approximation. Also, note that the errors in approximation of the mass flow rate average only 3.6% for the model using the geometric mean

approximation, while the average error for the model using the arithmetic mean

approximation is about 19%.

In addition, it should be noted that the model based on the arithmetic mean

approximation in a number of instances at lower loads results in physically impossible return temperatures, i.e., ones lower than the room air temperature of 20°C.

Note that at lower loads (lower q/q0 values), the errors of approximation tend to be

larger. This is predicted by our error analysis carried out earlier, since the approach

factor is lower in these cases. Also note that, like our basic geometric mean

approximation, our model based on it is conservative and under-predicts the heat

transfer on the average. An exception sometimes occurs at supply temperatures

lower than the design condition of 90°C. Since our model for the radiator predicts

the return temperature based on the return temperature at the design 0 condition,

the error in approximation of the return temperature at the design condition has an

effect on the error in our model at conditions other than the design condition. At

lower supply temperatures and high loads, for example Ts

 = 85°C and q/q0 = 1.0, the

approach factor is 0.83 for the geometric mean approximation and is thus higher

than the approach factor of 0.71 encountered at the design condition. Because the

error in the geometric mean approximation decreases with increasing approach

factor, our model for the prediction of the return temperature actually underpredicts it slightly at that point. This under-prediction is not, however, a cause for

concern, since it is so slight and in addition it exists at a load–temperature condition

that would not normally be encountered because it would require mass flow rates

greater than the design condition.

DESIGN OF A SINGLE PIPE SEGMENT WITH A CONSUMER MODEL

In Chapter 2 we developed a methodology to determine the optimal pipe

diameter for a single pipe segment. In the example given, it was assumed that both

the supply temperature and return temperature were constant over the entire yearly

cycle. This is of course not the case, and now that we have a simple model for the

consumer’s space heating substation, we can examine what the effect is of coupling

this model with our design methodology.

First, let’s consider what the effect is of assuming a constant supply temperature,

as we had done earlier, but rather than assuming a constant return temperature as

well, let this be determined by our consumer model. The varying return temperature

will affect the heat losses by altering the A1 parameter (eq 2-5) to some degree; this

will be addressed later. The primary effect, however, will be on the mass flow rate.

The heat load will be assumed to vary sinusoidally, as before, except now the

variation in the mass flow rate will not be sinusoidal itself, but will be determined

by the load and the return temperature from the consumer model. The relationship

between the mass flow rate, load, design supply–return temperatures and the actual

supply–return temperatures was given by eq 3-21. If we substitute our expression

for the return temperature as determined using the geometric mean approximation

27
(eq 3-17) into eq 3-21, we have

 

m ̇ / m ̇ d = (q/qd) (Ts – Tr

)d/{Ts – 

 

T T T T qq n

a s a mg 0 0 − − 



 



 − ( ) ( )(/ ) / 1 2 2 1 } . (3-22)

Notice that we have used both q/qd and q/q0 in this expression. In both cases q

represents the actual load on the system. The quantity q0 represents the maximum

load for which the consumers’ radiators were designed, while the quantity qd

represents the maximum load for which the piping network was designed. These

two “design” loads will in most cases not be equal, as discussed earlier. They will,

however, be related by some “over-design factor,” which will be a constant

(q/qd) = A13 (q/q0) (3-23)

where A13 is the over-design factor for the consumer’s radiators (dimensionless).

In the example of Chapter 2, we assumed a sinusoidal form for the variation of the

load over the yearly cycle. In general terms this can be written as

q/qd = A14 + [A15 cos(2πt/8760)] (3-24)

where A14 is the midpoint of the load curve (dimensionless) and A15 is the amplitude

of the load curve (dimensionless). The midpoint of the load curve A14 is simply the

average of the maximum and minimum loads. The amplitude of the load curve A15

is the maximum load minus the minimum load divided by two.

Now if we combine eq 3-22–3-24, we have the following expression for the

normalized mass flow rate over the yearly cycle

 

 ̇ /  ̇ ( ) (/ )

( ) ( ) ( / )/

/ m m

TT A A t

TT TT T A A t A

n d s r d

s a s a mg

cos

cos

= − + ( )

−− − + ( ) [ ] 















−

14 15

1

0

2

14 15 13

2

2 8760

2 8760 1

π

π

 .

(3-25)

If we select the same values for A14 and A15 as we used in the Chapter 2 example

(0.575 and 0.425 respectively), we can compare the resulting mass flow rate function

of eq 3-25 to the mass flow rate function without the consumer model (eq 2-28). We

Figure 6. Mass flow rate function with and without consumer

models, Tr

 = 60°C.

1.0

0.8

0.6

0.4

0.2

0 2,000 4,000 6,000 8,000

Time (hr)

Normalized Mass Flow Rate 

m/m

w/o Consumer Model

d

m/m

w/Consumer Model

d

28
have done so in Figure 6. Notice, for the case where the consumer model has been

included, that the maximum value of the normalized mass flow rate is less than

unity. This results from an inconsistency in the design conditions chosen for the

original problem of Chapter 2 and the consumer’s radiator model. Our model for the

consumer’s radiators assumes design temperatures of 90 and 70°C for the supply

and return respectively. With the supply temperature of 120°C, the radiator model

will predict a return temperature of 55°C, rather than the value of 60°C we assumed

earlier. In the case where the return temperature is assumed rather than determined

by a consumer model, the choice of 60°C is entirely appropriate, as are the normal

design temperatures of 90 and 70°C for supply and return radiator temperatures

respectively. Normal design practice would be to assign these temperatures independently. Thus, the lower flow rate that results when we include the effect of the

consumer model illustrates one of the inaccuracies encountered when normal

design practice is followed. If we make the design return temperature for the piping

network equal to that which results from the consumer’s radiator model, i.e., 55°C,

it results in the normalized mass flow rates shown in Figure 7. From Figure 7 it is

clear that the effect of including the consumer model on the mass flow rate is still

significant once the load condition drops slightly from its maximum value. Averaged over all load conditions, the normalized mass flow rate is 20% less for the case

that includes the consumer model.

To determine what effect this change in normalized mass flow rate will have on

our optimal pipe diameter determined in the example of Chapter 2, we need to

recompute I3, which is the only parameter affected by the changes in the mass flow

rate. We have modified the computer program used in Chapter 2 to calculate I3 by

including eq 3-25 in place of the original normalized mass flow rate as given by eq

2-29. The modified program is included in Appendix B as Program I2-C-GMT. All

constants were assumed to have the same values as before, with A13 taken as unity.

This results in the value of the I3 parameter decreasing by 15% from 44.1 to 37.5 ($

m5.095) because of the effect of the consumer model.

Now we need to determine what other parameters in the solution of the Chapter

2 example would be affected by our consumer model. The only additional effect will

be on the I

1 parameter. Since this parameter arises out of heat loss considerations, the

varying return temperature caused by the consumer model will change it somewhat. In the example of Chapter 2, since the supply and return temperature were

1.0

0.8

0.6

0.4

0.2

0 2,000 4,000 6,000 8,000

Time (hr)

Normalized Mass Flow Rate 

m/m

w/o Consumer Model

d

m/m

w/Consumer Model

d

Figure 7. Mass flow rate function with and without consumer

models, Tr = 55°C.

29
both constant over the yearly cycle, the integration needed to find I1 (see eq 2-5) is

unnecessary. With the varying return temperature produced by our consumer

model, however, we will need to carry out this integration. Assuming that the cost

of heat Ch is constant, our new equation for I1 becomes

 

I A

T T

T A

T

T T

AA t A t

n

1 16

0

2

14 15 13

2 1

0

8760

2 2 = 2 8760 ( ) + − 









 + − ( + ] 



 





















 ∫ s a m t

mg

s a

d ( )

( )

cos( / ) /

/

π

(3-26)

where A16 = PVFh L Ch 4πki ($/[°C hr]).

Because the cosine function is raised to a non-integer power, it is not possible to

carry out the integration in eq 3-26 analytically. Once again we have used the

Romberg method of numerical integration to evaluate the integral. The calculation

of I

1 was done using the FORTRAN program I1EQ3-26, which is included in

Appendix B. Using the parameter values assumed earlier in this section, we obtain

I

1 = $7.33 × 104. Thus, we find that including the consumer model reduces the value

of the I

1 parameter by 14.4% from 8.56 × 104 to 7.33 × 104.

Now that we have new values for the parameters that are affected by the

consumer model, we can recompute the optimal diameter for the sample application

given in Chapter 2. We proceed as before, i.e., before solving eq 2-24 to determine the

optimum diameter, we first find an approximate solution using eq 2-20, which

neglects the heat losses. From eq 2-20 we solve for the diameter directly, obtaining

d = 0.210 m. Using this value of d as an initial estimate, we can proceed to solve eq

2-24. We know that the solution to eq 2-24, which includes heat losses, will be a

smaller diameter than the solution to eq 2-20, which does not include heat losses,

since heat losses are an increasing function of the diameter. Guided by the value

obtained above, a simple trial-and-error method was once again used here. This

method yielded a solution to three significant digits with only four function

evaluations. The optimal diameter d was found to be 0.203 m. The total cost for this

design is found to be Ct

 = $1.064 × 106 using eq 2-19. By coincidence, the optimal

diameter we have found also is one of the standard discrete diameter pipes

available; thus, it is not necessary for us to compute total costs for other discrete

diameters as before.

The addition of the consumer model has changed the optimal diameter from

0.208 to 0.203 m, a decrease of only 2.5%. The optimal discrete diameter remains

unchanged. While in this particular case, the inclusion of the consumer model had

no net effect on the choice of optimal discrete diameter, this obviously will not

always be the case.

The cost predicted for any pipe diameter is also changed slightly by the addition

of the consumer. The total cost with and without the consumer model is given in

Table 4. Pressure drops and costs for discrete pipe sizes under maximum flow

conditions with and without the consumer model (pipe data from Marks 1978).

Nominal Inside diameter Ct , w/o Ct, with

pipe size schedule 40 ∆Pd consumer model consumer model

(in.) (in.) (m) (Pa/m) ($ × 106) ($ × 106)

— 8.187 0.208 340 1.111 1.065

8 7.981 0.203 384 1.112 1.064

10 10.020 0.255 120 1.178 1.140

12 11.938 0.303 50 1.305 1.267

30
Table 4 for each optimal diameter and the discrete diameters found when using the

rule of thumb based design method. Pressure drops at maximum flow conditions

are also given in Table 4. Note that these are unchanged from those in Table 1, since

the maximum flow condition remains the same. Thus, the rule of thumb based

design would remain the same and a 12-in. nominal diameter pipe would be

required. The cost saving of the optimal discrete design increases slightly once the

consumer model is added. Now the rule of thumb based design is 19% more costly

than the optimal discrete design. Also, note that the total life cycle costs are reduced

in all cases when the consumer model is added. Since it is important to have accurate

cost predictions when comparing district heating to alternatives, these seemingly

minor changes in total life cycle cost can be significant. For instance, the total life

cycle cost of our optimal discrete diameter design decreases 4% with the addition of

the consumer model. This is a very significant cost reduction. In our example 1-kmlong pipe segment with a design capacity of 25 MW, this refinement in predicted life

cycle cost amounts to $48,000. Note that since our optimal discrete diameter is

unchanged by the addition of the consumer model, the capital cost of this design is

unchanged as well. Thus, the optimal discrete design still represents a 30% reduction in capital costs from the rule of thumb based design.

HEAT CONSUMER CONSTRAINTS

Before leaving the topic of the consumers, let’s consider the constraints that they

place on the design. The consumers of heat place two very basic requirements on the

heat supply system:

1. That the delivered temperature of the heat be high enough to meet their

requirements.

2. That their heat demand be met at all times.

The first requirement will simply result in the following inequality constraint

Ts,i ≥ Tsmin,i (3-27)

where Ts,i is the supply temperature at the heat consumer i (°C) and Tsmin,i is the

minimum supply temperature required by heat consumer i (°C).

Satisfaction of the second requirement will result in an equality constraint that

must be obeyed at each heat consumer. This constraint will be based on the model

developed in the previous section. The load placed on the system by the consumer

will be known, expressed as a fraction of the load under the design condition, i.e.,

q/q0. The supply temperature will also be known. The model for the consumer heat

exchanger then becomes our constraint on the return temperature. Equation 3-17 is

modified slightly by removing the 2 subscript, the unsubscripted values now

representing the actual operating condition

 

T T T T T qq n

r a s a mg 0 = − 



 



 − ( ) ( )(/ ) / 1

0

2 2 1 . (3-28)

There is also an additional equality constraint on the mass flow rate that results

from eq 3-21

 

m ̇ / m ̇ d = (q/qd)(Ts

 – Tr )d/(Ts – Tr). (3-29)

In the next chapter, we will examine how these and other constraints interact

when multiple consumer designs are considered.

31
CHAPTER 4: CONSTRAINTS ON SYSTEMS WITH

MULTIPLE CONSUMERS AND PIPES

All district heating systems, with the exception of pure transmission systems, will

have multiple consumers and pipes. If each pipe were independent of the others, it

would be possible to apply the procedure developed earlier to each pipe and create

a complete design in that way. Of course each pipe segment does not operate

independently of the others and the system can thus not be designed completely in

that way for all but the most trivial cases. Many constraints are imposed on the

design by the physical process involved in the network, the consumer’s requirements and physical limitations of the piping. As we shall see, many of these

constraints will be inactive at the optimum system design and thus they can be

relaxed. Our task is then to formulate these constraints into mathematical expressions, and then identify those that must be active and use this information to develop

a solution methodology. All of this must be done with the minimum amount of

computational effort so as not to render the method intractable for large networks,

which often have hundreds or thousands of piping segments.

SYSTEM CONSTRAINTS

Constraints on the design of a heat distribution system originate from limitations

imposed in several distinct areas. Before we begin to formulate constraints into a

form suitable for inclusion in our problem, let’s consider where and why these

constraints arise. The source of constraints can be grouped into three basic categories:

1. Physical limitations of the piping systems.

2. Fluid dynamic and thermodynamic considerations for the network, consumers, and heat source.

3. Requirements dictated by the consumer’s equipment or processes.

In some instances considerations from each of these categories are coupled

together into a single constraint or set of constraints. Thus, as we formulate the

constraints below, we will address considerations from each of the categories above

and their interaction.

DIFFERENTIAL PRESSURE CONSTRAINTS

A very important set of constraints on the system arises from requirements for the

pressure difference between supply and return. At the consumer this differential

pressure must maintain a minimum level to ensure adequate flow through the

consumer’s heat exchanger. This pressure differential is consumed in both the heat

exchangers and control valves. In the heat exchanger, the pressure losses are caused

by fluid dynamic friction. In the control valve, the pressure losses are introduced by

a throttling process used to control the flow rate through the heat exchanger and

thus control its output. In the supply piping between the heat source and the

consumer, pressure losses occur due to friction. Similarly, in the return line from the

consumer back to the heat source, pressure losses also occur. There is then, in effect,

a requirement at each point in the piping network for a given differential between

supply and return pressure necessary to overcome downstream losses, including

those in the return system. Ultimately, at the heating plant pumps must be used to

provide the total differential pressure needed downstream of that point. In theory

it’s possible that pumps can be placed anywhere in the system or even dispersed

32
throughout. In practice this is not done very frequently, owing to the practical

considerations of monitoring, controlling and maintaining the pumps as well as

availability of power for them. Here, we will assume that all pumps are located at

one central heating plant, although some very interesting optimal design questions

arise if this limitation is removed, as we will see later.

We can write the constraint that arises from all of these differential pressure requirements easily by summing the pressure drops and increases around the system.

Since the entire district heating system, consisting of the heating plant, the piping

system and the consumer, forms an essentially closed loop, the pressure losses and

increases around this loop must sum to zero. Thus, we have the following result

 

∆ ∆∆ ∆ ∆ P PP P P

pipes hp = + ∑( ) s r cv + + he (4-2)

where ∆Php = pressure increase across the pump (N/m2)

∆Ps = pressure drop in the supply piping (N/m2)

∆Pr = pressure drop in the return piping (N/m2)

∆Pcv = pressure drop in the consumer control valves (N/m2)

∆Phe = pressure drop in the consumer heat exchangers (N/m2).

Each consumer will have at least one segment of the piping system that is not

shared with any other consumers. In addition, all consumers will have their own

control valve and heat exchanger. Thus, we will have one of these equations for each

consumer, each one representing a constraint on the design. Therefore, the summation

in eq 4-2 must be conducted over only the pipes that serve the consumer in question.

The pressure losses given above will vary with the flow rate in the system. In

many cases, flow rates in district heating networks are modulated over the course

of the year as a means of meeting varying loads. Flow can be modulated either by

using variable speed pumps or using what is called a “shunt” at the heating plant.

The shunt simply diverts a fraction of the flow from the pump back to its inlet. The

pressure increase across the pump is reduced as is the flow rate into the network.

Regardless of how it is done, if flow modulation is used, we must ensure that the

constraint of eq 4-2 is not only satisfied for each consumer, but in addition we must

also determine that this will be the case for all load (i.e., flow) conditions encountered. However, we will show later, after some other necessary constraints have

been introduced, that satisfying this set of constraints for only one load condition

will be sufficient, if we assume that all consumers have loads that vary in the same

manner. This is a reasonable assumption as long as the primary loads are space and

hot water heating, as is the case for most hot water based systems. Steam systems,

which are not addressed by this work, often have much larger fractions of industrial

and absorption air conditioning loads; thus, this assumption might not be reasonable for them.

In addition to the above equality constraint (eq 4-2), we have both equality and

inequality constraints on each of the quantities appearing within that constraint. At

this point we will formulate each of these additional constraints.

At the heating plant, the pressure increase by the pump must be related to the

pumping power attributable to the consumer in question. This results in the

following expression

∆Php= PPf,i ρ/ 

m ̇ i (4-3)

where the i subscript is the consumer index. The combined pressure loss of the

supply and return piping is simply given by

 

∆∆∆ P P P a A m Ld b c c bc

s r s&r j j j j += = ( ) + + − ++ ε π 4 2

6

2 5 /  ̇ , ( ) (4-4)

33
where the j subscript is the pipe segment index and ∆Ps&r is the combined pressure

loss of supply and return (N/m2).

The control valve will have varying amounts of pressure drop across it, depending on the consumer’s load. The minimum pressure drop for any given flow

rate condition will occur when the control valve is completely open and the flow rate

through the consumer’s heat exchanger is at its maximum value—what we have

called the design condition. Thus, we have the following simple constraint for each

control valve in the system

∆Pcv,i ≥ ∆Pcvm,i (4-5)

where ∆Pcvm is the minimum pressure drop in the control valve (N/m2).

And finally, for the heat exchanger the pressure drop will be related to the flow

rate in a manner very similar to that for the pipes as given by eq 4-4 above. First, let’s

assume the following simple form

 

∆P Am he,i he,i i =  ̇ β (4-6)

where Ahe relates the fluid properties and physical properties of the heat exchanger

to the pressure drop and flow rate (kg1–β/m s2–β) and β is an exponent yielding the

appropriate mass flow rate dependency for the heat exchanger (dimensionless).

In most cases Ahe and β would probably be empirically determined coefficients

and would depend on the type of heat exchanger and its specific design. The basic

form of Ahe would probably be dictated by the heat exchanger geometry. For

example, if a straight section of pipe formed the hydraulic passageway for the heat

exchanger, the form of Ahe based on our previous analysis for pressure loss in pipes

would be

 

A a Ld cb c b c

he,sp he he he he he = + −− − ++ ( / )( / ) ( ) 2 4 2 1 5 π ε ρμ (4-7)

where the sp subscript denotes straight pipe heat exchanger and the he subscripts

denote conditions within the heat exchanger or physical parameters of the heat

exchanger.

Equation 4-7 assumes the same form of approximation for the friction factor as

was derived for flow in the district heating pipes earlier in Chapter 2. Using this

approximation for the friction factor also determines β from eq 4-6 to be

βsp = 2 + c. (4-8)

We can now substitute the results for the pressure losses around the district

heating system loop (eq 4-3, 4-4 and 4-6) into our original pressure loss constraint (eq

4-2) to obtain the following constraint

 

PP m a A m L d P A m

j

b c c bc c

f,i i j jj j j cv,i he,i i

2 ρ επ /  ̇ (/)  ̇  ̇ , ( ) = ∑ ( ) + + + + − ++ + 4 2

6

2 5 ∆ (4-9)

Again, in eq 4-9 the summation over the j pipes only includes those pipes that serve

consumer i.

MAXIMUM ABSOLUTE PRESSURE CONSTRAINTS

Several constraints on the absolute pressure of the water within the system must

be considered. First, we consider the upper limit on pressure that results from the

absolute pressure limits of the piping. This limit will be established by the prevailing

34
piping code. In the case where all points in the distribution system are at or above

the level of the heating plant, the maximum absolute pressure will occur in the

supply pipe at the heating plant. In the general case, however, this will not always

be true and it will be necessary to determine that this constraint is not violated at any

point within the system. However, several heuristics will allow us to forgo computation of the absolute pressure level at many of the points. At any point in the supply

side of the system the absolute pressure is given by

 

P P P gz

j s = − hp,s ∑ ∆ ρ s,j s − (4-10)

where Ps = absolute pressure in supply pipe at point in question (N/m2)

Php,s = absolute pressure in supply pipe at heating plant (N/m2)

∆Ps,j = pressure loss in supply pipe j (N/m2)

z = elevation at point in question relative to heating plant (m).

Again, we have assumed that no intermediate pumping is employed and that the

summation over j includes only those pipe segments between the heating plant and

the point in question along the supply line. If Pmax is the maximum absolute pressure

for the piping system being used (N/m2), our constraint arising from it is then

 

P P P gz

j max ≥ − hp,s ∑ ∆ ρ s,j s − (4-11)

where Pmax is the maximum absolute pressure for the piping system being used

(N/m2).

We can easily eliminate the need to verify that the upper limit on absolute

pressure is not exceeded for the return side of the system. First, we assume that the

supply and return line are at the same elevation at any given point, certainly a

reasonable assumption. Since there will always be a finite pressure drop across the

consumer’s heat exchanger and control valve, the absolute pressure in the return

line will always be less than that in the supply line at any point, with the difference

being the smallest at the consumer. Thus, we need not verify that the maximum

absolute pressure constraint (eq 4-11) is satisfied for the return system.

We can also easily show that eq 4-11 only needs to be satisfied at certain points

along the supply line. For pipelines that are laid at a constant slope between junction

points, the hydrostatic component of the pressure gradient along the pipe will be

constant as well. This gradient is given by the following equation

(dP/dx)h = –ρg(∂z/∂x) (4-12)

where (dP/dx)

h = hydrostatic pressure gradient (N/m3)

(∂z/∂x) = partial derivative of the elevation of the pipe with respect to

its position (dimensionless)

x = position along the pipe in the direction of flow with x = 0

being defined as the inlet end to the pipe segment in question

(m).

Using our approximation for the friction factor previously determined (eq 2-12),

we can find the pressure gradient attributable to frictional losses in the flowing fluid

from the following equation

 

d d d Px a md cb c c bc / ( / )( / )  ̇ ( ) ( ) = + − − + − ++ 2 4 2 125 π ερ μ (4-13)

35
where (dP/dx)d is the hydrodynamic pressure gradient (N/m3).

The pressure at any point along a segment of the piping system is simply the

pressure at the inlet to the pipe segment plus the sum of the hydrostatic and

hydrodynamic gradients multiplied by the distance

Px = PI + [(dP/dx)

h + dP/dx)d]x (4-14)

where Px is the pressure at point x (N/m2) and PI is the pressure at the inlet to the

pipe segment (x = 0) (N/m2).

Since (dP/dx)

h and (dP/dx)d are both independent of x, we only need to know

if their sum is positive or negative to determine if the pressure will be higher or lower

than the inlet pressure at the outlet of the pipe segment. We can also easily show by

using monotonicity analysis (Papalambros and Wilde 1988) that the maximum

pressure must occur at either the inlet or the outlet of the pipe section and cannot

occur at an intermediate point. To do so we convert the maximization problem to its

equivalent minimization problem

min. Px = –{PI

 + [(dP/dx)h + (dP/dx)d]x} (4-15)

subject to

g1 = –x ≤ 0 g1(x –) (4-16)

g2 = –x – L ≤ 0 g2(x+) . (4-17)

The conventions for constraints and the labeling of monotonicity are from

Papalambros and Wilde (1988). If the sum of the gradients is positive, the objective

function is monotonically decreasing in x and must therefore be bounded from

above. Constraint g2 is the only constraint that bounds the objective from above, so

it must be critical and thus x = L. If the opposite is true, the sum of the gradients is

negative, the objective will be monotonically increasing in x and must therefore be

bounded below. Constraint g1 is the only constraint that bounds the objective from

below, so it must be critical and x = 0. If the sum of the gradients is zero, the pressure

will be the same at all points along the pipe segment. Thus, we have shown that the

maximum pressure must always be at one end of the pipe segment and it will only

be necessary for us to ensure that our absolute pressure constraint (eq 4-11) is

satisfied at these points. Remember that in arriving at this result, we assumed that

the pipe segment had a constant slope between end points. If in reality this is not the

case, the pipe segment in question can be broken up into two or more equivalent pipe

segments for applying this constraint.

The number of points at which the maximum absolute pressure constraint must

be checked for satisfaction may possibly be reduced even further if we proceed as

follows.

1. Starting at the heating plant, we proceed along the supply line checking only

the “junction” points as discussed earlier.

2. For any point that is at the same elevation or higher than the upstream point

previously identified as having the maximum pressure, we need not compute the

pressure.

3. When a point is identified that does not meet the above criteria, we proceed by

first computing the sum of the hydrostatic and hydrodynamic gradients. If this

quantity is negative, we do not need to compute the pressure.

4. If the sum of the gradients is positive, we will need to compute the pressure at

this point. To compute the pressure, we first find the elevation difference and

resulting hydrostatic pressure difference between the point in question and the

36
previously identified point of highest pressure. We then must calculate the pressure

losses attributable to the hydrodynamic gradients in each of the pipe segments between the previously identified point of highest pressure and the point in question.

5. If we maintain a running total of pressure losses computed as we complete the

above step, we can stop the calculation procedure as soon as this total pressure loss

exceeds the hydrostatic pressure difference computed above. Otherwise, we must

continue the calculation, in which case we will have identified a new maximum

pressure location.

6. Steps 2 through 5 are repeated until we reach the end of the piping network. At

branching points we will need to proceed out each branch following the procedure

as outlined.

MINIMUM ABSOLUTE PRESSURE CONSTRAINTS

Minimum allowable pressure constraints arise from three distinct considerations.

1. Net Positive Suction Head (NPSH) requirements of the pump.

2. Minimum pressure over atmospheric necessary to preclude the infusion of air

into the system.

3. Pressure necessary to prevent flashing of the liquid.

The constraint resulting from minimum NPSH requirements necessary to preclude pump cavitation only needs to be satisfied at the inlet to the pump. This

pressure requirement will be a function of the saturation pressure and hence the

temperature of the liquid at that point. The NPSH requirement is usually specified

by the manufacturer of the pump. Thus, this constraint is simply

Php,r ≥ PNPSH (4-18)

where Php,r is the pressure in the return line at the inlet to pump (N/m2) and PNPSH

is the minimum allowable pressure at the pump inlet from NPSH requirements

(N/m2).

The amount of pressure over atmospheric necessary to prevent infusion of air into

the system will be another area where engineering judgment will be required. This

will be an issue primarily in portions of the system that are operating at temperatures below 100°C, since the saturation pressure constraint (eq 4-20) will dominate

it at higher temperatures, given equal safety margins. If, as we assumed earlier, no

intermediate pumping is being used, then the minimum pressure level will be at the

inlet to the pump for a system that is at or below the level of the heating plant at all

points. For other systems, we must check for dominance of this constraint or the

saturation pressure constraint derived below, and then constraint satisfaction must

be verified at all points within the system. At the heating plant, this constraint can

be written as

Php,r ≥ Pa + Pasa (4-19)

where Pa is atmospheric pressure (≈ 105 N/m2) and Pasa is the minimum safety

margin above atmospheric pressure (N/m2).

The second constraint on minimum allowable absolute pressure results from the

requirement that the fluid must be maintained above its saturation pressure some

finite amount to preclude flashing to the vapor phase. The amount of excess pressure

above the saturation pressure of the fluid is a matter of engineering judgment.

Because localized areas of pressure lower than the “bulk” pressure of the fluid may

occur because of hydrodynamic effects, a safety margin above the saturation

37
pressure is prudent. The resulting constraint is

Px ≥ Px,sat + Psaf (4-20)

where Px,sat is the saturation pressure of the liquid at point x within the pipe segment

(N/m2) and Psaf is the minimum allowable safety margin on saturation pressure

requirements (N/m2).

The saturation pressure is a function of the fluid temperature, which will vary

between supply and return portions of the system, as well as within each portion.

Thus, it will be necessary to verify the satisfaction of this constraint at all points

within the system. Again, some simple rules will allow us to forgo the calculation at

many points, as with the maximum absolute pressure constraint described earlier.

As noted earlier, in some cases when the temperature is below 100°C, the air infusion

constraint above (eq 4-19) will dominate. The concept of constraint dominance is

illustrated later in Chapter 5.

The pressure level at any point in the supply side can be calculated with eq 4-10.

For the return side, the absolute pressure is given by

 

P P P gz P P P

j j r = − hp,s ∑ ∆ ∆∆ ∆ s,j s cv,i −− − − ρ he,i ∑ r,j (4-21)

where Pr,j is the pressure loss in the servicing return line j. The j subscript on the

return line summation indicates only those return pipes servicing consumer i

between consumer i and the point in question.

The evaluation of pressures in the return pipes using this expression requires

some care and forethought to avoid errors and redundant calculations. Errors can

result if the summations include pipes other than the appropriate ones, which will

be different in the case of supply and return. Equation 4-20, as written, could be

evaluated for each consumer at all locations in the piping system. However, all that

is required is to find the pressure at each location once for any consumer served

through that point. The evaluation of the equation for all remaining consumers

served through that point would yield the same result and thus is not required. Some

simple rules will allow us to reduce the number of locations where calculation of the

pressure will be necessary. For example, consider the case where the entire system

is at or below the elevation of the heating plant. In this case, the minimum pressure

in the return line would be at the heating plant. In the supply line, however, the

lowest pressure could be at any location, dependent on the relative magnitude of the

hydrodynamic gradient from friction and the hydrostatic gradient from elevation

differences. If the entire system was at the elevation of the heating plant, then the

lowest pressure in the return line would be at the consumer, who is, in a hydraulic

sense, the most distant from the heating plant.

It is important to note that in the case of this absolute pressure constraint,

the supply and return piping must be considered separately, since the temperature,

and thus saturation pressure, will usually be quite different in each. Strictly

speaking, it would be necessary to determine the actual temperature at each

location in the system and compare the saturation pressure requirement with

the other applicable low pressure constraints to determine which one is dominant

there.

Once the pressure has been calculated at the locations where the minimum

pressure constraint could be active, these pressures would be compared to the

minimum allowed pressure for that location determined from the dominant constraint of the applicable ones given above. Thus, our minimum pressure constraint

for the supply pipe becomes

 

P P gz P P

j hp,s − ∑ ∆ ρ s,j s x,sat −≥ + saf . (4-22)

38
For the return pipe our minimum pressure constraint is

 

P P gz P P P P P

j j hp,s − ∑ ∆ ∆∆ ∆ s,j s −− − − ρ he,i cv,i r,j x,sa ∑ ≥ +t saf. (4-23)

And at the inlet to the pump located at the heating plant, we have two constraints

Php,r ≥ PNPSH (4-24)

Php,r ≥ Pa+ Pasa. (4-25)

Now we have identified all of the constraints on the multiple consumer–multiple

pipe solution and derived mathematical expressions for each one. In the next

chapter, we will examine a method of finding an optimal solution that satisfies all

of the constraints.

39
CHAPTER 5: OPTIMAL DESIGN OF SYSTEMS WITH

MULTIPLE PIPES AND CONSUMERS

As noted in the Chapter 4, all district heating networks, with the exception of pure

transmission systems, will have multiple consumers and pipes. If each pipe were

independent of the others, it would be possible to apply the procedure developed

earlier in Chapters 2 and 3 to each pipe independently and develop a complete

design in that way. Unfortunately, as the constraints introduced in Chapter 4 show,

each pipe segment does not operate independently of the others. Thus, the system

can not be designed completely in that way for all but the most trivial cases.

However, the “optimal independent design,” as we will call it, is very useful, even

though we can not guarantee that it would be feasible, for we can use it to form a

lower bound on any other designs that we might propose. We know that it is not

possible for us to achieve a lower cost design for any one of the pipe segments than

the one we have determined independently. Thus, it follows that we also are assured

that no design for the entire system of pipe segments can be lower in cost than the

sum of the costs for the optimal independent designs. While this may appear to be

of little significance to the designer who is subject to system constraints, it’s actually

a very useful result. It will serve two very important functions for us. First, it will give

us a lower bound on total system cost to which we can compare other designs to see

if they are sufficiently close to render additional effort at achieving better designs

impractical or unnecessary. The second function of this “optimal independent

design” will be as a starting point for a solution strategy that will move towards an

optimal solution that satisfies all the system constraints. Both of these attributes of

the optimal independent design will be exploited in this chapter.

Our objective here will be to develop methodologies that will help us find the

optimal discrete pipe diameters for systems with multiple pipes and consumers,

while minimizing the computational effort necessary, such that large networks

often encountered in practice may be treated with an acceptable degree of effort.

Because of excessive computational effort, many of the methods that have been

previously applied to problems of this type are felt to be unsuitable. Several of the

more common approaches are discussed very briefly below.

SOME POSSIBLE APPROACHES

The classical approach to a constrained optimization problem like this one is to

include all the constraints in the problem solution and find a solution that satisfies

all the constraints. Here, we have many constraints that would need to be included.

For our problem from eq 4-2 and 4-5, we would have two constraints for every

consumer. Equations 4-11 and 4-22 would result in two constraints for every node

in the supply pipe. Equation 4-23 gives us one constraint for each node in the return

pipe, and at the heating plant there are two additional constraints (eq 4-24 and 4-25).

If, for example, we considered a moderate sized system, with 100 consumers and 125

nodal points, we would have 577 constraints. Several of the more common methods

for handling such problems are discussed below.

The method of linear programming (Wilde and Beightler 1967) is very efficient at

solving optimization problems with large numbers of constraints. However, as the

name of the method implies, the objective and constraints must be linear. Here, we

have a highly nonlinear problem because of the pressure losses being proportional

to the pipe diameters raised to approximately the negative five power. Methods

have been devised (Reklaitis et al. 1983) to use linear programming algorithms on

nonlinear problems by making linear approximations about a point using a Taylor

series expansion about that point. The computational effort involved in the use of

40
such methods can be considerable for problems with many variables and constraints, as is the case here. Additionally, if an “optimum” is found by such a method,

there is no guarantee that it is a global optimum.

Methods have also been developed for general nonlinear problems. Perhaps the

method most commonly referred to is Lagrange’s method of undetermined multipliers (Wilde and Beightler 1967). This method requires that the solution to a set of

nonlinear equations be found. The number of unknowns is equal to the number of

variables (pipe diameters and consumer control valve pressure losses) plus the

number of constraints. In our case this would result in very large systems of

nonlinear equations for all but the most trivial problems. The solution of large

systems of nonlinear equations can be a very difficult task, usually done by adapting

methods for the solution of linear equations. For this reason, Lagrange’s method is

felt to be impractical for this problem.

The Generalized Reduced Gradient (GRG) method is a popular one used for

nonlinear constrained problems (Reklaitis et al. 1983). It is based on extending

methods used for linear problems to nonlinear problems. The basic concept of the

GRG method is to follow along the direction of a constraint subset while seeking

improvement in the objective function. By requiring some subset of the constraints

to be satisfied, the number of degrees of freedom of the problem can be effectively

reduced. When inequality constraints are present in the problem, as is the case here,

either an active set strategy must be adopted or slack variables must be introduced

for each constraint. Gill et al. (1981) indicate that GRG methods can encounter

difficulties when highly nonlinear constraints are involved, as is the case here.

Because methods developed for linear problems are being used for nonlinear

problems, it is necessary to iterate at each step to achieve a feasible design. The

Newton-Rapson method is used for this iteration and it becomes the main computation burden of the GRG method (Arora 1989). Other quasi-Newton methods have

been proposed, but they can cause other problems with this method (Arora 1989).

Vanderplaats (1984) indicates that convergence of the Newton-Rapson method may

be a problem when using the GRG method for highly nonlinear problems.

Another class of methods for general nonlinear constrained problems is the

penalty function methods (Rao 1984). These methods reduce the constrained

problem to an unconstrained problem that can then be solved using any of the

various methods suitable for such problems. With many variables, as we have here,

the multidimensional optimization problem that results can be quite time consuming to solve. In addition, it’s usually necessary to solve the problem repeatedly for

different values of the penalty parameter until some convergence criterion has been

met. A feasible starting point is required as is an initial value for the penalty

parameter and the multiplication factor that is used to adjust the penalty parameter.

In this problem the diameters of the pipe segments must take on discrete values

in the final solution, while other variables such as the consumer control valve

pressure losses are continuous. Such a problem, which has both discrete and

continuous variables, can be formulated as what’s called a “mixed integer” problem

(Reklaitis et al. 1983). The methods described above can not be applied directly to

integer or mixed integer problems. They must be used in combination with another

technique, most notably the “branch-and-bound” approach, to find the solution for

the discrete variables. The branch-and-bound approach will be discussed later.

In search of a simpler and more efficient method than those described above, we

will proceed by starting with our optimal independent (unconstrained) design and

identifying methods to move from this design to one that satisfies all the constraints.

We will attempt to conduct this process of modifying the solution so that it satisfies

the constraints in a manner that will keep us as close as possible to the true globally

optimal design.

41
PIPE AND PIPE JUNCTION LABELING SCHEME

Before proceeding further with the development of our solution technique, we

first need to develop a methodology for identifying consumers, pipe junctions and

pipe segments. We would like this system to be as simple and intuitive as possible,

yet sufficiently general so as to be easily extendible to much larger networks. A

method that meets these requirements is a simple identification number for each

“node.” A node can be any one of the following items within the pipe network.

1. A source node where a net inflow of heat occurs, i.e., the heating plant.

2. A sink node where a net outflow of heat occurs, i.e., a consumer.

3. A pipe junction node where no net inflow or outflow of heat occurs.

Note that there are at least a couple of special cases of the pipe junction node that

might be of interest: a storage node where heat could be stored for release at later

times, and a “junction” node with only two pipe segments connected. The latter

could be simply a transition in pipe size or an intermediate pumping station for

instance. These special cases would be of interest for advanced system optimization

studies but are beyond the scope here.

The number of a node does not necessarily need to be assigned in any particular

fashion. They could be assigned sequentially from the plant or some consumer, or

in no particular sequence at all. In fact, alphanumeric characters could be used for

identification. The point is that the assigned identification characters have no

significance relative to one another, other than being unique to the node in question.

With an identification system established for our nodes, we need to establish the

identity of the pipes connecting these. The simple convention we will adopt is to use

the node numbers on either end of the pipe segment to identify the pipe segment that

connects them. For example, the pressure loss in the pipe segment between nodes

1 and 2 would be written as ∆P1,2. We will establish the convention of letting the first

node number in the pair be the upstream node in the supply line, with the second

node being the downstream node, again in the supply line. For the return line of the

same pipe segment, the convention will be established by the supply line, i.e., the

first node number in the pair will be the downstream node in the return line and the

second node will be the upstream node. Note that a system segment, as we have

currently defined it, can not have any intermediate nodes within it.

Now we are ready to begin the development of our solution method. As always,

we start by determining our objective function.

SYSTEM OBJECTIVE FUNCTION

The objective function for an entire system of pipes will include the sum of the

individual objective functions for each pipe segment. We must also include the cost

of pumping energy dissipated at the consumer and the capital cost of the pumps

needed to generate this pumping energy. At first it might seem unnecessary to

include costs associated with the consumer in our objective function when in fact

there are no decisions to be made about the consumer’s equipment. However,

constraints that the consumer places on the system will require that these costs be

included in order to achieve an optimal design that does not violate these constraints.

Additional costs would also need to be included if we were to expand the

objective of our design. For example, if we wished to determine an optimal

operational strategy for the system, as well as a design, it would be necessary to

include some additional costs in the objective function. These would be the costs of

generating the heat ultimately supplied to the consumer. Another example of an

42
expanded objective would be if we also wanted to determine the optimal size for the

consumer’s heat exchanger equipment, which would require including these costs

in the objective as well. Here, however, we will not address these additional issues.

With these limitations in scope our objective function becomes

 

min. 

 d 

st fixed hl pev pv pvc

e

p pm

i r

CC CC C C

C P Pm t

j

i yr

= + ++ ∑( ) + +

∑ ∫( ) + ( ) η η

∆ ∆ cv,i he,i  ̇ /ρ

(5-1)

where Cst = total system cost ($)

Cfixed = fixed cost of pipes and pumps, and the maintenance and repair on

this portion of their costs ($)

Cpv = diameter variable cost of pipes and the maintenance and repair on

that portion of pipe cost ($)

Cpev = diameter variable cost of pumps and pumping energy attributable

to piping pressure losses, and the maintenance and repair on that

portion of pump costs ($)

Cpvc = variable cost of pumps attributable to the pressure losses at the consumer ($).

Notice that the density used in the last term of this equation is the taken at the

return condition. This is done because the pumps are usually located on the return

side of the system at the heating plant. The cost of pumps, which was previously

lumped with the piping cost, has been broken out as a separate cost since the number

of pumps will be discrete for the system.

In general, the mass flow rate for any consumer 

m ̇ i and the pressure losses at the

consumer (∆Pcv + ∆Phe)i will be functions of time. Previously, we assumed that the

mass flow rate over the yearly cycle was given by eq 3-25. Since the pressure loss in

the consumer’s heat exchanger ∆Phe,1 is a function of mass flow rate, as given by eq

4-6, it will also be a function of time. As we will show later, the pressure loss in the

consumer’s control valve ∆Pcv,i will be used to “balance” the network. Hence, it will

become a function of time in most all cases as well. We will have some choices as to

the best way to balance the network using the consumer’s control valve, as will also

be shown later.

In eq 5-1 we have separated the cost of the pumps into the fixed costs, that portion

which does not depend on pump capacity, and the variable costs, which are

attributable to either pressure losses at the consumer or in the piping network. We

have also separated the fixed portion of the pipe cost as well from that portion that

depends on pipe diameter. Effectively, this does not change our objective function

as far as terms that contain the pipe diameters are concerned, since the fixed costs

of the pipes and pumps are not considered in determining the optimum pipe

diameter, as can be seen from eq 2-19. For a multiple pipe system, these fixed costs

are

 

C PVF A A n A L

j

fixed m&r m&r 1 p 3 j = ( + )( + ) 1 ∑ . (5-2)

The variable cost of pumps attributable to pressure losses at the consumer Cpvc

will be determined by the pressure losses and flow rate at the design condition. It is

this condition for which the pressure difference between the supply and return at the

heating plant, as well as the mass flow rate, are greatest. Thus, the pumps must be

sized for this condition. This portion of the pump cost will be given by

Cpvc = A2 (1+PVFm&rAm&r) 

 

∆ ∆ P Pm

i

cv,i he,i i r

d

∑[ ] ( ) + ( )  ̇ /ρ . (5-3)

43
SOLUTION STRATEGY

Now we are ready to formulate our solution strategy. Inspecting the objective

function, eq 5-1, we see that the costs have been grouped with respect to their source.

All of the costs that are dependent on our decision variables, the pipe diameters,

have been included in the first summation over j, the pipe segment index. The

summations in the third and fourth terms are those that arise from the pumping

energy expended at the consumer. The decision variables in the terms of these

summations are the ∆Pcv,i values.

We notice immediately from the form of the objective that it is a separable

function with regards to the pipe diameter for each of the system segments j. A

separable function is a function of more than one variable that may be written as a

combination of functions, one independent function for each variable. Thus, from

examining the objective function, it appears that we can consider each pipe diameter

function independent of the other pipe diameters and find its optimum. We will

proceed as if this is the case, although later we will see that the constraints will not

allow the diameters to be considered completely independent of one another in all

cases.

We begin by inspecting the objective function for monotonicity, since this will

help us simplify the solution as much as possible. Looking first at the terms in the

summation over the pipe segment index j, we look at each term in the summation

separately

Chl, j = I1/ln(A10/d) 

dj

+ ( )

Cpv, j = A9d 

dj

+ ( )

Cpev, j = I3 d–(5+b+c)

 

dj

− ( ) .

The monotonicities with respect to dj of each term are given and we see that we

have both increasing and decreasing terms, so we are unable to use monotonicity

analysis on these at the outset. This is consistent with our findings in Chapter 2,

where we first neglected the Chl term and then used geometric programming theory

to find a solution to the lower bounding problem thus formed. This result was used

as a starting point for a simple search to find the solution to the problem without

neglecting Chl. Since the objective function is separable for each of the dj values, we

will proceed with the same methodology and find the “optimal independent”

values for each dj in the same way.

The other remaining decisions variables in the objective are the ∆Pcv,i variables,

of which there is one for each consumer. The ∆Pcv,i variables appear in both of the

last two terms of the objective function once eq 5-3 has been substituted for Cpvc

 

A PVF A P P m

i

2 ( ) 1+ m&r m&r cv,i ∑ [ ] ( ) + he,i ( ) i r

d

∆ ∆  ̇ /ρ 

( ) ∆Pcv,i

+

 

C P Pm t

i yr

e

p pm

cv,i he,i i r d 

η η

∑ ∫( ) ∆ ∆ + ( )  ̇ /ρ 

( ) ∆Pcv,i

+ .

Both of these terms are monotonically increasing in each ∆Pcv,i and thus the objective

is monotonically increasing in each ∆Pcv,i. The First Monotonicity Principle, MP1

(see Papalambros and Wilde 1988), therefore tells us that each ∆Pcv,i must be

bounded below by at least one active constraint. We will examine the issue of

determining the constraint activity for these decision variables.

44
Constraint activity for consumer control valve pressure losses

We have two possible constraints for each ∆Pcv,i

 

h P PP P P

j

1 =− + ∑( ) + + 0 









 ∆ ∆∆ ∆ ∆ hp s r = j cv he

i

 

( ) ∆Pcv,i

m (4-2)

g1 = ∆Pcvm,i – ∆Pcv,i ≤ 0 

( ) ∆Pcv,i

− . (4-5)

The inequality constraint g1 is monotonically decreasing in ∆Pcv,i so it bounds ∆Pcv,i

in the proper sense. The equality constraint h1 could also be “directed” (see

Papalambros and Wilde [1988] for procedure for directing equality constraints) such

that it would bound ∆Pcv,i in the proper sense. At this point it is not clear which of

these two constraints would bound each of the ∆Pcv,i variables in the proper sense.

In fact it’s entirely possible that the active constraint may vary depending on the

particular consumer in question. If no other decision variables appeared in these two

constraints, they would form a conditionally critical set for each ∆Pcv,i (see

Papalambros and Wilde [1988] for definition of conditional criticality).

However, we see that pressure increase across the pump at the heating plant ∆Php

has not yet been fixed. ∆Php is monotonically increasing in constraint h1; thus, it

becomes a monotonic nonobjective variable in our problem. The Second Monotonicity Principle, MP2 (Papalambros and Wilde 1988) tells us that either ∆Php is

irrelevant and can be deleted from the problem together with all the constraints in

which it appears, or it is relevant and bounded by two active constraints, one

bounding it from above and one bounding it from below. If for just one consumer

i the constraint h1 is critical for ∆Pcv,i, then ∆Php becomes relevant. A critical

constraint is an active constraint whose deletion would cause the problem to become

unbounded. This active constraint would have the following monotonicities when

directed to bound ∆Pcv,i in the proper sense

 

h P PP P P

j

1 =− + ∑( ) + + 0 









 ∆ ∆∆ ∆ ∆ hp s r ≡ j cv he

i

< 

∆ ∆ P P cv,i hp − + ( ) , .

If one constraint is critical for more than one variable in the problem, it is said to

be “multiply critical” (Papalambros and Wilde 1988) and this would be the case for

h1 above. Papalambros and Wilde (1988) warn that multiply critical constraints

should be eliminated from the problem whenever possible. Here, we have that

option since we can combine the h1 constraints for any two consumers and eliminate

∆Php from the problem. Before we do so, let’s consider briefly what is physically

happening in our problem.

First, we note that the pressure increase required across the pump at the heating

plant ∆Php appears in eq 4-2 for each consumer. Since, physically, we know that ∆Php

can only assume one value, it must be the greatest value that results from consideration of all the consumers. For the remaining consumers, ∆Pcv,i must be greater than

the minimum value ∆Pcvm,i. The consumer who requires the greatest ∆Php will be

called the “critical” consumer. Notice that the equality constraint of eq 4-2 can be

satisfied for the remaining consumers by letting ∆Pcv,i > ∆Pcvm,i as allowed by eq 4-

5. This is, in fact, how it is done in practice in most cases; the ultimate balancing of

the pipe network flows is done by the consumer’s control valves. For the case of the

critical consumer, eq 4-5 will be satisfied as a strict equality, i.e., ∆Pcv,i > ∆Pcvm,i.

While these arguments of constraint activity would appear to be completely

intuitive, since we would not want to supply any more pumping energy than

necessary, they can also be shown analytically as follows.

45
The pressure difference across the pump at the heating plant ∆Php must be equal

in all of the constraints of eq 4-2 at any instance in time, so we can write this constraint

set in the form

 

s r cv he

k

s r cv he

i

∆∆∆ ∆ ∆∆∆ ∆ PPP P PPP P

j j

∑( ) + + +



















 ∑( ) + + +



















 – 0 = (5-4)

where i ≠ k.

If we have n consumers, there will be (n – 1) such equality constraints containing

∆Pcv,i that apparently could be “directed” such that they would bound ∆Pcv,i from

below as required. However, these constraints are not all independent. Since we

started with n independent equations and then eliminated ∆Php, we will have at

most (n – 1) independent equations remaining. Below, we will show that these

(n – 1) independent equations may bound at most (n – 1) of the ∆Pcv,i objective

variables in the proper sense. The arguments made will be for an instant in time but

must hold for any time during the yearly cycle.

We begin by examining the ∆Pcv,i term for the consumer arbitrarily chosen to be

consumer “1.” Now, ∆Pcv,1 for consumer 1 can be bounded from below as required

by any one of the (n – 1) constraints in which it appears with another consumer.

Suppose we let the constraint with consumer “2” bound ∆Pcv,1. Now we have

(n – 2) constraints remaining that can bound ∆Pcv,2 in the proper sense, since the

equality constraint with consumer 1 has been directed such that it would bound

∆Pcv,2 in the improper sense

Similarly, let ∆Pcv,2 be bounded by properly directing its equality constraint with

consumer “3.” Now, at first it would appear that ∆Pcv,3 could be bounded by (n – 2)

constraints as well, since we have only directed the constraint involving ∆Pcv,2 in the

improper sense and any one of the remaining (n – 2) constraints can be directed as

needed. However, since we directed the constraint between ∆Pcv,1 and ∆Pcv,2 such

that it bounded ∆Pcv,1 below, we are not free to direct the constraint between ∆Pcv,1

and ∆Pcv,3 as needed; in fact, it must be directed in the opposite sense of that

required, that is, if

 

∆∆ ∆ ∆ ∆∆ ∆ ∆ PP P P PP P P

j j s r j cv he s r j cv he ∑( ) + + + < 



















 − + ∑( ) + +



















 ≡

2 1

0

 

∆Pcv,1 ∆P – cv,2 , + ( )

and

 

∆∆ ∆ ∆ ∆∆ ∆ ∆ PP P P PP P P

j j s r j cv he s r j cv he ∑( ) + + + < 



















 − + ∑( ) + +



















 ≡

3 2

0

 

∆Pcv,2 ∆Pcv,3 – , + ( )

then

 

∆∆ ∆ ∆ ∆∆ ∆ ∆ PP P P PP P P

j j s r j cv he s r j cv he ∑( ) + + + < 



















 − + ∑( ) + +



















 ≡

3 1

0

 

∆Pcv,1 ∆P – cv,3 , + ( ).

Thus, ∆Pcv,3 has only (n – 3) ∆Pcv,3 ∆P ∆P ∆P ∆P ∆P – cv,4 cv,3

– ,,,, , cv,5 cv,3 cv,n ++ + – { } ( )( ) ( ) L

constraints that could be directed to bound it in the proper sense.

If we continue to follow this line of reasoning, we find that when we reach

46
consumer i = n, we have no constraints remaining that can be directed in the proper

sense to bound the ∆Pcv,n term in the objective. Note that because the assignment of

the numerical value for the consumer index i is arbitrary, its assignment will have

no impact on what we have shown here and that this result would hold for any set

of consumer indices.

Since we have n variables ∆Pcv,i in the objective, one for each of the n consumers,

and the objective is monotonically increasing in each of these variables, we must

have n constraints bounding the ∆Pcv,i from below. Above we have shown that at

most (n – 1) of these constraints could result from eq 4-2. The only remaining

constraints on the ∆Pcv,i are the set formed by eq 4-5. Let’s assume for the moment

that all consumers have the same minimum pressure differential requirement for

their control valves, i.e., that ∆Pcvm,i is the same for all i. Now we see that it must be

the consumer with the minimum value of ∆Pcv,i whose constraint from eq 4-5 is

active. This is true since a consumer with any greater value would cause at least one

of the other constraints from the set of eq 4-5 to be violated, i.e., if consumer i has the

minimum control value pressure loss

 

∆ ∆ P P cv,k > cv,i for all k ≠ i.

And if for all i and k the minimum allowable control valve pressure losses are equal

 

∆ ∆ P P cvm,i = cvm,k.

Now, if the constraint for consumer i is active

 

∆ ∆ P P cv,i cvm,i >

then the constraint for consumer k can not be active

 

∆ ∆ P P cv,k cvm,k > .

We have already shown that (n – 1) of the constraints from the set of eq 4-2

are active and can thus be treated as equalities. These (n – 1) constraints

force the pressure loss summations to be equal for all n values of the consumer’s

index i

 

∆∆ ∆ ∆ PP P P

j s r j cv he

i

∑( ) + + +



















 = ∆∆ ∆ ∆ PP P P

j s r j cv he

k

∑( ) + + +





















for all i ≠ k.

This means that we identify the consumer whose value of ∆Pcv,i is at its minimum allowed ∆Pcvm,i by finding the consumer with the maximum value for the

quantity

 

∆∆ ∆ ∆ PP P P

j s r j cvm,i he

i

∑( ) + + +





















.

Once the pipe sizes are known, this quantity is easily calculated. This consumer who

has the minimum value of ∆Pcv,i is our so called “critical” consumer.

It follows then that there are only two cases where more than one of the

constraints from the set of eq 4-5 may be active:

1. In the case where the 

∆Pcvm,i values are identical for all consumers, if the

maximum value of

47
 

∆∆ ∆ PP P

j s r j he

i

∑( ) + +





















is identical for more than one consumer.

2. If the consumers have different values of 

∆Pcvm,i such that the maximum value

of the quantity

 

∆∆ ∆ ∆ PP P P

j s r j he

i

∑( ) + + cvm,i



















 +

is identical for more than one consumer.

Under either of these two conditions, the corresponding constraints from the eq

4-5 set for the consumers with identical values as described above will be active.

However, under neither condition will more than one be critical, since the deletion

of any additional constraints from the problem will not make the objective unbounded. This is true because the constraints from the set of eq 4-2 directed as

described earlier will bound the ∆Pcv,i of all the consumers but one in the proper sense.

What is not immediately apparent is why all (n – 1) of the constraints from the set

of eq 4-2 must be critical for some ∆Pcv,i objective variable. To illustrate why this is

so, consider the case where for an arbitrary consumer k when we let ∆Pcv,k be

bounded below by the constraint from the set of eq 4-5, thus ∆Pcv,k = ∆Pcvm,k. Now

suppose that

 

∆∆ ∆ ∆∆ ∆ PP P PP P

j j s r j he

k

s r j he

i

∑( ) + +



















 = + ∑( ) +





















and that ∆Pcvm,k > ∆Pcvm,i, where consumer i is the consumer with the minimum

value of ∆Pcv,i. The constraint from the set of eq 4-2 that states that

 

∆∆ ∆ ∆ ∆∆ ∆ ∆ PP P P PP P P

j j s r j cv he

i

s r j cv he

k

∑( ) + + + = 0



















 − + ∑( ) + +





















would be violated. Notice that this constraint has not been directed and shown to be

active, yet it still must not be violated by any feasible solution. Thus, this is not a

feasible solution and we see that only one of the constraints from the set of eq 4-5 may

be critical, since (n – 1) of the critical constraints must come from the set of eq 4-2.

Thus, our general result is that we must have (n – 1) of the constraints from the set

of eq 4-2 active and critical for (n – 1) of the ∆Pcv,i values. Thus, we have (n – 1) critical

constraints of the form

 

∆∆ ∆ ∆ ∆∆ ∆ ∆ PP P P PP P P

j j s r j cv,k he,k

k

s r j cv,i he,i

i

∑( ) + + + < 



















 − + ∑( ) + +



















 ≡ 0

(5-5)

Where the consumer index i does not equal the consumer index k and the monotonicities are

 

∆ ∆ P P cv,i

– cv,k+ ( ) ,

we also have one critical constraint on the remaining ∆Pcv,i not bounded below by

one of the (n – 1) constraints of eq 5-5 of the form

48
 

∆ ∆ P P cvm,i

– cv,i+ − < 0 

∆Pcv,i

– ( ). (5-6)

The assessment of which consumer will be the “critical” consumer having his ∆Pcv,i

bounded below by eq 5-6 was found above to be determined by finding the

consumer who has the maximum value of the sum of the non-control-valve pressure

losses ∆Pncv,i and the minimum control valve pressure loss ∆Pcvm,i. The non-controlvalve pressure losses ∆Pncv,i are given by

 

∆ ∆∆ ∆ P PP P

j

ncv,i s r j he

i

= + ∑( ) +





















(5-7)

where the non-control-valve pressure losses (N/m2) are easily computed once the

pipe sizes are determined using the procedure discussed below.

Initial pipe size determination

As noted earlier, our objective function is separable in each pipe diameter. The

pipe diameter function for each pipe segment j is increasing in some terms while

decreasing in others. Thus, we should be able to find a minimum cost for each

diameter by proceeding exactly as we did in Chapter 2 if we at first ignore the

constraints. Therefore, we first find the optimal “independent” discrete diameters

using the methods developed in Chapter 2.

System constraint satisfaction

Once our pipe sizes are determined, we need to ensure that the constraints are

satisfied and, if not, determine a methodology for achieving this. Below is listed the

various constraints that were developed in Chapters 3 and 4, categorized by the

portion of the system in which they originate.

At each consumer

 

∆ ∆∆ ∆ ∆ P PP P P

j hp s r j = + ∑( ) + + cv he (4-2)

 

∆ ∆ P P cv,i cvm,i ≥ . (4-5)

In the supply pipe

 

P P P gz

j max ≥− − hp,s ∑ ∆ ρ s,j s (4-11)

 

P P gz P P

j hp,s − −≥ + ∑ ∆ ρ s,j s x,sat saf . (4-22)

In the return pipe

 

∆ ∆ ∆∆ ∆ P P gz P P P P P

j j hp,s − −− − − ∑ ∑ s,j s ρ he,i cv,i r,j x,sat ≥ + saf. (4-23)

In the return pipe at the heating plant

 

P P hp,r ≥ NPSH (4-24)

 

P PP hp,r ≥ + a asa. (4-25)

All of these constraints deal with pressure levels at various points within the

system. Note that in eq 4-2 we have expressed the total piping pressure loss as its

supply and return components because it will be necessary to compute these

49
independently for the constraints given by eq 4-11, 4-22 and 4-23. Since verification

of satisfaction for all of these constraints requires either directly or indirectly the

calculation of the pressure losses in the supply and return pipes, we begin by doing

so for each of the pipe segments. The pressure loss in either the supply or return line

is calculated by modifying eq 2-15 slightly so that it applies to each pipe independently. The results are

 

∆P a m Ld b c c c b c

d,s j d,s d

j

( ) = ( ) ( ) 







+ + ( ) + + / /  ̇ – – – 2 4 2 1 2 5 ε π ρμ (5-8)

 

∆P a m Ld b c c c b c

d,r j d,r d

j

( ) = ( ) ( ) 







+ + ( ) + + / /  ̇ – – – 2 4 2 1 2 5 ε π ρμ . (5-9)

Once the piping pressure losses are known, we can calculate the non-controlvalve pressure losses ∆Pncv,i for each consumer using eq 5-7 and sum this with the

minimum control valve pressure loss ∆Pcvm,i to find the consumer with the highest

value of this sum, our critical consumer. The sum of the pressure losses for this will

become our pressure increase across the pump at the heating plant ∆Php, as given by

eq 4-2. For this consumer the constraint of eq 4-5 will be active, as shown earlier.

Using the value of ∆Php calculated for the critical consumer, we can then calculate

the control valve pressure losses for all of the other consumers using eq 4-2.

With the piping and consumer pressure losses known, we can calculate the absolute pressure level at all nodes in the pipe network with either a maximum absolute

pressure assigned to the supply pipe at the heating plant, or a minimum absolute

pressure assigned to the return pipe at the heating plant. If we set the minimum

pressure level in the return pipe at the heating plant, we can use the constraints of

eq 4-23, 4-24 and 4-25 to guide our choice. Note that when eq 4-23 is evaluated at the

heating plant, the entire left-hand side of the equation reduces to the value ∆Php,r.

Dependent on the particular parameter values for the problem at hand, one of these

constraints will “dominate” (see Papalambros and Wilde [1988] for concept of

constraint dominance). The cases for constraint dominance are simply as follows. If

 

PP P P P a asa NPSH x,sat +≤ ≤ + saf

eq 4-23 dominates. If

 

PP P P P a asa x,sat +≤ +≤ saf NPSH

eq 4-24 dominates. If

 

P P P PP x,sat + ≤ ≤+ saf NPSH a asa

eq 4-25 dominates.

Alternately, as noted above, we can also assign the maximum absolute pressure

in the supply pipe at the heating plant and use that value to find the other absolute

pressures in the network. The logical choice for the maximum absolute pressure

value in the supply pipe at the heating plant would be the maximum absolute

pressure allowable for the piping system being used Pmax. In most cases the

maximum absolute pressure in the system will occur at the heating plant in the

supply pipe; thus, this is a logical choice. It is possible that this will not be the case,

however. Using eq 4-15 we have shown earlier that the maximum pressure must be

at a nodal point location. In the discussion after eq 4-15, we also developed a

procedure that can be used to minimize the number of nodes at which the absolute

pressure must be calculated. If this procedure is used, we can quickly determine if

the heating plant will be the location of the maximum absolute pressure.

50
Pmax

P @ Design Flow s

P @ Low Flow r

P @ Design Flow r

P Return min

P Supply min

P @ Low Flow s

Distance from Plant

Pressure

Figure 8. Hypothetical pressure distribution under high and low flow

conditions and absolute pressure constraints.

All the calculations made to check for absolute pressure constraint satisfaction

should use the maximum (design) mass flow rate. This will assure that constraint

satisfaction will be possible at all flow rate conditions. Since the pressure losses in

the piping and the consumer’s heat exchanger will be greatest under this load

condition, the difference between supply and return pressure at the heating plant

will also be greatest under this load condition. Thus, under this condition the least

flexibility exists to adjust the supply or return pressure at the plant without violating

either the maximum pressure constraint in the supply, eq 4-11, or one of the

minimum return line pressure constraints at the plant, eq 4-24 and 4-25. If the

various maximum and minimum pressure constraints are satisfied for all points in

the network at the higher flow rate condition, it will always be possible to satisfy

them at the lower flow rates. This is easily shown graphically by considering the

pressures in the system along the piping route out to a consumer and back, as shown

for a hypothetical consumer in Figure 8.

In Figure 8 the horizontal lines are the constraints on the absolute pressures that

must be satisfied at all points along the route to the consumer. The solid lines that

have both positive and negative slopes are the supply and return pressures under

maximum load conditions. The magnitude and the sign of the slope of these lines are

determined by the sum of the hydrodynamic and hydrostatic pressure gradients as

given by eq 4-14. The dotted lines that behave in a similar fashion are the supply and

return pressures under some mass flow condition that is lower than the maximum.

In the extreme case where there is no flow, the pressure losses in the piping and

consumer equipment all vanish and the absolute pressure level is identical in the

supply and return lines for any point along the route. Also notice that we have

shown the pressure drop at the consumer as being lower at the reduced flow

condition. This results from lower pressure losses in the consumer’s heat exchanger

at the reduced flow rate (see eq 4-6) as well as lower losses in the consumer’s control

valve. If the network were ideally balanced and this consumer were the critical

consumer, his control valve would be completely open at all levels of load (i.e., flow

rate) and the pressure losses would always be the minimum possible.

By studying Figure 8, we can see that if we are able to “fit” the supply and return

absolute pressure lines within the constraints at the maximum flow condition, then

we can always do so for any lower flow condition simply by adjusting the absolute

pressure of either the supply or return at the heating plant. This results from the

51
hydrodynamic losses in the piping always being lower for the lower flow rates and

thus the total pressure increase necessary at the plant is reduced.

A simplified objective function

If all the constraints are shown to be satisfied, then we have found the optimal

solution to the multiple-consumer–multiple-pipe problem and we need not do any

further calculations. If, however, we find that constraints have been violated, we

need to refine our design. The first reaction of the designer when faced with this

result should be to closely examine the nature of the constraint that has been

violated. Often these constraints are “soft” and may be changed if the optimum

design indicates so. For example, in this problem one such constraint would be the

maximum allowed pressure. The designer has the option of using a higher pressure

class of piping if he or she feels it is warranted. This, of course, will most likely add

to the cost and, if this additional cost is significant, the designer may choose to

evaluate the design subject to the original constraint and the revised constraint, as

well to determine which one yields the lowest cost when the additional cost of the

higher pressure class piping is included.

Now that we have shown activity for some of the constraints, let’s consider the

problem again with a reduced objective and determine if solution is possible. Thus,

in our reduced objective, we are only interested in the terms in the objective function

that relate to variable piping costs, since we have shown that constraint activity

determines the values of the other decision variables. Thus, our problem can be

restated as

 

min pt hl pev pv j.C C C C

j

′ = ++ ∑( ) (5-10)

where Cpt ′ is the total diameter variable pipe costs for the system ($).

The constraints to which this solution is subject are that the absolute pressure

levels not be exceeded. The activity of the constraints of eq 4-2 and 4-5 fixes the

pressure increase at the heating plant and therefore the pumping power for the

system. Thus, at this point we no longer need to include the pumping power

consumed in the piping in our reduced objective function, eq 5-10. If we remove the

pumping power from eq 5-10, it becomes a monotonically increasing function of

pipe diameter. Thus, for the problem to be bounded, we must have monotonically

decreasing constraints on the pipe diameters. The sum of eq 5-8 and 5-9 forms one

such equality constraint for each pipe diameter that can be directed to bound it

properly. This constraint is

 

a m Ld P P b c c cc b c ε π ρμ ρμ 4 0 2 1 12 5 /  ̇ –– –– – ( ) ( ) + ( ) 











 



 − + ( ) ≡ < + + ( ) + +

d,s d,r d

j

s r j ∆ ∆ (5-11)

with the monotonicities being 

dP P j s,j r,j

−− − ( ) , , ∆ ∆ .

Notice that the constraints in the set of eq 5-5 will bound ∆Ps,j and ∆Pr,j in the

opposite sense to this constraint, so these nonobjective variables are bounded above

and below as required by MP2.

Thus, our problem is now to use these active constraints to solve for the diameters.

We will have one constraint from eq 5-11 for each pipe segment in the system. In

addition, we have already shown that we have (n – 1) active constraints from eq 4-

2, where n is the total number of consumers. We have one additional active

constraint from eq 4-5 for the critical consumer. However, at this point we are unsure

as to which consumer is the critical consumer; thus, we must include all n of the

constraints from the set of eq 4-2. In addition, we would still need to include all of

52
the constraints from the sets of eq 4-11, 4-22 and 4-23, since at this point we have no

way of knowing which of these constraints will be active. For each of these latter sets

of constraints, we have one member for each node in the system. We can use the

constraint of either eq 4-23, 4-24 or 4-25 coupled with the fact that ∆Php must be the

difference between Php,s and Php,r to eliminate Php,s from eq 4-11, 4-22 and 4-23. Once

the parameter values for the problem are known, the choice of which equation to use

will be determined by the dominance argument given earlier. Thus, if we have nn

nodes in the system and n consumers, our result would be 3nn + 2n simultaneous

equations.

To introduce the inequality constraints of eq 4-5, 4-11, 4-22 and 4-23 directly into

a set of simultaneous equations, we would need to introduce a “slack” variable for

each inequality constraint. The slack variable allows us to convert the inequality

constraint to an equality constraint. For example, eq 4-5 would be converted into an

equality constraint of the form

Pcv,i = Pcvm,i + Pcvs,i (5-12)

where ∆Pcvs,i is equal to the slack variable for consumer control valve pressure losses

(N/m2).

Equation 5-12 would introduce n slack variables into the problem. It could,

however, be used immediately to eliminate the ∆Pcv,i unknowns in eq 4-2, thus

reducing both the number of unknowns and equations by n to 3nn + n.

The constraints of eq 4-11 and 4-22 would each introduce nn slack variables as

well. The constraint of eq 4-23 would introduce one less slack variable than these

constraints, since there will be no slack variable at the heating plant where the

pressure level is determined by the constraint dominance arguments discussed

earlier. Thus, eq 4-23 will result in (nn – 1) additional slack variables.

In addition to the slack variables, we would still have the diameters of our pipe

segments as unknowns as well. The number of pipe segments will always be one less

than the number of nodes. This result is easily shown if we consider the process of

building the network from one node to the next. The first two nodes introduced into

the system will require one pipe. Any subsequent nodes introduced will require one

pipe for each node, since one node will already be an existing connected node. The

only case where this would not be true is if we had a looped network, rather than the

pure branched networks to which we will limit our discussion. We have one

additional unknown ∆Php that appears in all of the constraints of eq 4-2. So, then our

total number of unknowns would be (n + 4nn – 1).

With (n + 4nn – 1) unknowns and only (n + 3nn) equations, we have no unique

solution. Recall, however, that (nn – 1) of the unknowns are the pipe diameters,

which must take on discrete values. If the pipe diameters were to be considered as

continuous, we would have an infinite number of solutions. It’s actually fortunate

that they are discrete because this limits the number of possible solutions. The

number of possibilities can still be quite large for a system of any significant size. For

example, if we were to consider only 3 possible pipe sizes for each pipe segment we

would have 3(nn – 1) possible solutions. For our system discussed earlier with 125

nodes (124 pipe segments), we would have 3124 ≅ 1.46 × 1059, a combinatorial

problem of staggering proportions by any measure.

Notice that by applying monotonicity analysis to this problem we were able

reduce it to one of solving for the variables using the constraint set, which has been

reduced somewhat. The constraint set is linear in all the variables except the pipe

diameters and the pipe diameters only appear in one set of constraints. We could

make the problem linear by making the transformation for pipe diameters of

 

ˆ – d d b c

j j = ( ) 5+ + .

53
We would then substitute this into the constraints of eq 4-2 and solve the resulting

problem linearly in the 

ˆ

dj values. However, we would still have a significant task in

the solution of the system of equations. For this reason we will abandon the

possibility of achieving a solution by this approach.

Constraint resolution by pipe size refinement

We have an infeasible solution from the unconstrained problem. For each of the

consumers whose absolute pressure distribution of Figure 8 exceeds the constraints,

we need to reduce the piping pressure losses by increasing pipe diameters enough

to allow for constraint satisfaction. Since the system constraints all deal with

pressure levels in the network, we need to find a strategy to resolve these constraint

violations. Let’s attempt to find a solution by starting with our optimal independent

design and identifying methods to move from this design to one that satisfies all the

constraints. We will attempt to conduct this process of modifying the solution so that

it satisfies the constraints in a manner that will keep us as close as possible to the true

globally optimal design. We have the distinct advantage of knowing that our optimal independent design will form an absolute lower bound on system cost. At any

point we can compare the cost of our feasible design to the cost of the optimal independent design and determine if further attempts at improvement are warranted.

Examining Figure 8, we see that to bring excessive pressure differences within the

bounds of the constraints, we will need to reduce the slope of the pressure vs.

distance lines. The slope of these lines is the pressure loss per unit length of pipe.

Equations 5-8 and 5-9 tell us that if we are to reduce the slope we must do so by

increasing the pipe size. We would like to identify a method of determining which

pipe sizes to increase and by how much to satisfy constraints with minimum cost

increase.

At first it might seem that the best procedure would be to start by increasing pipe

sizes at the consumer’s end of the system, where the sizes are smallest and the pipes

tend to be shorter. In the smaller pipe sizes, the incremental increases in diameters

are in general less than for the larger pipe sizes. Thus, we could make smaller moves

away from the lower bounding cost. Starting at the consumer appears to be the most

logical way to proceed if the critical consumer is the only consumer who has

exceeded the absolute pressure constraint. In the more general case, however, more

than one consumer will have violated the absolute pressure constraints; thus, we

will examine that case first.

If more than one consumer has violated the absolute pressure constraints, we

could achieve constraint satisfaction by increasing pipe sizes that serve each

consumer individually until all the constraints are satisfied. Alternately, we could

increase pipe diameters in pipes that serve all of the consumers with violated

constraints. Because the pipe sizes are discontinuous (discrete) and the incremental

differences between adjacent diameters are nonuniform, it’s not possible for us to

predict a priori which diameters would be the best candidates for increasing. Thus,

we need to identify a method that will guide our search for a feasible and acceptable

solution expediently. In deciding when to stop our search, we always have the

benefit of knowing our lower bounding cost.

If we refer back to Figure 8, we see that satisfaction of the absolute pressure

constraints relies on keeping the pressure in the supply and return lines within the

bounds prescribed by the maximum absolute pressure constraint and the two

minimum absolute pressure constraints. We can adjust pressures at the plant to

achieve a state that satisfies all the constraints, as long as the maximum pressure

difference within the system does not exceed the absolute pressure constraints

discussed above. Since the critical consumer previously identified will be the

consumer who requires the largest pressure differential within the system, we will

examine this consumer’s requirements first and attempt to resolve the constraint

54
1. For the set of consumers whose constraints are violated, find the pipe segments that they all share

in common. Identify those pipe segment within this group that are shared with no other consumers.

In the event that there are no pipe segments shared with no other consumers, choose those pipe

segments shared with the minimum number of other consumers.

2. Increase the pipe segment diameters within this set enough so that the consumer whose constraint

from eq 4-2 was closest to satisfaction in the original solution now has his constraint satisfied. The

first pipe segment to have its diameter increased should be the one that this consumer shares with

the largest number of other consumers within the set of consumers with violated constraints.

3. Now remove this consumer whose constraint is satisfied from the set and find the new set of

common pipe segments for the remaining set of consumers with violated constraints, again including

the minimum number of other consumers in this set.

4. Again increase the pipe segment diameters within the remaining set enough so that the consumer

whose constraint was closest to satisfaction in the new set now has his constraint satisfied. The choice

of the pipe segment diameter to increase first is done in the same way as in step 2 above.

5. Repeat the above steps until all consumers have their constraints satisfied.

Figure 9. Method A.

violations that result. In the process we will consider the other consumers whose

constraints have also been violated.

Starting from an infeasible point, which is at the lowest possible cost for any

design, we want to move in the direction that will satisfy all of the constraints that

are violated by this solution. Since the critical consumer is the consumer whose

constraint has been violated by the greatest amount, we will have to travel the

“furthest” from our infeasible point to the boundary of his constraint. Thus, it would

seem tempting to try to resolve this constraint first and then look and see what other

violated constraints remain. However, it’s possible that we can plot a course that will

take us straight to a point that will resolve all constraints rather than handling them

one at a time. To do so we might consider the algorithm given in Figure 9.

Note that the last consumer to have his constraint satisfied is the consumer who

was identified as the critical consumer in the original solution. Now, however, all of

the consumers whose pressure constraints were violated in the original solution are

“critical” consumers as well, having pressure levels just meeting the constraints,

within the tolerance achievable with the discrete pipe diameters available.

As an alternate to the above methodology, we could proceed by adjusting

diameters of the critical consumer first, but only enough to bring his piping pressure

loss to the level of the next highest consumer, i.e., using the method in Figure 10.

Since in many cases consumers will share more than one pipe segment, we still

may be left with a number of alternatives that must be evaluated at each of the steps

above. If in each case we choose the alternative that produces the minimum amount

of increase in cost over the previous design, we should be able to move to an ultimate

solution that satisfies all the constraints while reducing the cost as much as practical.

Because we may be faced with many possible alternatives when a number of pipe

segments are shared by two or more consumers, we may decide to stop the process

after finding an alternative whose cost is within some reasonable tolerance of the

55
1. Proceed by first finding the set of servicing pipe segments unique to the critical consumer; this will

be his final service pipe only. Increase the diameter of that pipe segment until it reduces his pressure

loss to the same level as the consumer with the next highest pressure losses.

2. Now identify the pipe segments that these two consumers alone share and increase those pipe

diameters enough to reduce their pressure losses to the level of the next highest consumer. Note that

it may be that there are no shared pipe segments for these two consumers alone. In that event proceed

to the next step directly.

3. Again look for pipe segments shared by the three consumers with the highest pressure losses and

increase the diameters of those pipe segments enough to bring the pressure losses of these three

consumers to the level of the consumer with the fourth highest pressure loss. Once again, in the event

that no shared pipe segments exist, proceed directly to the next step.

4. Repeat this procedure until no consumers remain with pressure losses exceeding the constraints.

Figure 10. Method B.

lowest cost up to that point in the process.

To address the instances where more than one alternative is available at a

particular step in either of the processes outlined in Figures 9 and 10, we would like

a strategy that minimizes cost. Let’s investigate the effect of pipe diameter to see if

it would be to our advantage to choose smaller or large pipes as candidates for the

diameter increase.

First, we note that the capital cost Cpv,j is a linearly increasing function of pipe

diameter. Thus, an incremental increase in pipe diameter would have the same effect

regardless of the absolute value of the pipe diameter.

The cost of heat loss Chl,j is a somewhat complicated function of the pipe diameter.

It also includes an approximation introduced in Chapter 2. Within the range of

validity of the approximation (0.025 m ≤ d ≤ 1.0 m), we can see how the heat loss cost

behaves by examining its slope as shown in Figure 11.

The slope of the heat loss cost as plotted below in Figure 11 is essentially the first

term of eq 2-24 with the values of the parameters taken from the example of Chapter

3. From Figure 11 we see that the slope of the heat loss curve is always positive within

our range of interest. This tells us that whenever we increase the pipe diameter we

will increase heat losses, as we would expect. We also see that the slope is a

decreasing function of the diameter, except for pipe diameters over about 0.75 m,

where it becomes a slightly increasing function. For the portion of the range where

the slope is decreasing, we know that an incremental change in pipe diameter will

result in less increase in heat loss for larger diameters than for smaller ones.

The pressure loss as a function of pipe diameter is given by the sum of eq 5-8 and

5-9, which is our former eq 4-4

∆Ps&r = (a εb (4/π)

2+c A6 

m ̇ d

2+c L d –(5+b+c)

)j. (4-4)

If we take the partial derivative of this pressure loss with respect to diameter, we

have

∂∆Ps&r/∂d = –(5 + b + c) (a εb (4/π)2+c A6 

m ̇ d

2+c L d –(6+b+c)

)j.

56
1200

800

400

0 0.2 0.4 0.6 0.8 1.0

Pipe Diameter (m)

Heat Loss Cost Slope ($/m)

0

– 2000

– 4000

0.2 0.4 0.6 0.8 1.0

Pipe Diameter (m)

Slope of Pressure Loss (N/m ) 3

Figure 11. Slope of the heat loss cost term as a function of pipe diameter.

Figure 12. Slope of the pressure loss as a function of pipe diameter.

This result is plotted in Figure 12, where we have arbitrarily set the aggregate of all

the coefficients of d equal to –1. We see that the slope is everywhere negative and that

it is increasingly negative for decreasing values of the pipe diameter. This tells us

that an incremental increase in the pipe diameter will have a greater effect on

reducing the pressure loss at smaller pipe diameters than at larger ones. Thus, if we

have a choice of several possible pipe segments whose diameters we can increase,

we can achieve a larger pressure loss reduction for a given increase in pipe size by

choosing the pipe segment with the smallest diameter.

Unfortunately, our results tell us we should favor the smaller diameters from a

pressure losses reduction standpoint, but from the standpoint of heat loss costs,

57
larger diameters would be better. Thus, it is unclear from this analysis whether we

should try to choose smaller or larger pipes for increased diameters to satisfy

constraints. We do know that incremental increases in the discrete sizes of pipes will

be smaller in general for the smaller pipes. However, this does not necessarily mean

that we will be able to get closer to just satisfying the constraint in all cases by

increasing the diameter of the smaller pipes. Fortunately, the branch-and-bound

method mentioned earlier will provide us with a general solution strategy despite

our inability to better characterize the nature of the path to the solution.

Branch-and-bound method

The objective of the branch-and-bound method is to use what is known about

designs already explored to reduce the number of remaining ones that must be

examined in detail. An additional caveat is that we would like to do so without

dismissing any designs superior to the best feasible ones identified. According to

Reklaitis et al. (1983), the branch-and-bound algorithm is the most widely used

method for mixed integer problems and is the basis for most commercial computer

codes for solving such problems. The essence of the branch-and-bound method is

that it breaks the problem down into branches, each of which corresponds to some

particular choice of a single discrete decision variable.

The first step in the method is to compute the optimal solution to the problem with

all the variables assumed to be continuous. This becomes our global lower bound.

The branching usually starts by choosing what is felt to be the most fruitful branch

(variable) to explore; some criteria are given in Reklaitis et al. (1983). A discrete value

for the variable of this branch is chosen; this will be one of the discrete values

bracketing the optimum continuous value. The lower bounding cost for this branch

is found by finding the optimum continuous values of the remaining decision

variables with the branching variable fixed at its discrete value. At this point we can

explore the discrete designs within this branch by branching on the other variables

or we can go directly to the next branch. If we continue to explore this branch and

we find a discrete design sufficiently close to our lower bounding continuous

design, we can stop searching and accept this design. Otherwise, we move on to the

next branch. If we decide to search this branch further, we do so by comparing costs

obtained for designs within the branch by making permutations of the other

decision variables in turn to their bracketing discrete values.

This process continues until a feasible design with discrete values for all those

variables requiring such values is found. This is our first candidate design and its

cost becomes our upper bounding cost. Through the remainder of the process all

solutions will be compared in cost to this one until another feasible discrete design

with a lower cost is found. Any designs with higher cost will immediately be

rejected, and if continuous variables are still included in these designs for variables

that must ultimately take on discrete values, then all other designs within that

branch will be rejected as well. This can be done since we know that restricting any

of the continuous variables to discrete values will only increase the cost. The process

of rejecting these other designs for which the cost is never computed is called

“implicit fathoming,” as opposed to “explicit fathoming” where the cost is computed and the design is rejected because its cost exceeds our upper bounding cost.

Any remaining feasible discrete designs are also compared to the lower bounding

cost for the branch found with continuous values of the non-branching variables. If

any are sufficiently close to this lower bound, the search of the branch is concluded

and the design found is accepted as an adequate design representative of what can

be expected within the branch.

The next branch to be explored will use the other bracketing value of the first

branching variable. Its lower bounding continuous cost design is compared to our

current upper bounding cost from the best design of the previous branch. If the

58
lower bounding cost of this branch is less, then we continue the search of this branch,

proceeding as we did in the first branch, with one exception. We now have an upper

bounding cost and if at any point we find a design, either fully discrete of not, higher

in cost than our upper bounding design, we fathom this node and implicitly fathom

any branches of this design. However, as long as improvement appears possible, the

branch is searched until a discrete design acceptably close to the lower bounding

cost is found or all alternatives are exhausted. If a discrete design with a lower cost

than the lowest cost discrete design of the previous branch is found, it becomes the

new upper bound on cost. If this cost is sufficiently close to the continuous lower

bounding cost, then the search is concluded. Otherwise, another variable is chosen

to branch on. When exploring alternatives within any main branch, the same basic

branch-and-bound approach is applied within the “sub-branch.” Below we will

show how the branch-and-bound method is applied to our problem.

Solution by the branch-and-bound method

Before we can apply the branch-and-bound method as described above, we must

first have a feasible solution point from which to start. To find such a solution, we

use one of the methods described earlier in the section entitled Constraint Resolution

by Pipe Size Refinement. For example, if we use the method A step 2 (Fig. 9), we

calculate the continuous pipe size necessary to reduce the pressure loss to the next

level as described there. This is done by combining eq 4-2 and 4-4 to obtain

 

d

P PP P P

a Am L

j j

bcc

b c

j

hp,i s r j cv he

i

+ d j

(++)

1

1 ( ( / ) ) =

− + ∑ ( ) + +

































































≠

+

( ) −

∆ ∆∆ ∆ ∆

1

4 2

6

2

1 5

ε π  ̇

/

(5-13)

.

In this case the pipe segment index j

1 is the pipe segment whose diameter we have

chosen to increase and the consumer index i is for the consumer with the second

highest pressure loss at the heating plant as determined by eq 4-2. We continue to

use method A, which essentially repeats steps 1 and 2 until all the constraints that

were previously violated are now satisfied. When we have finished, we can check

our application of the method by taking the pipe segments that were in the sets of

those constraints previously violated and decreasing them to the next lowest pipe

diameter one at a time. If our application of method A was correct, each pipe

diameter that is decreased should result in the violation of at least one constraint.

At the conclusion of the use of method A, we will have a number of pipe segments

whose diameters are continuous as a result of the refinement process used by the

method. To obtain discrete diameter values for these pipe segments, we use the

branch-and-bound method as described above. Note that the unconstrained discrete

solution, what we have called our “optimal independent design,” will form a greater

lower bound than the unconstrained continuous solution. It may not, however, form

a greater lower bound than the constrained continuous solution, which we have

chosen not to attempt to find owing to the computational effort involved, as

discussed earlier. What we have found using method A is a feasible solution whose

cost is greater than the unconstrained discrete solution. The solution we have must

also be higher in cost than the constrained continuous solution, since the deletion of

the requirement of discrete sizes for those pipe segments that already have discrete

sizes would allow us to find a feasible solution at some lower cost. We can not be

certain, however, that the solution of method A is lower in cost than the constrained

discrete solution that we seek. For this reason, when we use the cost found in method

59
Figure 13. Hypothetical branch-and-bound problem.

A as our lower bounding cost, we can not be certain that this is the greatest lower

bound. If the application of the branch-and-bound method finds feasible solutions

with costs lower than our solution found by method A, then we must from that point

on use the unconstrained discrete solution as the greatest lower bound.

Note that, similar to method A (Fig. 9), method B (Fig. 10) would also provide a

number of pipe segments whose diameters are continuous as a result of the

refinement process used. Thus, we would proceed in the same manner with the

branch-and-bound method regardless of which of the two methods had been used

for the pipe size refinement process.

The application of the branch-and-bound method to our feasible solution of

method A or B is straightforward and proceeds as described earlier. Many of the

branches will have infeasible designs, since we know decreasing any of the continuous pipe diameters will result in the violation of at least one constraint, unless other

pipe diameters have also been increased. Thus, any branch that only decreases pipe

diameters in any one of the combinations in which they appear in the previously

violated constraints will be infeasible and need not be explored.

Consider a hypothetical case where we have four pipe segments with continuous

diameters after the application of either method A or B. If we limit ourselves to only

the two discrete pipe sizes that bracket the continuous values found, we will have

24 = 16 possible discrete solutions. A “tree” diagram of the problem is shown in

Figure 13. The symbols within each “node” of the tree represent the particular case

being evaluated. The “0” symbol indicates the continuous pipe size as found by

either method A or B. A “+” symbol indicates the next larger discrete pipe size and

a “–” symbol indicates the next smaller discrete pipe size. The sequence of the

symbols represents the order of the four pipe segments under investigation. As

noted above, any possibilities that only decrease the size of one or more pipe

segments are immediately known to be infeasible. These nodes have been shown

with dotted outlines in Figure 13.

Applying the branch-and-bound method, we would proceed initially by choosing the first pipe segment to branch on. Starting at the top of Figure 13, we branch

on the first pipe segment by computing the value of the objective function with the

branching diameter, rounded both up and down to the adjacent discrete diameter

values. However, we notice that one of the options, the (– 000) case, is infeasible.

Thus, there is no use in computing the value of the objective at that point since it’s

of no use as a lower bound on the constrained problem. In this case we then proceed

further down this branch in search of a feasible case. Suppose that we find that the

(– + + 0) is the first feasible case in this branch. We then branch on the last pipe

segment by examining the cases (– + + +) and (– + + –). Assume that we find both of

these to be feasible, but the (– + + –) case is lower in cost. This is our first completely

discrete and feasible design. It now becomes our upper bounding cost. If this cost is

sufficiently close to our current lower bounding cost of the (0000) case, then we can

(0000)

(+-++) (+-+-) (+--+) (+---)

(+-+0) (+--0)

(+-00)

(++++) (++++) (++-+) (++--)

(+++0) (++-0)

(++00)

(+000)

(-+00)

(-+++) (-++-) (-+-+) (-+--)

(-++0) (-+-0)

(--++) (--+-) (---+) (----)

(--+0) (---0)

(--00)

(-000)

60
accept this solution and terminate the search. Let’s assume, however, that this is not

the case.

Thus, we now need to start the process of “backtracking“ and examining the

delayed cases that we passed over. We begin by examining the case (– + – 0). Since

we found that the (– + 00) case was infeasible, then we immediately know that the

(– + – 0) case must be infeasible as well, since it only decreases diameters from an

infeasible case. Consequently, the (– + – –) case must be infeasible as well for the same

reason. The (– + – +) case might be feasible and let’s assume that we find this to be

true. Assume that the cost of this completely discrete design is less than the cost of

the (– + + –) design and thus it becomes our new upper bound on cost. Again, we

would compare the cost of this upper bound with the lower bounding cost of the

(0000) design and decide if further searching is warranted. Let’s assume that we are

still not satisfied with the gap between the upper and lower bounding costs, so we

decide to continue our search.

Backtracking further, we explore the alternatives (– – 0 +), (– – + +), (– – + –) and

(– – – +). Assume that we find these to all be infeasible. We then return back to the

other side of the first branch explored, the alternatives in the (+ 000) branch. Assume

that we find that both the (+ + 00) and the (+ – 00) cases are feasible, but the (+ – 00)

case has a lower cost. Thus, we defer any further exploration of the (+ + 00) branch.

We continue in the (+ – 00) branch by exploring the (+ – + 0) case and the (+ – – 0) case.

Suppose that both cases are feasible, but the (+ – – 0) case has a lower cost. Thus, we

continue by checking the cost and feasibility of the (+ – – +) and the (+ – – –) cases.

Assume that the (+ – – –) case is infeasible, but the (+ – – +) case is not only feasible

but lower in cost than our former upper bounding cost of the (– + – +) design. Assume

that this cost is indistinguishable from the lower bounding cost of the (0000) case,

and thus we accept the (+ – – +) design and terminate the search.

As illustrated here, one of the techniques used by the branch-and-bound method

is to continuously move the upper and lower bounding costs closer together until

the remaining possible improvement (i.e., reduction in cost) does not justify further

effort. At this point, the search can be stopped regardless of how many alternatives

have actually been explored. This process is accomplished by moving the upper

bound down as low as possible, i.e., finding the “least upper bound.” The least upper

bounding cost is always the lowest cost feasible discrete design found up to that

point in the process. The difference between the lower and upper bounding costs is

also refined by finding the “greatest lower bound.” The lower bounding cost is

initially determined by the continuous feasible design found by either method A or

B. As we proceed with the branching process, we will find lower bounding costs

within that branch for designs that have some, but not all, of their diameters at

discrete values. These greater lower bounding costs allow us to refine the difference

between the upper and lower bounding costs for that branch only. They also,

however, may tell us that further exploration of the branch is unwarranted if they

exceed our current upper bounding cost.

Once we have reached a solution by the method described above, there is an

additional area where we can seek further cost reductions. Note first that those

consumers who did not have their constraint of eq 4-2 violated will have pressure

losses in their control valves greater than the minimum allowed values. We then

observe that if these excessive pressure losses were absorbed by decreasing the sizes

of the pipes servicing only these consumers, our constraints would still be satisfied,

but the cost of the piping network would be reduced by use of the smaller pipe sizes.

This possibility is explored in the next chapter.

The procedure we have developed in this chapter to solve for the pipe diameters

in a general multiple-pipe–multiple-consumer system may be used with any

number of pipes and consumers. In the next chapter we will illustrate the use of the

method presented here on a simple example.

61
CHAPTER 6: A SIMPLE MULTIPLE-PIPE–MULTIPLECONSUMER EXAMPLE

To demonstrate our solution strategy, it’s illustrative to use an example. The

example should be as simple as possible, but must at the same time include the

salient features of a realistic system. Thus, our simplistic system must include as a

minimum:

1. A heat generating plant.

2. At least one heat consumer.

3. A section of pipe feeding only other sections of pipe but no consumers.

4. A pipe junction feeding at least two consumers with unequal hydraulic

characteristics.

5. A pipe junction feeding only one consumer.

With these requirements in mind we will examine the system shown in Figure 14.

First, we define the physical parameters for the network. For each of the nodes,

we assign an arbitrary elevation and, in the case of the consumers, an arbitrary

maximum heat demand expressed as maximum flow mass rate at the design

condition. The maximum demands assigned to the consumers are representative of

multiple residential consumers or large commercial loads. The assigned values are

given in Table 5.

We will also need to assign lengths to the pipe segments between the nodes. In

Table 6 some arbitrary pipe segment lengths have been assigned. In addition, in

Table 6 the flows in each pipe segment and the elevation change for that segment

have been determined based on the data in Table 5.

Note that the maximum flow in each pipe segment has been assumed to be the

sum of all consumer flows downstream (again, in the supply line sense) from it.

Although this must be true at any instance during the operation of the system

because of simple conservation of mass principles, for design it may not be the most

appropriate assumption. Since the demand for heat by each consumer will most

likely not be completely coincident in time, the maximum aggregate demand of the

consumers will always be somewhat less than the total of all consumers’ maximum

demands. This concept, called “demand diversity” is recognized by the district

heating industry and is sometimes accommodated to some extent in design calculations. If it’s not included directly in the design, it has the effect of providing an

additional safety factor. While at this point we will make no effort to include the

effect of demand diversity in our solution methodology, it is important to make note

of it, since it’s here that it would be introduced.

Figure 14. Simple multiple-pipe–multiple-consumer system.

5 4

6 1

7 3

2

Heating

Plant

8

62
Table 5. Assigned values for nodes of example in Figure 14.

Node Elevation Maximum demand

number Node type (m) (kg/s)

1 heat consumer 40 10

2 heat consumer 30 10

3 heat consumer 20 10

4 heat consumer 10 10

5 pipe junction 0 —

6 pipe junction 0 —

7 pipe junction 0 —

8 heating plant 0 40

Now we proceed to find the optimal pipe diameters for our example system of

Figure 14. As we found in the last chapter, a convenient starting point is what we will

call the “optimal independent design.” This is the design that we would arrive at if

we use the procedure developed in Chapters 2 and 3 for each pipe segment as if it

were independent of the others and its design were unconstrained by the system

constraints identified in Chapter 4. We will make use of all of the assumptions for

parameter values and operating strategy that we used in the examples of Chapters

2 and 3. For clarity these are repeated below:

Am&r = 2%/yr

Aη = 0.90 (dimensionless)

A1 = $1060/pump

A2 = $0.242/W

A3 = $218

A4 = $2180/m

A13 = 1.0 (dimensionless)

A14 = 0.575 (dimensionless)

A15 = 0.425 (dimensionless)

Tm = 6.4°C

ki = 0.030 W/m °C

ks = 1.3 W/m °C

Hp = 1.0 m

xi = 0.050 m

Table 6. Assigned values for the

pipe segments in the example of

Figure 14.

Elevation Maximum

Pipe Length change flow

segment (m) (m) (kg/s)

6,1 100 40 10

7,2 25 30 10

7,3 50 20 10

5,4 100 10 10

6,7 50 0 20

5,6 100 0 30

8,5 200 0 40

63
ε = 5 × 10–5 m

a = 0.119 (dimensionless)

b = 0.152 (dimensionless)

c = –0.0568 (dimensionless)

Ts = 120°C

Tr = determined by consumer model (eq 3-18) (°C)

Ce = $7.0 × 10–5/Wh

Ch = $3.4 × 1–5/Wh

PVFe = PVFh = PVFm&r = 9.08 (dimensionless)

Load control by flow modulation with consumer model (eq 3-25).

There are a number of additional parameters that were introduced in developing

the multiconsumer constraints for which we have not yet assigned any typical

values. They are:

∆Pcvm,i = 5 × 104 N/m2 (for all consumers, i = 1,4)

∆Phe,i = 1.0 × 105 N/m2 (for all consumers, i = 1,4)

Pmax = 1.0 × 106 N/m2

Psaf = 1.0 × 105 N/m2

PNPSH = 2.0 × 105 N/m2

Pa = 1.0 × 105 N/m2

Pasa = 0.5 × 105 N/m2.

Before we can find the optimal independent diameters for the pipe segments, we

need to calculate the remaining parameters that are determined by the assumptions

above. Because the optimal pipe diameter for a single pipe segment is independent

of the pipe length and elevation (see eq 2-23), the optimal independent diameter will

be the same for pipe segments (6,1), (7,2), (7,3) and (5,4). Thus, we construct Table

7 with the parameter values needed and the resulting optimal independent diameters. In each case, we have proceeded as before by solving the Lower Bounding

Problem (LBP) (eq 2-20), which neglects heat losses first and, subsequently, using

that as a starting point for finding the solution to the complete problem including

heat losses (eq 2-24). Also, as earlier, FORTRAN programs I1EQ3-26 and I2-C-GMT

were used to compute I1 and I

3 respectively.

The optimal diameters found above do not necessarily correspond to actual

discrete pipe diameters available, so before we check this solution to see if it satisfies

the constraint set, we first need to determine what the optimal discrete diameters

would be. Table 8 contains pipe size data for standard metric pipe sizes in the range

needed for our example.

To find the optimal discrete diameters, we proceed as before in the example of

Chapter 2 by simply examining the total cost of the discrete pipe diameters that

Table 7. Parameter values and optimal independent diameters for example of Figure 14.

Pipe I1/L I3/L d by LBP d by eq (2-24)

segment ($/m) ($ m4.095) (m) (m)

(6,1), (7,2), 73.3 4.276 × 10–5 0.0691 0.0666

(7,3), (5,4)

(6,7) 73.3 3.289 × 10–4 0.0966 0.0932

(5,6) 73.3 1.085 × 10–3 0.1175 0.1134

(8,5) 73.3 2.529 × 10–3 0.1350 0.1304

64
Table 9. Discrete bounding diameters and variable costs for the example (optimal non–

discrete diameters shown in boldface, optimal discrete diameters shown in italic).

Variable costs ($/m)

Pipe d Discrete cost

segment (m) Heat loss Capital Pumping Total premium (%)

0.0545 64.76 140.61 117.31 322.68 12.44

(6,1), (7,2), 0.0666 72.90 171.83 42.23 286.96 —

(7,3), (5,4) 0.0703 75.34 181.37 32.06 288.78 0.63

0.0825 83.24 212.85 109.11 405.20 4.14

(6,7) 0.0932 90.02 240.46 58.62 389.09 —

0.1071 98.65 276.32 28.87 403.84 3.79

0.1071 98.65 276.32 95.21 470.18 0.85

(5,6) 0.1134 102.51 292.57 71.16 466.24 —

0.1325 114.03 341.85 32.19 488.07 4.68

0.1071 98.65 276.32 222.03 597.00 12.51

(8,5) 0.1304 112.77 336.43 81.44 530.64 —

0.1325 114.03 341.85 75.07 530.95 0.058

Table 8. Standard metric steel pipe sizes

(data from DFF 1985).

Nominal Outer Wall Inner

diameter diameter thickness diameter

(mm) (mm) (mm) (mm)

50 60.3 2.9 54.5

65 76.1 2.9 70.3

80 88.9 3.2 82.5

100 114.3 3.6 107.1

125 139.7 3.6 132.5

bound the optimal diameters we have found. The bounding diameter that has the

lowest total cost will become our optimal discrete diameter. Obviously, it is only

necessary for us to compute the portion of the total cost that is dependent on pipe

diameter in making this decision. Table 9 gives the cost data for the bounding

discrete diameter for each of the pipe segments of our example. The costs are

calculated on a unit length basis using eq 2-23 divided through by L, the pipe length.

The portions of the total variable cost ascribable to each of the major component costs

are also given in Table 9. The parameter values of I

1/L and I3/L used for each pipe

segment are the same as those given in Table 7.

Now we need to consider the constraints on our multiconsumer system as

derived earlier in Chapters 3 and 4. These are summarized in Chapter 5 in the System

Constraint Satisfaction subsection. All of these constraints deal with pressure levels

at various points within the system. Since verification of satisfaction for these

constraints requires calculation of the pressure losses in the supply and return line,

we begin by doing so for each of the pipe segments. The pressure losses in the supply

or return pipes are calculated with eq 5-8 and 5-9 using the optimal discrete

diameters we have determined independently. The results are given in Table 10.

Satisfaction of the constraint of eq 4-2 at each of the consumers requires that we

65
Table 10. Pressure losses for the pipe segments in the example.

Pipe Discrete Length Elevation Flow rate Pd ,s Pd ,r

segment diam. (m) (m) diff. (m) (kg/s) (N/m2) (N/m2)

6,1 0.0703 100 40 10 91,585 91,652

7,2 0.0703 25 30 10 22,896 22,913

7,3 0.0703 50 20 10 45,793 45,826

5,4 0.0703 100 10 10 91,585 91,652

6,7 0.1071 50 0 20 20,615 20,630

5,6 0.1071 100 0 30 90,655 90,721

8,5 0.1325 200 0 40 107,219 107,297

sum the pressure losses in the portions of the piping system servicing each consumer. Recall that earlier we assigned values for the minimum pressure drop in the

consumer’s control valve and the pressure drop in the consumer’s heat exchanger,

∆Pcvm,i and ∆Phe,i, respectively, as follows:

∆Pcvm,i = 5 × 104 N/m2 (for all consumers, i = 1,4)

∆Phe,i = 1.0 × 105 N/m2 (for all consumers, i = 1,4).

We will also assume initially that for each consumer ∆Pcv,i = ∆Pcvm,i. The process of

calculating the ∆Php for each of the consumers is summarized in Table 11.

The pressure increase required across the pump at the heating plant ∆Php that we

have calculated is different for each consumer. As shown in Chapter 5, since ∆Php

can only assume one value, it must be the greatest value that results from consideration of all the consumers, and for the other consumers ∆Pcv,i will be greater than

∆Pcvm,i. The consumer who requires the greatest ∆Php is called the critical consumer.

Referring to Table 11, we see that consumer 1 is our critical consumer and thus

∆Pcv,1 = Pcvm,1. This has been determined for the maximum load condition, but will

hold for all load conditions since we have assumed that all consumers have loads

varying in the same manner over the yearly cycle. For the other consumers, eq 4-2

will require that ∆Pcv,i > ∆Pcvm,i. By using the ∆Php calculated for consumer 1, our

critical consumer, we find the following values for the ∆Pcv,i of the other consumers

at the maximum load condition:

∆Pcv,2 = 146,183 N/m2

∆Pcv,3 = 100,374 N/m2

∆Pcv,4 = 231,376 N/m2.

Notice that the control valve pressure drops for all these consumers are high

compared to the minimum value of 5 × 104 N/m2. This situation sometimes makes

it difficult for the control valve to function properly. It also represents a wasteful

Table 11. Heating plant pressure increase required by eq 4-2 for each consumer.

Consumer Servicing Pcvm,i + Phe,i 

∆Php

index, i pipes, j (N/m2) (N/m2) (N/m2) (N/m2)

1 (6,1), (5,6), (8,5) 289,459 289,670 150,000 729,129

2 (7,2), (6,7), (5,6), (8,5) 241,385 241,561 150,000 632,946

3 (7,3), (6,7), (5,6), (8,5) 264,281 264,474 150,000 678,755

4 (5,4), (8,5) 198,804 198,949 150,000 547,753

 

∆P

j

∑ r,j

 

∆P

j

∑ s,j

66
practice, since pumping energy, inherently more expensive than heat energy, is

being converted into frictional heating of the fluid. As noted at the close of the last

chapter, there is an alternative to these high control valve pressure losses: reduce the

pipe sizes further such that the pressure differential at the consumer’s control valve

is reduced. Such a practice was proposed by DFF (1985), where they suggest

reducing the size of the “service pipes,” those ultimately connecting the consumer

to the network. It may also be possible to reduce some of the pipes sizes within the

network as well. For our example problem, we see that we have limited options.

Consumer 1 is our critical consumer, so we cannot reduce any of the pipe sizes

servicing this customer; this rules out the pipe segments (6,1), (5,6) and (8,5). The

remaining pipe segments are (7,3), (6,7), (5,4) and (7,2). Thus, we investigate the

possibility of reducing the size of these pipes.

First, let’s look at pipe segment (5,4). This is the only pipe segment serving

consumer 4; thus, this is the only option for reducing the pressure loss in this

consumer’s control valve. What we would like to do is find the minimum pipe size

that will not violate the constraint of eq 4-2 for consumer 4. Effectively, what has

happened here is that we have removed the pumping power term from the objective

function so it now becomes monotonically increasing in d(5,4). We need to find the

constraint that will bound d(5,4) below. While not immediately obvious, eq 4-2 forms

such a constraint on d(5,4) when directed as follows

 

h PP P P P

j

1 = + ∑( ) ∆∆ ∆ ∆ ∆ ++− ≡ 0 s r j cv he hp < 

∆ ∆ P P s, 5,4 r, 5,4 ( )

+

( )

+ ( ) , .

From eq 4-4, we see that ∆Ps,(5,4) and ∆Pr,(5,4) are related to d(5,4) by

h2 = {[(ρ–1 μ–c)d,s + (ρ–1 μ–c)d,r]

[(a/2) εb (4/π)

2+c

 

m ̇ c d

2+ L d–(5+b+c)

](5,4)} – (∆Ps)(5,4) – (∆Pr)(5,4) ≡< 0

with the monotonicities being 

hd P P 2 5 4 = ( ) () () () ,

–– – , , ∆ ∆ s, 5,4 r, 5,4 .

So, we see that d(5,4) is bounded below by h2 and that the non-objective variables

∆Ps,(5,4) and ∆Pr,(5,4) are bounded below by this constraint and above by h1,, as

required by the second monotonicity principle (see Papalambros and Wilde 1988).

Now we can use constraints h1 and h2 to find the optimal value of d(5,4). To do so we

treat h1 as a strict equality and solve for (Ps)(5,4) + (Pr)

(5,4). We then substitute the

result into h2, again treating it as a strict equality, and solve for d(5,4). The result is

d(5,4) = 0.0614 (m) .

The discrete diameters that bracket this value are 0.0545 and 0.0703 m. The lower

bracketing discrete diameter will cause constraint h1 to be violated since a decrease

in d(5,4) will increase 

 

∆ ∆ P P

j s r j ∑( ) + . The optimal discrete diameter determined previously was 0.0703 m; thus, we are unable to improve on this result.

Let’s look at the remaining pipe segments (7,3), (6,7) and (7,2). These pipe

segments serve both consumers 2 and 3. Consumer 2 is served by pipe segments (6,7)

and (7,2) and consumer 3 is served by pipe segments (7,3) and (6,7). Both consumers

are served by pipe segment (6,7); thus, any decisions we make about this pipe

segment must be checked to ensure that both consumer constraints (eq 4-2) are

obeyed. Also, notice that if we decrease one of the pipe sizes and this violates

constraint h1, we may be able to increase the other pipe size in the pair serving that

consumer such that the total costs for the pipes and heat losses are reduced but

constraint h1 is still satisfied. It is also possible that a pipe size could be reduced or

67
Table 12. Pipe size combinations for the example.

Pipe segment

Combination

number (6,7) (7,2) (7,3)

1 0* 0 0

200+

300–

40+0

50++

60+–

70–0

80–+

90––

10 + 0 0

11 + 0 +

12 + 0 –

13 + + 0

14 + + +

15 + + –

16 + – 0

17 + – +

18 + – –

19 – 0 0

20 – 0 +

21 – 0 –

22 – + 0

23 – + +

24 – + –

25 – – 0

26 – – +

27 – – –

*0 = pipe size unchanged; + = pipe size increased; – = pipe size

decreased.

increased by more than one discrete size. We will ignore this possibility for the

moment and return to it later, since it would result in many more combinations to

be checked, most of which would violate h1.

If we first look at all the possible combinations of increasing or decreasing the

three pipe sizes without regard to the constraints, we have 33 = 27 independent

possibilities; they are enumerated in Table 12. Combination number 1 is our design

as it now stands, the “do nothing” option. A number of these combinations are

known not to yield improvement in our design, however, and may be immediately

dismissed without further evaluation.

Specifically, any combination that increases any pipe sizes while decreasing none

will only result in additional pipe capital and heat loss costs and thus will be worse

than our design as is. Thus, the combinations 2, 4, 5, 10, 11, 13 and 14 can be

dismissed.

In addition, we know that any combination that increases the diameter of either

the final pipe servicing consumer 2 [(7,2)] or consumer 3 [(7,3)], while decreasing the

other and leaving pipe segment (6,7) unchanged, would be more costly than doing

the same yet not increasing the diameter of the one pipe; thus, we eliminate

combinations 6 and 8. As we proceed to explore the various combinations remaining, we will discover that many other possible combinations will immediately be

shown to be infeasible by the infeasibility of related combinations.

In Table 13 we have listed the remaining combinations. Table 13 also gives the

68
status of the two consumer constraints that must be satisfied and the total of the

variable portions of the capital costs and heat loss costs for each combination. We see

by examining the constraint satisfaction that only two combinations are feasible, i.e.,

they satisfy the h1 constraint for both consumers 2 and 3. However, when we

calculate the cost of these feasible combinations, we find that both cost more than our

original design. Thus, we are left with the result that none of the alternatives

investigated so far are better than our original design. There are some additional

designs that we have not investigated, however. Recall that earlier we dismissed the

possible designs that would increase or decrease pipe sizes by more than one

discrete size from the original design. Depending on how many pipe sizes we are

willing to deviate from our original design, there are many alternate designs. Of

course, there is no guarantee that these designs will be feasible, let alone lower in cost

than the original design. To explore these designs without resorting to “exhaustive

enumeration,” i.e., calculating the constraint satisfaction and cost of each, we can use

the branch-and-bound technique described in detail in Chapter 5. Below we apply

this technique to our example problem. In the process of doing so, we will not only

explore additional designs not considered yet, but we will show how the technique

would have allowed us to dismiss some of the alternatives in Table 13 without

computing the constraint satisfaction or total variable cost.

As noted in the previous chapter, the objective of the branch-and-bound technique is to use what is known about designs already explored to reduce the number

of remaining ones that must be examined in detail. In addition, we would like to do

so without dismissing any designs superior to the best feasible ones identified. We

have effectively already used the technique above to dismiss nine of the possible

combinations of Table 12. In that case, we used the fact that the variable portions of

the heat losses and capital pipe costs were monotonically increasing in pipe

diameter. This allowed us to dismiss cases that only increased pipe size.

After our initial elimination of nine combinations, as discussed above, we see that

half of our remaining combinations involve the case where d(6,7) is reduced; thus, we

69

Table 13. Constraint satisfaction and costs for the remaining combinations. Our original design

(combination no. 1) is shown in bold, and the other feasible designs are shown in italic.

Cost

Comb. no. d(6,7) d(7,2) d(7,3) (h1)2 (h1)3 Variable premium

and type (m) (m) (m) (N/m2) (N/m2) costs ($) (%)

1 (0,0,0) 0.1071 0.0703 0.0703 –96,183 –50,374 38,002 0

3 (0,0,–) 0.1071 0.0703 0.0545 –96,183 193,207 35,435 –6.76

7(0,–,0) 0.1071 0.0545 0.0703 25,607 –50,374 36,718 –3.38

9(0,–,–) 0.1071 0.0545 0.0545 25,607 193,207 34,151 –10.13

12(+,0,–) 0.1325 0.0703 0.0545 –123,482 165,908 39,480 3.89

15(+,+,–) 0.1325 0.0825 0.0545 –149,022 165,908 40,464 6.48

16(+,–,0) 0.1325 0.0545 0.0703 –1,692 –77,673 40,764 7.27

17(+,–,+) 0.1325 0.0545 0.0825 –1,692 –128,753 42,732 12.45

18(+,–,–) 0.1325 0.0545 0.0545 –1,692 165,908 38,196 0.51

19(–,0,0) 0.0825 0.0703 0.0703 18,468 64,278 34,058 –10.38

20(–,0,+) 0.0825 0.0703 0.0825 18,468 13,198 36,027 –5.20

21(–,0,–) 0.0825 0.0703 0.0545 18,468 307,859 31,491 –17.13

22(–,+,0) 0.0825 0.0825 0.0703 –7,071 64,278 35,042 –7.79

23(–,+,+) 0.0825 0.0825 0.0825 –7,071 13,198 37,011 –2.61

24(–,+,–) 0.0825 0.0825 0.0545 –7,071 307,859 32,475 –14.54

25(–,–,0) 0.0825 0.0545 0.0703 140,259 64,278 32,774 –13.76

26(–,–,+) 0.0825 0.0545 0.0825 140,259 13,198 34,743 –8.58

27(–,–,–) 0.0825 0.0545 0.0545 140,259 307,859 30,207 –20.51
will explore that “branch” first. With d(6,7) assigned a discrete diameter one size

lower than our original design, we can use the constraint h1 for consumers 2 and 3

to find the lower bounding continuous values for d(7,2) and d(7,3). We obtain

d(7,2) = 0.0778

d(7,3) = 0.0892.

Thus, any combinations with discrete pipe diameters less than these need not be

considered, since they would violate the h1 constraint. This rules out combinations

19, 20, 21, 25, 26 and 27 because these would violate the h1 constraint for both

consumers. Also, we see that combinations 22, 23 and 24 would all violate the h1

constraint for consumer 3, so they are infeasible as well. Thus, we have eliminated

all the combinations in this branch as originally proposed. As noted earlier, there are

combinations that deviate by more than one pipe size from our original design that

were not considered. Before exploring any of these, we compute the cost of the

design above with continuous diameters to see if it is an improvement on our

original design. When doing so we find that the variable cost portion of the heat loss

and pipe capital costs is slightly less than our original design: a 0.77% reduction. At

this point we could decide not to further explore this branch, since it offers such a

small potential for improvement; however, we will continue since it illustrates the

method to be used. From combinations 22, 23 and 24, we know that if we increase

the pipe size of d(7,3) to the next discrete pipe size greater than 0.0892, the h1

constraint for consumer 3 will be satisfied as well. Thus, we propose the discrete

design

d(6,7) = 0.0825

d(7,2) = 0.0825

d(7,3) = 0.1071.

We know that this design is feasible, so now we need to compute its cost to see if

it’s an improvement over our original design. When the variable portion of the heat

loss and pipe capital cost is computed, we see that it’s 7.77% greater than the original

design. Thus, we dismiss this design as well as any other feasible designs in this

branch, since all other feasible designs would need to have larger discrete diameters

and would thus be more costly yet.

We have two other major branches yet to explore: one where d(6,7) remains the

same as in the original design and one where it is increased one discrete pipe size.

The latter branch has four combinations remaining, one more than the other branch,

so we will explore it first. We proceed as before by using the h1 constraint for

consumers 2 and 3 to find the lower bounding values for the continuous diameters

of d(7,2) and d(7,3), obtaining

d(7,2) = 0.0544

d(7,3) = 0.0624.

As before, we also compute the total variable cost portion of the heat losses and pipe

capital costs for this design. We find that this cost is 3.88% greater than our original

design. Thus, we need not look at any discrete designs in this branch, since all will

require larger discrete diameters than those continuous diameters found above and

thus they will be more costly. Note that the two feasible combinations 16 and 17 in

this branch identified in Table 13 do in fact have costs in excess of 3.88% above the

original design.

Now we explore the remaining branch, where d(6,7) is the same discrete pipe size

as found in our original design. As before, we compute the minimum continuous

70
diameters using the h1 constraint for consumers 2 and 3 to find the lower bounding

values for the continuous diameters of d(7,2) and d(7,3), obtaining

d(7,2) = 0.0564

d(7,3) = 0.0646.

We know that this design must have a lower cost than the original since it has the

same pipe size for segment (6,7) and smaller pipe sizes for the other two pipe

segments. However, we go ahead and compute this cost saving to see if it justifies

exploring this branch further. We find that the saving is a significant 5.40%. We have

three combinations (3, 7 and 9) from Table 13 that have not been previously

eliminated from this branch. We see, however, that each of these will violate at least

one of the h1 constraints, since at least one of the pipe sizes is smaller than the

continuous minimums found above. Thus, we can dismiss all of these combinations.

In addition, we can dismiss any other designs in this branch as well, since they would

have pipe sizes greater than those of our original and would thus be more costly.

Note that our original design is actually in this branch, using the first discrete pipe

sizes greater than those found above for d(7,2) and d(7,3).

By using the branch-and-bound technique, we have eliminated all of the combinations of Table 13 and have only computed the cost four times. In addition, we

have computed diameters using the h1 constraint six times. These computations

compare favorably with those required for total “exhaustive enumeration” of the

possibilities (27 cost and 54 constraint calculations) and favorably to the computations of Table 13, which eliminated nine possibilities based on monotonicity considerations. We have also shown that no other discrete designs in the branches

explored, i.e., even those deviating by more than one discrete pipe size, could be both

feasible and less costly that the original discrete design. What remains to be shown

is that other branches that allow d(6,7) to deviate by more than one discrete pipe size

are either infeasible or not cost effective, or both.

To explore the branches where d(6,7) is more than one discrete diameter away from

our original design, we once again look to the constraint h1 for consumers 2 and 3.

We notice that there is a limit on how much we can decrease either d(7,2) or d(7,3) and

still find values of d(6,7) that will satisfy the constraints. Physically, what has

occurred is that we have decreased the pipe sizes of d(7,2) or d(7,3) to the point where

all of the pressure loss available between the pipe junction node 6 and either

consumer 2 or 3 is being absorbed in the pipe segment (7,2) or (7,3). To utilize this

condition, we first ignore the pressure loss of pipe segment (6,7) and calculate the

minimum continuous values for d(7,2) and d(7,3) that will satisfy h1 for consumers 2

and 3 respectively. We then find the next largest discrete diameters in each case,

since any actual design would be bounded by these. The results are

d(7,2) = 0.0536 (continuous); 0.0545 (discrete)

d(7,3) = 0.0614 (continuous); 0.0703 (discrete).

Now, with these discrete diameters, we calculate the minimum continuous value

of d(6,7) that would satisfy the constraint h1 for both consumers 2 and 3. This value

is determined to be 0.1296. The minimum discrete value of d(6,7) is then 0.1325. We

see that these discrete diameters are identical to those of combination 16 in Table 13.

It was shown earlier that for this combination the cost exceeded our original design.

Thus, any larger discrete diameters would also exceed the cost of our original

design. Since this result is for the minimum possible discrete diameters for pipe

segments (7,2) and (7,3), regardless of the size of pipe segment (6,7), no lower cost

alternatives can exist since their diameters for pipe segments (7,2) and (7,3) would

be greater and thus the designs more costly.

71
Now the remaining branches yet to be explored are those with discrete pipe sizes

for pipe segment (6,7) more than one size below our original design. To explore these

branches, we neglect the pressure losses of the pipe segments (7,2) and (7,3) and then

calculate the minimum value of d(6,7) that will satisfy constraint h1 for both consumers 2 and 3. We find that the minimum continuous diameter of d(6,7) is 0.0800. The

next largest discrete diameter is d(6,7) = 0.0825 and we see that this branch has already

been searched. Thus, there are no other feasible branches with discrete values of d(6,7)

less than that of our original design. We have now exhausted all the alternatives and

have found our original discrete design to be the optimal discrete design.

We still have several constraints remaining that must be checked for satisfaction.

The remaining constraints are eq 4-11, 4-22, 4-23, 4-24 and 4-25. These constraints all

deal with the absolute pressure level in the piping. Before we can compute the

absolute pressure at any point, we must first assign an absolute pressure in the

supply pipe at the heating plant. Since we suspect that this will be the place of highest

pressure in the network, we let the absolute pressure at that point be equal to the

maximum allowed, i.e.

 

P P hp,s max 6 2 = =× 1 10 N/m .

We start with eq 4-11, which is a constraint on the maximum pressure in the

supply pipe. The right-hand side of eq 4-11 equals the pressure level in the supply

pipe. As we have shown earlier in this chapter (see eq 4-15), the maximum pressure

must occur at a pipe node and not at an intermediate point. In Table 14, we have

calculated the pressure in the supply pipe at each of the nodes. We see that the

constraint of eq 4-11 is satisfied, since the pressure level does not exceed the

maximum allowed at any point in the supply piping.

Equation 4-22 requires that the pressure at each point in the supply pipe (the lefthand side of the equation) exceed the sum of the saturation pressure Psat and a

safety margin Psaf. For the supply pipe temperature of 120°C, the saturation pressure

is 1.985 × 105 N/m2 (Reynolds and Perkins 1977). Thus, the sum of these two

becomes 2.985 × 105 N/m2. We see by examining Table 14 that this constraint is also

satisfied at all nodes.

We have a similar constraint for the return pipe, eq 4-23. The left-hand side of this

equation equals the pressure in the return pipe, which has also been computed and

is given in Table 14. The temperature and thus the saturation pressure in the return

pipe are different from those in the supply pipe, of course. The return temperature

will vary with load as per our consumer model, eq 3-25. The maximum temperature

will occur at the design condition of maximum load, as can be seen from Figure 5,

and its value is 55°C as determined in Chapter 3 for our supply temperature of 120°C

and our assumptions regarding the radiator design conditions. The saturation

pressure will be greatest at the highest temperature, so if our constraint is satisfied

Table 14. Pressure levels in the piping network.

Node ∆Ps Ps ∆Pr Pr Ps – Pr

number (N/m2) (N/m2) (N/m2) (N/m2) (N/m2)

1 659,272 340,728 150,000 190,728 150,000

2 518,744 481,256 246,183 235,073 246,183

3 449,187 550,813 200,374 350,439 200,374

4 291,257 708,743 231,376 477,367 231,376

5 107,219 892,781 –385,267 394,554 498,228

6 197,874 802,126 –294,546 485,275 316,852

7 218,488 781,512 –273,916 505,905 275,607

8 — 1,000,000 — 287,256 712,744

72
at this condition, it will be satisfied at the lower temperature conditions.

For a return temperature of 55°C, the saturation pressure is 1.576 × 104 N/m2

(Reynolds and Perkins 1977). The sum of this and our safety margin is 1.157 ×

105 N/m2. By examining the return pressures in Table 14, we see the constraint of

eq 4-23 is satisfied at all points in the return piping.

Now we look at the two constraints on the pressure in the return line at the heating

plant, eq 4-24 and 4-25. We see that, for the parameter values chosen, eq 4-24 will

dominate. Equation 4-24 requires that the pressure in the return line at the heating

plant be greater than 2.0 × 105 N/m2. The pressure in the return line at node number

8, which is our heating plant node, is given as 287,256 N/m2 in Table 14, so we see

this constraint is satisfied as well. Thus, our design has satisfied all the constraints

specified.

73
CHAPTER 7: CONCLUSIONS AND RECOMMENDATIONS

SUMMARY

In Chapter 2 we found a suitable method for determining the optimal size for a

single pipe, independent of any others. This method was developed to be as simple

as possible yet complete and accurate enough for design calculations. The method

is general enough to allow for any set of economic or physical parameter values. In

addition, any form of load management, i.e., temperature or flow modulation, or

both, can be accommodated by the integral form of the coefficients in the cost

equation. A new approximation was developed for the friction factor. The form of

this expression was a simple power function of the Reynolds number and the

relative pipe roughness. This form allowed us to easily incorporate it into the head

loss equation without additional complication or rendering the result implicit. We

made use of geometric programming theory to identify a lower bounding problem

that can be used to provide us with a very good first estimate of our solution and a

global lower bound on cost. At the end of Chapter 2, an example is presented that

shows a 17% saving in life cycle cost over a design based on a common rule of thumb.

In Chapter 3 we looked at the heat consumers and the effect that they have on the

piping system. We developed a new model for the consumer’s heat exchanger that

uses the geometric mean temperature difference as an approximation for the

logarithmic mean temperature difference. This allowed us to develop an explicit

expression for return temperature, a result not possible when using the logarithmic

mean temperature difference. We conducted a complete error analysis on the

geometric mean approximation and our new consumer model based on it. This

analysis confirmed that the resulting error from this model was acceptable for

design purposes and much less than the error resulting from using the arithmetic

mean temperature difference as an approximation of the logarithmic mean temperature difference, as has been suggested by others. We integrated our new consumer

model into our single pipe model of Chapter 2 and for a sample case determined

what effect the addition of the consumer has on the integral coefficients of the cost

equation. At the end of Chapter 3, we reworked the example of Chapter 2, including

the effects of the new consumer model.

In Chapter 4 we developed the constraints for systems with multiple pipes and

consumers. Both absolute and differential pressure constraints were derived. By

using the monotonicity of the hydrodynamic and hydrostatic pressure gradients,

we were able to easily show that the maximum pressure within a pipe segment must

occur at one of the end points. We then developed a strategy to allow for constraint

satisfaction at all points implicitly without considering every point in the system.

In Chapter 5 we briefly reviewed general methods for constrained nonlinear

optimization. For various reasons these alternatives are all abandoned in favor of the

approach taken. Subsequently, our general solution strategy is developed for

systems with multiple pipes and consumers. The method makes use of the solution

to the problem, unconstrained by the network constraint requirements, as a starting

point for the constrained solution. Monotonicity analysis was then used to prove

activity of some of the constraints and thus simplify the problem. In addition, the

concept of constraint dominance is used to reduce the number of constraints that

must be considered. Before proceeding with the problem solution, a brief graphical

analysis verified that we only needed to provide for constraint satisfaction at the

maximum load condition to ensure satisfaction at all other load conditions. The

resulting reduced problem was then used as a starting point for two methods

proposed to find a solution to the constrained problem with continuous values for

some of the pipe diameters. Finally, the branch-and-bound technique is explained

and then shown to be suitable for finding a design with discrete values for all the

pipe diameters.

74
In Chapter 6 we worked a simple example with only four consumers and seven

pipe segments. The example illustrated the use of our method and also showed how

the branch-and-bound technique can be used to quickly eliminate candidate designs. A method is also demonstrated for further refinement of the pipe network to

eliminate excessive throttling losses in the consumer’s control valves.

CONCLUSIONS

The method developed here should be feasible for designing the piping networks

for district heating systems of moderate size, and in particular the systems used on

military facilities or college campuses, which tend to be smaller and less complex

than those of large cities. For very large systems, the branch-and-bound method

used for finding the discrete diameters may become cumbersome and

computationally too expensive. However, this remains to be shown and it may be

that, with the commercially available software and the enormous power of today’s

computers, this perceived problem is quite manageable.

The major advantage of the method developed here is its flexibility to accommodate any set of economic and physical parameters and operating strategy. In

addition, the approximations, where used, are much more suitable than some made

in the past: for instance, linearization of the equations, neglecting heat losses, and

oversimplification of the effect of varying load. It is felt that a significant contribution has been made by the derivation of mathematical expressions for all of the major

constraints. Perhaps the most significant contribution of this work has been the

analysis of constraint activity and the development of a method to exploit that

knowledge to arrive at a solution. In addition, we have shown what bounds can be

put on the solution such that the designer can be reasonably assured of whether or

not further significant cost reduction is possible. This not only gives the designer

some comfort in knowing what possible improvement remains, but it also avoids

excessive calculations that often result when no such knowledge is available.

Another possible use of the methodology developed here is for studies of the

relative merits of various operating strategies and what effect they have on the

design of the system. The general form of the cost coefficients can be useful for such

studies and can not only be used to develop designs based on the methods presented

here, but they may also be used to evaluate the effect of economic or physical

parameter changes, including operating strategy changes, on existing designs.

With many systems already in use in Europe, the issue of optimal design is of

lesser importance there. Currently, however, the interest in optimal operation of

district heating systems is significant in Europe, as evidenced by several recent

conferences devoted almost entirely to this topic alone (Nordic Council of Ministers

1989, 1994). Most of the efforts seem to be centered on real-time simulations of

operation and subsequent forecasting of short-term operating strategy. It seems that

a method such as the one developed here would be useful for studies at a higher level

to determine optimal overall operating strategies for the yearly load cycle.

RECOMMENDATIONS

It is recommended that the methodology developed here be field tested on the

design of a moderately sized system, such as would be found on a military base or

college campus. The design should be compared with a completely independent

design, as would be achieved by methods normally used by the district heating

design profession.

Under the assumption that the results of the field test were positive, it is

recommended that the method be coded for computer execution to the maximum

extent possible. The resulting CAD program could then be incorporated into one or

75
more of the currently available CAD systems explicitly created for district heating

system design and feasibility studies. If the economic benefits are as great as

indicated by the examples worked here, the incentive for doing so is significant.

Perhaps the most troubling aspect of using optimal design methods for sizing

district heating piping systems is the level of pressure losses resulting. These

pressure losses are rather high when compared to those encountered in operating

systems and designs based on other methods. This result is apparent when we

examine the pressure losses given in Table 1 for both the method presented here and

a common rule of thumb based design. This result has also been observed by others

(Bøhm 1986, Koskelainen 1980). It seems that current design practice and the

systems that result are ill-suited to the application of optimal design techniques.

Several possible solutions to this conflict exist. The first is simply to increase the

maximum pressure capability of the piping system used. The logic of this approach

can be seen in European practice where small district heating systems use piping

rated for only 6-bar (600-kPa) maximum pressure and often the connections to the

consumers are direct, i.e., without heat exchangers. For larger systems, piping rated

for 15-bar (1500-kPa) or greater maximum pressure is often used and heat exchangers are used to isolate the consumers’ equipment from the high system pressures.

The designer should always find the unconstrained optimum pressure level for the

network first before making a decision of which pressure class of piping to use. It’s

quite possible that pressures higher than those used in current practice may be

justified in some instances.

The advent of friction reducing additives (Nordic Council of Ministers 1991),

which are currently being field tested, offers some relief for the problem of excessive

system pressures. Such additives promise to reduce friction and thus pressure losses

by 50% or more. Such a change in something always assumed to be a basic given in

design renders much of what has been learned to date about district heating system

design, optimal or otherwise, nearly worthless. The ability of the method developed

here to rapidly reevaluate either existing or proposed designs clearly illustrates the

value and necessity of such a design tool. For instance, one possibility is that some

of the frictional reducing additives will be relatively short lived when compared to

the life of the district heating system, or even the annual operating cycle. However,

since the maximum flow rate and thus maximum pressure drop are only encountered over a short period of the yearly cycle, it may be that such friction reducing

additives can be very effective. The method developed here would allow for rapid

evaluation of any possibility to see if they are worthy of further study or consideration for field testing.

76
LITERATURE CITED

Aamot, H. and G. Phetteplace (1976) Long distance heat transmission with steam

and hot water. In Proceedings of the 1st International Total Energy Congress, 4–8 October,

Copenhagen, p. 517–549.

Arora, J.S. (1989) Introduction to Optimum Design. New York: McGraw-Hill Book

Company.

Bøhm, B. (1986) On the optimal temperature level in new district heating networks.

Fernwarme International, 15(5): 301–306.

Bøhm, B. (1988) Energy-economy of Danish district heating systems: A technical

and economic analysis. Lyngby, Denmark: Laboratory of Heating and Air Conditioning, Technical University of Denmark.

COWIconsult (1985) Computerized planning and design of district heating networks. Virum, Denmark: COWIconsult Consulting Engineers and Planners.

CRC Press, Inc. (1987) Standard Mathematical Tables, 28th Edition. Boca Raton,

Florida.

DFF (1985) Dimensionering af stikledninger for fjernvarme. Fjernvarme projectreport

23, DFF (Danske Fjernvarmeværkers Forening), Kolding, Denmark (in Danish with

English summary), ISBN 87-7367-038-3.

Ferziger, J.H. (1981) Numerical Methods for Engineering Application. New York: John

Wiley & Sons.

Forsythe, G.E., M.A. Malcolm and C.B. Moler (1977) Computer Methods for Mathematical Computations. Englewood Cliffs, New Jersey: Prentice-Hall, Inc.

Frederiksen, S. (1982) A thermodynamic analysis of district heating. Doctoral

Dissertation. Lund, Sweden: Lund Institute of Technology (unpublished).

Gartman, H. (1970) DeLaval Engineering Handbook, Third Edition. New York: McGraw-Hill Book Company.

Gill, P.E., W. Murray and M.H. Wright (1981) Practical Optimization. London:

Academic Press, Inc.

Hart, B. and J. Ponsness (1992) HeatMap—A computerized tool for analyzing

district heating and cooling system feasibility. Olympia, Washington: Washington

State Energy Office.

Jeppson, R.W. (1976) Analysis of Flow in Pipe Networks. Ann Arbor, Michigan: Ann

Arbor Science.

Kays, W.M. and A.L. London (1964) Compact Heat Exchangers, Second Edition. New

York: McGraw-Hill Book Company.

Koskelainen, L. (1980) Optimal dimensioning of district heating networks. Fernwarme International, 9(2): 84–90.

Marks, L.S. (1978) Marks’ Standard Handbook for Mechanical Engineers, Eighth Edition

(T. Baumeister, Ed.). New York: McGraw-Hill Book Company.

Marconcini, R. and G. Neri (1979) Numerical simulation of a steam pipeline

network. Geothermics, 7: 17–27.

McDonald, C.L. and C.H. Bloomster (1977) The geocity model: Description and

application. Battelle Pacific Northwest Laboratories, BNWL-SA-6343.

Morofsky, E.L. and V. Verma (1979) DESA analytical manual. Publics Works

Canada, Guidelines Series No. 29.

Munser, H. (1980) Fernwärmeeversorgung. Liepzig: VEB Deutscher für Grundstoffindustrie.

Nordic Council of Ministers (1989) Proceedings, International Symposium on District

Heat Simulation, 13–16 April, Reykjavík, Iceland.

Nordic Council of Ministers (1991) Proceedings, International Symposium on Fluids for

District Heating, 10–11 April, Copenhagen, Denmark.

Nordic Council of Ministers (1994) Proceedings, Symposium and Workshop on Simulation and Operational Optimization of District Heat Systems, 3–4 March, Lyngby,

77
Denmark (B. Bøhm, Ed.).

NRC (National Research Council) (1985) District Heating and Cooling in the United

States, Prospects and Issues. Committee on District Heating and Cooling, Energy

Engineering Board and Building Research Board, Commission on Engineering and

Technical Systems. Washington, D.C.: National Academy Press.

Papalambros, P.Y. and D.J. Wilde (1988) Principles of Optimal Design, Modeling and

Computation. Cambridge: Cambridge University Press.

Phetteplace, G.E. (1981) Transient analysis of heat transmission systems. USA Cold

Regions Research and Engineering Laboratory, CRREL Report 81-24.

Phetteplace, G.E., W. Willey and M. Novick (1981) Heat losses from the Fort

Wainwright heat distribution system. USA Cold Regions Research and Engineering

Laboratory, Special Report 81-14

Phetteplace, G.E. and V.P. Meyer (1990) Piping for thermal distribution systems.

USA Cold Regions Research and Engineering Laboratory, Internal Report 1059

(unpublished).

Rao, S.S. (1984) Optimization: Theory and Applications, Second Edition. New Delhi:

Wiley Eastern Limited.

Rasmussen, C.H. and J.E. Lund (undated) Computer aided design of district

heating systems. From District Heating Research and Technological Development

in Denmark, Danish Ministry of Energy.

Reisman, A.W. (1985) District Heating Handbook. Vol. 2, A Handbook of District

Heating and Cooling Models. Washington, D.C.: International District Heating and

Cooling Association.

Reklaitis, G.V., A. Ravindran and K.M. Ragsdell (1983) Engineering Optimization:

Methods and Applications. New York: Wiley-Interscience, John Wiley & Sons, Inc.

Reynolds, W.C. and H.C. Perkins (1977) Engineering Thermodynamics, Second

Edition. New York: McGraw-Hill Book Company.

R. S. Means Co., Inc. (1987) Mechanical Cost Data, 11th Annual Edition. Kingston,

Massachusetts.

Segan, E.G. and C. Chen (1984) Investigation of tri-service heat distribution

systems. U.S. Army Construction Engineering Research Laboratory, Technical

Report M-347.

Soumerai, H. (1987) Practical Thermodynamic Tools for Heat Exchanger Design Engineers. New York: John Wiley and Sons.

Stoner, M. (1974) Modeling the steady state pressure-flow response of steam and

water systems. In Proceedings of the Annual Conference of the International District

Heating Association, p. 53–69.

Szepe, S. and J. Calm (1979) Thermal fluid selection for long-distance heat transmission. In Proceedings of the 14th Intersociety Energy Conversion Engineering Conference,

5–10 August, Boston, Massachusetts, p. 1985–1990.

Vanderplaats, G.N. (1984) Numerical Optimization Techniques for Engineering Design:

With Applications. New York: McGraw-Hill Book Company.

Werner, S.E. (1984) The Heat Load in District Heating Systems. Goteborg, Sweden:

Chalmers University of Technology.

Wilde, D.J. and C.S. Beightler (1967) Foundations of Optimization. Englewood Cliffs,

New Jersey: Prentice-Hall, Inc.

Wood, D. J. (1966) An explicit friction factor relationship. Civil Engineering, 36(12).

Zinger, N.M., K.S. Andreeva, E.M. Fedorov and T.B. Oshchepkova (1976) Hydraulic design of a branched heating network, employing an electronic digital computer.

Teploenergetika, 23(11): 20–24.

78
APPENDIX A: APPROXIMATION OF THE FRICTION FACTOR

For incompressible flow in circular conduits (pipes), the head losses can be

calculated using the Darcy-Weisbach equation

hf = f L v2/2 g d (A-1)

where f = friction factor (dimensionless)

L = pipe length (m)

v = flow velocity (m/s)

g = acceleration of gravity (9.8 m/s2)

d = inside diameter of the pipe (m).

Numerous other expressions have been proposed for calculating frictional head

losses. The Darcy-Weisbach equation is, however, the most fundamentally appealing as it can be derived analytically while the other relationships are empirical

in nature (Jeppson 1976). For laminar flow it is possible to show that the friction

factor f is a function of the Reynolds number alone. Unfortunately, the flow in heat

distribution piping is seldom laminar. For turbulent flow the friction factor has been

determined empirically to be a function of the Reynolds number and the relative

roughness of the pipe. A number of correlations have been proposed for the friction

factor. Those correlations that give the best agreement with the experimental data

are implicit in the friction factor. This renders them impractical for analyses such as

this one. For this analysis, and other applications, it would be desirable to have a

simple expression that would provide sufficient accuracy over limited ranges of

interest. To keep the expression as simple as possible, while allowing it to be an

accurate approximation, a method is developed here that yields a one-term power

function.

To approximate friction factor information in the form of implicit equations or

empirical data, we can develop our approximation using the least-squares method.

First, we assume a desired form for our expression for the friction factor

 

ˆ

f a RR Re b c = (A-2)

where 

ˆ

f = predicted friction factor (dimensionless)

a, b, and c = coefficients determined by the least-squares method (dimensionless)

RR = ε/d = relative roughness of the pipe (dimensionless)

ε = absolute roughness of the piping (m)

Re = Reynolds number for the pipe flow (dimensionless).

If we assume that for any set of values for RR and Re we have an observed friction

factor f, we would like to minimize the sum of the squares between the f ’s predicted

by our equation and all the observed f ’s within the range of interest

 

min ˆ ∑( ) f f − 2

. (A-3)

The summation is taken over all the observations available within the range of

interest for the parameters RR and Re. In the event that we are trying to approximate

an implicit empirical expression, we would choose incremental values of RR and Re

over the range of interest and use these to calculate a corresponding f value. This

approach will be illustrated later in this appendix. To accomplish the minimization,

we first convert to a linear form by making the following substitutions

Y = ln f (A-4)

79
ˆ

Y = ln ˆf = ln (a RRbRec

) = ln a + b ln RR + c ln Re

= β0+ β1X1+ β2X2 (A-5)

where β0 = ln a

β1 = b

β2 = c

X1 = ln RR

X2 = ln Re.

Now we can restate the problem in a linear form

 

min YY Y X X ˆ

i

n

i

n

∑( ) i i i ,i ,i − = −− − ∑( ) = = 1

2

0 11 22

1

2

ββ β . (A-6)

The summation index i has now been added. The summation occurs over the total

number of observations n. To find the minimum for this expression with respect to

the parameters β0, β1 and β2, we take the partial derivative of the expression with

respect to each of the parameters and set the result to zero in each case

 

∂ ∂β β β β / 0 0 11 22

1

= −− − ∑( ) = 0 =Y XX

i

n

i ,i ,i (A-7)

 

∂ ∂β β β β / 1 1

1

= ∑ ( ) −− − 0 11 22 = 0 =XY X X

i

n

,i i ,i ,i

(A-8)

 

∂ ∂β β β β / 2

1

= ∑ ( ) −− − 0 11 22 = 0 =XY X X

i

n

2,i i ,i ,i .

(A-9)

We now have three linear equations in the three unknown parameters by rearranging as follows

 

ββ β 0 + + 1 ,i 2,i i ∑ ∑ = ∑ = ==

X XY

i

n

i

n

i

n

1

1

2

1 1

(A-10)

 

ββ β 0 1,i 1,i 2 2,i 1,i i 1,i X X X X YX

i

n

i

n

i

n

i

n

∑ ∑ + ∑ + = ∑ = == = 1

1

2

1 11

( )

(A-11)

 

ββ β 0 2,i 1,i 2,i 2 2,i i 2,i X X X X YX

i

n

i

n

i

n

i

n

∑ ∑ + ∑ + = ∑ = == = 1

1 1

2

1 1

( ) () .

(A-12)

These equations can be written in matrix form as

 

AAA

AAA

AAA

C

C

C

11 12 13

21 22 23

31 32 33

0

1

2

1

2

3

β

β

β

= (A-13)

where A11 = n

 

AA X

i

n

12 21 1,i = = ∑

=1

 

AA X

i

n

13 31 2,i = = ∑

=1

 

A X

i

n

22 1,i = ∑

=

( )

1

2

80
Table A1. Constants for the friction factor equation.

Water Flow Pipe Reynolds

temp. velocity diameter number Max. Avg.

(°C) (m/s) (m) × 10–6 error error

min/max min/max min/max min/max abc (%) (%)

50/130 0.5/3.3 0.05/0.77 0.04/11 0.123 0.146 –0.0626 6.2 1.0

50/130 0.5/4.5 0.05/0.77 0.04/15 0.119 0.152 –0.0568 6.9 1.1

70/150 0.3/6.3 0.03/0.93 0.02/29 0.129 0.156 –0.0589 10.8 2.0

50/90 0.5/2.9 0.05/0.41 0.04/3.7 0.140 0.141 –0.0762 4.1 0.8

90/130 0.5/2.9 0.10/0.46 0.16/5.9 0.116 0.150 –0.0563 2.5 0.6

50/90 0.5/2.9 0.10/0.46 0.09/4.1 0.128 0.132 –0.0751 3.4 0.7

50/90 0.5/3.7 0.10/0.46 0.09/5.3 0.125 0.137 –0.0698 3.9 0.8

90/130 0.5/3.7 0.10/0.46 0.16/7.5 0.113 0.154 –0.0520 2.8 0.6

 

A A XX

i

n

23 32 2,i 1,i = = ∑

=1

 

A X

i

n

33 2,i = ∑

=

( )

1

2

 

C Y

i

n

1 i = ∑

=1

 

C YX

i

n

2 i 1,i = ∑

=1

 

C YX

i

n

3 i 2,i = ∑

=1

.

This system of linear equations can be solved by forward elimination and

subsequent back solution. The resulting expressions for the parameters are

 

β2

3 1 31 11 2 1 21 11

32 12 31 11

22 12

2

11

33 13

2

11 23 13 21 11

32 12 31

=

{ } − ( ) − − [ ] ( ) ( ) − ( )

( ) − ( )









































{ } − ( ) − − [ ] ( ) −

C CA A C CA A

A AA A

A AA

A A A A AA A

A AA A

/ / /

/

/ / / 11

22 12

2

11

( ) ( )

( ) − ( )





































 A AA/ 

 

(A14)

 

β β

1

2 1 21 11 23 13 21 11 2

22 12

2

11

= { } − ( ) − − { } [ ] ( ) [ ]

{ } − ( )

C CA A A A A A

A AA

/ /

/

(A15)

 

β ββ 0 1 12 1 13 2 11 = − { } CA A A [ ] − [ ] / . (A16)

A FORTRAN program FFCONST was written to evaluate the A’s and C’s in the

above expressions and then solve eq A14, A15 and A16 to find the parameters β0, β1

and β2. The parameters in our original expression, eq A2, for the predicted friction

factor 

ˆ

f can then be found. For this program, the “observed” friction factor f is found

using the Colebrook-White equation (Jeppson 1976)

f = [1.14 – 0.869 ln(RR + 9.35/Re f )]−2 . (A-17)

The Colebrook-White equation is implicit in the friction factor f and thus it cannot

be solved directly. A number of methods can be used to solve implicit equations such

81
as this one. The Newton-Raphson method has been used in SUBROUTINE CWFF.

This method uses knowledge of the first derivative of the function to find the

solution. A description of the method can be found in nearly any reference on

numerical methods, such as Forsythe et al. (1977).

To use the Newton-Raphson method to solve the implicit Colebrook-White

equation, an initial estimate of the f value is needed. An explicit equation for the

friction factor f can be used for this. The explicit equation does not need to be

extremely accurate to yield a suitable first estimate. The equation given by Wood

(1966) is a good explicit relationship for turbulent flow and can be used. Wood’s

equation is

 

f = ++ ( ) 



























0 094 0 53 88 0 225 0 44 1 62 0 134

.. / . . . .RR RR RR Re

RR . (A-18)

To calculate the friction factor using either eq A17 or A18 requires that we know

the Reynolds number Re and the relative roughness RR. To calculate these parameters, we need to specify the fluid density and dynamic viscosity as well as the pipe

diameter and absolute roughness and the flow velocity. The fluid properties are a

function of the temperature of the fluid and to a lesser extent its pressure as well.

Here, we will assume that the fluid is at its saturation pressure for the temperature

specified. Two FORTRAN subroutines were written to determine the fluid properties. The first, SUBROUTINE SATLN, calculates the saturation pressure for water

given the temperature. The second, SUBROUTINE WTRTBL, calculates the density

and dynamic viscosity given the temperature and pressure. The main program

FFCONST and each of the subroutines mentioned above are included in Appendix

B.

Using the program FFCONST, a number of the constants a, b and c were

determined for several sets of parameters. Table A1 summarizes the results. In each

of the examples of Table A1 the absolute roughness of the pipe was taken as 4.6 ×

10–5 m.

82
APPENDIX B: COMPUTER PROGRAM LISTINGS

Program FFCONST

PROGRAM FFCONST

DIMENSION VI(10),DI(10),FI(10,10,10),FFCALC(10,10,10),

*ERROR(10,10,10),RNI(10,10,10),T(10)

100 FORMAT(15X,3E12.5)

200 FORMAT(////,15X,4E15.6,///,15X,4E15.6)

300 FORMAT(////,20X,3F18.6,/////)

400 FORMAT(2X,1F10.0,1F10.2, 1F10.3,1E15.4,3F10.4)

500 FORMAT(2X,4E12.4)

600 FORMAT(////,15X,3F18.6)

R=.046E-3

ESUM=0.

A1=0.

A2=0.

A3=0.

A4=0.

A5=0.

A6=0.

A7=0.

A8=0.

A9=0.

A10=0.

A11=0.

A12=0.

A13=0.

A14=0.

A15=0.

A16=0.

EMAX=0.

EMIN=1.

JV=9

ID=9

DELTAD=.025

DL=.025

VL=.25

VDELTA=.75

VI(1)=VL

KT=9

TDELTA=10.

TL=70.

T(1)=TL

N=JV*ID*KT

DO 4 K=1,KT

CALL SATLN(T(K),P)

CALL WTRTBL(T(K),P,RHO,XH,DV)

* PRINT 500,T(K),P,RHO,DV

83
DO 3 J=1,JV

DI(1)=DL

DO 2 I=1,ID

RR=R/DI(I)

RNI(I,J,K)=DI(I)*VI(J)*RHO/DV

RN=RNI(I,J,K)

F=.094*(RR**.225)+(.53*RR)+(88.*(RR**.44))/(RN**(1.62*(RR**

*.134)))

FI(I,J,K)=F

A1=A1+LOG(F)

A2=A2+LOG(RR)

A3=A3+LOG(RN)

A4=A4+LOG(F)*LOG(RR)

A5=A5+((LOG(RR))*(LOG(RR)))

A6=A6+LOG(RN)*LOG(RR)

A7=A7+LOG(RN)*LOG(F)

A8=A8+((LOG(RN))**2.)

* PRINT 100,RNI(J),DI(I),FI(I,J)

DI(I+1)=DI(I)+(DELTAD*I)

2 CONTINUE

VI(J+1)=VI(J)+VDELTA

3 CONTINUE

T(K+1)=T(K)+TDELTA

4 CONTINUE

A9=A6-(A2*A3/N)

A10=A5-(A2*A2/N)

A11=A6-(A2*A3/N)

A12=A4-(A1*A2/N)

A13=A8-(A3*A3/N)

A14=A7-(A1*A3/N)

A15=A13-(A11*A9/A10)

A16=A14-(A12*A9/A10)

C=A16/A15

B=A12/A10-C*A11/A10

ALN=A1/N-C*A3/N-B*A2/N

A=EXP(ALN)

* PRINT 200,A1,A2,A3,A4,A5,A6,A7,A8

PRINT 300,A,B,C

RN=RNL

DO 7 K=1,KT

DO 6 J=1,JV

D=DL

DO 5 I=1,ID

RN=RNI(I,J,K)

FFCALC(I,J,K)=A*((R/D)**B)*(RN* *C)

ERROR(I,J,K)=(FFCALC(I,J,K)-FI(I,J,K))/FI(I,J,K)

* PRINT 400,T(K),VI(J),D,RN,FI(I,J,K),FFCALC(I,J,K),ERROR(I,J,K)

D=D+(DELTAD*I)

EABS=ABS(ERROR(I,J,K))

84
ESUM=ESUM+EABS

EMAX=MAX(EMAX,EABS)

EMIN=MIN(EMIN,EABS)

5 CONTINUE

6 CONTINUE

7 CONTINUE

EAVG=ESUM/N

PRINT 600,EMIN,EMAX,EAVG

STOP

END

Subroutine WTRTBL(T,P,RHO,XH,DV)

SUBROUTINE WTRTBL(T,P,RHO,XH,DV)

* THIS SUBROUTINE CALCULATES THE THERMODYNAMIC AND

* TRANSPORT PROPERTIES OF WATER GIVEN THE TEMPERATURE AND

* PRESSURE CONDITIONS.

* INPUT VALUES ARE T(C) AND P(BARS). THE OUTPUTS ARE DENSITY

* RHO (KG/M3), THE ENTHALPY XH(KJ/KG), AND THE DYNAMIC

* VISCOSITY DV IN KG/M-SEC. THE EQUATIONS FOR THE VISCOSITY ARE

* NOT VALID FOR TEMPERATURES GREATER THAN 300 C.

T=T+273.16

* CALL SATLN(T,P)

DATA RI1,B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,Ul ,Wl/22129.,-37444.8692,

*466453.368,-2666876.77,9030271.53,- 19769400.2,28949239.9,

*-28309932.7,17808942.6,-6534676.01,1065198.53,.58620689,

*.41666667/

DATA G,H,RK,RL,RM,RN,Fl,Gl ,Hl,RKl,RLl,RMl,RNl,Ql,Rl,Zl/

*.417,1.139706E-4,9.949927E-5,7.241165E-5,.7676621,1.052358E- 11,

*3.7E8,3.122199E8,199985.,1.72,1.362926E16,1.500705,.6537154,

*62.5,13.10268,1.5108E-5/

DATA A2,A3,A4,A5/.3828209486,.2162830218,.1498693949,.4711880117/

P1=221.287

T1=647.3

T=T/T1

P=P/P1

U=F1-(G1*T*T)-(H1*(T**(-6.)))

W=U+SORT((RK1*U*U)+(RL1*(P-(RM1*T))))

V1=(G/(W**(1./3.4)))-H+(RK*T)

V2=((ABS(RN1-T))**2.)*(RL+(((ABS(RN1-T))**8.)*RM))

V3=(RN*(R1+(R1*P)+(P*P)))/(Z1+(T**11.))

RHO=1./(V1+V2-V3)

H0=B0+(B1*T)+(B2*T*T)+(B3*T**3.)+(B4*T**4.)+(B5*T**5.)+

*(B6*T**6.)+(B7*T**7.)+(B8*T**8.)+(B9*T**9.)

Y1=RL1*RM1/2.

Q=2.*G/RL1

V=(-2.*G1*T*T)+(6.*H1/(T**6.))

H1=((U1 * W)-(W1*((3.4*U)-V)))*W

H1=Q/(W**(1./3.4))*(H1+(Y1*T)-(.72*V*U))

H2=(((RN1-T)*((RL*(RN1+T))+(RM*((ABS(RN1-T))* *8.)*(RN1+(9.*T)))))

85
*-H)*P

H3=(RN*(Z1+(12.*(T**11.))))/((Z1+(T**11.))**2.)

H3=H3*(Q1+(((R1/2.)+(P/3.))*P))*P

XH=H0+(RI1*(H1+H2-H3))

TC=(T*T1)-273.16

CALL SATLN(TC,PS)

PS=PS/P1

DV=.02414+(10.**(A2/(T-A3)))*(1.+((P-PS)*A4*(T-A5)))

T=T*T1

P=P*P1

DV=DV*1.E-3

T=T-273.16

RETURN

END

Subroutine SATLN(T,P)

SUBROUTINE SATLN(T,P)

* THIS SUBROUTINE CALCULATES THE SATURATION PRESSURE FOR A

* GIVEN WATER TEMPERATURE. THE WATER TEMPERATURE IS IN C

AND

* THE SATURATION PRESSURE IS RETURNED IN BARS.

DATA RK,A,B,C,D,E,F/2.937E5,5.426651,-2005.1,1.3869E-4,1.1965E-11,

*-.0044,-.0057148/

T=T+273.16

T2=T+.01

Y=647.26-T

X=(T2*T2)-RK

T=T/647.3

A1=A+(B/T2)+(C*X/T2*(10.**(D*X*X)-1.))+(E*(10.**(F*(Y* *1.25))))

P=(1.01325*(10.**A1))+((T-.422)*(.577-T)*EXP(-12.+(T**4.))*

*9.80665E-3)

T=T*647.3

T=T-273.16

RETURN

END

Subroutine ROMBRG(FUN,A,B,C,ERR,RES)

SUBROUTINE ROMBRG(FUN,A,B,C,ERR,RES)

C THIS SUBROUTINE COMPUTES INTEGRALS OF A USER SUPPLIED

C FUNCTION USING ROMBERG’S METHOD. THIS SUBROUTINE IS FROM

C “NUMERICAL METHODS FOR ENGINEERING APPLICATION”,J.H.

C FERZIGER, JOHN WILEY AND SONS, 1981. THE ARGUMENT ARE:

C FUN = THE FUNCTION TO BE INTEGRATED

C A = LOWER LIMIT

C B = UPPER LIMIT

C C = ARRAY OF FUNCTION DEFINITION PARAMETERS IF REQUIRED

C ERR = THE DESIRED ACCURACY

C RES = THE RESULTING VALUE FOR THE INTEGRAL

86
EXTERNAL FUN

C Z IS THE ARRAY OF APPROXIMATIONS

DIMENSION Z(10,10)

C INITIALIZE AND COMPUTE THE FIRST APPROXIMATION.

I=1

DEL=B-A

Z(1, 1)=.5*DEL*(FUN(A,C)+FUN(B,C))

C THE MAIN LOOP. THE FIRST PART COMPUTES THE INTEGRAL USING A

C 2J+1 POINT TRAPEZOID RULE. THE METHODS MAKES MAXIMAL USE

C OF THE VALUES ALREADY COMPUTED.

10 J=2**(I-1)

DEL=DEL/2.

I=I+1

Z(I,1)=.5*Z(I-l,l)

DO 1 K=1,J

X=A+(2.*K-1)*DEL

Z(I,1)=Z(I,1)+DEL*FUN(X,C)

1 CONTINUE

C NOW WE DO THE RICHARDSION EXTRAPOLATION

DO 2 K=2,I

Z(I,K)=(4.**(K-1)*Z(I,K-1)-Z(I-1,K-1))/(4.**(K-1)-1.)

2 CONTINUE

C ERROR CONTROL

DIFF=ABS(Z(I,I)-Z(I,I-1))

PRINT 108, DIFF,Z(I,I)

108 FORMAT(5X,2E15.5)

IF(DIFF.LT.ERR) GO TO 20

C THE MAXIMUM NUMBER OF ITERATIONS ALLOWED IS 10.

IF(I.LT.10) GO TO 10

PRINT 100

100 FORMAT(‘ MORE THAN lO ITERATIONS REQUIRED, CHECK

PARAMETERS.’)

STOP

20 RES=Z(I,I)

RETURN

END

Function FUN(T,Z)

FUNCTION FUN(T,Z)

C THIS FUNCTION CONTAINS THE INTEGRAND OF THE I2 PARAMETER.

A=0.135

B=0.161

C=-0.0555

PVFE=9.077

PVFH=9.077

EPS=5.E-5

PL=1000.

FRD=100.

87
CE=7.E-5

CH=3.4E-5

AEDA=0.9

TSD=120.

TRD=60.

CALL WTRTBL(TRD,PRD,RHORD,RHD,RMUD)

CALL WTRTBL(TSD,PSD,RHOSD,SHD,SMUD)

RHOD=(RHOSD+RHORD)/2.

TS=120.

TA=20.

A12=A*PVFE*PL*(EPS**B)*(1.273240**(2+C))

QF=(0.425*COS(6.283*T/8760.))+0.575

CALL CLMTD(TS,TA,QF,FRF,TR,TRG,FRFG)

FT=FRF

C CALCULATE QUANITIES WHICH MAY BE A FUNCTIONS OF TIME.

CALL WTRTBL(TS,PS,RHOS,SH,SMU)

CALL WTRTBL(TR,PR,RHOR,RH,RMU)

A7=(((1./(SMU**C))/(RHOS*RHOS))+((1./(RMU**C))/(RHOR*RHOR)))/2.

RHOA=(RHOS+RHOR)/2.

FUN=A7*((((CE*RHOA)/(AEDA*RHOD))*(FT* *(2+C)))-((PVFH/PVFE)

&*CH*(FT**(3+C))))

FUN=FUN*A12*(FRD**(3+C))

PRINT 102, T,FT,FUN

102 FORMAT(5X,3E 15.4)

RETURN

END

Program I2

PROGRAM I2

EXTERNAL FUN

CALL ROMBRG(FUN,0.,4380.,0.,1.E-7,PI2)

PRINT 101,PI2+2.

101 FORMAT(l0X,F15.5)

STOP

END

Program I2-C-GMT

PROGRAM I2-C-GMT

SUBROUTINE ROMBRG(FUN,A,B,C,ERR,RES)

C THIS SUBROUTINE COMPUTES INTEGRALS OF A USER SUPPLIED

FUNCTION

C USING ROMBERG’S METHOD. THIS SUBROUTINE IS FROM

“NUMERICAL

C METHODS FOR ENGINEERING APPLICATION”,J.H. FERZIGER, JOHN

WILEY

C AND SONS, 1981. THE ARGUMENTS ARE:

C FUN = THE FUNCTION TO BE INTEGRATED

C A = LOWER LIMIT

88
C B = UPPER LIMIT

C C = ARRAY OF FUNCTION DEFINITION PARAMETERS IF REQUIRED

C ERR = THE DESIRED ACCURACY

C RES = THE RESULTING VALUE FOR THE INTEGRAL

EXTERNAL FUN

C Z IS THE ARRAY OF APPROXIMATIONS

DIMENSION Z(10,10)

C INITIALIZE AND COMPUTE THE FIRST APPROXIMATION.

I=1

DEL=B-A

Z(1,1)=.5*DEL*(FUN(A,C)+FUN(B,C))

C THE MAIN LOOP. THE FIRST PART COMPUTES THE INTEGRAL USING A

C 2J+1 POINT TRAPEZOID RULE. THE METHODS MAKES MAXIMAL USE

OF THE

C VALUES ALREADY COMPUTED.

10 J=2**(I-1)

DEL=DEL/2.

I=I+1

Z(I,1)=.5*Z(I-l,l)

DO 1 K=1,J

X=A+(2.*K-1)*DEL

Z(I,l)=Z(I,l)+DEL*FUN(X,C)

1 CONTINUE

C NOW WE DO THE RICHARDSION EXTRAPOLATION

DO 2 K=2,I

Z(I,K)=(4.**(K-1)*Z(I,K-l)-Z(I-l,K-1))/(4.**(K-1)-1.)

2 CONTINUE

C ERROR CONTROL

DIFF=ABS(Z(I,I)-Z(I,I-1))

PRINT 108, DIFF,Z(I,I)

108 FORMAT(5X,2E15.5)

IF(DIFF.LT.ERR) GO TO 20

C THE MAXIMUM NUMBER OF ITERATIONS ALLOWED IS 10.

IF(I.LT.10) GO TO 10

PRINT 100

100 FORMAT(‘ MORE THAN 10 ITERATIONS REQUIRED, CHECK

PARAMETERS.’)

STOP

20 RES=Z(I,I)

RETURN

END

Function FUN(T,C)

FUNCTION FUN(T,C)

C THIS FUNCTION CONTAINS THE INTEGRAND OF THE I2 PARAMETER.

C It has been modified to include the effect of the consumers model

C using the GMTD model. See EQ. 3.25 for explanation of symbols

C used in flow rate equation.

89
TSD=120.

TRD=55.

TA=20.

TS=120.

TGMTD=59.1608

RN1=1.3

A13=1.0

PL=1000.

DMF=4.

TF=(0.425*COS(6.283 *T/8760.))+0.575

TOP=(TSD-TRD)*TF

BOTTOM=TS-TA-((1/(TS-TA))*(TGMTD*TGMTD)*((TF/A13)**(2/RN1)))

A=TOP/BOTTOM

C “A” is the normalized mass flow rate, m/md.

FUN=(2.043E-8*DMF*((DMF*A)**1.9432))-(8.924E-9*((A*DMF)**2.9432))

FUN=FUN*PL/1000

PRINT 102, T,A,FUN

102 FORMAT(5X,3E15.4)

RETURN

END

Program I2

PROGRAM I2

EXTERNAL FUN

CALL ROMBRG(FUN,0.,4380.,0.,1.E-9,PI2)

PRINT 101,PI2*2.

101 FORMAT(l0X,E15.5)

STOP

END

Program I1EQ3-26

PROGRAM I1EQ3-26

FUNCTION FUN1(T,Z)

C THIS FUNCTION CONTAINS THE INTEGRAND OF THE I1 PARAMETER.

C It uses equation 3.26 which was derived using the GMTD appx.

C Modified from “IlFT.for” on 1/12/94.

RN1=1.3

A14=0.575

A15=0.425

A13=1.0

FUN1=(((A15*COS(6.283*T/8760.))+A14)/A13)**(2/RN1)

PRINT 102, T,FUN1

102 FORMAT(5X,2E15.4)

RETURN

END

90
Program I1

PROGRAM I1

EXTERNAL FUN1

CALL ROMBRG(FUN1,0.,4380.,0.,1.E-7,PI1)

PVFH=9.077

PL=1000.

CH=3.4E-5

TS=120.

TM=6.4

TA=20.

TCINS=0.030

AT=8760

TGMTD=59.1608

A16=PVFH*PL*12.56637*TCINS*CH

RI1=A16*(((PI1*2.)*(TGMTD *2/(2*(TS-TA))))+((((TS+TA)/2)-TM)*AT))

PRINT 101,RI1

101 FORMAT(10X,E15.7)

STOP

END

91
1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHORS

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

 REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS 15. NUMBER OF PAGES

16. PRICE CODE

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT

 OF REPORT OF THIS PAGE OF ABSTRACT

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and

maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,

including suggestion for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington,

VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

August 1995

Optimal Design of Piping Systems for District Heating

Gary Phetteplace

U.S. Army Cold Regions Research and Engineering Laboratory

72 Lyme Road

Hanover, N.H. 03755-1290

CRREL Report 95-17

Approved for public release; distribution is unlimited.

Available from NTIS, Springfield, Virginia 22161

103

UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UL

Applied optimization Heat distribution systems

Branch-and-bound method Heat transfer

District heating Optimal design

First, a method for determining the optimal size for a single pipe segment in a district heating system is developed.

The method is general enough to allow for any set of economic or physical parameter values. In addition, any form

of load management, i.e., temperature or flow modulation, or both, can be accommodated by the integral form of

the coefficients in the cost equation. An example is presented that shows a 17% savings in life cycle costs over a design

based on a common rule of thumb. Next the heat consumer and his effects on the piping system are studied. A new

model is developed for the consumer’s heat exchanger that uses the geometric mean temperature difference as an

approximation for the logarithmic mean temperature difference. The new consumer model is integrated into the

previous single pipe model and, for a sample case, its effect is determined. For systems having multiple pipes and

consumers, the constraints are first developed and then the general solution strategy. The method makes use of the

solution to the unconstrained problem as a starting point for the constrained solution. Monotonicity analysis is then

used to prove activity of some of the constraints, and thus simplify the problem. Finally, the branch-and-bound

technique is shown to be suitable for finding a design with discrete values for all the pipe diameters. A simple

example is provided. In addition, a method is also demonstrated for further refinement of the pipe network to

eliminate excessive throttling losses in the consumer’s control valves. The method developed here should be feasible

for designing the piping networks for district heating systems of moderate size, and its major advantage is its

flexibility.