Determining the optimum economic insulation thickness of double pipes buried in the soil for district heating systems

Abstract

The insulation thickness (IT) of double pipes buried in the soil (DPBIS) for district heating (DH) systems was optimized to minimize the annual total cost of DPBIS for DH systems. An optimization model to obtain the optimum insulation thickness (OIT) and minimum annual total cost (MATC) of DPBIS for DH systems was established. The zero point theorem and fsolve function were used to solve the optimization model. Three types of heat sources, four operating strategies, three kinds of insulation materials, seven nominal pipe size (NPS) values, and three buried depth (BD) values were considered in the calculation of the OIT and MATC of DPBIS for DH systems, respectively. The optimization results for the above factors were compared. The results show that the OIT and MATC of DPBIS for DH systems can be obtained by using the optimization model. Sensitivity analysis was conducted to investigate the impact of some economic parameters, i.e., unit heating cost, insulation material price, interest rate, and insulation material lifetime, on optimization results. It is found out that the impact of sensitivity factors on the OIT and MATC of DPBIS for DH systems is different.

Keywords

double pipes / optimization model / optimum insulation thickness / minimum annual total cost

Introduction

DH systems have been developed widely because their application is considered as an effective method to improve the energy efficiency of space heating in buildings. DH systems are very common in many countries, such as China, Russia, European countries, etc. In Sweden, the energy used in the building sector accounts for 40% of the total energy. 55% of the heating demand of Swedish buildings is met by DH systems [¹,²]. In northern heating regions of China, the area of buildings connected to DH systems increased from 5 billion m² in 2001 to 12 billion m² in 2013 [³]. Energy conservation is an important part of China’s national energy strategy. In general, energy consumption can be examined in four primary sectors, i.e., industry, construction, transportation, and agriculture. In the EU countries and the US, the building sector is one of the largest energy consumers, accounting for about 40% of the total energy consumption [⁴–⁸]. In northern heating regions of China, the heating energy consumption was 181 million tons of standard coal equivalent in 2013, accounting for 24% of the building energy consumption [³]. However, the application of insulated heating pipes is an effective way to reduce the energy consumption caused by DH systems [⁹–¹⁸].

DH systems mainly consist of heat sources, piping networks, and consumers. Choosing the reasonable laying method of the DH piping networks is of great significance for saving investment, ensuring the reliable operation of DH systems, facilitating the construction and maintenance of the DH piping networks, etc. [¹⁹]. There are two laying methods of the DH piping networks, overhead laying and underground laying. Underground laying does not affect the appearance and traffic of cities. Therefore, it is widely used for the DH piping networks in towns. Besides, underground laying is divided into trench laying and buried laying. Buried laying is quite common for the DH piping networks, especially in China [¹⁹]. Generally speaking, the heat loss accounts for about 10%–30% of total heat supply [²⁰,²¹]. However, for some old DH piping networks in China, the proportion of heat loss to heat supply is up to 30% [³]. In Western and Northern European countries, the typical heat loss accounts for 8%–15% of the heat supply [²²,²³]. Due to the lower linear heat density and lower insulation quality, the corresponding proportion in Eastern European countries is about 15%–25% [²²].

In Refs. [²⁴,²⁵], the status of the 4th generation district heating (4GDH) was investigated. The 4GDH systems have advantages over previous DH systems in heat loss of heat distribution pipes, utilization of low-temperature heat sources, efficiency in the heat production, etc. The focus of transformation of the current DH systems to 4GDH systems is the improvement of DH piping networks. One of the measures to improve the DH piping networks is to reduce the heat loss. As adding insulation materials to heat distribution pipes is essential to reducing the heat loss, it is increasingly necessary to study the IT of DH pipes.

Heat loss cost can be reduced when the IT of DH pipes increases. However, the investment cost of insulation materials increases with the increase in IT. It can be seen that there must be an OIT to appropriately balance the heat loss cost and investment cost of the insulation materials. Therefore, it is absolutely necessary to study the OIT of buried DH pipes. In China, the OIT of heat distribution pipes for DH systems is recommended in Ref. [²⁶]. However, these values are always obtained based on the empirical data. Besides, with the development of economy, these data become unreasonable. In Ref. [²²], an optimization model on insulated DH pipes above ground was presented. The OIT was determined by minimizing the annual total cost. Zhang et al. [²⁷] investigated the economic and energetic performance of insulated DH pipes buried in the soil. The life cycle cost analysis (LCCA) method and the heating degree-days (HDD) method were used to perform the thermo-economic evaluation for optimizing the IT. In some studies [²⁸–³⁹], the HDD method and the LCCA method were used to calculate the annual heat loss of insulated DH pipes above ground. In Refs. [²⁹–³¹], environmental impact was also considered in optimizing the IT of insulated DH pipes above ground. Daşdemir et al. [³²] reported the effect of air gap on the IT and economic cost for DH pipes with different nominal pipe sizes. In Ref. [³³], the artificial neural network was used to evaluate the life cycle cost and OIT of insulated DH pipes above ground. Daşdemir et al. [³⁴] aimed to optimize the IT of insulated DH pipes made of different materials (i.e., steel, plastic, and copper). In Ref. [³⁵], the OIT of insulated DH pipes above ground, energy savings over a lifetime of 10 years, and payback periods were calculated for five different insulation materials (i.e., foam board, XPS, rock wool, EPS, and fiberglass), respectively.

In recent years, exergy analysis method has become more popular [²⁸,³⁶–³⁹]. The combination of exergy analysis [²⁸,³⁶,³⁹], LCCA [²⁸,³⁶], and life cycle assessment [³⁹] was used to optimize the IT of insulated DH pipes above ground. However, few studies are available concerning the use of the thermo-economic method to optimize the IT of DH pipes based on heat loss, insulation cost, total cost, and exergy efficiency, such as Refs. [³⁷,³⁸]. Keçebaş et al. [³⁷] established an optimization model depending on LCCA via the P1–P2 method. The results showed that the OIT of insulated DH pipes above ground varied from 0.085 m to 0.228 m in Afyonkarahisar, Turkey. Ucar [³⁸] used thermo-economic analysis to calculate the OIT of insulated DH pipes above ground for four cities (i.e., Aydın, Tekirdağ, Elazığ, and Kars).

Generally speaking, there remain some problems in the above studies. To begin with, the HDD method based on outdoor temperature cannot be used to calculate the heat loss of buried pipes for DH systems. Instead, soil temperature should be used in such calculations. Besides, the method used to calculate the heat loss of buried pipes for DH systems is different from that used to calculate the heat loss of DH pipes above ground [²²]. Next, double pipes, i.e., supply and return pipes, are always used in the DH piping networks. The complementary heat resistance from coinciding temperature fields should be considered in the calculation of the total heat resistance of DH pipes. Therefore, the heat loss of supply and return pipes should be calculated, respectively. Finally, the supply and return temperature were always regarded as constant parameters in the above analysis. Design values of such parameters were used to calculate the heat loss of DH pipes. However, such parameters vary because of different heating demand, operating strategies, heating equipment, etc. Therefore, the heat loss of buried pipes for DH systems should be calculated based on the operating supply and return temperature instead of design values.

The objective of this paper is to obtain the MATC by selecting the optimum economic IT of DPBIS for DH systems. Therefore, the heat loss of the supply and return pipes buried in the soil are calculated, respectively. The temperature interaction between the supply and return pipes is considered in the calculation of the heat loss. Besides, the mathematical model on the OIT of DPBIS for DH systems is established. In addition, the zero point theorem and fsolve function are introduced in the program to solve the optimization model by using MATLAB. Moreover, the impact of related factors, such as heat source, operating strategy, insulation material, NPS, and BD, on the OIT and MATC is analyzed, respectively. Furthermore, sensitivity analysis is conducted to investigate the impact of some economic parameters (i.e., unit heating cost, insulation material price, interest rate, and insulation material lifetime) on optimization results.

Methodology

The DPBIS for DH systems used in this paper are described in Fig. 1 in which the left and right pipes stand for the supply and return pipes, respectively. Hot water is used as the heat medium. The operating supply and return temperature vary during a heating season. Soil temperature depends on the meteorological conditions of the city where DH pipes are installed.

Equivalent annual pipe insulation cost

The pipe insulation cost is proportional to the cross-section area of the insulation materials added to the DH pipes, which can be expressed as Eq. (1).

(1)C_{\mathrm{p}}={\displaystyle \frac{\pi }{\mathrm{4}}}[{(d_{\mathrm{w}}+\mathrm{2}\delta )}^{\mathrm{2}}-d_{\mathrm{w}}^{\mathrm{2}}]c_{\mathrm{ins}}.

The dynamic method is used to calculate the annual capital recovery factor, which is expressed as Eq. (2) [⁴⁰,⁴¹].

(2)\phi ={\displaystyle \frac{I{(\mathrm{1}+I)}^y}{{(\mathrm{1}+I)}^y-\mathrm{1}}}.

The equivalent annual pipe insulation cost can be obtained by using Eq. (3).

(3)C_{\mathrm{pa}}={\displaystyle \frac{\pi }{\mathrm{4}}}[{(d_{\mathrm{w}}+\mathrm{2}\delta )}^{\mathrm{2}}-d_{\mathrm{w}}^{\mathrm{2}}]c_{\mathrm{ins}}\phi .

Annual heat loss cost

As for buried DH pipes, the heat resistance of insulation material, the soil heat resistance, and the complementary heat resistance from coinciding temperature fields are considered in the calculation of the heat loss. It should be noted that the heat resistance from coinciding temperature fields is a result of the temperature interaction between supply and return pipes.

The heat loss of the supply pipe per unit length can be obtained by using Eq. (4) [¹⁹,²²,⁴²].

(4)q_{\mathrm{s}}={\displaystyle \frac{(t_{\mathrm{osp}}-t_{\mathrm{so}})(R_{\mathrm{ins}}+R_{\mathrm{so}})-(t_{\mathrm{orp}}-t_{\mathrm{so}})R_{\mathrm{c}}}{{(R_{\mathrm{ins}}+R_{\mathrm{so}})}^{\mathrm{2}}-R_{\mathrm{c}}{}^{\mathrm{2}}}}.

The heat loss of the return pipe per unit length can be obtained by using Eq. (5) [¹⁹,²²,⁴²].

(5)q_{\mathrm{r}}={\displaystyle \frac{(t_{\mathrm{orp}}-t_{\mathrm{so}})(R_{\mathrm{ins}}+R_{\mathrm{so}})-(t_{\mathrm{osp}}-t_{\mathrm{so}})R_{\mathrm{c}}}{{(R_{\mathrm{ins}}+R_{\mathrm{so}})}^{\mathrm{2}}-R_{\mathrm{c}}{}^{\mathrm{2}}}}.

Therefore, the total heat loss of supply and return pipes can be calculated by using Eq. (6).

(6)q_{\mathrm{t}}=q_{\mathrm{s}}+q_{\mathrm{r}}={\displaystyle \frac{t_{\mathrm{osp}}+t_{\mathrm{orp}}-\mathrm{2}{t}_{\mathrm{so}}}{R_{\mathrm{ins}}+R_{\mathrm{so}}+R_{\mathrm{c}}}}.

The operating supply and return temperature in the primary side (t_osp and t_orp) can be calculated based on that in the secondary side (t_oss and t_ors). When the mass flow rate of the secondary side keeps the design value, the operating supply and return temperature in the secondary side can be calculated by using Eqs. (7) and (8), respectively [¹⁹,⁴²].

(7)t_{\mathrm{oss}}=t_{\mathrm{nd}}+\mathrm{0.5}(t_{\mathrm{dss}}+t_{\mathrm{drs}}-\mathrm{2}{t}_{\mathrm{nd}}){\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{1}/(\mathrm{1}+\alpha )}+\mathrm{0.5}(t_{\mathrm{dss}}-t_{\mathrm{drs}}){\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}.

(8)t_{\mathrm{ors}}=t_{\mathrm{nd}}+\mathrm{0.5}(t_{\mathrm{dss}}+t_{\mathrm{drs}}-\mathrm{2}{t}_{\mathrm{nd}}){\left({\displaystyle \frac{t_{\mathrm{n}\mathrm{d}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{1}/(\mathrm{1}+\alpha )}-\mathrm{0.5}(t_{\mathrm{dss}}-t_{\mathrm{drs}}){\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}.

When the mass flow rate of the secondary side varies, it is always assumed that the relative flow ratio of the secondary side is equal to the relative heating load ratio [¹⁹,⁴²]. The relative flow ratio is the ratio of the operating mass flow rate to the design mass flow rate [¹⁹,⁴²]. The relative heating load ratio means the ratio of the operating heating load at a certain outdoor temperature to the design heating load [¹⁹,⁴²]. Then, the operating supply and return temperature in the secondary side can be calculated by using Eqs. (9) and (10), respectively [¹⁹,⁴²].

(9)t_{\mathrm{oss}}=t_{\mathrm{nd}}+\mathrm{0.5}(t_{\mathrm{dss}}+t_{\mathrm{drs}}-\mathrm{2}{t}_{\mathrm{nd}}){\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{1}/(\mathrm{1}+\alpha )}+\mathrm{0.5}(t_{\mathrm{dss}}-t_{\mathrm{drs}}).

(10)t_{\mathrm{ors}}=t_{\mathrm{nd}}+\mathrm{0.5}(t_{\mathrm{dss}}+t_{\mathrm{drs}}-\mathrm{2}{t}_{\mathrm{nd}}){\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{1}/(\mathrm{1}+\alpha )}-\mathrm{0.5}(t_{\mathrm{dss}}-t_{\mathrm{drs}}).

When the mass flow rate of the primary side keeps the design value, the operating supply and return temperature in the primary side can be calculated by using Eqs. (11) and (12), respectively [¹⁹,⁴²].

(11)\begin{array}{l}{t}_{\mathrm{osp}}={\displaystyle \frac{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}})-t_{\mathrm{oss}}+t_{\mathrm{ors}}\exp \left[{\displaystyle \frac{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}})-t_{\mathrm{oss}}+t_{\mathrm{ors}}}{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\Delta t_{\mathrm{dm}}}}{\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}\right]}{\exp \left[{\displaystyle \frac{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}})-t_{\mathrm{oss}}+t_{\mathrm{ors}}}{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\Delta t_{\mathrm{dm}}}}{\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}\right]-\mathrm{1}}}\\ +{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}}).\end{array}

(12)t_{\mathrm{orp}}={\displaystyle \frac{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}})-t_{\mathrm{oss}}+t_{\mathrm{ors}}\exp \left[{\displaystyle \frac{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}})-t_{\mathrm{oss}}+t_{\mathrm{ors}}}{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\mathrm{\Delta }{t}_{\mathrm{dm}}}}{\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}\right]}{\exp \left[{\displaystyle \frac{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}(t_{\mathrm{dsp}}-t_{\mathrm{drp}})-t_{\mathrm{oss}}+t_{\mathrm{ors}}}{{\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\mathrm{\Delta }{t}_{\mathrm{dm}}}}{\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}\right]-\mathrm{1}}}.

The design logarithmic mean temperature difference of heat exchanger can be calculated by using Eq. (13) [¹⁹,⁴²].

(13)\mathrm{\Delta }{t}_{\mathrm{dm}}={\displaystyle \frac{t_{\mathrm{dsp}}-t_{\mathrm{dss}}-t_{\mathrm{drp}}+t_{\mathrm{drs}}}{\ln {\displaystyle \frac{t_{\mathrm{dsp}}-t_{\mathrm{dss}}}{t_{\mathrm{drp}}-t_{\mathrm{drs}}}}}}.

When the mass flow rate of the primary side varies, it is always assumed that the relative flow ratio of the primary side is equal to the relative heating load ratio [¹⁹,⁴²]. Then, the operating supply and return temperature in the primary side can be calculated by using Eqs. (14) and (15), respectively [¹⁹,⁴²].

(14)t_{\mathrm{osp}}={\displaystyle \frac{t_{\mathrm{dsp}}-t_{\mathrm{drp}}-t_{\mathrm{oss}}+t_{\mathrm{ors}}\exp \left[{\displaystyle \frac{(t_{\mathrm{dsp}}-t_{\mathrm{drp}}-t_{\mathrm{oss}}+t_{\mathrm{ors}}){\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}}{{\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{t}_{\mathrm{dm}}}}\right]}{\exp \left[{\displaystyle \frac{(t_{\mathrm{dsp}}-t_{\mathrm{drp}}-t_{\mathrm{oss}}+t_{\mathrm{ors}}){\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}}{{\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{t}_{\mathrm{dm}}}}\right]-\mathrm{1}}}+t_{\mathrm{dsp}}-t_{\mathrm{drp}}.

(15)t_{\mathrm{orp}}={\displaystyle \frac{t_{\mathrm{dsp}}-t_{\mathrm{drp}}-t_{\mathrm{oss}}+t_{\mathrm{ors}}\exp \left[{\displaystyle \frac{(t_{\mathrm{dsp}}-t_{\mathrm{drp}}-t_{\mathrm{oss}}+t_{\mathrm{ors}}){\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}}{{\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{t}_{\mathrm{dm}}}}\right]}{\exp \left[{\displaystyle \frac{(t_{\mathrm{dsp}}-t_{\mathrm{drp}}-t_{\mathrm{oss}}+t_{\mathrm{ors}}){\left({\displaystyle \frac{t_{\mathrm{oss}}-t_{\mathrm{ors}}}{t_{\mathrm{dss}}-t_{\mathrm{drs}}}}\right)}^{\mathrm{0.5}}}{{\left({\displaystyle \frac{t_{\mathrm{nd}}-t_{\mathrm{o}}}{t_{\mathrm{nd}}-t_{\mathrm{od}}}}\right)}^{\mathrm{0.5}}\mathrm{\Delta }{t}_{\mathrm{dm}}}}\right]-\mathrm{1}}}.

The complementary heat resistance from coinciding temperature fields can be calculated by using Eq. (16) [¹⁹,²²,⁴²].

(16)R_{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi \lambda { }_{\mathrm{so}}}}\ln \sqrt{{\left({\displaystyle \frac{\mathrm{2}h}{b}}\right)}^{\mathrm{2}}+\mathrm{1}}.

The IT of the supply and return pipes is assumed to be the same. Then, the heat resistance of insulation material can be calculated by using Eq. (17) [¹⁹,²²,⁴²].

(17)R_{\mathrm{ins}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{ins}}}}\ln {\displaystyle \frac{d_{\mathrm{w}}+\mathrm{2}\delta }{d_{\mathrm{w}}}}.

The soil heat resistance can be calculated by using Eq. (18) [¹⁹,²²,⁴²].

(18)R_{\mathrm{so}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{so}}}}\ln {\displaystyle \frac{\mathrm{4}h}{d_{\mathrm{w}}+\mathrm{2}\delta }}.

Considering the variations of the operating supply and return temperature, it is necessary to calculate the hourly heat loss. Then, the annual heat loss cost of supply and return pipes can be calculated by using Eq. (19).

(19)C_{\mathrm{hl}}={\mathrm{10}}^{-\mathrm{9}}{c}_{\mathrm{h}}{(R_{\mathrm{ins}}+R_{\mathrm{so}}+R_{\mathrm{c}})}^{-\mathrm{1}}{\displaystyle \sum _{i=\mathrm{1}}^n[\mathrm{3600}(t_{\mathrm{osp(}i)}+t_{\mathrm{orp(}i)}-\mathrm{2}{t}_{\mathrm{so(}i)})]}.

where n is the total number of operating hours during a heating season (h), i is the index for operating hours, t_osp(i) is the operating supply temperature in the primary side at time i (°C), t_orp(i) is the operating return temperature in the primary side at time i (°C), and t_so(i) is the soil temperature at time i (°C).

Optimization model

From the above analysis, the annual total cost includes the equivalent annual pipe insulation cost and the annual heat loss cost, which can be obtained by using Eq. (20).

(20)\begin{array}{l}{C}_t={\displaystyle \frac{\pi }{4}}[{(d_{\mathrm{w}}+\mathrm{2}\delta )}^{\mathrm{2}}-d_{\mathrm{w}}^{\mathrm{2}}]c_{\mathrm{ins}}\phi +{\mathrm{10}}^{-\mathrm{9}}{c}_{\mathrm{h}}{\displaystyle \sum _{i=\mathrm{1}}^n[\mathrm{3600}(t_{\mathrm{osp(}i)}+t_{\mathrm{orp(}i)}-\mathrm{2}{t}_{\mathrm{so(}i)})]}\\ \times {\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{ins}}}}\ln {\displaystyle \frac{d_{\mathrm{w}}+\mathrm{2}\delta }{d_{\mathrm{w}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{so}}}}\ln {\displaystyle \frac{\mathrm{4}h}{d_{\mathrm{w}}+\mathrm{2}\delta }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{so}}}}\ln \sqrt{{\left({\displaystyle \frac{\mathrm{2}h}{b}}\right)}^{\mathrm{2}}+\mathrm{1}}\right]}^{-\mathrm{1}}.\end{array}

When other variables are determined, δ is the only variable that influences the value of C_t. In Eq. (20), the minimum value of C_t can be obtained by selecting the OIT of DPBIS for DH systems. As for Eq. (20), the derivative on δ can be expressed as Eq. (21).

(21)\begin{array}{l}{\displaystyle \frac{d{C}_{\mathrm{t}}}{d\delta }}\mathrm{=\pi }{c}_{\mathrm{ins}}\phi (d_{\mathrm{w}}+\mathrm{2}\delta )-{\mathrm{10}}^{-\mathrm{9}}{c}_{\mathrm{h}}\left[\left({\displaystyle \frac{\mathrm{1}}{\pi {\lambda }_{ \mathrm{ins}}}}-{\displaystyle \frac{\mathrm{1}}{\pi {\lambda }_{ \mathrm{so}}}}\right){\displaystyle \frac{\mathrm{1}}{d_{\mathrm{w}}+\mathrm{2}\delta }}\right]{\displaystyle \sum _{i=\mathrm{1}}^n[\mathrm{3600}(t_{\mathrm{osp(}i)}+t_{\mathrm{orp(}i)}-\mathrm{2}{t}_{\mathrm{so(}i)})]}\\ \times {\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{ins}}}}\ln {\displaystyle \frac{d_{\mathrm{w}}+\mathrm{2}\delta }{d_{\mathrm{w}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{so}}}}\ln {\displaystyle \frac{\mathrm{4}h}{d_{\mathrm{w}}+\mathrm{2}\delta }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\pi {\lambda }_{ \mathrm{so}}}}\ln \sqrt{{\left({\displaystyle \frac{\mathrm{2}h}{b}}\right)}^{\mathrm{2}}+\mathrm{1}}\right]}^{-\mathrm{2}}. \end{array}

To facilitate construction and maintenance, the distance between the outer surfaces of DPBIS for DH systems should not be less than 0.2 m [¹⁹], which is expressed as Eq. (22).

(22)b-d_{\mathrm{w}}-\mathrm{2}\delta \ge \mathrm{0.2}

When {\displaystyle \frac{d{C}_{\mathrm{t}}}{d\delta }}=\mathrm{0}, the optimal value of δ can be obtained. Besides, the minimum value of C_t can be obtained by using Eq. (20). Thus, the optimization model used to obtain the OIT of DPBIS for DH systems is established.

Solution of the optimization model

MATLAB can be used to solve this kind of problem. The zero point theorem and fsolve function can be applied in the program [⁴³–⁴⁵]. The zero point theorem can be expressed as:

If f(x) is continuous in [A,B], and f(A)f(B)<0, then there exists M(A,B) such that f(M)=0.

When the fsolve function of MATLAB is used, it means that the least square method is used to solve nonlinear equations. The program flowchart used to solve the optimization model is shown in Fig. 2.

Results

The DPBIS for DH systems in Beijing, China was analyzed in this paper. The outdoor design temperature for heating is –7.6°C [⁴⁶]. The heating season lasts from November 12 to March 14 [⁴⁶]. Stainless steel pipes were used as the heat distribution pipes for DH systems. The properties of stainless steel pipes at various nominal pipe sizes (i.e., 50 mm, 100 mm, 200 mm, 400 mm, 600 mm, 800 mm, and 1000 mm) are presented in Table 1 [¹⁹]. Foam rubber, rigid phenolic foam (RPHF), and rigid polyurethane foam (RPOF) were assumed to be the insulation materials of DPBIS for DH systems, respectively. According to market survey, the properties of the insulation materials used in this paper are tabulated in Table 2. It was assumed that three types of heat sources were used, including coal-fired combined heat and power plant (CFCHPP), coal-fired boiler (CFB), and natural gas-fired boiler (NGFB). The unit heating cost of CFCHPP, CFB, and NGFB is 2.324 $/GJ [⁴⁷], 3.486 $/GJ [⁴⁷], and 12.933 $/GJ [⁴⁷,⁴⁸], respectively. As for DH systems, operating supply and return temperature depends on the operating strategy, which may influence the annual heat loss. So it is necessary to analyze the optimization results based on different operating strategies. Four operating strategies based on the variations of heating parameters are listed in Table 3 [²]. In Table 3, constant mass flow rate means that the operating mass flow rate is equal to the design mass flow rate, while the operating supply and return temperature varies in the process of the operation of DH systems. Variable mass flow rate means that the operating mass flow rate, operating supply temperature, and operating return temperature all vary in the process of the operation of DH systems. According to the design requirements specified in Refs. [⁴⁶,⁴⁹], the design values of the primary supply temperature, primary return temperature, secondary supply temperature, and secondary return temperature are 130°C, 70°C, 85°C, and 60°C, respectively. The operating values of the above parameters during a heating season were used to calculate the hourly heat loss (see Section 2. 2). Figure 3 depicts the variations of the operating supply and return temperature based on different operating strategies during a typical heating season. The distance between pipe center and ground surface is not determined by pipe insulation requirements [⁴⁹]. So the cost related to the BD of DH pipes was not considered in the annual total cost. However, the annual heat loss cost is still influenced by the BD of DH pipes. The BD of DPBIS for DH systems was assumed to be 1.5 m, 2.0 m, and 2.5 m, respectively. In addition, the soil temperature for different BD in a typical year is exhibited in Fig. 4. Some related parameters used in calculations are summarized in Table 4.

Annual total cost

The annual total cost should be calculated in order to obtain the OIT. Here, the annual total cost includes the equivalent annual pipe insulation cost and the annual heat loss cost (see Section 2). Figure 5 illustrates the impact of IT on the annual cost. As can be seen in Fig. 5, the annual heat loss cost decreases with the increase in the IT, while the equivalent annual pipe insulation cost increases approximately linearly when the IT increases. When the IT is smaller than the OIT, the annual total cost decreases with the increase in the IT. On the contrary, when the IT is larger than OIT, the annual total cost increases with the increase in the IT. The reason for this is that the downtrend of the annual heat loss cost becomes less obvious when the IT increases.

The annual total cost is affected by many factors, e.g., heat source, operating strategy, insulation material, BD, etc. The impact of the above factors on the annual total cost is demonstrated in Figs. (6)–(9). Figure 6 describes the impact of heat sources on the annual total cost. As can be observed in Fig. 6, the annual total cost for NGFB is much larger than that for other heat sources. The reason for this relates to the fact that the unit heating cost of NGFB is much higher than that of other heat sources. More annual heat loss cost is caused when NGFB is used compared with that when other heat sources are used. Figure 7 describes the impact of operating strategies on the annual total cost. As can be seen in Fig. 7, the annual total cost for Strategies 1 to 4 gradually increases. The reason for this is that the sum of the operating supply and return temperature for Strategies 1 to 4 increases (see Fig. 3). It can be also seen from Fig. 7 that the annual total cost for Strategies 1 and 2 is similar, while that for Strategies 3 and 4 is similar. The reason for this is that the operating supply and return temperature for Strategies 1 and 2 is similar, while the operating supply and return temperature for Strategies 3 and 4 is similar (see Fig. 3). Figure 8 describes the impact of insulation materials on the annual total cost. It can be seen from Fig. 8 that the annual total cost for foam rubber is higher than that for other insulation materials. The main reason for this is that the price of foam rubber is higher than that of other insulation materials. More equivalent annual pipe insulation cost is caused when foam rubber is used compared with that when other insulation materials are used. Figure 9 describes the impact of BD on the annual total cost. It can be seen from Fig. 9 that as the BD increases, the annual total cost decreases. The reason for this relates to the fact that the soil temperature during most of the heating season increases with the increase in BD (see Fig. 4). Moreover, the annual heat loss cost decreases when the BD increases.

Optimization results

The optimization model (see Section 2.3) was used to calculate the OIT and MATC of DPBIS for DH systems. The optimization results of DPBIS for DH systems are displayed in Figs. (10)–(12). It can be seen that the OIT and MATC in Figs. (10)–(12) are consistent with those described in Figs. (6)–(9), which indicates that the optimization model can be used to obtain the OIT and MATC of DPBIS for DH systems.

From Figs. (10)–(12), it can be obtained that the OIT varies from 0.055 m to 0.108 m, 0.065 m to 0.133 m, and 0.104 m to 0.250 m for CFCHPP, CFB, and NGFB, respectively. The MATC varies from 0.553 $/(m∙a) to 7.426 $/(m∙a), 0.738 $/(m∙a) to 9.468 $/(m∙a), and 1.998 $/(m∙a) to 20.775 $/(m∙a) for CFCHPP, CFB, and NGFB, respectively. The OIT and MATC for NGFB are much larger than those for other heat sources. The reason for this relates to the fact that the unit heating cost of NGFB is much higher than that of other heat sources. The unit heating cost of NGFB is 5.56 times and 3.71 times that of CHCHPP and CFB, respectively. As the unit heating cost increases, the annual heat loss cost increases. Therefore, more annual heat loss cost is needed when NGFB is used compared with that when other heat sources are used. Besides, the OIT and MATC for Strategies 1 to 4 gradually increase. The reason for this relates to the fact that the sum of the operating supply and return temperature for Strategies 1 to 4 increases (see Fig. 3), resulting in the increase in the annual heat loss cost. Therefore, it is necessary to consider the operating strategy in the calculation of OIT and MATC. In addition, the OIT for foam rubber, RPHF, and RPOF varies from 0.055 m to 0.191 m, 0.070 m to 0.250 m, and 0.065 m to 0.232 m, respectively. The MATC for foam rubber, RPHF, and RPOF varies from 0.820 $/(m∙a) to 20.755 $/(m∙a), 0.553 $/(m∙a) to 13.146 $/(m∙a), and 0.747 $/(m∙a) to 18.120 $/(m∙a), respectively. The OIT for RPHF is larger than that for other insulation materials. But the MATC for foam rubber is larger than that for other insulation materials. The reason for this relates to the fact that the price and heat conductivity of the insulation materials influence the OIT and MATC. Moreover, the MATC increases with the increase in NPS. This is because the equivalent annual pipe insulation cost and annual heat loss cost increase with the increase in NPS. However, the OIT does not always increases with the increase in NPS. The reason for this relates to the fact that the effect of NPS on the heat loss decreases with the increase in NPS. As the NPS increases, the outer surface area of the insulated DH pipes increases. Thus, the annual heat loss increases with the increase in NPS. But it should be noted that the above effect reduces when the NPS increases. Furthermore, the OIT and MATC decrease with the increase in BD. The reason for this relates to the fact that the soil temperature during a heating season increases when BD increases, and the annual heat loss decreases with the increase in BD. Therefore, as BD increases, the annual heat loss cost decreases. Compared with operating strategy, insulation material, NPS, and BD, the heat source has a greater impact on the OIT and MATC. This is because the annual heat loss cost accounts for a large proportion of the MATC. And the annual heat loss cost is determined by the unit heating cost related to heat sources. The average proportion of annual heat loss cost to the MATC for foam rubber, RPHF, and RPOF is 64.10%, 64.04%, and 64.64%, respectively.

Sensitivity analysis

The OIT and MATC are influenced by some economic parameters, e.g., unit heating cost, insulation material price, interest rate, insulation material lifetime, etc. The sensitivity analysis method was used to investigate the impact of these parameters on OIT and MATC. The one-factor sensitivity analysis method was used in this paper. The initial values of sensitivity factors can be seen in Tables 2 and 4. In the process of sensitivity analysis, one sensitivity factor varied while the other factors were kept constant. Then the OIT and MATC changed with the variation of the sensitivity factor. Because the results for different operating strategies, NPS, and BD are similar, it was assumed that the DH system was operated under Strategy 1, the NPS was 400 mm, and the BD was 1.5 m.

The sensitivity analysis results are shown in Fig. 13. When the sensitivity factors vary by 10%, the average variation rate of the OIT and MATC for different types of heat sources and insulation materials are shown in Table 5. It can be seen from Fig. 13 and Table 5 that the OIT and MATC are influenced by sensitivity factors. However, the impact of each sensitivity factor on the OIT and MATC is different. Then, it can be concluded that the OIT and MATC increase with the increase in the unit heating cost. This is related to the fact that the annual heat loss cost increases when the unit heating cost increases. The increase in insulation material price results in the decrease in the OIT and the increase in the MATC. The reason relates to the fact that the equivalent annual pipe insulation cost increases with the increase in the insulation material price. The increase in the interest rate results in the decrease in the OIT and the increase in the MATC. The reason lies in the fact that the equivalent annual pipe insulation cost increases with the increase in the interest rate. The increase in the insulation material lifetime results in the increase in the OIT and the decrease in the MATC. The reason relates to the fact that the equivalent annual pipe insulation cost decreases with the increase in the insulation material lifetime. Besides, compared with other parameters, the interest rate and insulation material lifetime have a smaller impact on the OIT and MATC. The reason relates to the fact that compared with other parameters, the interest rate and insulation material lifetime have a smaller impact on the equivalent annual pipe insulation cost. Moreover, compared with other parameters, the unit heating cost and insulation material price have a greater impact on the OIT. The reason lies in the fact that compared with other parameters, the unit heating cost and insulation material price have a greater impact on the annual heat loss cost and equivalent annual pipe insulation cost, respectively. Furthermore, of all parameters, the unit heating cost has the greatest impact on the MATC. This is due to the fact that the proportion of the annual heat loss cost to the annual total cost is larger than that of the equivalent annual pipe insulation cost to the annual total cost.

Conclusions

The IT of DPBIS for DH systems was optimized to minimize the annual total cost. An optimization model was established to calculate the OIT and MATC of DPBIS for DH systems. The results show that the OIT and MATC of DPBIS for DH systems can be obtained by using the optimization model. Based on the results, the following conclusions can be drawn:

The annual total cost for NGFB is larger than that for other heat sources. The annual total cost for Strategies 1 to 4 gradually increases. The annual total cost for foam rubber is larger than that for other insulation materials. The annual total cost decreases with the increases in BD.

The OIT and MATC for NGFB are larger than those for other heat sources. The OIT and MATC for Strategies 1 to 4 gradually increase. The OIT for RPHF is larger than that for other insulation materials. But the MATC for foam rubber is larger than that for other insulation materials. The MATC increases with the increase in NPS. However, the OIT does not always increases with the increase in NPS. The reason for this relates to the fact that the effect of NPS on the heat loss decreases with the increase in NPS. The OIT and MATC decrease with the increase in BD.

According to sensitivity analysis, the OIT and MATC are influenced by unit heating cost, insulation material price, interest rate, and insulation material lifetime. The increase in unit heating cost results in the increase in the OIT and MATC. However, the increase in insulation material price or interest rate results in the decrease in the OIT and the increase in the MATC. The increase in the insulation material lifetime results in the increase in the OIT and the decrease in the MATC. Compared with other parameters, the interest rate and insulation material lifetime have a smaller impact on the OIT and MATC. The unit heating cost and insulation material price have a greater impact on the OIT. The unit heating cost has the greatest impact on the MATC.

However, it should be noted that the objective of this paper is to obtain the OIT of DPBIS for DH systems from an economic perspective. The OIT of DPBIS for DH systems obtained from other perspectives may be different. Therefore, it is necessary to conduct a more comprehensive study in the future.

References

[1] Åberg M. Investigating the impact of heat demand reductions on Swedish district heating production using a set of typical system models. Applied Energy, 2014, 118: 246–257

[2] Jie P F, Zhu N, Li D Y. Operation optimization of existing district heating systems. Applied Thermal Engineering, 2015, 78: 278–288

[3] Building Energy Research Center of Tsinghua University. 2015 Annual Report on China Building Energy Efficiency. Beijing: China Architecture & Building Press, 2015

[4] Kaynakli O. A study on residential heating energy requirement and optimum insulation thickness. Renewable Energy, 2008, 33(6): 1164–1172

[5] Rocha P, Kaut M, Siddiqui A S. Energy-efficient building retrofits: an assessment of regulatory proposals under uncertainty. Energy, 2016, 101: 278–287

[6] Mikulić D, Bakarić I R, Slijepčević S. The economic impact of energy saving retrofits of residential and public buildings in Croatia. Energy Policy, 2016, 96: 630–644

[7] Walter T, Sohn M D. A regression-based approach to estimating retrofit savings using the Building Performance Database. Applied Energy, 2016, 179: 996–1005

[8] Ciulla G, Galatioto A, Ricciu R. Energy and economic analysis and feasibility of retrofit actions in Italian residential historical buildings. Energy and Buildings, 2016, 128: 649–659

[9] Lund R, Mohammadi S. Choice of insulation standard for pipe networks in 4th generation district heating systems. Applied Thermal Engineering, 2016, 98: 256–264

[10] Berge A, Hagentoft C, Adl-Zarrabi B. Field measurements on a district heating pipe with vacuum insulation panels. Renewable Energy, 2016, 87(Part 3): 1130–1138

[11] Masatin V, Volkova A, Hlebnikov A, Latosov E. Improvement of district heating network energy efficiency by pipe insulation renovation with PUR foam shells. Energy Procedia, 2017, 113: 265–269

[12] Bloomquist R G. Geothermal district energy system analysis, design, and development. Dissertation for the Doctoral Degree. Romania: University of Oradea, 2001

[13] Kalyon M, Sahin A Z. Application of optimal control theory in pipe insulation. Numerical Heat Transfer Part A: Applications, 2002, 41(4): 391–402

[14] Dalla Rosa A, Li H, Svendsen S. Method for optimal design of pipes for low-energy district heating, with focus on heat losses. Energy, 2011, 36(5): 2407–2418

[15] Wojdyga K, Chorzelski M. Chances for Polish district heating systems. Energy Procedia, 2017, 116: 106–118

[16] Berge A, Adl-Zarrabi B, Hagentoft C. Assessing the thermal performance of district heating twin pipes with vacuum insulation Panels. Energy Procedia, 2015, 78: 382–387

[17] Yarahmadi N, Vega A, Jakubowicz I. Determination of essential parameters influencing service life time of polyurethane insulation in district heating pipes. Energy Procedia, 2017, 116: 320–323

[18] Wang X L, Li H, Wang Q L. Influence of insulation layer position on short-time intermittent heating room. Procedia Engineering, 2017, 205: 1484–1493

[19] He P, Sun G, Wang F, Wu H X. Heating Engineering. Beijing: China Architecture & Building Press, 2009

[20] Köfinger M, Basciotti D, Schmidt R R, Meissner E, Doczekal C, Giovannini A. Low temperature district heating in Austria: energetic, ecologic and economic comparison of four case studies. Energy, 2016, 110: 95–104

[21] Tunzi M, Boukhanouf R, Li H W, Svendsen S, Ianakiev A. Improving thermal performance of an existing UK district heat network: a case for temperature optimization. Energy and Buildings, 2018, 158: 1576–1585

[22] Frederiksen S, Werner S. District Heating and Cooling. Lund: Studentlitteratur, 2013

[23] Nord N, Nielsen E K L, Kauko H, Tereshchenko T. Challenges and potentials for low-temperature district heating implementation in Norway. Energy, 2018, 151: 889–902

[24] Lund H, Østergaard P A, Chang M, Werner S, Svendsen S, Sorknæs P, Thorsen J E, Hvelplund F, Mortensen B O G, Mathiesen B V, Bojesen C, Duic N, Zhang X L, Möller B. The status of 4th generation district heating: research and results. Energy, 2018, 164: 147–159

[25] Lund H, Duic N, Østergaard P A, Mathiesen B V. Future district heating systems and technologies: on the role of smart energy systems and 4th generation district heating. Energy, 2018, 165 (Part A): 614–619

[26] China Institute of Building Standard Design & Research. National Technical Measures for Design of Civil Construction (Heating, Ventilation and Air Conditioning). Beijing: China Planning Press, 2009

[27] Zhang L Y, Wang Z N, Yang X H, Jin L W, Zhang Q L, Hu W J. Thermo-economic analysis for directly-buried pipes insulation of district heating piping systems. Energy Procedia, 2017, 105: 3369–3376

[28] Keçebaş A. Determination of insulation thickness by means of exergy analysis in pipe insulation. Energy Conversion and Management, 2012, 58: 76–83

[29] Başoğul Y, Demircan C, Keçebaş A. Determination of optimum insulation thickness for environmental impact reduction of pipe insulation. Applied Thermal Engineering, 2016, 101: 121–130

[30] Özel G, Açıkkalp E, Görgün B, Yamık H, Caner N. Optimum insulation thickness for piping system using exergy and environmental methods. International Journal of Global Warming, 2017, 11(1): 107–123

[31] Başoğul Y, Keçebaş A. Economic and environmental impacts of insulation in district heating pipelines. Energy, 2011, 36(10): 6156–6164

[32] Daşdemir A, Ertürk M, Keçebaş A, Demircan C. Effects of air gap on insulation thickness and life cycle costs for different pipe diameters in pipeline. Energy, 2017, 122: 492–504

[33] Kayfeci M, Yabanova İ, Keçebaş A. The use of artificial neural network to evaluate insulation thickness and life cycle costs: pipe insulation application. Applied Thermal Engineering, 2014, 63(1): 370–378

[34] Daşdemir A, Ural T, Ertürk M, Keçebaş A. Optimal economic thickness of pipe insulation considering different pipe materials for HVAC pipe applications. Applied Thermal Engineering, 2017, 121: 242–254

[35] Kayfeci M. Determination of energy saving and optimum insulation thicknesses of the heating piping systems for different insulation materials. Energy and Buildings, 2014, 69: 278–284

[36] Keçebaş A. Variation of insulation thickness and exergetic cost saving with outdoor temperature in pipe insulation. Environmental Progress & Sustainable Energy, 2013, 32(3): 784–789

[37] Keçebaş A, Ali Alkan M, Bayhan M. Thermo-economic analysis of pipe insulation for district heating piping systems. Applied Thermal Engineering, 2011, 31(17–18): 3929–3937

[38] Ucar A. A thermo-economic approach for determination of optimum thickness of thermal insulation applied to piping systems. Journal of Renewable and Sustainable Energy, 2018, 10(3): 035102

[39] Keçebaş A. Determination of optimum insulation thickness in pipe for exergetic life cycle assessment. Energy Conversion and Management, 2015, 105: 826–835

[40] Jie P F, Kong X F, Rong X, Xie S Q. Selecting the optimum pressure drop per unit length of district heating piping network based on operating strategies. Applied Energy, 2016, 177: 341–353

[41] Liu X M. Engineering Economic Analysis. Beijing: Peking University Press, 2014

[42] Li D Y. Heating Engineering. Beijing: China Architecture & Building Press, 2018

[43] Yuan S C, Zhang J H, Zheng J X. MATLAB Language and Applied Technology. Beijing: National Defense Industry Press, 2011

[44] Su J M, Ruan S Y. MATLAB Practical Tutorial. Beijing: Electronic Industry Press, 2008

[45] Zhang L J, Xiao G, Yang D, Xu S B. MATLAB Data Analysis and Data Mining. Beijing: China Machine Press, 2016

[46] China Academy of Building Research. Design Code for Heating Ventilation and Air Conditioning of Civil Buildings (GB 50736–2012). Beijing: China Architecture & Building Press, 2012

[47] Jie P F, Zhang F H, Fang Z, Wang H B, Zhao Y F. Optimizing the insulation thickness of walls and roofs of existing buildings based on primary energy consumption, global cost and pollutant emissions. Energy, 2018, 159: 1132–1147

[48] Beijing Municipal Commission of Development and Reform. Inquiry about price. 2020, available at the website of beijing.gov.cn

[49] Beijing Gas and Heating Engineering Design Institute. Design Code for City Heating Network (CJJ 34–2010). Beijing: China Architecture & Building Press, 2010

RIGHTS & PERMISSIONS

Higher Education Press