Time-Dependent Integration of Solar Thermal Technology in Industrial Processes

Abstract

Solar energy is currently an underutilized renewable energy source that could fulfill low-temperature industrial heat demands with significant potential in high solar irradiance counties such as Malaysia. This study proposes a new systematic method for optimization of solar heat integration for different process options to minimize the levelized cost of heat by combining different methods from the literature. A case study from the literature is presented to demonstrate the proposed method combined with meteorological data in Malaysia. The method estimates capital cost and levelized cost of solar heating considering important physical constraints (e.g., available space) and recovery of waste heat. The method determines and optimizes important physical dimensions, including collector area, storage size, and control design. As the result of the case study, the solar thermal integration with Clean-In-Place streams (hot water) gives the lowest levelized cost of heat with RM 0.63/kWh (0.13 EUR/kWh) due to its lowest process temperature requirement. The sensitivity analysis indicates that collector price and collector efficiency are the critical parameters of solar thermal integration.

1. Introduction

A driver behind the rapid technological growth in renewable energy generation is the harmful environmental effects of using fossil fuels. The generation of energy by the consumption of fossil fuel resources has led to global warming and climate change. The overall temperature of the atmosphere will be increased and it will negatively affect not only the environment but also mankind itself. This deeper awareness of the environment has fast-tracked a transition towards cleaner, safer, and more sustainable energy systems [1]. The declined price in solar energy equipment, coupled with greenhouse gas (GHG) pricing, peak energy considerations, and renewed global economic growth, have accelerated most of the predictions of experts on the time when solar will achieve price parity with traditional fossil-fuel sources [2]. Minimal maintenance effort and high reliability make solar energy a favorable source of energy to be exploited increasingly in the future [3]. Solar energy can provide energy for cities and industry using large solar systems, while at the same time offering a stand-alone capability in the most isolated rural regions. There are two conventional ways to directly harness solar energy—photovoltaic (PV) and solar thermal. PV cells directly convert light into electricity via a solar panel and it requires direct solar radiation for maximal function. Solar thermal technology, on the other hand, captures solar radiation by heating a fluid, which is then used to generate electricity or process heating. A solar thermal energy system is often used in the industrial sector to supply low-grade thermal heat (80 to 100 °C) [4], and it becomes more favorable in situations of high solar irradiances and high ambient air temperature to guarantee high efficiency of a solar collector [5]. The solar thermal system design is essential to assure its optimal performance and maximize economic benefits. However, changes in one of the design parameters could cause a significant variation in the performance of the solar system. Therefore, a method is required to design solar thermal systems to optimally integrate into industrial processes.

There are several case studies on the viability and optimization of solar thermal integration. For example, Schnitzer et al. [6] presented an example in Austria and introduced appropriate tools for designing solar thermal integration in the supply energy structure of plants. They simulated the system and calculated return of investment and environmental performance, but the only independent variable in their simulation was the area of the solar collector. Atkins et al. [7] focused on the utilization of solar thermal technology by conducting a pinch analysis of a dairy processing plant in New Zealand. Their solar thermal integration was designed for a specific process, given a 1000 m² solar collector and different integration approaches, and reported the energy saving potential for each integration option. Nemet et al. [8] developed a procedure that combined time slices of the solar thermal supply profile and process heat load profile to handle the unpredictability of solar irradiance. Their study, however, only considered a combination of solar energy supply and energy demand by processes in time slices, but not the temperature of the solar energy supply and the temperature required by processes. Walmsley et al. [9] investigated the effects of adding solar as an additional heat source to into a heat recovery loop (HRL) for supplying heat energy to all the processes in a New Zealand dairy factory case study. By a similar method to [9], Walmsley et al. [10] investigated three general methods for integrating solar thermal heating into an HRL by using a New Zealand dairy factory case study and concluded that solar heat should feed into a heat storage tank to dampen temperature fluctuations. Allouhi et al. [11] examined the utilization of a centralized solar water-heating system in a milk processing industry located in Casablanca based on economic feasibility. This study mainly focused on the effect of solar collector area, tilt angle, and storage tank volume on the economic feasibility. Schweiger et al. [12] analyzed the viability of solar energy technology in Spain and Portugal with some examples in the milk, brewing, malting, textile, and paper industries, in which solar energy was integrated into a centralized utility system with the simulation of different types of solar collectors. Hassine et al. [13] developed a general procedure for the identifying of suitable integration points and integration concepts for solar thermal process heat systems, and the integration concepts were classified based on the variety of thermal process operation concepts. Their guideline, however, only considered the selection of a solar thermal integration concept but not the detailed sizing of solar heat equipment.

The above literature review highlighted several research gaps. Firstly, stagnation control and availability of space for solar collector installation have not been considered in the literature on solar integration. Stagnation control is essential to prevent possible damage to the solar collector from overheating that occurs during stagnant flow. The availability of space analysis is essential for the determination of maximum solar collector area that can be installed in the process plant precinct. Many studies propose designs for integrating solar thermal with multiple processes. However, different supply and return temperatures may be required depending on the requirement of different process unit operations [13]. Atkins et al. [7], for example, offered possible configurations for integrating solar heating with dairy processes but did not evaluate the effectiveness of solar heat integration for a wide range of process demands as the efficiency and sizing of solar heat equipment will be different depending on the process demands. Their study also did not conduct an economic assessment. Integration studies also have not fully considered the time profile of solar and process fluctuations, control design, storage design, system stability, collector area, and storage size. Nemet et al. [8] considered the time profile of all the processes in the plant, but they did not attempt to optimize collector area to meet the time profile of process demand without any extra utility. Lastly, there are few studies on solar thermal integration in the context of Malaysia. Since local meteorological data affects the energy saving potential of solar integration, it is critical to study an installation for a specific country and location. To summarize, the main gap in the literature is a comprehensive method that can optimize the design of solar heat integration in terms of economic performance and compare the economic performance while considering a wide range of control and design issues that add important constraints and costs to the system.

This study aims to develop a robust method for the optimization of solar heat integration for different process options to minimize the levelized cost of heat by combining different methods from the literature. A case study from the literature is presented to demonstrate the proposed method by using meteorological data in Malaysia.

The novel contributions of this study to the literature are:

• A new systematic method for designing solar thermal integration by combining different methods from the literature for the design of solar integration.
• The inclusion of critical constraints into the method, including stagnation control and availability of space for solar collector installation.
• A method to optimize the design of solar thermal integration, minimizing the levelized cost of heat.
• A detailed case study of solar thermal integration using meteorological data from Malaysia, being a first case study of its kind.

2. Methods

The procedure of the proposed method is shown in Figure 1. The method determines and optimizes critical physical dimensions, including collector area, storage size, and control design while considering significant physical constraints (e.g., available space) and recovery of low-grade heat. After that, the method estimates the capital cost and levelized cost of heat based on the optimized integration and physical dimensions. The method assumes the output temperature of the solar thermal system was set to the exact process temperature required and remained constant, as well as the input and output temperature in the solar collector and heat storage.

The subsequent subsections then describe the various elements of the overall method. Section 2.1 describes the determination of potential stream for solar integration based on stream data and a grand composite curve. Section 2.2 explains the determination of solar collector heat output from the classified solar irradiance and ambient temperature profile. Section 2.3 gives an insight into the determination of solar collector area and storage size based on a process heat demand profile. Section 2.4 is dedicated to the determination on the availability of space for solar collector, which is one of the important physical constraints. Section 2.5 gives an overview of the stagnation control design. Section 2.6 describes the estimation of capital cost and levelized cost of heat of the solar integration, followed by Section 2.7, which explains the optimization of solar collector area based on the levelized cost of heat. The collector area that results in the lowest levelized cost of heat is selected, and the final storage volume is determined by the method described in Section 2.3.

2.1. Identification of Potential Streams

Stream data can be established by measuring heating requirements (cold streams) and cooling requirements (hot streams) in the process plant. Each stream is defined by its mass flow rate, specific heat capacity, supply temperature, and target temperature. The grand composite curve (GCC) can be drawn from stream data to obtain the possible solar collector operation temperatures and the potential stream within that temperature range. As stated by Atkins et al. [7], an energy penalty will occur when the solar thermal integration is integrated below the pinch temperature, since it adds to heat rejection below the Pinch rather than substitutes hot utility above the Pinch. The GCC is an important tool to identify some potential streams that are above the pinch temperature. Figure 2 shows the system configuration for the solar thermal integration in this study.

The configuration with two heat exchangers (HEX) is recommended when the process medium, storage medium, and the solar heat transfer fluid are different. Similar configurations have been used in the brewery industry in Germany [14] and proposed for the milk processing industry in Casablanca, Morocco [11]. Sensible heat storage is used in this study, as it is currently the most commonly integrated storage type in solar thermal systems [15]. A single stratified tank system is used for the integration because it is cheaper than two tank systems [16], reduces storage heat loss, and requires a simpler storage control strategy. The hot and cold fluids in the same tank are separated by a thermocline where the hot fluid is stored on top (lower density), and the cold fluid is stored on the bottom (higher density) [16].

Since there are two heat exchangers in the system, inlet and outlet temperatures of the solar collector should be higher than the target process temperature to create a sufficient temperature driving force for heat transfer. The temperatures can be calculated by adding the desired output temperature with the multiplication of selected ΔT_min and the number of heat exchangers in the system.

2.2. Determination of Solar Collector Heat Output

Determination of heat output requires the efficiency of the solar collector. The solar irradiance profiles are classified and the most common profiles are used for the calculation. The efficiency of the solar collector is given in Equation (1) [5]:

$${\mathsf{\eta }}_{\mathrm{coll}}={\mathrm{a}}_0-\frac{{\mathrm{a}}_1}{\mathrm{G}}\left({\mathrm{T}}_{\mathrm{m},\mathrm{coll}}-{\mathrm{T}}_{\mathrm{a}}\right)-\frac{{\mathrm{a}}_2}{\mathrm{G}}{\left({\mathrm{T}}_{\mathrm{m},\mathrm{coll}}-{\mathrm{T}}_{\mathrm{a}}\right)}^2$$

(1)

The average collector internal fluid temperature (T_m,coll) can be assumed to be the average of the collector inlet and outlet temperatures. G is the solar irradiance, and T_a is ambient air temperature. Optical efficiency (a₀), linear heat loss coefficient (a₁), and the quadratic heat loss coefficient (a₂) depend on the solar collector type and are determined experimentally.

Solar collector heat output per area, Q_out,coll, can be obtained by Equation (2):

$${\mathrm{Q}}_{\mathrm{out},\mathrm{coll}}={\mathsf{\eta }}_{\mathrm{coll}}\mathrm{G}$$

(2)

The Q_out,coll can be used in sizing the solar loop HEX, and also in the calculation of the initial solar collector area required. This will be explained further in Section 2.3.

Solar irradiance is inconsistent throughout a year, and a method is required to classify the solar irradiance profiles, and also identify the most common type of solar irradiance profiles. The identified profiles can be used to obtain the typical daily collector performance for the sizing of the collector area. The other type of profiles can also be used in targeting the possible excess heat produced and extra heat required from a backup utility. This will be explained further in Section 2.3.

Solar irradiance profiles can be classified and grouped based on their irradiance level, the variability, and the distribution over the time of day [17]. Each feature can be characterized by its specific index. The transmittance index (k) can be used to characterize the irradiance level, and it is given by [18]:

$$\mathrm{k}=\frac{{\mathrm{G}}_{\mathrm{d}}}{{\mathrm{G}}_{\mathrm{d},\mathrm{c}\mathrm{s}}}$$

(3)

where G_d is the measured daily solar irradiance and G_d,cs is the daily solar irradiance in clear sky conditions.

The variability index (VI) can be used to characterize the variability of solar irradiance. The VI is shown in Equation (4) and is defined as the ratio between the length of the solar irradiance profile and the length of the maximum enveloping clear sky solar irradiance profile [17]:

$$\mathrm{VI}=\frac{{{\displaystyle \sum }}_{\mathrm{t}=2}^{\mathrm{n}}\sqrt{{\left({\mathrm{G}}_{\mathrm{t}}-{\mathrm{G}}_{\mathrm{t}-1}\right)}^2+\Delta {\mathrm{t}}^2}}{{{\displaystyle \sum }}_{\mathrm{t}=2}^{\mathrm{n}}\sqrt{{\left({\mathrm{G}}_{\mathrm{c}{\mathrm{s}}_{\mathrm{t}}}-{\mathrm{G}}_{\mathrm{c}{\mathrm{s}}_{\mathrm{t}-1}}\right)}^2+\Delta {\mathrm{t}}^2}}$$

(4)

where Δt is the time interval in minutes, n is the number of time intervals of the day, and G_cs is the clear sky solar irradiance. The solar irradiance in clear sky condition in Equations (3) and (4) can be calculated by the method of Hottel [19].

The morning fraction (f_m) can be used to characterize the time distribution of solar irradiance profile in the day. The f_m is defined as the ratio between the total solar irradiance from sunrise to solar noon and the total solar irradiance in the whole day, and it is determined using [17]:

$${\mathrm{f}}_{\mathrm{m}}=\frac{{\mathrm{G}}_{\mathrm{n}}}{{\mathrm{G}}_{\mathrm{d}}}$$

(5)

where G_d,m is the measured solar irradiance from sunrise to the solar noon. The time fraction (f_t) can be used to segment morning fraction (f_m) into three groups for the clustering/grouping of solar irradiance profile, and the f_t is given by [17]:

$$\begin{gathered} {\mathrm{f}}_{\mathrm{t}}=1\hspace{1em}\mathrm{i}\mathrm{f}{\mathrm{f}}_{\mathrm{m}}\le 0.3 \\ {\mathrm{f}}_{\mathrm{t}}=0.5\hspace{1em}\mathrm{i}\mathrm{f}0.3\le {\mathrm{f}}_{\mathrm{m}}\le 0.7 \\ {\mathrm{f}}_{\mathrm{t}}=0\hspace{1em}\mathrm{i}\mathrm{f}{\mathrm{f}}_{\mathrm{m}}\ge 0.7 \end{gathered}$$

(6)

The calculated value of k and VI can also be segmented into three groups based on the similar criteria definition of f_t in Equation (6), and it is given by Equation (7):

$$\begin{array}{cc}\mathrm{Low}\left(\mathrm{L}\right)\hfill & \mathrm{i}\mathrm{f}\mathrm{k}\le 0.3\left({\mathrm{k}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{t}}-{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\right)+{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\\ \mathrm{Medium}\left(\mathrm{M}\right)\hfill & \begin{array}{ll}\mathrm{i}\mathrm{f}0.3({\mathrm{k}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{t}}& -{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}})+{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\le \mathrm{k}\\ & \le 0.7\left({\mathrm{k}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{t}}-{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\right)+{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\end{array}\\ \mathrm{High}\left(\mathrm{H}\right)\hfill & \mathrm{i}\mathrm{f}\mathrm{k}\ge 0.7\left({\mathrm{k}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{t}}-{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\right)+{\mathrm{k}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}}\end{array}$$

(7)

The k in Equation (7) can be replaced by VI for segmentation of the calculated VI value.

To help in the identifying of solar irradiance profile with its features, a letter is given to the profile for each type of the features. For k and VI, the letters are High (H), Medium (M), and Low (L) as shown in Equation (7). For f_t, the letters are Morning (M) when f_t = 0, Balanced (B) when f_t = 0.5, and Afternoon (A) when f_t = 1.

Table 1. k, VI and ft values corresponding to the group.

k-VI-ft Group	Transmittance Index, k	Variability Index, VI	Time Fraction, ft
Low-Medium-Morning (LMM)	0.278	1.355	0
Medium-Low-Balanced (MLB)	0.530	1.036	0.5
High-High-Afternoon (HHA)	0.739	1.753	1

shows some example for k, VI, and f_t values of the solar irradiance profile corresponding to the group; the first letter refers to the k value, the second letter refers to the VI value, and the last letter refers to the f_t value.

The frequency of each type of group can be recorded, and the group with the highest frequency is identified as the most common type of solar irradiance profile. The arithmetic mean of the most common solar irradiance profile can be calculated and used for the calculation in Equations (1) and (2). Figure 3 summarizes the main inputs and the connections between the equations used in this sub-section.

2.3. Calculation of Solar Collector Area and Storage Volume

After identifying some potential streams for solar thermal integration, the initial solar collector area required can be determined based on the daily process heat demand, H_d, of the stream and the Q_out,coll from Section 2.2. Equation (8) illustrates the calculation of the initial solar collector area:

$${\mathrm{A}}_{{\mathrm{coll}}_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}}=\frac{{\mathrm{H}}_{\mathrm{d}}}{{{\displaystyle \sum }}_{\mathrm{t}}{\mathrm{Q}}_{\mathrm{o}\mathrm{u}\mathrm{t},{\mathrm{coll}}_{\mathrm{t}}}}$$

(8)

The initial solar collector area can be used to form a heat storage profile based on the cascade analysis proposed by Ho et al. [20]. The heat storage profile shows the load of solar heat supply and discharging of energy from storage in a time period, which in this study is 24 hours.

As seen in the configuration, the solar heat will first charge the storage. The heat will then discharge from the storage to the process. Equations (9) and (10) can be used to determine heat for charging (E_{in_t}) and discharging (E_{out_t}) of the storage:

$${\mathrm{E}}_{\mathrm{i}{\mathrm{n}}_{\mathrm{t}}}={\mathrm{A}}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{\mathrm{Q}}_{\mathrm{o}\mathrm{u}\mathrm{t},\mathrm{c}\mathrm{o}\mathrm{l}{\mathrm{l}}_{\mathrm{t}}}{\mathsf{\eta }}_{\mathrm{sto}}$$

(9)

$${\mathrm{E}}_{\mathrm{o}\mathrm{u}{\mathrm{t}}_{\mathrm{t}}}=\frac{{\mathrm{H}}_{\mathrm{t}}}{{\mathsf{\eta }}_{\mathrm{s}\mathrm{t}\mathrm{o}}}$$

(10)

where h_sto is the storage charging and discharging efficiency, including heat loss during the storage period and the charging/discharging cycle. H_t is the heat demand by the process, and subscript t refers to a specific time interval.

The heat surpluses and deficits (N_t) of the system can be estimated by Equation (11) based on heat charging and discharging of the storage [20]:

$${\mathrm{N}}_{\mathrm{t}}={\mathrm{E}}_{\mathrm{i}{\mathrm{n}}_{\mathrm{t}}}-{\mathrm{E}}_{\mathrm{o}\mathrm{u}{\mathrm{t}}_{\mathrm{t}}}$$

(11)

Heat surpluses (positive number) indicate the excess heat available, and heat deficits (negative number) indicate that the demand is higher than the heat generated.

Heat in the storage throughout the time interval (E_t) is given by Equation (12) [20]:

$${\mathrm{E}}_{\mathrm{t}}={\mathrm{E}}_{\mathrm{t}-1}+{\mathrm{N}}_{\mathrm{t}}$$

(12)

The results must be analyzed to ensure there is no negative number, as negative numbers indicate that the system is infeasible. If there is a presence of a negative number, the largest negative value of E_t is added to the beginning of the analysis (t = 0).

Table 2. Example of the initial heat storage profile.

Time	Heat Charging, ${\mathbf{E}}_{\mathbf{i}{\mathbf{n}}_{\mathbf{t}}}$ [kWh]	Heat Discharging, ${\mathbf{E}}_{\mathbf{o}\mathbf{u}{\mathbf{t}}_{\mathbf{t}}}$ [kWh]	Heat Surpluses and Deficits, Nt [kWh]	Infeasible Heat Storage Cascade, Et [kWh]	Feasible Heat Storage Cascade, Et [kWh]
00:00					895
00:00–01:00	0	100	−100	−100	795
01:00–02:00	0	100	−100	−200	695
02:00–03:00	0	100	−100	−300	595
03:00–04:00	0	100	−100	−400	495
04:00–05:00	0	100	−100	−500	395
05:00–06:00	0	100	−100	−600	295
06:00–07:00	0	100	−100	−700	195
07:00–08:00	0	100	−100	−800	95
08:00–09:00	5	100	−95	−895	0
09:00–10:00	139	100	39	−856	39
10:00–11:00	256	100	156	−701	194
11:00–12:00	346	100	246	−454	441
12:00–13:00	352	100	252	−202	693
13:00–14:00	352	100	252	50	945
14:00–15:00	267	100	167	216	1111
15:00–16:00	161	100	61	277	1172
16:00–17:00	67	100	−33	244	1139
17:00–18:00	0	100	−100	144	1039
18:00–19:00	0	100	−100	44	939
19:00–20:00	0	100	−100	−56	839
20:00–21:00	0	100	−100	−156	739
21:00–22:00	0	100	−100	−256	639
22:00–23:00	0	100	−100	−356	539
23:00–00:00	0	100	−100	−456	439

shows an example of the initial heat storage profile; the largest negative value of infeasible heat (895) is added to the initial value in the column of feasible heat storage cascade. The charging and discharging efficiencies of storage used in this example are assumed to be 90%.

The next criterion to be fulfilled is that the initial heat storage (895) must be equal to the final heat storage (439). If the heat storage at the end is less than at the start, there will be insufficient heat in some future time period, whereas if the heat storage at the end is more than at the beginning, the storage will store up excess heat every day and make the system unstable if the excess heat is not disposed. Goal seek analysis in Microsoft Word Excel can be used to adjust the solar collector area until the initial heat storage is equal to the final heat storage.

After there is no negative number and the heat storage is the same at the beginning and the end of the analysis, a new profile can be established with the new solar collector area. The final heat storage profile is presented in

Table 3. Example of the new heat storage profile.

Time	Heat Charging, ${\mathbf{E}}_{\mathbf{i}{\mathbf{n}}_{\mathbf{t}}}$ [kWh]	Heat Discharging, ${\mathbf{E}}_{\mathbf{o}\mathbf{u}{\mathbf{t}}_{\mathbf{t}}}$ [kWh]	Heat Surpluses and Deficits, Nt [kWh]	Infeasible Heat Storage Cascade, Et [kWh]	Feasible Heat Storage Cascade, Et [kWh]
00:00					894
00:00–01:00	0	794	−100	−100	795
01:00–02:00	0	694	−100	−200	695
02:00–03:00	0	594	−100	−300	595
03:00–04:00	0	494	−100	−400	495
04:00–05:00	0	394	−100	−500	395
05:00–06:00	0	294	−100	−600	295
06:00–07:00	0	194	−100	−700	195
07:00–08:00	0	94	−100	−800	95
08:00–09:00	6	0	−94	−894	0
09:00–10:00	171	71	71	−823	39
10:00–11:00	316	287	216	−607	194
11:00–12:00	428	615	328	−279	441
12:00–13:00	435	949	335	55	693
13:00–14:00	435	1,284	335	390	945
14:00–15:00	329	1,513	229	619	1111
15:00–16:00	198	1,611	98	717	1172
16:00–17:00	83	1,594	−17	700	1139
17:00–18:00	0	1,494	−100	600	1039
18:00–19:00	0	1,394	−100	500	939
19:00–20:00	0	1,294	−100	400	839
20:00–21:00	0	1,194	−100	300	739
21:00–22:00	0	1,094	−100	200	639
22:00–23:00	0	994	−100	100	539
23:00–00:00	0	894	−100	0	439

. The required storage capacity is the maximum value in the column of “Feasible Heat Storage Cascade”.

The solar collector area must also be checked against the available space in the process plant, which will be explained in Section 2.4. If there is insufficient space for the new solar collector, a project can still proceed, but extra heat from a backup utility system will be required. The extra heat from a backup utility plant can be targeted using Equation (13).

$${\mathrm{E}}_{\mathrm{u}}{}_{\mathrm{t}}={\mathrm{E}}_{\mathrm{t}-1}+{\mathrm{N}}_{\mathrm{t}}\hspace{1em}\mathrm{i}\mathrm{f} {\mathrm{E}}_{\mathrm{t}-1}+{\mathrm{N}}_{\mathrm{t}}<0$$

(13)

The effect of different solar irradiance profiles to the heat storage profile can be determined by applying different types of solar irradiance profiles based on the classification method presented in Section 2.2. This will help in the targeting of possible extra heat required (Equation (13)) or excess heat (Equation (14)) from each time period, and also help in the establishment of a predictive control scheme:

$${\mathrm{E}}_{\mathrm{e}\mathrm{x}}{}_{\mathrm{t}}={\mathrm{E}}_{\mathrm{t}-1}+{\mathrm{N}}_{\mathrm{t}}-{\mathrm{c}}_{\mathrm{sto}}\hspace{1em}\mathrm{i}\mathrm{f} {\mathrm{E}}_{\mathrm{t}-1}+{\mathrm{N}}_{\mathrm{t}}>{\mathrm{c}}_{\mathrm{sto}}$$

(14)

where c_sto is the storage capacity.

Figure 4 summarizes the main inputs and the connections between the equations used in this sub-section.

2.4. Space Constraints for Solar Collector Installation

The new solar collector area must be checked with the available empty area in the processing plant to determine whether the available area is sufficient with the installation of the required collector area, as mentioned in Section 2.3. In an actual solar collector installation, some distance (d) must be left between the collectors to avoid overshadowing of one solar collector to another. The concept from trigonometry [21] can be used to determine the availability of space for collector installation, as shown in Figure 5.

The dimensions of the solar collector, the inclination angle (β) of the solar collector, and the incident solar radiation angle (ε) will affect the distance, d, between the solar collectors required, and hence, limiting the area of solar collector that can be installed in the process plant. The incident solar radiation angle (ε) can be determined by using the Duffe and Beckman formula, which is presented by Elhab et al. as Equation (3) in their study [22]. D and d in Figure 5 are given by Equation (15) and Equation (16), respectively [21]:

$$\mathrm{D}=\frac{{\mathrm{W}}_{\mathrm{coll}}\sin \left(180°-\mathsf{\beta }-\mathsf{\epsilon }\right)}{\sin \mathsf{\epsilon }}$$

(15)

$$\mathrm{d}=\frac{{\mathrm{h}}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}}{\mathrm{Tan}\mathsf{\epsilon }}$$

(16)

There are two possible ways to install solar collectors, which will affect the total collector area and the total solar irradiation on the collectors. Figure 6 shows the possible orientations, where the blue lines represent the solar collector. L_b and W_b are the length and width of the available space in the building for solar collector installation.

The two possible ways in Figure 6 should be evaluated to obtain the maximum area utilization from the available space in the process plant. Within the available space, the maximum number of rows allowed for a collector can be determined by Equation (17):

$${\mathrm{Row}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{{\mathrm{L}}_{\mathrm{b}}-2{\mathrm{L}}_{\mathrm{e}}}{\mathrm{D}}+1$$

(17)

where L_e is the space left between the collector and the edge of the available space.

Next, the maximum number of collectors allowed per row is:

$${\mathrm{Collector}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{{\mathrm{W}}_{\mathrm{b}}-2{\mathrm{W}}_{\mathrm{e}}}{{\mathrm{L}}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}}$$

(18)

where L_coll is the length of a solar collector. The Row_max and Collector_max must be rounded down to the nearest whole number.

The actual area available for the installation of a collector, A_max, is:

$${\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\left(\left(\mathrm{D}-\mathrm{d}\right)\mathrm{R}\mathrm{o}{\mathrm{w}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)\left({\mathrm{L}}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}\mathrm{C}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{o}{\mathrm{r}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(19)

Note that Equations (17)–(19) are feasible for Case 1 only. For Case 2, L_b and L_e are switched to W_b and W_e. On the other hand, W_b and W_e are switched to L_b and L_e.

An area utilization ratio, r_Area, can be calculated by Equation (20) for both cases to identify the case with more utilized space:

$${\mathrm{r}}_{\mathrm{A}\mathrm{r}\mathrm{e}\mathrm{a}}=\frac{{\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{\mathrm{L}}_{\mathrm{b}}{\mathrm{W}}_{\mathrm{b}}}$$

(20)

Case 1 and Case 2 in Figure 6 show a different orientation of the solar collector surface, and a different orientation of the solar collector surface will result in a different solar yield. Solar yield with different orientations must be determined and compared to choose the case of installation with the most solar yield. Equation (21) can be used to determine the average daily solar irradiation on a titled surface, ${\mathrm{i}}_{\mathrm{t}\mathrm{i}}$:

$${\mathrm{i}}_{\mathrm{t}\mathrm{i}}={\mathrm{r}}_{\mathrm{t}\mathrm{i}}\mathrm{i}$$

(21)

where i is the average daily irradiation on a horizontal surface, and r_ti is defined as the ratio of average daily irradiation on a titled surface to that on a horizontal surface. r_ti can be estimated by Equation (22) assuming the diffuse and reflected radiation is uniformly distributed over the sky (isotropic) [23]:

$${\mathrm{r}}_{\mathrm{t}\mathrm{i}}=\left(1-\frac{{\mathrm{i}}_{\mathrm{D}}}{\mathrm{i}}\right){\mathrm{r}}_{\mathrm{t}\mathrm{i},\mathrm{b}\mathrm{e}}+\frac{{\mathrm{i}}_{\mathrm{D}}}{\mathrm{i}}\left(\frac{1+\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }}{2}\right)+\mathsf{\rho }\left(\frac{1-\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }}{2}\right)$$

(22)

where i_D is the average daily diffuse irradiation, r_ti,be is the ratio of average beam irradiation onto the titled surface to that onto a horizontal surface, and ρ is the ground reflectance. If the measurement of i_D is not available, the Equation (3) developed by Page [24] can be used to determine the value of i_D. r_ti,be can be determined for any tilt and orientation of a surface using [25]:

$$\begin{array}{ll}{\mathrm{r}}_{\mathrm{t}\mathrm{i},\mathrm{b}\mathrm{e}}& =\left[\right(\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\delta }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\varphi })\left({\mathsf{\omega }}_{\mathrm{s}\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{s}\mathrm{r}}\right)\left(\mathsf{\pi }/180\right)\\ & \hspace{1em}\hspace{1em}\hspace{1em} -\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\delta }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\varphi }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\beta }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\gamma }\right)\left({\mathsf{\omega }}_{\mathrm{s}\mathrm{s}}-{\mathsf{\omega }}_{\mathrm{s}\mathrm{r}}\right)\left(\mathsf{\pi }/180\right)\\ & \hspace{1em}\hspace{1em}\hspace{1em} +\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\varphi }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\delta }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\beta }\right)\left(\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\omega }}_{\mathrm{s}\mathrm{s}}+\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\omega }}_{\mathrm{s}\mathrm{r}}\right)\\ & \hspace{1em}\hspace{1em}\hspace{1em} +\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\delta }\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\gamma }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\varphi }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\beta }\right)\left(\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\omega }}_{\mathrm{s}\mathrm{s}}-\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\omega }}_{\mathrm{s}\mathrm{r}}\right)\\ & \hspace{1em}\hspace{1em}\hspace{1em} -\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\delta }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\beta }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\gamma }\right)\left(\mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\omega }}_{\mathrm{s}\mathrm{s}}-\mathrm{c}\mathrm{o}\mathrm{s}{\mathsf{\omega }}_{\mathrm{s}\mathrm{r}}\right)]\\ & \hspace{1em}\hspace{1em}\hspace{1em} /\left[2\right(\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\varphi }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\delta }\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\omega }}_{\mathrm{s}}+\left(\mathsf{\pi }/180\right){\mathsf{\omega }}_{\mathrm{s}}\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\varphi }\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\delta }\left)\right]\end{array}$$

(23)

where γ is the surface azimuth angle or orientation, ϕ is the latitude of the location, and δ is the solar declination, which can be approximately estimated by Equation (8) established by Cooper [26]. ω_sr and ω_ss are sunrise and sunset hour angles on the tilted surface and are given by Equations (12)–(15) presented by Klein [25], and ω_s is the sunset hour angle which can be determined by Equation (24).

$${\mathsf{\omega }}_{\mathrm{s}}={\cos }^{-1}(\mathrm{t}\mathrm{a}\mathrm{n}\mathsf{\varphi }\mathrm{t}\mathrm{a}\mathrm{n}\mathsf{\delta })$$

(24)

Total solar irradiation on the titled solar collector can be calculated from the multiplication of ${\mathrm{i}}_{\mathrm{t}\mathrm{i}}$ and A_max from Equation (19).

The case with more solar irradiation on the titled solar collector should be selected, and the A_max of that case can be used to check whether the available space in the process plant is sufficient for the proposed new solar collector area required from Section 2.3. Figure 7 summarizes the main inputs and the connections between the equations used in this sub-section.

2.5. Stagnation Control of Solar Collectors

Stagnation occurs when the solar energy absorbed by a collector exceeds the capability of its heat transfer fluid circuit to cool it, resulting in excessive collector temperatures. During this period, the absorber plate of the solar collector will increase in temperature until it is high enough to reject heat to the surrounding ambient through the housing. The magnitude of absorber temperature is referred to as the stagnation temperature. This stagnation temperature is the highest temperature that can be handled by solar collector units [27]. The magnitude of stagnation temperature can be estimated by solving the quadratic Equation (25) from Harrison and Cruickshank [28] and is highly depends on the heat loss coefficient of the solar collector. The stagnation temperature can be used for further design of the cooling system or so-called stagnation control. Fluid stagnation and overheating can potentially damage a solar collector as well as system components.

$${\mathrm{a}}_2{\left({\mathrm{T}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}-{\mathrm{T}}_{\mathrm{a}}\right)}^2+{\mathrm{a}}_1\left({\mathrm{T}}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}-{\mathrm{T}}_{\mathrm{a}}\right)-{\mathrm{a}}_0\mathrm{G}=0$$

(25)

T_stag is stagnation temperature, G is the amount of solar irradiance, and T_a is ambient air temperature. The optical efficiency (a₀), linear heat loss coefficient (a₁), and quadratic heat loss coefficient (a₂) depend on the solar collector type and are determined experimentally.

Maximum solar irradiance and ambient temperature data should be used to determine the maximum stagnation temperature that can be reached by a solar collector. The temperature can be used to design the cooler HEX for the stagnation control. Figure 8 shows the concept of the cooler HEX design. The cooling water loop is activated by the control system in case the solar collector loop is overheating, and it is one of the methods of Frank et al. [27] for reducing the temperature of solar collector loop.

Equation (1) from Section 2.2 shows that the efficiency of the solar collector also depends on the average collector internal fluid temperature (T_m,coll). Therefore, the efficiency of the same solar collector will be different during stagnation condition. A new heat output per area of solar collector, Q_out,coll, should be calculated. The heat transfer rate of the cooler unit can be determined by multiplying the Q_out,coll with the area of solar collector installed.

The required cooler area can be determined and it contributes to the techno-economic analysis for their cost estimation. The techno-economic analysis will be explained further in Section 2.6. Figure 9 summarizes the main inputs and the connections between the equations used in this sub-section.

2.6. Techno-Economic Analysis

Based on the known solar collector area and storage volume from Section 2.3, the purchase cost of solar collectors and storage can be calculated by multiplying with their corresponding unit price.

2.6.1. Heat Exchanger (HEX)

The purchase cost for a HEX is given by Equation (26) [29]:

$${\mathrm{C}}_{\mathrm{P}}={\mathrm{F}}_{\mathrm{p}}{\mathrm{F}}_{\mathrm{M}}{\mathrm{F}}_{\mathsf{\lambda }}\left(\frac{{\mathrm{I}}_{\mathrm{i}}}{{\mathrm{I}}_{{\mathrm{B}}_{\mathrm{i}}}}\right){\mathrm{C}}_{\mathrm{B}}$$

(26)

where I_i is the cost index at a later year and I_{B_i} is the base cost index corresponding to the year that applies to the purchase cost. The index used in this study is the Chemical Engineering Index as it is one of the commonly used indices [30]. F_p is the pressure factor based on the design pressure of HEX. F_M is a material factor for various combinations of shell and tube materials. F_λ is the tube-length correction of the HEX. The base cost, C_B for the HEX is based on the cost correlations of Corripio et al. [31].

2.6.2. Centrifugal Pump

The purchase cost for a centrifugal pump can be estimated by Equation (27) [29]:

$${\mathrm{C}}_{\mathrm{P}}={\mathrm{F}}_{\mathrm{T}}{\mathrm{F}}_{\mathrm{M}}\left(\frac{{\mathrm{I}}_{\mathrm{i}}}{{\mathrm{I}}_{{\mathrm{B}}_{\mathrm{i}}}}\right){\mathrm{C}}_{\mathrm{B}}$$

(27)

where F_T is the pump-type factor and F_M is the material factor of the centrifugal pump. The base cost for the centrifugal pump, C_B, is based on the cost correlations of Seider et al. [29], which were adopted from Corripo et al. [32].

2.6.3. Electric Motor

An electric motor is usually used to drive a centrifugal pump, and its purchase cost is given by Equation (28) [29]:

$${\mathrm{C}}_{\mathrm{P}}={\mathrm{F}}_{\mathrm{T}}\left(\frac{{\mathrm{I}}_{\mathrm{i}}}{{\mathrm{I}}_{{\mathrm{B}}_{\mathrm{i}}}}\right){\mathrm{C}}_{\mathrm{B}}$$

(28)

where F_T is the motor-type factor, and the base cost for an electric motor is based on the electric motor cost correlations of Corripio et al. [32].

2.6.4. Total Capital Cost

The total capital cost of equipment can be calculated by Equation (29), which according to the method developed by Lang [33,34,35], requires the sum of the purchase cost of the equipment, C_P:

$${\mathrm{C}}_{\mathrm{T}\mathrm{C}\mathrm{I}}=1.05{\mathrm{F}}_{\mathrm{L}}{\displaystyle \sum }_{\mathrm{i}}{\mathrm{C}}_{{\mathrm{P}}_{\mathrm{i}}}$$

(29)

where 1.05 is the estimation for the delivery cost of equipment and F_L is the recommended Lang factor that was detailed by Peters et al. [36]. The detailed breakdown of the Lang factor by Peters et al. assumes that major plant additions are made to an existing site. Although the detailed breakdown by Peters et al. includes all the equipment cost, direct plant cost, indirect plant cost, and working capital, a first-order capital cost estimate for solar thermal integration can be made by using equipment-specific Lang factors.

2.6.5. Levelized Cost of Heat

To measure the economic performance of the proposed solar thermal integration, the levelized cost of heat (LCoH) proposed by Task 54 [37] can be used, as presented in Equation (30):

$$\mathrm{LCoH}=\frac{{\mathrm{C}}_{\mathrm{T}\mathrm{C}\mathrm{I}}+{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}\frac{\mathrm{O}\&{\mathrm{M}}_{\mathrm{t}}}{{\left(1+\mathrm{R}\right)}^{\mathrm{t}}}}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}\frac{{\mathrm{E}}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{}_{\mathrm{t}}}{{\left(1+\mathrm{R}\right)}^{\mathrm{t}}}}$$

(30)

where O&M_t is the annual operation and maintenance cost, R is the discount rate, and E_{coll_t} is the annual heat generated by a solar collector. The E_{coll_t} is given by Equation (31):

$${\mathrm{E}}_{\mathrm{c}\mathrm{o}\mathrm{l}{\mathrm{l}}_{\mathrm{t}}}=365{\mathrm{A}}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{\displaystyle \sum }_{\mathrm{t}}{\mathrm{Q}}_{\mathrm{c}\mathrm{o}\mathrm{l}{\mathrm{l}}_{\mathrm{t}}}$$

(31)

Figure 10 summarizes the main inputs and the connections between the equations used in this sub-section.

2.7. Solar Collector Area Optimization

The solar collector area (A_coll) can be adjusted to indicate its effect on the LCoH in Equation (30). Increasing the solar collector area will increase the annual heat generated (E_{coll_t}) and heat savings, but it will also increase the capital cost (C_TCI) and the annual operation and maintenance cost (O&M_t). This trade-off needs to be balanced to minimize the total life cycle cost.

There are some constraints when selecting the solar collector area. First, there is a minimum value of solar collector area (Section 2.3) to make the system feasible and stable. Second, there is a maximum value of solar collector area that is governed by the available space (Section 2.4).

Increasing the solar collector area will increase the amount of heat, likely increasing the required heat storage volume. To stabilize the system, the net daily heat gain can be determined by deducting the final heat storage with the initial heat storage. The calculated value of the net heat gain can help establish a predictive control scheme to utilize the excess heat for other process heat demands. Another way to stabilize the system is to reject the excess heat from the storage by running the solar collector loop during nighttime to dissipate heat to the environment [38] or by using a cooling tower. The optimization is followed by a sensitivity analysis of various parameters for a range of LCoH outcomes.

The collector area that results in the lowest levelized cost of heat is selected, and the final storage volume is determined based on the cascade analysis proposed by Ho et al. [20] and described previously in Section 2.3. The required storage capacity is the maximum value in the column of “Feasible Heat Storage Cascade” in the new heat storage profile.

3. Case Study

Malaysia is located in the equatorial region and has a tropical rainforest climate where the weather is hot and humid all year. That means Malaysia receives higher solar radiation than most of the other countries in the world; the monthly solar irradiation received in Malaysia is around 400–600 MJ/m² [39]. Generating heat energy from solar thermal energy system in Malaysia is highly attractive. The industrial sector of Malaysia contributed to 26.2% of the total final energy demand in 2013, which was approximately 2160 PJ [40]. Among the total energy used in industry, 67% is contributed by heat energy. Suresh and Rao [41] reported that with current collector technologies, low-temperature heat demands are most suitable for the integration of solar heat. Solar energy has a significant industrial market in Malaysia given that 30% of the industrial process heat requires temperatures below 100 °C [42].

As a demonstration of the proposed method, a case study of solar thermal integration in a milk powder plant was investigated with the solar irradiance profile in Malaysia. The selected case study was first presented by Atkins et al. [7].

Table 4. Stream data for the milk powder plant.

Stream Name	Stream Type	Supply Temperature [°C]	Target Temperature [°C]	Heat Flow [kW/°C]
Raw milk	Cold	10	43	115.8
Skim milk	Hot	45	10	110.6
Cream	Cold	45	80	4.8
Cream A	Hot	80	10	4.8
Cleaning (CIP)	Cold	15	80	25.1
Dryer feed (product)	Cold	48	80	21.0
Cow water	Hot	45	15	118.5

shows the stream list for the relevant hot and cold streams. The Clean-In-Place (CIP) stream is the water used for cleaning the process equipment and the cow water is the water that evaporated from the milk during the process.

Solar irradiance and ambient temperature data from May 2018 until April 2019 were recorded at the weather station in the MJIIT building (UTM KL). For simplicity, the processes were considered to operate continuously for 8000 h/y. The output temperature of the solar thermal system was set to the exact process temperature required. Here, the aim is to optimize the design of solar heat integration for different process options by minimizing LCoH using the proposed method.

4. Results and Discussion

This section presents and discusses some important results obtained from the integrated design method, and also the input variables. After the optimization of LCoH, the ranges of LCoH results from the sensitivity analysis are presented.

4.1. Determination of Potential Stream

Figure 11 shows the Grand Composite Curve (GCC) for all the streams in Table 4, with a ΔT_min of 5 °C. Based on the GCC, a pinch temperature of 42.5 °C, and hot and cold utility targets of 2023 and 3493 kW, respectively. From the GCC, the potential solar operating temperature is between 42.5 and 82.5 °C, and the heat demand on this temperature range is 2023 kW. This heat demand is composed of heating of the cream, cleaning (CIP), and dryer feed. There is one HEX network design with solar integration on the CIP stream by [7], as shown in Figure 12.

The GCC is drawn with an assumption of full exploitation of the thermodynamic heat recovery potential with all other streams based on a 5 °C minimum approach temperature. In the first instance, one should design an efficient HEX network to minimize utility use and then consider the most effective use of solar heating. However, the primary purpose of this case study is to demonstrate the solar integration design method, and it is assumed that the system already achieves the maximum potential heat recovery. Therefore, the design in Figure 12 was not modified, and further analysis was done on this solar integration and two additional potential streams of cream and dryer feed, despite their possible constraint in solar integration.

Table 5. Stream data for the potential streams.

Option	Stream Name	Supply Temperature [°C]	Target Temperature [°C]	Mass Flow [kg/s]	Heat Flow [kW/°C]	Heat Load [kW]
1	Cream	45	80	1.2	4.8	168
2	Dryer feed (product)	48	80	5.3	21	672
3	Cleaning (CIP)	40	80	6.0	25.1	1,004

shows the stream data of the three potential streams. In the remainder of the article, the streams will be referred to by their stream number.

4.2. Determination of Solar Collector Area

Indices k, VI, and f_t from Equations (3–7) were calculated for every solar irradiance profile at 1 h intervals. All of the solar irradiance profiles were grouped, and the top three types of solar irradiance profiles are shown in

Table 6. Top three types of solar irradiance profile in terms of their frequency of appearance.

k-VI-ft Group	Appearance Percentage (%)
Medium-Medium-Balanced (MMB)	50.5
High-Medium-Balanced (HMB)	12.9
Medium-High-Balanced (MHB)	7.8

, in terms of their frequency of appearance.

However, to demonstrate the effect of different solar irradiance level on the heat output of solar collector and heat storage profile, LMB solar irradiance profiles were used instead of MHB. LMB solar irradiance profiles have an appearance percentage of 4.1%. The arithmetic means of solar irradiance and ambient temperature data of MMB, HMB, and LMB are shown in Figure 13.

A flat plate solar collector (FPC) was selected for this case study, as its low design temperature range is suitable for the selected stream [43]. The FPC used in this study has an optical efficiency (a₀) of 0.817, linear heat loss coefficient (a₁) of 2.205, and 0.014 as its quadratic heat loss coefficient (a₂) [44]. The charging and discharging efficiencies of sensible heat storage for this study are 90% [45].

Table 7. Results of solar collector area for each option.

Option	Daily Heat Demand, Hd [kWd]	Solar Collector Heat Output per Area, Qout,coll [kWh/m2]	Initial Solar Collector Area, ${\mathbf{A}}_{{\mathbf{coll}}_{\mathbf{initial}}}$ [m2]	New Solar Collector Area [m2]
1	4032	2.09	1930	2382
2	16,128	2.05	7869	9715
3	24,096	2.16	11,157	13,774

shows the results of solar collector area for all options based on the MMB profiles which has the highest appearance percentage.

4.3. Space Availability for Solar Collector Installation

Since there is no information about the available space in the case study from [7], it is assumed to be a rectangular shape with a length of 140 m and a width of 100 m.

Table 8. Data for space availability analysis.

Collector length (m) [44]	5.96
Collector width (m) [44]	2.27
Collector height (m) [44]	0.14
Collector inclination angle (°) [46]	3.2
Incident solar radiation angle (°) [22]	25.95
Ground reflective [47]	0.2

shows the data input to the space availability analysis. The results of the maximum area utilization ratio from both arrangements shown in Figure 6 are shown in

Table 9. Space availability analysis.

	Case 1	Case 2
Maximum number of rows allowed for collector, Rowmax	55	39
Maximum number of collectors per row, Collectormax	16	23
Length of space allowed for collector [m]	124.66	137.08
Width of space allowed for collector [m]	95.36	88.39
Area that can be installed with collector, Amax [m2]	11,887	12,117
Area utilization ratio, rArea	0.85	0.87

.

Table 10. Calculation of solar irradiation with different orientations.

	Case 1		Case 2
Surface azimuth angle, ϕ [°]	0	180	90	−90
Ratio of average beam irradiation onto the titled surface to that onto a horizontal surface, rti,b	1.0403	0.9566	0.9936	0.9984
Ratio of average daily irradiation on a titled surface to that on a horizontal surface, rti	1.0010	0.9977	0.9992	0.9993
Average daily irradiation on a horizontal surface, i [MJ/m2]	1.205
Average daily solar irradiation on an inclined surface, iti [MJ/m2]	1.2062	1.2022	1.2040	1.2042
Area that can be installed with collector, Amax [m2]	11,887		12,117
Total solar irradiation on the titled solar collector [MJ]	14,338	14,291	14,588	14,591

shows the results of the calculation of solar irradiation with different orientations. The orientations are assumed 0 ° (south) and 180 ° (north) for Case 1, and 90 ° (west) and −90 ° (east) for Case 2.

Although Case 1 has the most r_ti with a surface azimuth angle of 0 °, the larger A_max in Case 2 still results in more solar irradiation. Since Case 2 has the most solar irradiation with a surface azimuth angle of −90 °, the A_max of Case 2 (12,117 m²) was selected as the available area for solar collector installation. After checking with the new solar collector area for all options, only the area for Option 3 is more than A_max, which means that the area for Option 3 will be 12,117 m² instead of 13,943 m². Figure 14 shows the feasible heat storage cascades in storage for Options 1, 2, and 3. The new heat storage profiles for all the options are presented in Appendix A.

Figure 14 shows that the final heat storage for Option 3 (6533 kWh) is less than the initial heat storage (10,040 kWh) over a day due to insufficient solar collector area. In some future time period, the extra heat from the backup utility system (Equation (13)) will be required to make up the difference. Required storage volumes for Options 1, 2, and 3 are 78 m³, 341 m³, and 383 m³, respectively. The heat storage profiles established from the different types of solar irradiance profiles are presented in Appendix B.

4.4. Stagnation Control of Solar Collector

The maximum solar irradiance and ambient temperature were obtained from data with a time step of 1 minute. Water is assumed as the heat transfer fluid for the preliminary calculation of cooler HEX.

Table 11. Stagnation temperature of solar collector and design parameters of cooler HEX.

Option	1	2	3
Max solar irradiance [Wh/m2]	1360
Max ambient temperature [°C]	36.7
Max stagnation temperature, Tstag [°C]	253.3
Design parameters of cooler HEX
Hot stream inlet temperature, Th_in [°C]	253.3
Hot stream outlet temperature, Th_out [°C]	223.35
Cold stream inlet temperature, Tc_in [°C]	27
Cold stream outlet temperature, Tc_out [°C]	57
Hot and Cold stream heat capacity [kJ/kgK]	4.2
Qout,coll of collector when stagnation [kWh/m2]	0.087
Heat transfer rate of HEX [kW]	281	1144	1622
Mass flow rate of hot and cold stream [kg/s]	2.23	9.08	12.87
Area of HEX, AHEX [m2]	220	897.3	1272

shows the result of stagnation temperature and the design parameters of the cooler HEX. Note that the stagnation temperature is the same for all options due to their usage of the same solar collector for integration.

4.5. Techno-Economic Analysis

The data inputs to the cost estimation are given in

Table 12. Data for cost estimation.

Unit price of flat plate solar collector (RM/m2) [43]	405
Unit price of sensible energy storage (RM/m3) [45]	593
Heat Exchanger
Material factor, FM (carbon steel/carbon steel)	1.00
Tube-length correction, Fλ (6.096 m)	1.00
Pressure factor, FP (1000 kPa [44])	1.01
Chemical Engineering Index of 2006, IB_i [29]	500
Chemical Engineering Index of September 2017, Ii [49]	606.9
Pump
Pump-type factor, FT	1.00
Material factor, FM (cast iron)	1.00
Chemical Engineering Index of 2006, IB_i [29]	500
Chemical Engineering Index of September 2017, Ii [49]	958.0
Electric Motor
Motor-type factor, FT	1.00
Chemical Engineering Index of 2006, IB_i [29]	500
Chemical Engineering Index of September 2017, Ii [49]	521.9
Levelized cost of heat
Number of operators required per day	6
Operator salary rate (RM/h) [50]	8.50
Discount rate (%) [51,52]	8

, and the Lang factors are presented in

Table 13. Lang factor, F_L [36].

Solids–Fluids Processing Plant
Delivered cost of equipment	1.00
Installation	0.39
Instrumentation and control	0.26
Piping	0.31
Total Lang factor, FL	1.96

. The HEX assumes counter-current flow and plate HEX. For the calculation of HEX, a heat transfer coefficient, U, is assumed the same in all processes. A centrifugal pump type is selected with a rate of rotation of 3600 rpm and a vertical split case. The electric motor type for the pump is an open, drip-proof enclosure motor. The O&M_t was calculated based on the cost sheet prepared by Busche [48], which was then modified by Seider et al. [29]. The annual maintenance cost is given by 10.35% of the total capital cost for the solids–fluids handling process [29]. The annual operations cost is given by 121% of the annual operator salary [29]. Techno-economic analysis results before optimization are presented in

Table 14. Techno-economic analysis results before optimization.

Option	1	2	3
Solar collector purchase cost [RM]	965,019	3,935,331	4,908,235
Storage purchase cost [RM]	46,143	202,446	227,085
Total HEX purchase cost [RM]	444,144	1,353,918	1,764,959
Total pump purchase cost [RM]	98,769	110,809	116,748
Stagnation control purchase cost [RM]	198,852	431,183	506,412
Total purchase cost, CP [RM]	1,752,926	6,033,687	7,523,439
Total capital cost, CTCI [RM]	3,607,522	12,417,327	15,483,238
Total annual O&M cost [RM/y]	867,059	1,778,873	2,096,195
Levelized cost of heat, LCoH [RM/kWh]	0.98	0.66	0.62
Levelized cost of heat, LCoH [EUR/kWh]	0.21	0.14	0.13

based on a project lifetime of 15 years.

4.6. Solar Collector Area Optimization

Equation (30) determines the effect of increasing the solar collector area to the LCoH. Figure 15 shows the graph of LCoH against the solar collector area for the three options.

Note that the graphs for Options 1 and 2 were plotted from the minimum value of collector area from Section 4.2 until the maximum value of 12,117 m². Although the collector area for Option 3 before optimization has reached the maximum available area, Option 3 is still included in Figure 15 to shows the effect of increasing collector area on the LCoH.

The LCoH values for all three options are decreased when the collector area increased, indicating that increasing the collector area is economically beneficial. Despite using the same type of solar collector and equipment, there is still a difference of LCoH for all options with the same solar collector area. This is caused by the difference in process temperature required by streams. Table 5 shows that Option 2 has the highest temperature required, followed by Option 1, and Option 3 has the lowest temperature required. Equation (1) shows that the efficiency of a solar collector is higher with the lower process temperature; hence the annual heat generated by the solar collector, E_{coll_t}, is higher for Option 3. However, the amount of difference in LCoH shown in Figure 15 is insignificant, as there is only a small difference in the process temperature required.

As a result, the collector area for Options 1 and 2 increased up to 12,117 m², since it shows a trend of decreasing LCoH and has not reached the maximum available area.

Table 15. Techno-economic analysis results after optimization.

Option	1	2	3
Solar collector purchase cost [RM]	4,908,235 (66.7%)	4,908,235 (66.0%)	4,908,235 (65.2%)
Storage purchase cost [RM]	321,341 (4.4%)	273,825 (3.7%)	227,085 (3.0%)
Total HEX purchase cost [RM]	1,509,433 (20.5%)	1,630,043 (21.9%)	1,764,959 (23.5%)
Total Pump purchase cost [RM]	116,033 (1.6%)	114,603 (1.5%)	116,748 (1.6%)
Stagnation control purchase cost [RM]	506,412 (6.9%)	506,412 (6.8%)	506,412 (6.7%)
Total purchase cost, CP [RM]	7,361,454	7,433,118	7,523,439
Total capital cost, CTCI [RM]	15,149,872	15,297,356	15,483,238
Total annual O&M cost [RM/y]	2,061,692	2,076,956	2,096,195
Levelized cost of heat, LCoH [RM/kWh]	0.62	0.64	0.62
Levelized cost of heat, LCoH [EUR/kWh]	0.135	0.138	0.133

shows the techno-economic analysis results after optimization; the percentages inside the brackets indicate the percentage of the purchase cost for each piece of equipment. As mentioned in Section 2.7, increasing the collector area will increase the daily net heat gain. Figure 16 shows the feasible heat storage cascade of Options 1 and 2 based on their new heat storage profiles after the optimization.

From Figure 16, the net daily heat gain for Options 1 and 2 are 18,306 kWh and 4430 kWh, respectively. The storage size for Options 1 and 2 are 542 m³ and 462 m³, respectively.

The solar thermal integration for Option 3 gives the lowest LCoH among all of the options. However, Figure 14 shows that Option 3 has less final heat storage than initial heat storage. Thus backup utility heat is required in some future time period. On the other hand, solar thermal integrations for Options 1 and 2 do not require backup utility, but they have some net daily heat gain (Figure 16) that requires a control scheme for utilizing the excess heat for other usages in the process plant. The cost and convenience need to be compared between setting up a backup utility and a control scheme in order to decide the best potential stream for solar thermal integration, but this will not be discussed in this paper.

The results of the sensitivity analysis for a range of LCoH is presented in Figure 17. All the parameters are increased and decreased by a 15% for the sensitivity analysis and they are represented by green and red bar. The LCoH of Option 1 (0.62) was selected for the base value of the sensitivity analysis.

As shown in Figure 17, the LCoH is most sensitive to the collector price, collector efficiency, and the discount rate. The sensitivity of LCoH to the collector price and efficiency can be justified by the highest cost percentage of collector shown in Table 15. This indicates that the LCoH of solar integration can be greatly reduced if the solar collector technology can be improved in terms of cost and efficiency. The ranges of the LCoH values are 0.57 to 0.68 RM/kWh due to collector price sensitivity, 0.47 to 0.93 RM/kWh for collector efficiency sensitivity, and 0.57 to 0.69 RM/kWh for discount rate sensitivity.

The conclusions made from the results of optimization and sensitivity analysis are based on the configuration as shown in Figure 2 only, different configurations that integrate solar thermal without HEX or heat storage will affect the energy savings and costs of the solar thermal system, and hence affect the strategy of optimization. However, the proposed method can be adjusted slightly by changing the calculation of T_m,coll, heat storage profile, and cost to encounter different integration configurations. Other limitations of this study are that only 24 hours per day of demand requirements and the most common weather profile in Malaysia are used in this case study; as different feasible collector areas and storage volumes will be required, it can be investigated further by using different demand and weather profiles.

The proposed method is suitable to estimate the feasible physical dimensions, including collector area, storage size, and control design, while considering significant physical constraints (e.g., available space) and recovery of low-grade heat based on the common scheme of solar thermal integration as shown in Figure 2. After that, the method can be used to optimize the physical dimensions for minimal LCoH. Therefore, the proposed method has a potential in the industrial sector, as the industrial applications typically require higher temperatures and greater volumes than both residential and commercial applications, and it makes the stagnation control and availability of space analysis more significant in industrial applications.

The LCoH of a solar thermal system estimated by Wahed et al. [51] is 0.33 RM/kWh, which is lower than the LCoH values shown in Table 15, but the amount of heat demand encountered by Wahed et al. is much larger than the demand in this study. Based on the concept of economies of scale, this will cause a lower LCoH due to a larger heat production scale.

Companies in Malaysia can apply for import duty and sales tax exemptions on imported machinery, equipment, materials, spare parts, and consumables used directly in the generation process of solar energy. For locally purchased machinery, equipment, materials, spare parts, and consumables, a full exemption is given on sales tax. An exemption is given for a period of one year. Under the National Budget of 2009, this incentive has been extended to grant exemption of import duty and sales tax to third party distributors of the solar thermal systems [53].

5. Conclusions

A new systematic method was used to optimize the solar heat integration for different process options by minimizing the levelized cost of heat through synergistically combining different techniques from the literature. Consequently, a case study of solar thermal integration in a milk powder plant from the literature was investigated with multiple classified solar irradiance profiles in Malaysia. This method is useful in the estimation of the capital cost and levelized cost of heat of the solar integration considering important physical constraints (e.g., available space) and recovery of waste heat. As a result, the solar thermal integration with the CIP stream gives the lowest levelized cost of heat of RM 0.63/kWh (0.13 EUR/kWh), makes the CIP the best potential stream for solar thermal integration. However, backup utility heat is required for the CIP stream in some future time period. On the other hand, solar thermal integration for cream and dryer does not require backup utility, but they have some excess final heat storage that requires a control scheme that can utilize the excess heat for other usages in the process plant. The cost and convenience need to be compared between setting up a backup utility and a control scheme in order to determine the best potential stream for solar thermal integration. The process temperature required causes a difference in the levelized cost of heat of solar thermal integration, but the difference is insignificant. From the sensitivity analysis, it has been shown that collector price, collector efficiency, and discount rate are the critical parameters of solar thermal integration. In future research, more case studies with different process temperature ranges, heat loads, solar collector types, discount rates, and also different weather data durations can be further investigated for different optimization strategies. Besides that, the proposed method can be adjusted slightly by changing the calculation of T_m,coll, the heat storage profile, and the cost of encountering different integration configurations.