Model Based Optimization of Energy Consumption in Milk Evaporators

1. Introduction

Evaporation involves the conversion of a solvent into vapor, which is then removed from a solution or slurry. In most cases, water is the solvent used in evaporation systems. This process entails vaporizing a portion of the solvent to create a concentrated solution, thick liquor, or slurry [1].

Dairy manufacturers commonly employ concentration techniques to produce dairy products with higher levels of dry matter, increased value, reduced volume, and extended shelf-life [2]. Lowering the water activity and reducing transportation and storage costs are key benefits of dehydrating dairy products. This process involves converting a liquid product into a dry powder by removing nearly all the available water. However, dairy products are sensitive to heat, and their functional properties and digestibility can be negatively affected by excessive heat during the dehydration process. Therefore, a single water removal method cannot consistently achieve optimal performance. It thus, becomes necessary to employ multiple processing steps tailored to the specific properties of the material being processed, while considering both product quality considerations and processing costs [3].

Typically, the main steps involved in the production of milk powder include standardisation, homogenisation, pasteurisation, evaporation and drying [4]. Evaporation is a significant step in milk-powder production plants, serving not only to concentrate milk to the desired viscosity for subsequent spray drying but also to reduce the energy requirements during the spray drying process. In the evaporation stage, sterilized milk is concentrated under vacuum conditions at temperatures ranging from 40 to 70°C. This process leads to a significant increase in the total solids content, which typically rises from around 13% to 50% [5]. Vacuum conditions are used to mitigate the negative effects of heat on heat-sensitive milk components like fats and to prevent the degradation of essential nutrients, such as vitamins.

Milk powder production consists of many thermal processes, including evaporation and drying, and is responsible for 15% of the total energy use in the dairy industry [4]. In France, approximately 25% of the total energy consumption in food industry is attributed to the dairy sector [3]. Therefore, energy consumption attributed to the process of evaporation has a substantial impact on the cost of milk powder production.

To this end, there are several energy consumption optimization strategies in the evaporation of milk that have been discussed in the literature. The ones that are examined in this work involve the most commonly used in the industrial practice, such as the Mechanical Vapor Recompression (MVR) and the Thermal Vapor Recompression (TVR) technologies.

• In the evaporator unit, the excess heat generated by the secondary steam from the evaporator is typically released as waste heat. However, this waste heat can be effectively utilized to preheat the feed. An important feature of MVR technology is the utilization of the secondary steam cycle. MVR employs a mechanical fan, typically powered by electricity, to recompress low-pressure vapor to a slightly higher pressure and temperature.
• On the other hand, TVR utilizes a thermos-compressor, which employs high-pressure vapor to recompress low-pressure vapor to a slightly higher pressure and temperature. Numerous studies have demonstrated that multi-effect evaporation reduces energy consumption costs by enhancing steam economy. This is achieved by utilizing the secondary steam generated by the preceding effect as the heat source for the subsequent effect [6].

Jebson and Chen [7] assessed the effectiveness of falling film evaporators used in the New Zealand dairy sector for concentrating whole milk by calculating the kg steam utilized to kg water evaporated ratio and the heat transfer coefficient of each evaporator pass. The steam consumption of full and skim milk was comparable. Schuck et al. [3] introduced a methodology for assessing and comparing the energy consumption involved in the production of dairy and feed powders at various stages of the dehydration process. The findings of the study revealed that the energy consumption for fat-filled and demineralized whey powders were 9.072 kJ/kg and 15.120 kJ/kg, respectively.

Walmsley et al. [8] conducted a study on applying Pinch Analysis to an industrial milk evaporator case study to quantify the potential energy savings. By the appropriate placement of mechanical vapor recompression in a new improved two-effect milk evaporation system design, a 78% reduction to steam (6397 kW) at the expense of 16% (364 kW) more electricity use was achieved and the emissions reduction was 3416 t CO₂-e/y. Srinivasan et al. [9] studied the energy efficiency at India’s largest milk processing plant and proposed retrofits for improving the plant’s sustainability. The results reveal that exergy efficiency of certain units is very low (<20%) while significant improvements in energy efficiencies can be achieved through simple, low-cost retrofits to these units. Moejes [4] studied the possibilities of upcoming milk processing technologies such as membrane distillation, monodisperse-droplet drying, air dehumidification, radio frequency heating, and radio frequency heating paired with renewable energy sources such as solar thermal systems. It was illustrated that the combination of developing technologies has the potential to cut operational energy consumption for milk powder manufacturing by up to 60%.

Zhang et al. [10] conducted a study by simulating a "pseudo" milk composition using hypothetical components in a commercial process simulator. The purpose of their work was to model a falling film evaporator commonly (FFE) employed in milk powder production plants. The study demonstrated that commercial process simulators have the ability to accurately simulate dairy processes. Building upon this research, Munir et al. [11] further enhanced the capabilities of commercial process simulators, providing valuable insights for practicing engineers to identify potential process improvements in the dairy industry. Bojnourd et al. [12] developed two types of dynamic models for an industrial four-effect falling-film evaporator used to condense whole milk: lumped and distributed. The findings indicate that while the distributed model demonstrates slightly better predictive capabilities compared to the lumped model, the latter outshines in terms of performance due to its simpler structure and significantly reduced simulation time. Zhang et al. [5] developed models for two commonly used types of milk powder evaporators: a conventional five-effect FFE without MVR and a three-effect evaporator with MVR. Heat-recovery processes were incorporated into the models to enable a comparison of energy consumption between the two processes. The results revealed that a three-effect FFE with MVR could achieve a 60% reduction in energy consumption compared to a conventional five-effect evaporator.

Gourdon and Mura [13] created a modeling tool based on experimental correlations built at industrial scale circumstances. The complicated interaction between the generated vapor and the liquid flow is included in their model. The results show that pressure drop is important in evaporator performance because of its influence on saturation temperature. Daz-Ovalle et al. [14] provided a set of dynamic models to study fouling of falling-film evaporators by considering fouling thickness, film thickness, temperature, and solids mass percentage. Hu et al. [15] developed a model for a water-to-water FFE simulation, which was employed in water vapor heat pump systems. That study focused on an existing FFE with four working tubes.

Bouman et al. [16] conducted experiments with a one-tube evaporator using whole and skim milk to ascertain the heat transfer and pressure decrease in evaporator tubes. Based on the findings, a computer software was created to optimize the design of multistage falling-film evaporators for dairy products. Silveira et al. [17] investigated the evaporation of water and skim milk using a pilot-scale, single-stage falling-film evaporator. In comparison to water, a thicker and slower film was formed at the end of the skim milk concentration procedure. They concluded that the behavior of a product during the evaporation process cannot be predicted solely by the overall heat transfer coefficient, and that a wide range of information, such as residence time distribution, product viscosity, and surface tension, is required to understand the evaporation process. Gourdon et al. [18] examined the flow behavior of a dairy product under falling film evaporation. The effects of varying dry solids concentrations, flow rates, and driving temperature changes were investigated. All three factors were shown to have a significant impact on the flow characteristics. Mura et al. [19] investigated the absolute vapor flow pressure losses during dairy product falling film evaporation using an experimental internal tube evaporator setup, adjusting the co-flow input velocity and product dry solid content. They concluded that pressure losses are strongly dependent on the co-flowing vapor rate for dry solid content between 13 and 40%.

Wijck et al. [20] conducted an evaluation of tools used for dynamic modeling and supervisory multivariable control design of multiple-effect falling-film evaporators. They specifically focused on the NIZO four-effect evaporator as a case study. The research aimed to achieve improved process operation for industries, leading to economic benefits such as increased yield, enhanced product quality, reduced energy consumption, and minimized material waste. Sharma et al. [21] created an Excel-based multi-objective optimization tool based on the elitist non-dominated sorting genetic algorithm (NSGA-II) and tested it on benchmark tasks. It is then used for multi-objective optimization (MOO) of the design of a falling-film evaporator system for milk concentration, which includes a pre-heater, evaporator, vapor condenser, and steam jet ejector. Haasbroek et al. [22], conducted a study utilizing historical data from a FFE to develop models for control purposes, without requiring knowledge of the plant's physical dimensions. The results indicated that fuzzy predictive controllers and LQR control exhibited the highest performance, followed by cascade control, and finally, PI control. Galvản-Ángeles et al. [23] analyzed a thermo-compression evaporation method for milk. The suggested tool considers the cost optimization of the evaporation system while incorporating thermo-physical parameters of the foodstuff as a function of composition and temperature. The results revealed that the evaporation economy is proportional to the percentage of recycled steam and the location of the effect that recycles the steam, and inverse to the thermodynamic efficiency of the thermo-compressor.

This work aims at optimizing the milk evaporation process by minimizing steam consumption under the operational and product quality constraints. To this end, 5 different Cases of evaporator layouts are examined using a global system analysis (GSA) and an advanced optimization approach. GSA is employed to revise decisions to improve system robustness, validate the process and reduce parameter uncertainty. An uncertainty analysis allowed the investigation of the impact of design and operational decisions and environmental inputs on Key Performance Indicators (KPIs). Afterwards, Cases 1-4 are optimized under 5 different scenarios, either minimizing the cost of steam while maximizing product yield or minimizing the total annualized cost which is applicable to a new plant design or a capacity expansion option. Steam economy, energy consumption profile and heat transfer areas are assessed and compared. Moreover, this study investigates to what extend switching from milk powders to new products known as milk concentrates effects the energy consumption in the evaporation process. Thus, Cases 1-4 are optimized under different end-product specifications (30, 35 or 50% solid content). Finally, it is evaluated whether the use of MVR or TVR is more economical for the evaporation process, based on current steam and electricity prices, economic trends, and costs of steam generated from renewable energy sources.

The remaining of the article is structured as follows: Section 2 describes the material composition along with some of its properties. Section 3 provides a detailed presentation of the mathematical model that describes the operation of a falling film evaporator. Sectio 4 presents the studied evaporator layouts. The attained global system analysis and optimization results of each studied Case/evaporator layout are presented and thoroughly discussed in Section 5. Finally, Section 6 summarizes the concluded remarks emerging from the study.

2. Materials

This section provides information about milk properties and conditions. In the present study, the analysis of FFEs is applied to the milk evaporator process, while the thermo-physical properties of the food vary with composition. Two different composition profiles of whole milk before evaporation are examined (Table 1) while Table 2 summarizes equations for the calculation of properties in the mixture. Thermodynamic and physical properties such as specific heat capacity $c_p$, thermal conductivity $k$, and density, $\rho$, of ‘pseudo’ milk are modelled using the process simulator and average values are compared to the literature. The results are summarized in Table 3.

3. Process model

The concentration of milk in this study is focused on falling film evaporators. In this process, the milk, which is close to its boiling point, is introduced in a uniform manner at the upper section of the inner surface of a tube. These tubes are arranged side by side, fixed in place, and surrounded by a jacket. As the milk descends within each tube, it forms a thin film and undergoes boiling as a result of the heat exchanged with the steam. The concentrated liquid is collected at the lower part of the equipment, while the remaining portion is separated from the steam in a subsequent separator. In evaporators equipped with multiple effects, the concentrated liquid is pumped to the next stage, while the steam serves the purpose of heating the subsequent [2].

A dynamic model of a simple falling-film evaporator is developed in this Section. The model examines the flash calculations of liquid entering the evaporator, its distribution through distributor plates, and its evaporation as it flows downwards through pipes. The model has been developed based on the following assumptions:

• It is assumed that the product immediately reaches its boiling temperature once it passes above the distributor plate, with either steam flashing or condensing.
• All the heating steam condenses at its saturation temperature. The heat released during condensation is utilized to preheat and evaporate the feed stream.
• The pressure difference on the steam side is insignificant, thus resulting in a constant temperature on the steam side.
• The unit is assumed to be perfectly mixed at all times, implying there are no spatial variations in intensive properties within it.
• A steady state balance is applied to the tube flow, but the mass hold-up of liquid on the distributed plate is calculated.

3.1. Mass balance

The rate of change of the mass hold-up, $M$, of any given species in phase $p$ in the unit is given by Equation (4):

$$\frac{d{M}_{i,p}}{dt}=\sum _{j=1}^{N_{inlet}}{F}_{j,p}^{in}·w_{i,j,p}^{in}-F_p^{out}·w_{i,p}^{out}+\sum _{p_n\ne p\in P}{R}_{p/p_{n,}i} ,\forall i\in I,\forall p\in P$$

(4)

Where $F_{j,p}^{in}$is the mass flowrate of the j^th feed stream in phase p, $w_{i,j,p}^{in}$ is the mass fraction of species i in the j^th feed stream in phase p, $F_p^{out}$ is the mass flowrate of material phase p in the outlet stream, $w_{i,p}^{out}$ is the mass fraction of species i in the outlet stream of phase p and $R_{p/p_{n,}i}$ is the rate of mass transfer between the phase p and the phase p_n.

Assuming that the unit is perfectly mixed, the composition of any outlet stream is the same as the composition within the unit:

$$w_{i,p}^{out}=w_{i,p},\forall i\in I,\forall p\in P$$

(5)

Where $w_{i,p}$is the mass fraction of species i in phase p in the unit.

The composition of material within the unit is given by Equations (6),(7):

$$w_{i,p}=\frac{M_{i,p}}{M_{total,p}},\forall i\in I,\forall p\in P$$

(6)

$$M_{total,p}=\sum _{i\in I}{M}_{i,p},\forall i\in I,\forall p\in P$$

(7)

Where $M_{p,total}$is the total holdup of material in phase p within the unit.

3.2. Energy balance

The rate of accumulation of enthalpy, $H$, within the unit is given by Equations (8),(9):

$$\frac{dH}{dt}=\sum _{j=1}^{N_{inlet}}\sum _{p\in P}(F_{j,p}^{in}·h_{j,p}^{in})-\sum _{p\in P}(F_p^{out}·h_p^{out})-R_{p/p_{n,}i}{h}_{p/p_{n,}i}+Q_{trans}$$

(8)

$$H=\sum _{p\in P}(M_{total,p} ·h_p)$$

(9)

Where $h_p$is the specific enthalpy of the material in the unit in phase p, $h_{j,p}^{in}$ is the specific enthalpy of the j^th feed stream in phase p, $h_p^{out}$ is the specific enthalpy of the outlet stream in phase p, $h_{p/p_{n,}i}$ is the enthalpy of phase change between phase p and phase p_n and $Q_{trans}$ is the enthalpy transferred into the unit through boundary due to heat loss, heating and so forth.

Due to the assumption that the unit is perfectly mixed, the specific enthalpy of each of the outlet streams is equal to the specific enthalpy of the material within the unit:

$$h_p^{out}=h_p$$

(10)

The specific enthalpy is assumed to be a function of the composition and temperature, $T$, of the material within the unit:

$$h_p=h_p(T,T_{ref,}{w}_{i,p}),\forall p\in P$$

(11)

3.3. Heat transfer rates

The total heat transfer rate, $Q_{trans}$, is the sum of energy transfer by steam, $Q_{heating}$ and energy loss, $Q_{loss}$, to the environment:

$$Q_{trans}=Q_{heating}+Q_{loss}$$

(12)

$$Q_{heating}=UA\Delta {{\rm T}}_{lm}$$

(13)

$$A=\pi D_{i}NL$$

(14)

$${\Delta T}_{lM}=\frac{{\Delta {\rm T}}_1-{\Delta {\rm T}}_2}{\mathrm{l}\mathrm{n}\left(\frac{{\Delta {\rm T}}_1}{{\Delta {\rm T}}_2}\right)}$$

(15)

Where $U$ is the overall heat transfer coefficient, $A$ is the heat transfer area, $N$ is the number of tubes, $L$ is the length of the tubes and finally $\Delta {\rm T}$ and ${\Delta {\rm T}}_{lm}$ are the linear and logarithmic temperature differences, respectively.

3.4. Flow through distribute plate

The distribute plate is an inserted plate that distributes product stream evenly around the inner periphery of the evaporator tubes. Aside from the liquid holes that allow liquid to pass through, there are also openings to let steam flow through the plate. Two designs for vapor transport through the distribute plate are considered:

• Plate with upstanding vapor tubes.
• Plate with an upstanding rim at the edges of the circular plate and a ring-shaped gap between plate and evaporator body.

The vapor flowing downwards or upwards through the distribute plate causes a vapor pressure drop, while the liquid flowing downwards through the distribute plate and liquid holdup above the plate cause a liquid pressure drop.

The vapor flow through plate, ${\dot{F}}_{vap}$, for plate design with upstanding vapor tubes is expressed by the following relation.

$${\dot{F}}_{vap}=N_{vap tube}·\frac{\pi }{4}{d}_{vap tube}^2{\rho }_{vap}{u}_{vap}$$

(16)

$${\Delta P}_{vap}=\left(1+f_{f,vap}\frac{L_{vap tube}}{d_{ vap tube}}\right)·\frac{1}{2}{\rho }_{vap}{u}_{vap}^2$$

(17)

Where $N_{vap tube}$ is the number of vapour uprising tubes, $d_{vap tube}$ is the diameter of each uprising tube, $L_{vap tube}$ is the length of each tube, ${\rho }_{vap}$ is the vapor density above plate, $f_{f,vap}$is the vapor friction factor and $\Delta P_{vap}$ is the vapor pressure drop.

The vapor flow through plate for plate design with upstanding rim at the edges is expressed by Equations (18)-(19):

$${\dot{F}}_{vap}=\frac{\pi }{4}{(d}_{out}^2{-d_{in}^2)\rho }_{vap}{u}_{vap}$$

(18)

$${\Delta P}_{vap}=\left(1+f_{f,vap}\frac{h_{rim}}{d_{ out}-d_{ in}}\right)·\frac{1}{2}{\rho }_{vap}{u}_{vap}^2$$

(19)

Where $d_{out}, d_{in}$ are the diameter of the evaporator body and the plate respectively and $h_{rim}$ is the height of the outer rim.

The liquid flow through plate, ${\dot{F}}_{liq}$, is expressed by Equations (20)-(21):

$${\dot{F}}_{liq}=N_{liq hole}·\frac{\pi }{4}{d}_{liq hole}^2{\rho }_{liq}{u}_{liq}$$

(20)

$${\Delta P}_{liq}=\left(1+f_{f,liq}\frac{h_{plate}}{d_{liq hole}}\right)·\frac{1}{2}{\rho }_{liq}{u}_{liq}^2$$

(21)

Where $N_{liq hole}$ is the number of liquid holes, $d_{liq hole}$ is the diameter of liquid holes, $h_{plate}$ is the thickness of the distribute plate, ${\rho }_{liq}$ is the liquid density above plate, $v_{liq}$ is the liquid flow velocity through plate, $f_{f,ilq}$is the liquid friction factor and $\Delta P_{liq}$ is the liquid pressure drop.

3.5. Liquid flow through pipes

The wetting rate, $\Gamma$, of a falling film evaporator tube is correlated to the liquid mass flow of the tube, ${\Phi }_m$ and the tube diameter $D_{tube}$ through:

$$\Gamma =\frac{{\Phi }_m}{\pi D_{tube}}$$

(22)

The Reynolds number,$Re$, in the falling film evaporation used in empirical relations is defined by Equation (23):

$$Re=\frac{\Gamma }{\mu }$$

(23)

Where $\mu$is the liquid dynamic viscosity.

The liquid characteristic length, $l_c$, is given as:

$$l_c{=\left(\frac{{\mu }^2}{{\rho }^{2}g}\right)}^{\frac{1}{3}}$$

(24)

Where $\rho$is the liquid density and $g$ is the gravity acceleration.

The film thickness, $\delta$, is calculated by Equation (25):

$$\delta =\left\{\begin{array}{c}{\left(\frac{{2.4\mu }^2}{{\rho }^{2}g}\right)}^{\frac{1}{3}}{Re}^{\frac{1}{3}}=1.34l_c{Re}^{\frac{1}{3}} Re<400\\ 0.302{\left(\frac{{3\mu }^2}{{\rho }^{2}g}\right)}^{\frac{1}{3}}{Re}^{\frac{8}{15}}=0.436l_c{Re}^{\frac{8}{15}} Re \ge 400\end{array}\right\}$$

(25)

Assuming that the film thickness is small compared to the tube diameter, the liquid velocity downwards on the tube, $u$, and the residence time, $t$, for liquid flow to the bottom are calculated based on Equations (26)-(27):

$$u=\frac{\Gamma }{\rho \delta }$$

(26)

$$t=\frac{L_{tube}}{u}$$

(27)

3.6. Boiling point elevation

The boiling point elevation of a solution can be calculated by combining the laws of Raoult and Clausius Clapeyron:

$$\Delta {\rm T}=-\frac{{RT}_{water}^2}{H_{evap}{M}_{w,water}}{\mathrm{l}\mathrm{n}(X}_{water})$$

(28)

Where $R$ is the gas constant, $T_{water}$ is the boiling temperature of water, $H_{evap}$ is heat of evaporation, $M_w$ is molar weight and $X_{water}$ is the molar fraction of water in the solution.

It is assumed that fat does not dissolve in the solvent and therefore does not affect the boiling point. The molar fraction of water in the fat-free product can therefore be calculated:

$$X_{water}=\frac{\frac{f_{m,water}}{M_{w,water}}}{\sum _{i\in {comp}_{fatfree}}\frac{f_{m,i}}{M_{w,i}}}$$

(29)

Where $f_m$ denotes for the mass fraction and $com{p}_{fat free}$ refers to all but fat components.

3.7. Liquid friction factor

The most widely used correlations for liquid friction factor is the Wallis correlation for annular falling film and it is incorporated to the model for better prediction of the liquid friction factor [19].

$$f_{f liq}=0.05·\left(1+300·\frac{\delta }{D_i}\right)$$

(30)

3.8. Energy cost

In order to compare the energy consumption of TVR and MVR, the variable Energy cost is used, defined by Equation (31). When employing TVR technology, the annual cost is given as:

$$Energy cost=Stea{m}_{unit}cost·Stea{m}_{flowrate}·hour{s}_{peryear}$$

(31a)

For MVR technology, the annual energy cost is expressed as:

$$Energy cost={Electricity}_{unit}cost·Powe{r}_{consumption}·hour{s}_{peryear}$$

(31b)

The current price of steam based on natural gas is 20.01$/t [25]. According to Eurostat [26], the electricity unit cost for non-household users in 2023 is 0.21$/kWh. The operating hours per year are considered as 7,920h/year, thus variable Energy cost refers to annual energy cost.

4. Case Studies

Five different evaporators are considered, combined with different energy consumption minimization strategies, namely MVR and TVR. These process alternatives are often met in the food processing industry. A brief description of each studied Case is given below, while the corresponding falling film evaporator layouts is showed in Figure 1. The design characteristics of the evaporators of each layout are provided in Tables S1–S4 of the Supplementary material.

Cases 1 and 5 represent industrial scale evaporator systems while Cases 2,3 and 4 are related to pilot conditions (1,000 kg/h).

4.1. Case 1

Case 1 (Figure 1(a)) examines the operation of a local, South African plant that consists of two evaporation chambers (referred to as effects), a TVR system for vapor recycling from 1^st stage to the same stage and a condenser to remove heat and keep the unit under vacuum [22,27].

4.2. Case 2

Raw milk enters the preheater tubes of the 2^nd and the 1^st evaporator tubes before entering the 1^st evaporator chamber, as shown in Figure 1(b). As the evaporation process continues in FFE1, the generated steam works as the heating medium in the 2^nd evaporator effect where the semi-concentrated product is further evaporated resulting to end-product (concentrate). The steam generated in the 2^nd effect is recycled to TVR where fresh steam is added, increasing its pressure before entering the 1^st effect.

4.3. Case 3

In Case 3, 1,000 kg/h raw milk are fed to the preheater tubes of the 3^rd, 2^nd, and 1^st evaporator effect subsequently, before entering the 1^st evaporator chamber, as shown in Figure 1(c). As the evaporation process continues in FFE1 and FFE2, the generated steam at each effect works as the heating medium for the next one, and the whole process results to end-product (concentrate). The steam generated in the 3^rd effect is recycled to TVR where fresh steam is added, increasing its pressure before entering the 1^st effect.

4.4. Case 4

In Case 4 (Figure 1(d)), 1,000 kg/h raw milk are fed in a single-stage evaporator, with no addition of external steam. The evaporation occurs only by recompressing the vapor generated and recycling it to the unit by an MVR. The only source of energy needed in this process is electricity.

4.5. Case 5

Case 5, presented in Figure 1(e), examines the proposed evaporation system by Galvản-Ángeles et al. [23] and it has four effects with co-current flow, capable of processing 16,000 kg/h of whole milk. The arrangement employs TVR technology, recycling vapor from the 2^nd to the 1^st stage of evaporation.

5. Results

5.1. Solution method

In the present study, global system analysis and optimization were conducted in gPROMS^TM FormulatedProducts [28] version 2.3.0 from Siemens Process Systems Enterprise. A database containing physical and thermodynamic properties of milk components was developed using gPROMS FormulatedProducts Utilities.

5.2. Global system analysis (GSA)

GSA allows for a systematic exploration of the relative impact of input parameters on model output responses. The fundamental principle of GSA involves sampling from a defined particulate system. This entails extracting various values for individual parameters and calculating the corresponding output values.

Uncertainty analysis, can be conducted using the Monte Carlo method, is employed to determine the effects of uncertain input parameters on the system responses. The Monte Carlo method involves performing multiple model evaluations using either deterministic or probabilistic distributions. The resulting output values from these evaluations are then utilized to assess the impact of system's uncertainty. In this study, uncertainty analysis was conducted to assess key decision space optimal areas rather than single optimal values of the degrees of freedom. Uncertainty analysis was performed for all Cases, implementing standard Monte Carlo method. Quasi-random (Sobol) sampling was used as the sample generation method which fills the space uniformly and it is proved efficient for medium-sized problems and a high number of samples. The number of uncertainty scenarios is set to 1000.

5.2.1. Steam cost analysis

This study also considers the uncertainty of steam price. Steam can be generated from a wide range of sources, including fossil fuels such as natural gas (main source), as well as renewable energy sources like lignocellulosic biomass, solar and waste. The rise of renewable energy sources presents an opportunity to examine the prices of steam generated from these sources, taking into account their future values, which could potentially fluctuate by around ±40%. Pérez-Uresti et al. [25] evaluated the production cost of renewable-based steam including solar, lignocellulosic biomass and biogas. Results showed that, considering the current price of steam based on NG, lignocellulosic biomass offers better perspectives to become a competitive resource when compared to the other two sources. More specifically, the evaluated production cost of steam are as follows: NG-based steam costs 16.81 $/tn, biomass-based steam costs 20.37 $/tn, solar based steam costs 32.5$/tn and for biogas-based steam costs 27.42 $/tn.

The analysis conducted in this work takes into account future values of those described above, which could potentially fluctuate by a maximum of 40%, to determine to what extent and in what circumstances renewable-based steam can become a competitive resource compared to natural gas. The specific analysis involves only Case 1, because the same pattern of results would emerge when investigating the rest of the Cases.

Table 4 summarizes the prices examined, to determine to what extent and in what circumstances renewable-based steam can become a competitive resource compared to natural gas.

Uncertainty analysis results related to operating annual cost of evaporation depending on various steam sources and future steam unit costs are summarized in Table 5, which provide a detailed comparison and assessment of steam sources. Despite natural gas being currently the most economical choice, a future possible reduction of biomass-based steam cost by only 20% (or more) would lead to lower values of annual expenditure for the evaporation process of Case 1, compared to the nominal scenario of using NG-based steam. Solar-based steam and biogas-based steam do not seem economically attractive after the analysis as reduction of 40% for the former or over 50% for the latter would be needed in order for those sources to be competitive to natural gas. However, as predictions indicate a rise in natural gas’ price, renewable-based steam would become more and more competitive both for environmental and economic reasons.

5.2.2. GSA-Case 1

The factors used for the uncertainty analysis for Case 1 are the steam cost, the TVR discharge pressure, the TVR suction ratio, $SR$, (given as SR = Steam mass flowrate (suction)/ Steam mass flowrate (motive)), and the feed temperature. The lower and upper bounds along with the probability distribution of the factors are presented in Table S5 of the Supplementary material.

Figure 2(a) illustrates the variation of concentrate mass fraction of water versus feed temperature and TVR discharge pressure. It is clear that higher outlet pressures and higher feed temperatures result to a lower mass fraction of water in the end product, thus, leading to a product with higher dry mass fraction. Higher feed temperatures can be achieved by preheating the feed stream.

As shown in Figure 2(b), the TVR suction ratio illustrates the largest effect on the steam economy. The steam economy is defined as sum of the vapor mass flowrate of each evaporator effect i divided by the feed mass flowrate (Steam economy=Total mass flowrate of evaporated water(kg/h)/Fresh steam mass flowrate (kg/h)). Suction ratios over 1.3 results to higher values of yield, while high feed temperatures along with suction ratios over 1.6 can result to a steam economy close to 4. On the contrary, motive steam mass flowrate is mostly affected by TVR discharge pressure, as shown in Figure 2(c). Low values of discharge pressure result to minimum mass flowrates (around 600 kg/h steam) whereas high outlet pressures along with low suction ratios can result to even 3200 kg/h steam which quintuple cost of evaporation.

From Figure 2(d) and Figure 2(e) it can be seen that TVR discharge pressure affects both the pressure of concentrate product and the solid mass fraction of the FFE1 outlet product. Higher discharge pressure values can be mutually beneficial both for achieving greater vacuum (lower concentrate pressure) and resulting to more concentrated products (lower moisture content) which is ideal for downstream spray drying processes.

5.2.3. GSA-Case 2

Uncertainty analysis was performed for Case 2 to examine the effect of key factors on various key responses (model variables). The factors of this analysis are the TVR discharge pressure, the TVR suction ratio and the feed temperature. information about the factors of this Case is provided in Table S6 of the Supplementary material.

As responses of several variables were discussed in the GSA of Case 1 which also includes two evaporator effects with TVR in plant scale, in this section it is worth highlighting differences in responses related to laboratory scale or responses of other important variables. From Figure 3(a) it can be seen that in laboratory scale, higher values of steam economy can be achieved as the lower and upper bounds of this variable are higher than that of Figure 2(b) (1.4 as opposed to 1 and 5.2 as opposed to 5 respectively).

Figure 3(b) illustrates that although the lower bounds of solid mass fraction after the first evaporation effect are the same as Figure 2(e) (plant scale), the upper bounds are remarkably higher (almost 0.45) indicating that a single evaporation effect would be sufficient in some cases depending on the underlying product specifications.

Figure 3(c) illustrates the annual cost of steam versus feed temperature and TVR outlet pressure. The steam unit cost is fixed at 20$/t (NG – based steam) while the operating hours are 7,920 h. The results show that the operating conditions can greatly affect the annual cost ranging from 6,000 to 44,000 $/year. However, these values are related to different end product specifications, so the final optimization results are presented in Section 4.2.

5.2.4. GSA-Case 3

Uncertainty analysis was also performed for Case 3. The factors used in this analysis are the same as Case 2.

Figure 4(a) shows the variation of steam economy versus feed temperature and suction ratio. It is demonstrated that an arrangement of 3 evaporator effects with TVR leads to high steam economy values compared to Case 2. Specifically, lower and upper bounds of steam economy in Case 3 are 2.5 and 7.5 respectively as opposed to 1.4 and 5.2 of Case 2. Figure 4(b) and Figure 4(c) illustrate the variation of dry mass fraction at FFE1 and FFE2 evaporator tubes outlet respectively, while Figure 4(d) shows the water mass fraction of the concentrate (end product). Higher values of TVR discharge pressure result to more concentrated products while FFE1 liquid product can have 13% to 18% solid content in Case 3 as opposed to 14% to 46% of that of Case 2. In Case 3, the maximum solid concentration of FFE2 outlet product is 26% which implies that the evaporation process is mainly enhanced in the 3^rd effect, rather than the 1^st one in the case of two evaporator effects with TVR (Case 2).

Figure 4(e) illustrates the annual cost of steam versus TVR suction ratio and TVR outlet pressure. The steam unit cost is fixed at 20$/t (NG – based steam) and the operating hours are 7,920. The results show that higher values of outlet pressure combined with low values of suction ratios can greatly increase the annual cost. On the other hand, Case 3 results to lower steam annual cost than that of Case 2. Specifically, the analysis of Case 2 resulted to operating costs of 6,000 – 44,000 $/year (Figure 4(c)) as opposed to 4,000 – 32,000 $/year of Case 3.

5.2.5. GSA-Case 4

The factors consided for the uncertainty analysis for Case 4 are the MVR compression ratio and the feed temperature. The compression ratio is expressed as: CR = Steam pressure (outlet)/ Steam pressure (inlet). The lower and upper bounds along with the probability distribution of the factors is provided in Table S7 of the Supplementary material. The results are shown using various figures, while TVR and MVR operation can be assessed and compared.

Figure 5(a) illustrates the moisture mass fraction of concentrate versus feed temperature and MVR compression ratio. Higher values of compression ratio and feed temperature result to more concentrated products and especially compression ratios over 1.25 lead to an efficient evaporation (solid mass fraction over 0.3). Consequently, the operating cost of evaporation as expressed by Equation (31b), rises with higher values of compression ratios (Figure 5(b)).

It is worth mentioning that the electricity unit cost is fixed at 0.21$/kWh and the operating hours are 7,920. The results show that the operating conditions can greatly affect the annual cost ranging from 2,000 to 36,000 $/year. Comparing with Cases 2 and 3, Case 4 leads to the lower minimum cost (2,000$/year as opposed to 4,000$/year and 6,000$/year of Cases 3 and 2 respectively) and an intermediate upper bound (36,000$/year as opposed to 32,000$/year and 44,000$/year of Cases 3 and 2 respectively). However, these values are related to different end product specifications and a comprehensive comparison is made using dynamic optimization techniques and results are presented in Section 5.3.

Finally, Figure 5(c) depicts the variation of concentrate product pressure versus MVR compression ratio and feed temperature. The latter one illustrates the greatest impact on the achieved vacuum, as lower values of feed temperature result to greater vacuum (around 2500 Pa pressure). However, as shown in Figure 5(c), end-product dry mass fraction over 0.5 cannot be achieved for feed temperatures lower than 40°C.

5.2.6. GSA-Case 5

Uncertainty analysis was finally performed for Case 5 to examine the effect of the split fraction, TVR suction ratio and feed temperature (factors) on various KPIs. More information about the studied factors is provided in Table S8 of the supplementary material. As Case 5 involves a plant scale process with 4 evaporator effects by recycling steam from the 2^nd effect to TVR, different recycling fractions (split fractions) can be examined. Furthermore, both plant and laboratory processes can be evaluated.

Figure 6(a) illustrates the water mass fraction of concentrate product versus the feed temperature and split fraction. Empty spaces indicates that simulations cannot be provided in certain areas, thus, recycling over 60% of FFE2’s steam (providing FFE3 with the remainder steam) while keeping the feed temperature below 50°C results to failed simulations. However, higher feed temperature and lower split fractions enhance the evaporation process as the concentrated product has lower moisture content. This implies that by providing more steam to the next evaporator effect, the evaporation is more enhanced than recycling larger quantities of steam in the TVR unit.

Figure 6(b) shows the variation of yield versus the feed temperature and split fraction. The evaporation yield is defined as the mass flowrate of end product (concentrate) divided by the feed mass flowrate (Yield = Concentrate mass flowrate (kg/h)/Fresh steam mass flowrate (kg/h)). It follows the exact opposite pattern with Figure 6(a), as higher feed temperature and lower split fraction result to lower product yield. However, despite producing smaller amounts of concentrated product, the lower yield indicates products with reduced moisture content.

Figure 6(c) illustrates the annual steam cost (steady state values) versus TVR suction ratio and split fraction. Suction ratio influences steam cost the most, as higher values of suction ratio result to higher costs, ranging from 370,000 – 500,000 $/year. It is demonstrated that there a relatively linear relationship between laboratory and plant scale evaporation processes.

Table 6 summarizes minimum and maximum values of outlet dry mass fraction of each effect resulted from the uncertainty analysis. It is concluded that for each of the first 3 evaporator effects there is a narrow range of dry mass fraction. FFE1 liquid product outlet dry mass fraction is limited to 0.17 while the FFE2 liquid product outlet dry mass fraction ranges from 0.225 – 0.254. It is thus indicated that the evaporation process is mainly enhanced in the 3^rd and 4^th evaporator effects with maximum outlet product dry mass fraction of 0.36 and 0.61 respectively.

5.3. Dynamic optimization

Depping et al. [29] investigated whether switching from milk powders to new products known as milk concentrates diminishes the overall environmental impact along the supply chains of dairy-containing products. For relevant environmental indicators such as cumulative energy demand, global warming potential, eutrophication potential, and acidification potential, concentrates were found to have a lower environmental impact than powders, even if the former are trucked up to 1000 km. To this end, as drying milk into milk powders is a highly energy-intensive process and energy demand per kilogram of water removed rises nonlinearly with increasing dry-matter content [30], it is crucial to investigate the amount of energy saved when switching to milk concentrates. This can be done by optimizing the annual operating cost of steam or electricity used in the evaporation process with respect to different product specifications. The new products considered are shelf-stable milk concentrates, which are compared to the current benchmark product, milk powders. Specifically, 3 concentrate dry mass contents are considered: 30%, 35% and 50%. The first two are related to concentrates and the third one is the one necessary for milk powder production before drying.

The optimization of the process aims to maximize or minimize an objective function subject to design and operating constraints related to the nature of the process, the end-product specifications or the environment. In this work, advanced dynamic optimization technics are employed for Cases 1 – 4, aiming at increasing process performance while minimizing operating cost (in case of existing industries) or total cost (in case of new plant design) under different end-product specifications. Specifically, Cases 1–4 are optimized under 6 scenarios, described in Table 7 along with product quality specifications. Optimization results are compared with the nominal Cases expressed by simulation results to highlight the benefits of optimization.

In many cases, the optimization is applied to already existing industries the design of which is fixed. However, it is worth examining an objective function that includes annualized fixed cost along with the operating cost of steam, for the cases of new plant design, capacity expansion or even the examination of capital investment due to trade-offs between capital and operating costs. Equation (32) represents the total annualized cost [31].

$$Total cost=CCF·F_{BM}·\left(\sum _{i=1}^n{10}^{3.9119+0.8627\mathit{log}\left(A_i\right)-0.0088{\left[log{(A}_i\right)]}^2}\right)+stea{m}_{unit}cost·Stea{m}_{flowrate}·hour{s}_{peryear}$$

(32)

Where $CCF(=1/3)$ is the fixed capital charge coefficient, $FBM (=2.45)$ is the installed equipment coefficient and n the number of evaporator effects. Operating hours are fixed at 7,920 h.

The optimization variables represent the degrees of freedom of the optimization problem. For Cases 1-3 and optimization scenarios 1-4,6, the optimization variables are: (a) the TVR – Pressure (discharge), (b) TVR – Suction ratio, (c) Feed Temperature while for optimization scenario 5 the optimization variables are: (a) FFE1 tube inner diameter, (b) FFE1 tube length, (c) FFE2 tube inner diameter and (d) FFE2 tube length. The optimization variables for Case 4 and scenarios 1-4,6 are: (a) MVR-Adiabatic efficiency, (b) MVR-Compression ratio and (c) Feed temperature for scenarios. Finally, for Case 4 and scenario 5 the optimization variables are: (a) FFE1 tube inner diameter and (b) FFE1 tube length. More comprehensive information regarding the optimization variables, including their upper and lower limits, for every examined Case and scenario, is provided in Tables S9–S11 of the Supplementary materials.

The time horizon of interest is fixed at 1h in all Cases and scenarios, while the optimization variables are time-invariant, and the constraints (equality or inequality) must be satisfied at the end of the time horizon. In each optimization case, a full description of the constraints and optimization results are provided for each Scenario.

5.3.1. Dynamic optimization of Case 1

The constraints imposed on the process optimization and state variables are summarized in Table 8, along with their explanation.

The results of the various optimization Scenarios for Case 1 are summarized in Table 9. In Case 1, the minimum steam annual cost when producing product with 50% moisture content (necessary for the subsequent spray drying) is 232.118 $/year, while when the end-product specification is reduced to 35% dry mass content, the annual steam cost decreases by 14%, to 199,609 $/year (Scenario 2). Further reduction in concentrate dry mass fraction constraint to 30% leads to approximately 22% reduction in annual steam cost (181,563 $/year, Scenario 3). In all Cases, the optimization strategies were able to reduce the annual steam cost to over 37%. Product yield is also affected by the product specifications, as it increases with the rise of product water content. Consequently, higher production rates can be achieved when producing concentrates as opposed to milk powders, with the same feed flowrate. Feed temperature seems to be optimal at its upper bound in the first 3 scenarios, indicating that preheated feed streams should be considered.

Comparing scenarios 4 and 6, high values of product yield cannot be achieved when restricting the product dry mass fraction to 50% even when maximizing yield. However, when this constraint is relaxed (Scenario 6), the yield can take values up to 0.42 as opposed to 0.25 of Scenario 4.

Regarding the annualized total cost, it is worth considering the optimization Scenario 5 as opposed to the base Case 1. For a new plant design, optimization strategies can result to around 45% reduction in the total cost. The plant under consideration has 16 m tubes with 0.0508 m inner diameter in both evaporator effects. To minimize the total annualized cost, regarding an expansion of this plant or a new plant design, the optimal values of these dimensions would be 17 m and 5 m tube length and 0.0435m and 0.0254 m inner diameter for the 1^st and the 2^nd evaporator effect respectively. The total cost is reduced to 2.33·10⁶ $/year for a plant processing 10,600 kg/h raw milk (Case 1). From the aforementioned results it is clear that optimal strategies would enhance evaporation in the 1^st effect (longer tubes and tube diameter) rather than sticking to the identical design of the two evaporator chambers with moderate dimensions.

5.3.2. Dynamic optimization of Case 2

The constraints imposed on the process optimization and state variables are the same with Case 1, summarized in Table 8. The optimization results are presented in Table 10.

As Case 2 involves a laboratory/pilot scale evaporation process with TVR, various results can be concluded from Table 10 compared to Case 1 (plant scale). Regarding Case 2, the minimum annual steam cost when producing a product with 35% dry mass content (Scenario 2) is 18,887 $/year, indicating a 10.75% reduction from the minimum cost when producing milk powder restricted products (Scenario 1). Further reduction to the restriction in the concentrate dry mass fraction to 30% lead to approximately 18.8% saving in annual steam cost (17,184 $/year, Scenario 3). In all Cases, optimization strategies were able to reduce the annual steam cost to over 31%. Results regarding product yield (Scenario 4 and 6) follow the same trend as Case 1.

The total annualized cost is optimized in Scenario 5 as opposed to the base Case 2. The results indicate a 46% reduction in the total cost. The pilot plant under consideration has 6 m tubes with 0.0508 m inner diameter in both evaporator effects. To minimize the total annualized cost, a new laboratory scale process would select the optimal values of these dimensions to be 5.7 m and 0.6 m tube length and 0.0508m and 0.0254 m inner diameter for the 1^st and the 2^nd evaporator effect respectively. Thus, the total annualized cost is reduced to 274,768 $/year for processing 1,000 kg/h raw milk (Cases 2, 3 and 4).

5.3.3. Dynamic optimization of Case 3

The constraints imposed on the process optimization and state variables are summarized in Table 11, along with their explanation. The results of the various optimization Scenarios are summarized in Table 12.

$W_W\left(t_f\right)=0.5 (\mathrm{o}\mathrm{r} 0.65\mathrm{o}\mathrm{r} 0.7)$$283.15 K\le T_{FFE1}\left(t_f\right)\le 343.15 K$$283.15 K\le T_{FFE2}\left(t_f\right)\le 341.15 K$$283.15 K\le T_{FFE3}\left(t_f\right)\le 339.15 K$$0.5\le W_W\left(t_f\right)\le 0.7$

The results demonstrate that the minimum annual steam cost when producing product with 50% dry mass content (Scenario 1) is 15,078 $/year, while a 13.85% (or more) reduction can be achieved when producing milk concentrates (Scenarios 2 and 3). Optimization strategies were able to reduce annual steam cost to over 32% compared to Base Case.

The total annualized cost is optimized in Scenario 5 as opposed to the base Case 3. The results indicate an approximately 35% reduction in the total cost. The pilot plant under consideration has 6 m tubes with 0.0508 m inner diameter in all three evaporator effects. To minimize the total annualized cost, a new laboratory scale process would select the optimal values of these dimensions to be 4.25 m, 4.12 m and 4.13 m tube length and 0.0508 m, 0.0381 m and 0.0381 m inner diameter for effects 1, 2 and 3 respectively. The total annualized cost is therefore reduced to 449,236$/year for processing 1,000 kg/h raw milk (Cases 2, 3 and 4).

5.4.4. Dynamic optimization of Case 4

The constraints imposed on the process optimization and state variables are summarized in Table 13 and the attained results for each Case scenario are presented in Table 14.

Case 4 is compared to Cases 2 and 3 to investigate whether MVR or TVR is the most economically attractive option, and under which circumstances the opposite figure could be profitable. However, various results can be concluded from Table 14, solely for Case 4. The minimum annual steam cost when producing product with 50% moisture content is 34,319 $/year, while when the end-product specification is reduced to 35% dry mass content, the annual steam cost decreases by 30.7%, to 23,815 $/year (Scenario 2). Further reduction in the restriction of the concentrate dry mass fraction to 30% leads to approximately 44% reduction in annual steam cost representing 19,299 $/year (Scenario 3).

Regarding the total annualized cost, the optimization strategies resulted to a 9% reduction in the total cost (significantly less than in Cases 2 and 3). However, its value is the minimum among the 3 Cases due to the existence of a single evaporator effect. The evaporator examined has 6 m tubes with 0.0508 m inner diameter. To minimize the total annualized cost, the optimal values of these dimensions are 5.7 m tube length while the optimal inner diameter is the same. The total cost is therefore reduced to 251,240 $/year for a pilot plant processing 1,000 kg/h raw milk (Case 4).

6. Conclusion

In this work, five different and industrially relevant milk evaporator Cases have been studied using a model-based approach. For each Case study various conclusions were drawn regarding the simulation, model validation, global system analysis and process optimization. In this section, the results are discussed in terms of comparing TVR and MVR.

Comparing Cases 2, 3 and 4 for current steam and electricity prices, when processing 1,000 kg/h raw milk the most economical option includes 3 evaporator effects with TVR (Case 3) to meet the desired 50% product dry mass content. The same figure is reported for optimization Scenarios 2 and 3. However, Case 4 indicates the most significant reduction in the annual cost when reducing the product specification to 30 or 35% dry mass content.

It is worth mentioning that current high electricity prices (0.21$/kWh) lead to Case 4 being the most unprofitable choice. When producing product with 35% dry mass content, only a 11% reduction in the unit electricity price leads to Case 4 being more cost effective than Case 2 with only a single evaporator effect. Simultaneous reduction by 7% of electricity price along with a 5% increase of gas-based steam price would also lead to Case 4 being the most profitable option among these Cases.

Regarding the maximum values of product yield in Cases 2-4 (Scenario 4), Case 2 can achieve slightly higher values than Cases 3 and 4 (0.28 as opposed to 0.25 respectively). Moreover, for a new plant design, the minimum total annualized cost is achieved in Case 4 which includes a single evaporator effect with MVR, thus indicating that the capital fixed cost in such processes has the dominant contribution on the total annualized cost compared to the operating costs.

Cases 1 and 2 can be compared as they both include 2 evaporator effects with TVR, the former relevant to a plant scale and the latter to a pilot scale. The annual steam cost seems to have a relatively linear relationship with capacity while lower product yield values can be achieved in Case 1 when producing products with 50% dry mass content.

Overall, switching from milk powder production to milk concentrates results to a reduction in the annual cost from 10.75 to 44% depending on the Case under consideration. Furthermore, a forecasted reduction of biomass-based steam cost by only 20% (or more) leads to lower values of the annual expenditure in all Cases, than that of the currently used NG-based steam. As predictions indicate a rise in natural gas’ price, renewable-based steam would potentially become more and more competitive. Finally, assuming a simultaneous increase in the price of NG-based steam by 10% and a reduction of biomass-based steam by 10%, the former one is no longer the most economically attractive solution.