Mathematical modeling of LFR was conducted to compute the useful energy harnessed from solar radiation by the thermal fluid at the factory location. The modeling framework comprised two main components: optical and thermal models. The optical model evaluated the optical efficiency of the collector, which is influenced by the sun’s position, geometry of the collector, reflectivity of its components and incidence angle modifier. The critical solar angles that influence the redirection of reflected solar radiation towards the receiver were determined. Figure 1 illustrates these critical angles, including the longitudinal angle of incidence, Φ_L, and transve.rsal incidence angle, Φ_T, along with their determinants. A north-south (N-S) mirror orientation was chosen for this study since it optimizes annual irradiance utilization, particularly at solar noon [26], [27].

Figure 1.

LFR system angles (a) Mirror inclination angle and (b) Transverse solar altitude angle and (c) Longitudinal and transversal incidence angles [27], [28]

The mirror inclination angle, ψ_i was calculated using eq. (1) [27]:

To determine the angle between reflecting surface of the i^th mirror and horizontal ground plane, β_i , the dependence on the receiver height above the primary mirrors, h_fr, and the transversal location of the i^th mirror, d_i, were considered, as shown in eq. (2):

Correlations for the other critical parameters were adopted from reference models as detailed in Table 1.

The other two key parameters for computing optical efficiency: the incident power on the collector’s surface directly from the sun (P_Ap) and the incident power on the absorber tube, after reflection by primary mirrors (P_abs) are given by eq. (3) and eq. (4) [34]:

where DNI represents location’s direct normal irradiance (W/m²) and A_p denotes area of collector’s aperture (m²) and θ_z is the Zenith angle.

where N_ph represents the number of photons reaching the absorber’s surface, and P_ph denotes individual photon’s power (W) at a specific location. P_ph was determined using Tonatiuh software, a Monte Carlo Ray Tracing (MCRT) tool [35].

The geometry-dependent determinants of the number of photons reaching the absorber’s surface, N_ph, namely, intercept, shading loss and blocking loss factors were then determined. Appendix A1 provides correlations for computing these determinants. To account for additional optical losses from mirror surface imperfections, the nominal mirror reflectivity of the reference collector was reduced from 0.95 to 0.88 in the Tonatiuh simulations, consistent with reported reductions due to soiling [36].

The collector’s optical efficiency was then calculated using eq. (5):

The thermal model, on the other hand, was developed to calculate the energy absorbed by the heat transfer fluid and thermal efficiency while accounting for heat losses from the receiver to the atmosphere. Here, a one-dimensional heat transfer model of an evacuated tube receiver with secondary reflector was developed. This approach was chosen because one-dimensional models, when developed with sufficient physical detail have demonstrated a good agreement with experimental data [36]. Appendix A2 illustrates the heat flow in the receiver and outlines the correlations for heat transfer.

A thermal energy balance algorithm was developed in MATLAB to evaluate the heat gain and losses per unit length along the absorber tube. The model couples local energy balance with temperature-dependent thermophysical properties of water and steam obtained directly from the XSteam library, allowing continuous tracking of the fluid’s thermodynamic state from sub-cooled liquid at the inlet temperature, T_i, of 87 °C to saturated conditions at the outlet temperature, T_o, of 188.15 °C and a pressure of 12 bar.

The algorithm progresses axially along the receiver in uniform segments of 0.5 m. At each segment, N, the inlet and outlet heat transfer fluid temperatures, (T_HTF,i)_N and (T_HTF,o)_N are determined, computing absorbed solar energy, convective heat transfer to the working fluid and the external convective and radiative losses. The sections are shown in Figure 2.

Figure 2.

Axial portioning of absorber tube

The convective heat transfer coefficient was determined according to the local flow regime while the outer losses were evaluated by solving coupled energy balance equations for the absorber, glass envelope and secondary reflector using an iterative solver (fsolve). The wall temperature was updated iteratively until convergence, after which the local enthalpy was incremented based on the computed heat gain per unit length. The useful heat gain by the heat transfer fluid per unit length is equivalent to the convective heat transfer between the absorber tube’s inner surface and HTF.

The thermal efficiency was computed using eq. (6) [10]:

where, Qu_total denotes the total useful heat gain (W).

The model was validated using data reported in the literature, with key parameters, namely, mirror inclination, IAM and heat gain, benchmarked against the experimental results of Pino et al. [34] and Hafner et al. [37] for real LFR configurations. The validation process involved incorporating the geometric parameters and operating conditions of the experimental plant into the model and comparing the simulated outcomes with the experimental data. The discrepancies between the model results, experimental data and real LFR configuration were evaluated using the Coefficient of Determination (R²) statistical method, computed using eq. (7) [38]:

where A_i denotes measured values, P_i denotes the predicted model values and Ē_i denotes the mean of measured values.

Individual optical loss components (blocking, shading and spillage) were not independently validated due to the absence of published experimental datasets reporting detailed geometric loss decomposition for identical LFR configurations. However, the adopted loss formulations are consistent with established analytical LFR optical models reported in literature.

To aid the development of the model, several assumptions were made. The thermal model was based on steady-state conditions, ignoring potential transient effects. Additionally, a uniform circumferential heat flux distribution was assumed on the absorber’s surface.

Although one-dimensional steady state formulation was employed, previous studies have demonstrated that such approaches can reliably predict receiver thermal performance, producing results consistent with experimental data. Liang et al. [36] compared one-dimensional models assuming uniform circumferential heat flux with more detailed formulations accounting for non-uniform flux and demonstrated that, while circumferential non-uniformities may induce localized wall temperature gradients, the differences in predicted outlet fluid temperature and overall thermal efficiency remain small. In the present configuration, the secondary reflector further enhances circumferential redistribution of reflected radiation, thereby reducing peak flux concentrations and temperature asymmetry on the absorber surface. Thus, while transient behaviour and local flux non-uniformities may occur under real operating conditions, their influence on overall thermal input is expected to be limited. Future work may incorporate transient and multi-dimensional modeling to enable refined prediction of local temperature fields and detailed thermal stresses.

A MATLAB simulation modelled the hourly energy flows of the LFR system integrated with storage and fuelwood backup, evaluating performance across six months: high irradiation months (January, February, December) and low irradiation months (April, May, October). The simulation accounted for continuous 24-hour factory operation, with monthly operational days adjusted accordingly. For each month, the hourly thermal energy demand was kept constant, with its value varying across months depending on respective monthly demand patterns. Hourly direct normal irradiance values were sourced from National Solar Radiation Database (NSRD) [39], while hourly optical and thermal efficiencies were calculated using the model developed in this study. The collector area was held constant throughout the simulations. Also, the thermal storage capacity was not treated as an independent design variable but was coupled to the optimized solar field size within the system optimization framework. It, therefore, remained fixed during the analysis.

The algorithm operates through nested temporal loops, progressing from months to days to hours, to replicate energy dispatch priorities while accounting for system constraints. At the start of each monthly cycle, the storage state of charge (SOC) is reset to zero and hourly, location-specific DNI and efficiency data are processed. Solar generation is prioritized to meet process demand, with surplus energy directed to thermal storage subject to capacity constraints, and any excess dissipated. During deficits, the system first deploys stored energy, accounting for discharging losses; followed by dispatch of a fuelwood-fired boiler, treated as a load-following backup. This hierarchy maximizes renewable utilization while ensuring operational reliability under varying irradiation conditions. A flowchart of the simulation process is presented in Figure 3.

Figure 3.

Flowchart of hourly energy flow simulation

The hourly dispatch algorithm assumes idealized priority-based control without explicitly modeling actuator response times, valve dynamics, or ramp-rate constraints of the storage and back-up systems. Energy flows are resolved at an hourly time step, and switching between solar, storage and fuelwood sources treated as instantaneous within each time interval. This approach captures system-level energy balances while neglecting short duration transient control effects occurring at sub-hourly time scales. Incorporating detailed dynamic control modeling would increase computational complexity and require additional system-specific control parameters. Future work may integrate transient control simulation to quantify short term operational effects and refine dispatch predictions.

The simulation outputs include monthly solar energy generated, energy supplied from storage, and energy supplied from fuelwood boiler system. In addition, the solar fraction – a key performance metric – is calculated monthly to capture the system’s responsiveness to seasonal irradiation patterns. The hourly solar energy generated was computed using eq. (8):

where DNI denotes the hourly direct normal irradiance (MW/m²), h denotes the number of hours (set to 1), A_c denotes the total collector area in (m²), η_opt and η_th represent collector’s hourly optical and thermal efficiencies, respectively.

The charging and discharging efficiencies were set at 0.99 [40] while the distribution system efficiency was set at 0.85.

The reduction in CO₂ emissions resulting from partial displacement of fuelwood was quantified using an emission factor of 2.27 kg CO₂ per kg MT reported for Kenyan tea factories [16], [17]. Baseline factory data indicate that an annual production of 3,450,000 kg MT requires approximately 11,500 trees [18]. These values were used as reference case. The annual solar fraction was determined as the average of the representative high and low irradiation months obtained from the dynamic hourly simulations. The corresponding reduction in annual fuelwood contribution was then applied proportionally to the baseline tree consumption and associated emissions.

Table 4 presents the optimized geometric parameters of the solar field – including receiver height, mirror width and row spacing – along with the corresponding optical efficiency.

Figure 7 presents the two-dimensional contour plots that provide a sensitivity-based assessment of the optimized configuration, showing how optical efficiency varies with small deviations from the optimal design parameters. The contour plot of receiver height and mirror width in Figure 7a exhibits a single peak centered near a receiver height of 4.4 m and mirror width of 0.5 m, respectively. The contours form approximately elliptical patterns around the optimum, indicating a unique global solution. Below 4.4 m, the optical efficiency increases rapidly with receiver height due to reductions in shadowing, blocking of reflected radiation and spillage losses. Beyond 4.4 m, the contours become more widely spaced in the height direction, indicating reduced gradient magnitude and greater tolerance to further height increases. In contrast, tighter contour spacing for mirror widths exceeding 0.5 m reflects stronger sensitivity to aperture enlargement, where increased shading of adjacent reflectors and blocking of reflected radiation outweigh gains in reflective aperture area.

Similarly, the mirror width – row spacing contour in Figure 7b demonstrates a pronounced peak at approximately 0.5 m and 0.2 m. Optical efficiency increases rapidly as row spacing approaches 0.2 m from the lower values due to reduced inter-row shading but declines sharply beyond this point as reflected radiation increasingly misses the receiver. The comparatively dense contour spacing along the row spacing axis indicates higher sensitivity to spacing variations relative to mirror width.

The receiver height-row spacing contour in Figure 7c further confirms broader tolerance in receiver height compared to row spacing, particularly above the 4.4 m optimum. Optical efficiency rises steadily with height up to 4.4 m and decreases more gradually thereafter, whereas row spacing exhibits steeper gradients on both sides of 0.2 m.

Figure 7.

Optical efficiency contour plots: (a) Receiver Height vs Mirror Width, (b) Mirror Width vs Row Spacing, and (c) Receiver Height vs Row Spacing

Overall, the contour topology provides insight into the response surface of the objective function. Receiver height primarily controls interception efficiency and ground clearance effects, while mirror width and row spacing strongly influence inter-mirror shading and reflected beam divergence. The smooth contour regions surrounding the optimum indicate a stable global solution within the defined design space. A comparison with previous studies, summarized in Table 5, reveals both consistencies and variation in the optimization trends.

There are notable similarities between the results of this study and those reported in previous works with respect to mirror width and row spacing, which generally fall within the ranges of 0.4 – 0.5 m and 0.06 – 0.25 m, respectively. The consistencies reinforce the validity of the presented model. However, a significant discrepancy is observed in the optimized receiver height. While reported optimal heights reach up to 11.81 m, the present study yields a lower optimal height of 4.4 m for a medium-sized collector. This variation reflects both the collector scale and the adopted optimization methodology. In this study, the receiver height is optimized simultaneously with mirror width and row spacing, capturing the interplay among these parameters, which significantly influences optimal system performance.

In general, the identification of optimal parametric setpoints ensures a balanced system design that maximizes the optical efficiency by minimizing shading, blocking, and spillage.

The MATLAB-based simulation results for the evacuated tube receiver are presented in Figure 8. The inlet HTF and ambient temperatures were maintained at 87 °C and 18 °C, respectively. As the HTF outlet temperature increases from 100 °C to 300 °C, the thermal efficiency decreases gradually from 65.05% to 62.6%. However, beyond 300 °C, thermal efficiency drops steeply, reaching 58.51% at 400 °C. In contrast, the heat loss per unit length rises significantly with increasing outlet temperature, from 22.76 W/m at 100 °C to 219.76 W/m at 400 °C. This inverse relationship indicates that as the outlet temperature increases, convective and radiative heat losses become more pronounced, leading to a marked reduction in thermal efficiency. The sharp decline beyond 300 °C underscores the dominant influence of the temperature-dependent radiation losses at elevated fluid outlet temperatures.

Figure 8.

Influence of heat transfer fluid outlet temperature on LFR receiver’s heat loss per unit length and thermal efficiency

An average daily thermal energy demand of 56 MWh serves as the benchmark to evaluate the capability of each system size in meeting the energy requirements. Figure 9 illustrates the influence of solar multiple on solar energy generated and solar fraction. The solar field area ranges between 7680 m² and 23,040 m² as SM increases from 1 to 3. From Figure 9a, at a SM of 1, the system only meets the 26 MWh demand during sunny hours. For solar multiples between 1.25 and 2, the system meets the demand during sunny hours while storing 7.39 – 27.47 MWh of the 30 MWh nighttime demand. At SM 2.25, the system generates a total of 60.86 MWh, marginally exceeding the daily demand. The stored 34.86 MWh fully meets nighttime demand and offers a 4.86 MWh buffer for variability. Between SM 2.5 and 3, the system produces a surplus for storage, increasing from 41.2 MWh to 54.93 MWh, with equivalent storage hours varying between 17 and 22.

Secondly, achieving 100% solar fraction is crucial for maximizing the substitution of fuelwood. Figure 9b shows that at SM values of 2 or lower, the solar field area is insufficient to meet this target. SM 2.25 marginally exceeds this target. When SM exceeds 2.25, SF significantly exceeds 100%, indicating a system overdesign for fuelwood substitution. Therefore, optimal system features a SM of 2.25, 1.03 SF and an area of 17,280 m², striking a balance between maximizing fuelwood substitution and minimizing the solar field size.

Figure 9.

Influence of solar multiple on (a) thermal energy generated and stored thermal energy and (b) solar field size and solar fraction

The technical specifications of the LFR plant are detailed in Table 6.

Table 8 shows the sensitivity of economic indicators to the variation of performance and cost inputs for the optimal configuration (SM = 2.25). Results indicate that LCOH varies proportionally with capital expenditure (±5%). Variations in DNI and thermal efficiency produce an inversely proportional change in LCOH due to its direct influence on annual thermal energy yield. NPV is most sensitive to capital expenditure, exhibiting approximately ∓25% variation for ±5% change in capital expenditure, reflecting the influence of upfront investment on long-term cash flow. O&M cost variations exert moderate influence on NPV, while DNI and thermal efficiency variations produce smaller NPV shifts, with ±4% and ±6% variations occasioning ±5% and ±7% changes, respectively. IRR shows comparatively limited sensitivity, all within ±2%, remaining close to but above the prevailing interest rate under the tested uncertainty bounds. Overall, the system retains economic feasibility under realistic variations in solar resource, performance and cost parameters, although capital expenditure remains the dominant economic risk factor.

Table 9 shows the sensitivity of economic indicators to the variation of key financial parameters for the optimal configuration (SM = 2.25). The results indicate that NPV is particularly sensitive to variations in interest rate. Since the optimized configuration yields a marginal baseline NPV (1,288 USD), ±4.25% changes in the interest rate produce comparatively large relative shifts (∓30%) in NPV over the 25-year analysis period. Higher interest rates substantially reduce discounted cashflows, thereby compressing profitability margins, while lower financing costs improve economic attractiveness. The variations also occasion a ∓6.9% change in LCOH and ∓1.8% change in IRR, indicating that financing conditions materially influence both cost competitiveness and investment performance. At an elevated interest rate of 17%, the IRR decreases to 12.55%, falling below the baseline financing threshold and rendering the project economically unattractive under such borrowing conditions. VAT variations primarily affect capital expenditure and therefore induce modest shifts in LCOH and NPV. ±2% VAT variation occasions ∓1.97% and ∓4.96% changes in LCOH and NPV, respectively. However, given the marginal profitability of the optimized configuration, even small capital cost adjustments translate to measurable changes in IRR, reinforcing the sensitivity of this configuration to fiscal policy parameters. Inflation predominantly influences long-term operating cost escalation and revenue progression over the 25-year analysis period. Within the tested ±4.51% range, inflation variations result in approximately ∓4% change in LCOH and ∓10% in NPV(1,159 – 1,417 USD). The associated IRR varies by approximately ∓0.35%, yielding a range of 12.43% to 13.13%. While these variations are less pronounced than those occasioned by interest rates fluctuations, higher inflations slightly compress long-term profitability by reducing cash flow margins.

Overall, the analysis demonstrates that while optimized LFR configuration remains economically viable under moderate macroeconomic fluctuations, its marginal NPV makes it particularly sensitive to financing conditions. These findings underscore the importance of stable borrowing environments and supportive financial frameworks in enhancing the attractiveness of industrial solar thermal systems.