*Article* **Optimal Configuration with Capacity Analysis of a Hybrid Renewable Energy and Storage System for an Island Application**

#### **Chih-Ta Tsai, Teketay Mulu Beza, Wei-Bin Wu and Cheng-Chien Kuo \***

Department of Electrical Engineering, National Taiwan University of Science and Technology 43, Sec. 4, Keelung Rd., Taipei 10607, Taiwan; marcotsai3@gmail.com (C.-T.T.); mteke24@gmail.com (T.M.B.); andytouchwork@gmail.com (W.-B.W.)

**\*** Correspondence: cckuo@mail.ntust.edu.tw; Tel.: +886-920881490; Fax: +886-227376688

Received: 26 October 2019; Accepted: 16 December 2019; Published: 18 December 2019

**Abstract:** The Philippines consists of 7100 islands, many of which still use fossil fuel diesel generators as the main source of electricity. This supply can be complemented by the use of renewable energy sources. This study uses a Philippine offshore island to optimize the capacity configuration of a hybrid energy system (HES). A thorough investigation was performed to understand the operating status of existing diesel generator sets, load power consumption, and collect the statistics of meteorological data and economic data. Using the Hybrid Optimization Models for Energy Resources (HOMER) software we simulate and analyze the techno-economics of different power supply systems containing stand-alone diesel systems, photovoltaic (PV)-diesel HES, wind-diesel HES, PV-wind-diesel HES, PV-diesel-storage HES, wind-diesel-storage HES, PV-wind-diesel-storage HES. In addition to the lowest cost of energy (COE), capital cost, fuel saving and occupied area, the study also uses entropy weight and the Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) method to evaluate the optimal capacity configuration. The proposed method can also be applied to design hybrid renewable energy systems for other off-grid areas.

**Keywords:** renewable energy; hybrid energy systems; cost of energy; energy storage; distributed generation (DG); sensitivity analysis

#### **1. Introduction**

To achieve good economic life and growth, off-grid communities require an affordable and reliable energy supply. Besides, the growing concern about climate change and environmental pollution, especially since fossil fuels are the main source of energy on Earth, has pushed the power generation systems towards the use of renewable energy [1,2]. In addition to high transportation and fuel cost, energy delivered to isolated areas frequently used fossil fuel-based generators which threaten the anthropogenic and natural ecosystems [3–5]. The main factors pushing increased energy access are regular power interruptions, limited power grid accessibility, availability of renewable resources, increasing concern about the need to decrease fossil fuel dependency and high oil prices, and also the significance of escalating energy access to development and availability of climate- related financial support for low-carbon and climate-resilient development that prioritizes renewables [6–8].

Nowadays, the main challenges facing hybrid systems is to design an energy management strategy to meet the demand for loads, despite the intermittent nature of renewables [9,10], in addition to cost savings and total versatility and multi-faceted goals that can be accomplished [11–13].

PV and wind turbine (WT) have been considered the most promising renewable energy options for off-grid areas or islands to fulfill the energy demand [14–16]. Even stand-alone wind and solar energy may fulfill the low load requirements, while these systems need a significant energy storage for higher loads, resulting in high COE [17–19]. The other option to alleviate this problem is autonomous hybrid renewable energy systems (HRES) which combine two or more energy resources, to fulfill higher energy requirements of off-grid areas and resolve the inherent problem of single renewable energy (RE) resource [1,20–24]. Furthermore, hybridization of energy sources increases the reliability of the system as the shortcomings of any component are compensated for by the selection of other but appropriate components and their sizing is essential during design of such systems [25–27]. Some of the hybrid energy systems with different storage technologies and performance measure criteria found in literatures are presented in Table 1 [12,28,29].


**Table 1.** Hybrid power systems with various storage technologies.

Key: **CC**: Cyclic Charging; **CD**: Combined Dispatch; **DF**: Duty Factor; **EE**: Excess Energy; **COE**: Cost of Energy; **HOGA**: Hybrid Optimization by Genetic Algorithms; **LA**: Lead Acid; **Li-ion**: Lithium-ion; **RF**: Renewable Friction.

In most of the literature on hybridizing energy systems, cost reduction and CO2 emission reduction are the focal points of the research work. However, the number of scenarios considered to pinpoint the best trade-off configuration among system reliability, cost of energy and environmental sustainability; the magnitude of excess electricity and the way it is managed; the percentage RF, still needs further investigation, depending on where the HES is going to be installed.

The Batanes Islands are the northernmost provinces of the Philippines. They are composed of 10 islands including Batan, Sabtang, Itbayat and so on. They have a typical volcanic island hilly terrain and a rich natural ecology. The largest island is Batan Island with a longitude of 120.968◦ and a latitude of 20.445◦, as shown in Figure 1, about 161 km from Luzon, Philippines, and only about 190 km from Taiwan [30]. The fuel and other large resources needed on the islands are transported by ship from Luzon. At present, the population of Batan Island is about 12,000 and the main economic activities are sightseeing, agriculture and fishery.

**Figure 1.** Geographic location of Batan Island.

The power on Batan Island is offered 24 h a day by the Basco diesel power plant of the Philippine National Power Corporation (NPC), as shown in Figure 2. A total of five generators are mainly usied, of which four generators are 600 kW and the one is 450 kW. Because of the problems such as the difficulty of diesel transportation, high pollution and high power generation cost, the NPC is actively seeking to use renewable energy to solve these issues. With a subsidy of the Asian Productivity Organization (APO), the APO Center of Excellence on Green Productivity (APO COE GP) and NPC have conducted renewable energy research, and analyzed hybrid renewable energy systems based on the load demand and climatic conditions on Batan. The proposed capacity allocation scheme in different power supply systems can be as a reference of investment construction for NPC or private power plants.

**Figure 2.** The Basco diesel power plant.

The rest of this paper is arranged in the following way: In Section 2, the materials and methods followed are discussed. Hybrid energy system descriptions are discussed in Section 3. In Sections 4 and 5 component cost and financial assumptions, results and discussions are described, respectively. Conclusions are discussed in Section 6.

#### **2. Materials and Methods**

The study cases in the research contain seven kinds of power systems. Using entropy weight and TOPSIS method to analyze the techno-economic and evaluate optimal capacity configuration for hybrid energy system (HES). The discussions are as follows:


The required information for our simulation analysis include the status of the diesel generators, the load power consumption, the statistics of the oil prices, interest rates and inflation rates in the Philippines in recent years. In addition, the weather data on the island can be obtained on the NASA website, and e market surveys for equipment costs are also needed.

• Based on the load conditions and the specifications of existing diesel generator, the electrical and economic results of the stand-alone diesel system (Case 1) can be analyzed by HOMER (Pro.3.7.6.0, HOMER Energy, Boulder, CO, USA)) and the result is considered as a reference for discussing different hybrid energy systems.


#### *2.1. Simulation Software Description*

The HOMER software was developed by the U.S. National Renewable Energy Laboratory (NREL). The software is helpful in analyzing the power system's electricity and economy to model the optimal power grid. Furthermore, the user can define the input parameters and the constraints as a reference for modeling profitable power system [31]. The architecture diagram of the simulation analysis procedure of the HOMER software is as shown in Figure 3.

**Figure 3.** HOMER software simulation analysis process architecture diagram.

#### *2.2. Global Horizontal Irradiation and Wind Speed*

Entering the latitude and longitude of Batan into the HOMER software, NASA's global horizontal irradiation data (GHI) and wind speed statistics can be obtained through the Internet connection. The obtained information contains the annual GHI is 1876 kWh/m2/yr, and the daily average GHI is about 5.14 kWh/m2/d. The sunshine is abundant and suitable for the development of solar photovoltaic applications. As shown in Figure 4, the varies for the daily average GHI of each month from 3.31 kWh/m2/d to 6.41 kWh/m2/d, and the highest and lowest average GHI are in April and December, respectively.

**Figure 4.** Global horizontal irradiation-monthly data.

According to the NASA anemometer data measured 50 m above the surface of the Earth, the annual average wind speed is 7.22 m/s. As shown in Figure 5, the variation of average wind speed is from 4.95 m/s to 10.04 m/s in each month. The highest and lowest average wind speeds are in December and May, respectively. The monthly trend of the whole year is opposite to that of GHI. The solar energy and wind energy can be supplied power in turn and effectively use natural resources.

**Figure 5.** Wind speed monthly data.

#### *2.3. Load Profile*

Through the data collection of the frequency distribution of the load power consumption, the average of daily electricity consumption is 16,974 kWh/d, July is the month with higher electricity consumption, and January is the lower month for electricity consumption. The load profile for each month is as shown in Figure 6.

**Figure 6.** Daily load profile for each 12 months.

Figure 7 shows load power consumption frequency distribution. In the sampling period of hourly, the frequency of in 550 kW–599 kW and 600 kW–649 kW is more than 10%, the total is 20.35%; the frequency exceeding 1000 kW is only 7.02%, and the maximum 1442 kW.

**Figure 7.** Load power consumption frequency distribution.

#### *2.4. Diesel Price Data*

As shown in Figure 8, the oil price information regularly announced by the Philippines' Department of Energy in recent years, the Philippines is affected by fluctuations in international oil prices. From October 2015 to February 2019, domestic oil prices have risen sharply. The retail common price has reached 49.15 PHP/L in 15 October 2018 [32], converted to \$ is 0.91 \$/L (1 & = 52.26 PHP). The lowest oil price during the statistical period is 20.5 PHP/L, the maximum is 49.15 PHP/L, and average is 33.335 PHP/L. In the offshore islands, the transportation cost must be added to the oil price. Therefore, the actual oil price is according to the geographic location. In more remote areas, the oil price will be affected by the transportation cost. The actual oil price may be 1–1.5 times the estimated oil price [33], in this paper, the average price of the statistical period is 1.25 times the oil price (0.8 \$/L), which is used as the basic parameter of the simulation.

**Figure 8.** Diesel fuel price data.

#### *2.5. Assignment Indexes*

The evaluation indicators used in this paper include economic indicators (Sections 2.5.1–2.5.4), electrical indicators (Sections 2.5.5 and 2.5.6), and occupied area (Section 2.5.7). The calculation formulas and descriptions are as follows [3,34]:

#### 2.5.1. Annual Real Interest Rate

Annual real interest rate (*i*) is used to adjust between one-time cost and annualized costs. HOMER uses annual real interest rate to compute discount factor and to determine annualized costs from present costs:

$$i = \frac{i'-f}{1+f} \tag{1}$$

where *i* is annual real interest rate (%); *i '* is nominal interest rate (bank board rate) (%); *f* is expected inflation rate (%).

#### 2.5.2. Net Present Cost

The total of NPC value denotes the cost of system life cycle in HOMER. Equation (2) displays the amount of cash flow of t-year over the factor and initial capital cost. Costs include capital cost, operation cost, replacement cost, maintenance cost, fuel cost and etc. The income includes the selling of electricity and the residual price after the life cycle:

$$\begin{split} \dot{i} &= \frac{\dot{i}^{\prime} - f}{1 + \dot{f}} \text{NPC} = \text{CF}\_{0} + \left\{ \frac{\text{CF}\_{1}}{(1 + \dot{i})^{1}} + \frac{\text{CF}\_{2}}{(1 + \dot{i})^{2}} + \frac{\text{CF}\_{3}}{(1 + \dot{i})^{3}} + \dots + \frac{\text{CF}\_{N}}{(1 + \dot{i})^{N}} \right\} \\ \dot{i} &= \text{CF}\_{0} + \sum\_{t=1}^{N} \frac{\text{CF}\_{t}}{(1 + \dot{i})^{t}} \end{split} \tag{2}$$

where *CFt* is cash flow of t-year (according to HOMER software: expenditure is positive and income is negative.) (\$), *i* is the annual real interest rate (%), *N* is project life time (yr), *t* is number of year (yr), *CF*<sup>0</sup> is initial capital cost (\$).

#### 2.5.3. Capital Recovery Factor (CRF)

*CRF* is ratio used to determine present value of the annuity during the project lifetime:

$$\text{CRF}(i, t) = \frac{i(1+i)^t}{(1+i)^t - 1} \tag{3}$$

where *t* is the number of years, *i* is the annual real interest rate (%).

#### 2.5.4. Cost of Energy

HOMER defines levelized COE as average cost per kWh of useful electric energy generated by the system and computed by dividing total annualized cost (*TAC*) by total annualized useful electric energy generation. The unit of COE is \$/kWh. *TAC* is annualized value of NPC, and its unit is \$/yr. The relation is as follows:

$$TAC = \text{NPC} \times \text{CRF}(i, N) \tag{4}$$

$$\text{COE} = \frac{\text{TAC}}{E\_{\text{prin}}} \tag{5}$$

where *Eprim* is annualized primary served load (kWh/yr), *N* is project life time (yr).

#### 2.5.5. Renewable Fraction (RF)

RF is the fraction of energy delivered to the load that originated from renewable power sources. The equation is as follows:

$$\text{RF} = \left(1 - \frac{E\_{\text{non-ren}} + H\_{\text{non-ren}}}{E\_{\text{servad}} + H\_{\text{servad}}}\right) \times 100\% \tag{6}$$

where *Enon-ren* is the nonrenewable electrical production (kWh/yr), *Hnon-ren* is the nonrenewable thermal production (kWh/yr), *Eserved* is the total electrical load served (kWh/yr), *Hserved* is the total thermal load served (kWh/yr).

#### 2.5.6. Excess Electricity Fraction

The excess electricity fraction is the ratio of total excess electricity to total electrical production. The equation is as follows:

$$\text{excess electricity fraction} = \frac{E\_{\text{excess}}}{E\_{\text{prod}}} \times 100\% \tag{7}$$

where *Eexcess* is the total excess electricity (kWh/yr), *Eprod* is the total electrical production (kWh/yr).

#### 2.5.7. Occupied Area

In the calculation of the occupied area, the main considerations are PV array, wind turbine and energy storage [3]. The equation is as follows:

$$\text{Occupised area} = (Ap\nu \times \text{N}p\nu) + (A\mu\_{\text{WT}} \times \text{N}p\_{\text{WT}}) + (A\_{\text{Bat}} \times \text{N}\_{\text{Bat}}) \tag{8}$$

where *APV* is the PV array area per kWp (m2), *NPV* is the capacity of PV system (kWp), *AWT* is the occupied area of one wind turbine (m2), *NWT* is the number of wind turbines, *ABat* is occupied area of one energy storage rack (m2), *NBat* is the number of energy storage racks.

#### *2.6. Entropy Weight and TOPSIS Method*

Entropy weight method is derived from Shannon entropy, that was suggested by Shannon for the numerical measurement of uncertainty in information systems [35]. The weighting factor of this evaluation method depends entirely on the value of the indicator, not on the subjective evaluation. Therefore, it has been considered as an objective method of calculating weights and has been widely applied [3].

TOPSIS is a method for relatively evaluation in limited alternative. This method is evaluated by computing relative closeness of the alternative and the ideal solution, and sorting according to the relative closeness. The higher closeness means the solution is close to the positive ideal solution, that is, the farther away from negative ideal solution, which is the best solution [36].

This study combines these two methods, first using the entropy weight method to calculate the objective weight value of each indicator, and using the TOPSIS method to calculate the relative closeness with the positive ideal solution for the weighted alternative. The calculated results are sorted so that the best solution can be found [37]. The main calculation steps and formulas in this paper are as follows:

#### *Step 1: Initialize matrix*

Construct m alternatives, the initial data matrix V of *n* index, and the eigenvalue of the *j*th index of the *i*th alternatives is expressed as *vij*, so matrix V can be defined as follows:

$$\mathbf{V} = \begin{pmatrix} v\_{\bar{l}\bar{j}} \end{pmatrix}\_{\mathbf{m}\times\mathbf{n}} = \begin{bmatrix} v\_{11} & \cdots & v\_{1n} \\ \vdots & \ddots & \vdots \\ v\_{m1} & \cdots & v\_{mn} \end{bmatrix} \tag{9}$$

#### *Step 2: Normalize indexes*

The min-max normalization formula is used to normalize the index. If the index is positive, it is used for the benefit type. The higher value is the better. The formula is as follows:

$$r\_{ij} = \frac{v\_{ij} - \min\_{j}(v\_{ij})}{\max\_{j}(v\_{ij}) - \min\_{j}(v\_{ij})} \tag{10}$$

If the index is negative, it is used for cost type. The lower the value the better. The formula is:

$$r\_{ij = } \frac{\max\_j(v\_{ij}) - v\_{ij}}{\max\_j(v\_{ij}) - \min\_j(v\_{ij})} \tag{11}$$

The matrix after normalizing can be defined as:

$$\mathbf{R} = \left( r\_{\vec{i}\vec{j}} \right)\_{m \times n} \tag{12}$$

*Step 3: Calculate proportion and information entropy for each index*

Equation (13) is used to calculate the proportion of each *j* index in each *i* alternative:

$$P\_{ij=} \frac{r\_{ij}}{\sum\_{i=1}^{m} r\_{ij}} \; , j = 1, 2, \dots, n \tag{13}$$

Then calculate the j-index information entropy *ej* by using Equation (14), and use the value of *ej* to judge the discrete degree of the index. When *ej* is smaller, the discrete degree is greater, and the weight as well:

$$\sigma\_{\vec{j}} = -\mathbb{K} \times \sum\_{i=1}^{m} P\_{ij} \cdot \ln p\_{ij}, \vec{j} = 1, 2, \dots, n \tag{14}$$

where the constant K = <sup>1</sup> *lnm*

*Step 4: Calculate weight for each index*

Equation (15) is used to calculate the weight ω*<sup>j</sup>* of the *j*th index:

$$\omega\_{j} = \frac{1 - e\_{j}}{n - \sum\_{j=1}^{n} e\_{j}}, j = 1, 2, \dots, n \tag{15}$$

where *<sup>n</sup> j* ω*<sup>j</sup>* = 1 and 0 ≤ ω*<sup>j</sup>* ≤ 1. The larger ω*<sup>j</sup>* means the higher importance of the objective response of the j indicator.

*Step 5: Calculate the weight of normalized matrix R*

Multiply the *j* index corresponding to each *i* alternatives by ω*<sup>j</sup>* . The matrix Y can be calculated by Equation (16):

$$\mathbf{Y} = \begin{pmatrix} r\_{ij} \end{pmatrix}\_{m \times n} \times \boldsymbol{\omega}\_j = \begin{bmatrix} r\_{11} \cdot \boldsymbol{\omega}\_1 & \cdots & r\_{1n} \cdot \boldsymbol{\omega}\_n \\ \vdots & \ddots & \vdots \\ r\_{m1} \cdot \boldsymbol{\omega}\_1 & \cdots & r\_{mn} \cdot \boldsymbol{\omega}\_n \end{bmatrix} \tag{16}$$

*Step 6: Compute distance between each i alternatives and positive ideal solution D*<sup>+</sup> *<sup>i</sup>* and negative ideal solution *D*− *i .*

The Euclidean Distance can be calculated as follows:

$$D\_i^+ = \sqrt{\sum\_{j=1}^n \left( y\_{ij} - y\_j^+ \right)^2}, i = 1, 2, \dots, m \tag{17}$$

$$D\_i^- = \sqrt{\sum\_{j=1}^n \left( y\_{ij} - y\_j^- \right)^2}, i = 1, 2, \dots, m \tag{18}$$

$$\epsilon\_i = \epsilon\_j \tag{19}$$

where *y*<sup>+</sup> *<sup>j</sup>* = max 1≤*j*≤*n yij i* = 1, 2 ··· *m* and *y*− *<sup>j</sup>* = min 1≤*j*≤*n yij i* = 1, 2 ··· *m* .

*Step 7: Calculate relative closeness for each alternative*

Through the calculation of the relative closeness (*Ci*), it can be judged that alternative is close to ideal solution. If the *Ci* is close to 1, the *i*th alternative is close to the ideal solution:

$$\mathbf{C}\_{i} = \frac{D\_{i}^{-}}{D\_{i}^{+} + D\_{i}^{-}}, i = 1, 2, \dots, m \tag{19}$$

*Step 8: Sorting the Ci calculated by each i alternative, the alternative has the maximum Ci is the optimization alternative in this study.*

#### **3. Hybrid Energy System Description**

#### *3.1. Hybrid Energy System Schematic Diagram*

The schematic diagram of the hybrid energy system shown in Figure 9 consists of diesel generators, PV system (PVS), wind-power generation system (WGS), storage system, power conversion system (PCS) and load. PVS and WGS are operated in parallel with diesel generators and AC coupled. Electricity generated by PVS and WGS can be delivered to AC load to decline the output of diesel generator and to minimize fuel consumption [38,39].

**Figure 9.** Schematic diagram of the hybrid energy system.

#### *3.2. PV System*

The way of electricity generation by PV module is through the conversion of solar energy into DC power. Identification, adoption and utilization of dependable interconnection technology to assembly crystalline silicon solar cells in PV module are very important to guarantee the device functions persistently up to 20 years of design life span [33]. For this study, the PV module selected is the GTEC-G6S6A model mono-crystalline silicon solar cell type which specifications are presented in Table A1 of Appendix A. South installation direction, 10◦ incline, 5.42 m<sup>2</sup> per kWp, 25 years of serving

life and 80% derating factor are considered. Equation (20) is used by HOMER to calculate the power output of PV system [40,41]:

$$P\_{\rm pv} = P\_{\rm pv,STC} f\_{\rm pv} \frac{G\_T}{G\_{T,STC}} \left[ 1 + K\_P (T\_C - T\_{STC}) \right] \tag{20}$$

where *Ppv,STC* is PV system rated power (kWp), *fpv* is PV derating factor (%), *GT* is solar irradiance on the surface of the PV module (kW/m2), *GT,STC* is standard solar irradiance (1 kW/m2), *KP* is the PV module temperature coefficient (−0.4003%/ ◦C), *TC* is temperature of the PV module (◦C), *TSTC* is temperature of PV module under the standard test conditions (25 ◦C).

#### *3.3. WG System*

The WGs converts wind kinetic energy into electricity using wind turbine blades. The architecture of WG system in the study is AC coupled hybrid. The wind turbine used for simulation is the Bergey XL 10 Wind Turbine complete set [42] with AC output, rated power 10 kW at 12 m/s, minimum wind speed 2.5 m/s and swept area 38.5 m2. The hub height is 30m, power curve of wind turbine is displayed in Figure 10, and the lifetime is set to 20 years [43].

**Figure 10.** Wind turbine power curve.

The calculation of the power law profile in HOMER is based on Equation (21) to estimate the wind speed of different hub heights [44]:

$$\mathcal{U}I\_{\text{hub}} = \mathcal{U}I\_{\text{ref}} \times \left(\frac{H\_{\text{hub}}}{H\_{\text{ref}}}\right)^{\alpha} \tag{21}$$

where *Uref* is the reference wind speed (m/s), the value is provided by NASA measures at a height of 50 m from the surface. *Hhub* is the hub height (m). *Href* is the measured height of the reference wind speed (m). α is surface roughness, the value is set to 0.14 [44,45].

#### *3.4. Storage System*

A storage system will be used to regulate the power and minimize the influence the quality of energy on intermittent renewable energy sources in the proposed HES. It accumulates and transfers energy from PVs, WGs and diesel to a rechargeable Li-ion thin film during excess electricity generation. In this paper a high energy density Samsung SDI M2-R084 model lithium-ion battery [46] was selected as storage system with the following specifications: The lowest expected lifespan of each rack is five years and each rack holds 11 modules. Each module has 22 series-connected cells. The capacity is 94 Ah and the rated capacity of each storage rack is 84 kWh. The area occupied by one rack is 0.31 m2. The voltage rating is 774–1004 V and the lowest state of charge (SOC) is 20%. The C rate is 1 C and lifetime throughout is 337,123 kWh with 4000 cycles.

#### *3.5. Power Conversion System*

A dual-direction DC-AC converter is used for power conversion. In HES, when excess electricity is produced by diesel, PV and WG, the excess AC power will be changed into DC and stored in battery to supply it back to AC load when required. Conversion performance of 95% and 10 years of lifespan is assumed in the simulation.

#### *3.6. Diesel Generators*

The power supply system of the island discussed in this study uses diesel generators (DGs) to supply power 24 h a day. The main power supply, with a total capacity of 2850 kW, contains five DGs and the supply schedule is accordingly to meet the load demand.

To simplify the process of simulation analysis, set operation time of DG1 and DG2 to 0–11 a.m., operation time of DG3 and DG4 to 12–23 p.m., DG5 to backup, minimum load to 25%, and lifetime to 131,400 h. The simulation parameters of diesel generators are as shown in Table 2.


**Table 2.** Simulation parameters of the diesel generators.

#### *3.7. System Dispatch Strategy*

System dispatch approach in HOMER is primarily for managing the operational rules of the power generation equipment and energy storage. The simulation software has two options: cycle charging strategy and load following strategy and [34]. HES operation in this case study is based on a DG power supply. PV and WG are used as auxiliary power supplies. Increasing the penetration of the PVs and WGs reduces fuel consumption. The DGs only generate sufficient power to fulfill load demand. DG has the limited base load output, so when PVs and WGs have excess power, the storage will charge, but if the energy storage is full and there is still excess energy that cannot be stored or used, the software will consider that to be excess electricity, as shown in (7). As in the actual operation of the HES, in the above scenario, the output of the PVs or WGs power converter must be controlled to reduce their output to maintain balance and/or damping load is required. Based on the operation mode, the dispatch mode of the power supply system must select the load following strategy and so as to calculate various costs such as fuel, operation and maintenance, and replacement.

#### **4. Component Cost and Financial Assumption**

#### *4.1. System Component Cost*

The summary of component costs used in this study is from the Taiwan System Integration Company, and is presented in Table 3. The cost consists of capital cost, replacement cost, operation and maintenance (O&M). In order to simplify the analysis, the capital and replacement cost is the same.

Capital cost of PVs is 4000 \$/kW, including PV module, PV inverter, transportation, installation engineering. Capital cost of WGs is 5800 \$/kW, including wind turbine, wind inverter, 30 m tower and cement foundation, transportation, installation engineering and so on. PVs, WGs and storage operate in parallel with the existing five DGs. To simplify the analysis, the capital cost of DG is not included in the calculation. The estimates of the cost of remaining components include transportation, and installation engineering.


**Table 3.** Summary of component costs.

#### *4.2. Interest Rate and Inflation Rate*

Figures 11 and 12 depict the interest rate and inflation rates of Bangko Sentral NG Pilipina (BSP) in recent years. Interest rates have been adjusted to 4.75% in 7 February2019 [47]. Refering to the BSP announcement in February 2019, the inflation rate is 4.4% [48]. Since January 2013 to February 2019, the highest interest rate was 4.75%, the lowest was 3%, so the average was 3.579%. The highest inflation rate was 6.7%, the lowest was −2%, and the average was 2.716% [48,49]. In this paper, the real discount rate is the average value in the statistical period, and the value calculated by (1) is 0.84%. The lifetime is set to 20 years.

**Figure 11.** Interest rate information for Taiwan since January 2013 to February 2019.

**Figure 12.** Inflation rate information for Taiwan since January 2013 to February 2019.

#### **5. Simulation Results and Discussion**

The flowchart of systems operation for cases 1 to 7 are presented in Figure A1 of Appendix A.

#### *5.1. Case 1: Stand-Alone Diesel System*

Simulation results are as shown in Table 4. Total annual fuel intake of the five DGs is 1,623,100 L/y (4447 L/d), the total annual power supply is 6,195,877 kWh/y, average fuel intake per kWh is 0.26 L/kWh, and annual total carbon dioxide (CO2) emission is 4,274,156 kg/y.


**Table 4.** Electrical characteristics of stand-alone diesel system.

No capacity shortage problem throughout the year and it meets the annual load requirements. DG5 is an auxiliary generator, which starts when the load demand reaches 90% of the rated capacity of two DGs to avoid any power supply shortage. The load ratio of DG5 is lower than that of other diesel generator sets, so the fuel consumption required per kWh (liter/kWh) is also relatively higher than that of other diesel generator sets. As indicated in Table 5, total NPC of the five DGs is \$29,650,884 and the total COE is 0.2609 \$/kWh. Among costs, the cost of oil is the largest and it takes 80% of the total NPC.

**Table 5.** Economic characteristics of stand-alone diesel system.


#### *5.2. Hybrid Energy System without Storage*

#### 5.2.1. Case 2: PV-Diesel HES

The PV capacity configuration scheme for minimum COE of a PV-diesel HES is analyzed and the change of PV capacity is set from 50 to 2000 kWp and the step interval is 50 kWp in the simulation. The simulation results are shown in Figure 13. Because the analyzed system does not include an energy storage system, as the PV capacity rises, the percentage of renewable fraction (RF) and fuel saving will improve. The increment will reduce in accordance with the increase of PV capacity and the excess electricity fraction will increase relatively.

**Figure 13.** Effect of varying PVS capacity on fuel saving, RF and excess electricity fraction.

The varying PV capacity on COE and RF are as shown in Figure 14. When the capacity of PVs is 550 kWp and PV operates with the DG in parallel, the lowest COE is 0.2583 \$/kWh, which is slightly lower than the COE (0.2609 \$/kWh) in the DG-only case. As shown in Table 6, the RF is 12.4%, the excess electricity fraction is 0.2% and the fuel saving is −11.17%. When the PV capacity continuously increases, the investment cost will also increase. Because the power generated by PVs is not entirely used by loads to minimize the DG fuel consumption, the COE will increase relatively. Table 7 shows economic characteristics of the lowest COE PV-diesel HES. Cost of oil is still the largest proportion, which is 72% of the total NPC. The value of NPC of PVs with 550 kWp total capacity is only 8% of the total NPC. Figure 15 shows the annual PVS power production by the lowest COE PV-diesel HES. The output power of quarter 2 (April to June) and quarter 3 (July to September) can be observed, which is higher than quarter 1 (January to March) and quarter 4 (October to December).

**Figure 14.** Effect of varying PVS capacity on COE and RF.


**Table 6.** Electrical characteristics of lowest COE PV-diesel HES.


**Table 7.** Economic characteristics of the lowest COE PV-diesel HES.

**Figure 15.** Annual PVS power production by lowest COE PV-diesel HES.

#### 5.2.2. Case 3: Wind-Diesel HES

Using the analysis method mentioned in Section 5.2.1, we set the change of WG capacity between 200–800 kW and the interval at 10 kW to find the configuration with the lowest COE. As presented in Figure 16, trend of the curve is like the PV-diesel HES one. When the capacity of WGs is 380 kWp and with DG in parallel operation, the minimum COE calculated is 0.2541 \$/kWh and this result is slightly less than COE in Case 1 and Case 2, as shown in Figure 17.

**Figure 16.** Effect of varying WGS capacity on fuel saving, RF and excess electricity fraction.

**Figure 17.** Effect of varying WGS capacity on COE and RF.

As shown in Table 8, the capacity factor of WGs is 32.28%, which is higher than 16.32% of PVs. Therefore, the installation capacity of WG is lower than PV but the system still has lower COE and higher RF: RF is 13.59%, the excess electricity fraction is 2.2%, and the fuel saving is −13.59%.



Table 9 shows the cost simulation of the lowest COE wind-diesel HES in 380 kW total of capacity. The COE is lower than the COE of DG 1–5 and that of PVs. Figure 18 shows the distribution of annual WGs output power and displays that output power of quarter 1 and 4 is higher than that of quarter 2 and 3. The trend is opposite to the result of the curve of quarter output power of PVs.


**Table 9.** Economic characteristics of the lowest COE wind-diesel HES.

**Figure 18.** Annual WGS power production by the lowest COE wind-diesel HES.

#### 5.2.3. Case 4: PV-Wind-Diesel HES

As shown in Figures 15 and 18, PVs and WGs have different seasonal power generation characteristics, so the hybrid power supply can be complementary to improve RF. After simulation analysis, the hybrid power system has PVs of 200 kWp and WGs of 340 kW, which is the lowest capacity configuration. As depicted in Table 10 the analysis of electrical characteristics, RF is 18.1%, fuel saving is −16.18%, which is better than Case 1 and 2. In Table 11, the economic characteristics analysis indicates that COE is 0.2539 \$/kWh and result is slightly lower than Case 3. In the comparison, the result of Case 4 has the lowest COE in all of capacity configuration without the energy storage HES.


**Table 10.** Electrical characteristics of the lowest COE PV-wind-diesel HES.


**Table 11.** Economic characteristics of the lowest COE PV-wind-diesel HES.

#### *5.3. Hybrid Energy System with Storage*

#### 5.3.1. Case 5: PV/Diesel/Storage HES

As the capacity of PV or WG increases, the RF will increase, but the excess electricity fraction will also increase. Although it will waste the energy, the problem could be enhanced by installing an energy storage. Capacity configuration for PV-diesel-storage HES of the lowest COE in each RF level addressed in this section. The setting parameters are as follows: RF is 25%–50%, the interval is 5% and the excess electricity fraction must be less than 5% considering the utilization rate of the energy.

Table 12 lists the electrical characteristics of PV-diesel-storage HES. The result displays that higher RF can let the effect of fuel saving better.

**Table 12.** Electrical characteristics of PV-diesel-storage HES.


Table 13 presents the economic characteristics of PV-diesel-storage HES. In RF 25% hybrid system configuration: PVs capacity is 1200 kWp, energy storage capacity is 1260 kWh and operates in parallel with DG. Above all, the minimum COE can be obtained is 0.2864 \$/kWh. The cost and storage capacity have the significant difference in RF 45% and 50% because the simulation constraint should be satisfied: The capital cost difference is 1.35 times and the COE difference is 1.13 times.


**Table 13.** Economic characteristics of PV-diesel-storage HES.

#### 5.3.2. Case 6: Wind-Diesel-Storage HES

Simulation analysis shows that WG power gen ration was concentrated in quarters 1 and 4, and the load power consumption during this period is lower than that of quarters 3 and 4. When the installed capacity of WGs and storage is more than RF 35%, the excess electric fraction has exceeded 5% and there is no system capacity configuration meets the set simulation conditions.

Tables 14 and 15 are the simulation results of the electrical and economic characteristics, respectively. In RF 25% hybrid system configuration: The capacity of WGs is 680 kW, storage capacity is 2520 kWh and operates in parallel with DG. Above all, the lowest COE can be obtained is 0.2874 \$/kWh. The cost and storage capacity have the significant difference in RF 30% and 35% because the simulation constraint should be satisfied. The cost of the system configuration has increased significantly in RF 35%.


**Table 14.** Electrical characteristics of wind-diesel-storage HES.

**Table 15.** Economic characteristics of wind-diesel-storage HES.


#### 5.3.3. Case 7: PV-Wind-Diesel-Storage HES

The simulation results of the electrical and economic characteristics of the PV-wind-diesel-storage HES are shown in Tables 15 and 16, respectively. The capital cost, NPC and COE in case 7 are lower than in cases 5 and 6 when RF is 25–50%.

In RF 25% hybrid system configuration: The minimum capital cost is \$5,317,506, the NPC is \$31,100,367 and the COE is 0.2737 \$/kWh. Because solar energy and wind energy have characteristics of seasonal power complementary, the required capacity of PVs, WGs and storage are lower than that of Cases 5 and 6, are 400 kWp, 440 kW and 168 kWh, respectively.


**Table 16.** Electrical characteristics of PV-wind-diesel-storage HES.

5.3.4. Optimal Capacity Configuration Analysis of the PV-Wind-Diesel-Storage HES

As mentioned in Section 5.3.3, if the cost is the main indicator, Case 7 has the lowest COE HES. This section uses the entropy weight and TOPSIS method and regards the fuel saving, capital cost, and COE as the indicators. Because the available space on the island is limited, the area occupied by the equipment is also considered. Above all, we evaluate the optimal hybrid system configuration for Case 7 at RF 25%−50%. Substituting the values from Table 17 into Equations (9)–(15), the information entropy and weighting of each index can be obtained. The weight of fuel saving is the highest, as shown in Table 18. Through the Equations (16)–(19), the relative closeness value and sorting result of each alternative can be obtained, as shown in Table 19. Therefore, the best system capacity configuration for this stud: The maximum relative closeness is the system capacity configuration in RF 35%, the PVs capacity is 950 kWp, the WGs capacity is 410 kW, and the energy storage capacity is 1680 kWh.

**Table 17.** Economic characteristics of the PV-wind-diesel-storage HES.


**Table 18.** Information entropy and weighting factors of each index.


**Table 19.** The relative closeness and ranking results.


#### *5.4. Sensitivity Analysis*

According to the analysis result in Section 5.3.4, the GHI, wind speed, diesel fuel price, and load consumption for RF 35% PV-wind-diesel-storage HES are used to discuss the economic and electrical of the system.

#### 5.4.1. Global Horizontal Irradiation and Wind Speed

We set the value of GHI and wind speed 5.14 kW/m2/d and 7.22 m/s in <sup>±</sup>20% changes, respectively, i.e., GHI from 4.09 to 6.17 kW/m2/d, and wind speed from 5.78 m/s to 8.66 m/s. According to Figures 19 and 20, when GHI or wind speed increases, the RF will continue to increase and this will let the DG fuel consumption and COE reduce. When GHI is 6.17 kW/m2/d and wind speed is 8.66 m/s, RF can rise to 40.7% and COE can be reduced to 0.2747 \$/kWh.

**Figure 19.** Effect of varying GHI and wind speed on RF generated by optimal PV-wind-dieselstorage HES.

**Figure 20.** Effect of varying GHI and wind speed on COE generated by PV-wind-diesel-storage HES.

#### 5.4.2. Diesel Fuel Price

We set the diesel fuel price change from 0.49–1.18 \$/liter, which is 1.25 times the maximum and minimum values during the statistics period of oil price, as the change range of the oil price. Figure 21 show the impact on NPC and COE under different fuel price scenarios. When oil price is 1 \$/L, NPC and COE are as high as \$36,434,464 and 0.3207 \$/kWh, respectively.

**Figure 21.** Effect of varying diesel fuel price on NPC and COE.

#### 5.4.3. Load Consumption

We set the load consumption to vary from 16,974 to 20,000 kWh/d. As depicted in Figure 22, when the load consumption rises and the capacity configuration of PVs, WGs and storage do not change, the DGs will increase their power generation to meet the load demand, so the RF will decrease and the fuel usage will increase. As shown in Figure 23, when the load consumption increases, the fuel consumption, the fuel cost and the NPC rise because the power generation of DG increases. It can also be observed from (5) that the load consumption increases so the COE decreases.

**Figure 22.** Effect of varying load consumption on fuel consumption and RF.

**Figure 23.** Effect of varying load consumption on NPC and COE.

#### **6. Conclusions**

Promoting the infrastructure of renewable energy is a global trend, especially in areas with insufficient power supply and low electrification. Renewable energy is regarded as a key to improving economic and development. This paper uses the Philippine islands, and a renewable energy and hybrid power supply to supplement existing diesel generators, to be the study case.

We use HOMER to simulate and investigate the techno-economic of different power supply systems' containing stand-alone diesel, PV-diesel HES, wind-diesel HES, PV-wind-diesel HES, PV-diesel-storage HES, wind-diesel-storage HES, PV-wind-diesel-storage HES. The study also uses entropy weight and TOPSIS method to evaluate optimal capacity configuration. The conclusions of the analysis are as follows:


to fulfill the demand. This causes the RF to decline, the fuel consumption and NPC the rise, and COE conversely.

**Author Contributions:** Conceptualization, C.-C.K.; Data curation, W.-B.W.; Formal analysis, C.-T.T. and T.M.B.; Investigation, W.-B.W. and C.-C.K.; Methodology, C.-T.T. and T.M.B.; Resources, C.-T.T.; Supervision, C.-C.K.; Writing—original draft, W.-B.W.; Writing—review & editing, T.M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** Supported of this research by the Asian Productivity Organization Center of Excellence on Green Productivity (APO COE GP) Project and Ministry of Science and Technology of the Republic of China (Grant No. MOST108-3116-F-042A-004-CC2) are gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**


**Table A1.** Electrical parameters of the PV modules.

**Figure A1.** Flowchart of system operation in various cases 1–7.

#### **References**

1. Ma, T.; Javed, M.S. Integrated sizing of hybrid PV-wind-battery system for remote island considering the saturation of each renewable energy resource. *Energy Convers. Manag.* **2019**, *182*, 178–190. [CrossRef]


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Integrated Analysis of Influence of Multiple Factors on Transmission E**ffi**ciency of Loader Drive Axle**

#### **Jianying Li 1,\*, Tunglung Wu 1, Weimin Chi 1, Qingchun Hu <sup>2</sup> and Teenhang Meen 3,\***


Received: 2 October 2019; Accepted: 22 November 2019; Published: 28 November 2019

**Abstract:** In this study, a loader drive axle digital model was built using 3D commercial software. On the basis of this model, the transmission efficiency of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear of the loader drive axle were studied. The functional relationship of the transmission efficiency of the loader drive axle was obtained, including multiple factors: the mesh friction coefficient, the mesh power loss coefficient, the normal pressure angle, the helix angle, the offset amount, the speed ratio, the gear ratio, and the characteristic parameters. This revealed the influence law of the loader drive axle by the mesh friction coefficient, mesh power loss coefficient, and speed ratio. The research results showed that the transmission efficiency of the loader drive axle increased with the speed ratio, decreased when the mesh friction coefficient and the mesh power loss coefficient increased, and that there was a greater influence difference on the transmission efficiency of the loader drive axle.

**Keywords:** drive axle; differential planetary mechanism; power loss coefficient; friction coefficient

#### **1. Introduction**

The loader is widely used as a construction machine, and the drive axle is an assembly of the loader drive system. The loader drive axle is mainly composed of a main reducing gear, a differential planetary mechanism, and a wheel planetary reducing gear. The transmission efficiency of the loader drive axle is a topic of interest among scholars. Xu and Kahraman [1] proposed a model of hypoid gears for predicting friction power loss, and their research showed that the friction power loss of a hypoid gear is affected by the gear smoothness, lubrication characteristics, and other factors. Fanghella et al. [2] discussed the functional characteristics of a nutating gearbox with bevel gears, for which the expressions of the transmission ratio and efficiency were obtained, including the main design. Kakavas et al. [3] proposed a thermal coupled transient model of a hypoid gear for power loss calculation and analyzed and compared numerical simulations and experimental results. Paouris [4] presented a tribodynamic model of differential hypoid gears, for which the lubricated contact of the meshing gear pair was integrated with lubricants of varying rheological properties, and the influence of lubricant formulation and gear geometry on system efficiency was examined. The results showed that there is a trend of increasing power loss when using lubricants with higher pressure–viscosity coefficients. Simon [5] found that the efficiency of a hypoid gear pair was improved by reducing the pressure and temperature in the oil film; the results were obtained based on the conditions of the mixed elastohydrodynamic lubrication (EHL) characteristics affected by the manufacturing parameters. To improve the transmission efficiency of a loader, Lyu et al. [6] designed an output power-split

transmission system, in which a planetary gear set was applied as the dynamic coupling element to allow the output power of the loader to be split-transmitted. Xiong et al. [7] introduced a sizing approach for both hydrostatic–mechanical power-split transmission (PST) and hydrostatic transmission (HST) systems and presented a multidomain modeling approach for the powertrain of a wheel loader. The sizing method and simulation models should facilitate the development of powertrains for wheel loaders and other wheeled heavy vehicles. Barday et al. [8] analyzed the power losses of a truck drive axle, including the sliding and tooth friction of hypoid gear sets, rolling element bearings and oil churning, and then reviewed thermal exchanges and compared the bulk temperatures of the gear set.

The concept of virtual power was first proposed by Chen et al. [9] and was used to study the transmission efficiency of a simple planetary gear train. Then, on the basis of the virtual power theory, Chen et al. [10] analyzed the transmission efficiency of a dual-input, single-output planetary gear train. In Ref. [11–13], the virtual power theory was further developed to study the transmission efficiency of a compound planetary gear train. Li et al. [14] studied the power flow of a compound planet gear train with bevel gears via virtual power theory, and they analyzed the parameters' influence on transmission efficiency. In Ref. [15], Li et al. established a functional expression of transmission efficiency, including the characteristic parameters, the mesh friction coefficient, and the mesh angle of multistage microplanetary gearing transmission, and adopted the virtual power theory to study the influence of multiple factors on the transmission efficiency of a multistage microplanetary gearing transmission.

Although many studies have been carried out on the transmission efficiency of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear of the loader drive axle, the integrated analysis of the influence of multiple factors on the transmission efficiency of the loader drive axle assembly has not been reported. By establishing a three-dimensional digital model of a loader drive axle, the transmission efficiency of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear are studied, respectively. In doing so, the functional relationship between the transmission efficiency of the loader drive axle and the transmission efficiency of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear is obtained. This revealed the influence law between the transmission efficiency of the loader drive axle and the meshing friction coefficient, power loss coefficient, and speed ratio, which can provide a theoretical foundation for the design of heavy-duty automotive transmission systems.

#### **2. Experimental**

#### *2.1. Digital Model of the Loader Drive Axle*

According to the geometric characteristics of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear of the loader drive axle, the 3D commercial software parametric design method was used to establish the three-dimensional (3D) model of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear of the loader drive axle, as seen in Figure 1. The loader drive axle assembly is shown in Figure 2, in which the numbers 0–8 represent the loader drive axle frame, the main reducing driving gear, the main reducing driven gear, the differential planetary mechanism planet gear, the half shaft gear, the wheel planetary reducing planet gear, and the carrier, respectively. The basic parameters of the loader drive axle are listed in Table 1.

**Figure 1.** Three-dimensional (3D) digital model of the loader drive axle.

**Figure 2.** Schematic diagram of the loader drive axle.


#### **Table 1.** Basic parameters of the loader drive axle.

#### *2.2. Research on Transmission E*ffi*ciency of Loader Drive Axle*

2.2.1. Research on Transmission Efficiency of the Main Reducing Gear

The theoretical transmission efficiency of the main reducing gear spiral bevel gear of the loader axle is shown in Equation (1) [16]:

$$\eta = \frac{1}{1 + \sqrt{\frac{T\_{\text{max}}}{T} \left[\frac{\mu}{\cos \phi} \sqrt{\left(\tan \psi\_p - \tan \psi\_{\mathcal{S}}\right)^2 + \left(0.15 \left(1 - \frac{E}{D}\right)^{16}\right)^2}\right]}}\tag{1}$$

where η is the transmission efficiency; *T* is driven gear load torque; *T*max is the driven gear peak torque at 2.75 times the fatigue limit bending stress; μ is the mesh friction coefficient; φ is the normal pressure angle of the driving gear surface; ψ*<sup>p</sup>* is the spiral bevel gear driving gear nominal helix angle; ψ*<sup>g</sup>* is the spiral bevel gear driven gear nominal helix angle; *E* is the offset of the hypoid gear driving gear; and *D* is the pitch diameter of the hypoid gear.

As reported by Saiki [17], the modified transmission efficiency of the spiral bevel gear does not depend on the load torque, as shown in Equation (2):

$$\eta\_1 = \frac{1}{1 + \left[\frac{\mu}{\cos \phi} \sqrt{\left(\tan \psi\_p - \tan \psi\_{\mathcal{S}}\right)^2 + \left(0.15 \left(1 - \frac{\mu}{\mathcal{D}}\right)^{16}\right)^2}}\tag{2}$$

The nominal helix angle of the driving and driven spiral bevel gears is shown in Equations (3) and (4), respectively. *F* is the spiral bevel gear driven gear tooth width, and the rest of the symbol is in accordance with the meaning of Equation (1):

$$
\psi\_p = 25 + 5\sqrt{\frac{z\_2}{z\_1}} + 90(E/D) \tag{3}
$$

$$
\psi\_{\mathcal{S}} = \psi\_p - \arcsin[2\mathcal{E}/(\mathcal{D} - \mathcal{F})] \tag{4}
$$

2.2.2. Research on Transmission Efficiency of the Differential Planetary Mechanism

The diagram of the differential planetary mechanism of the loader drive axle shown in Figure 2 can be converted into the diagram shown in Figure 3 without considering the gear mesh power loss. The rules for the graphical conversion process are detailed in [14]. The virtual split power ratio β<sup>3</sup> and the split power ratio β <sup>3</sup> of component 3 of the differential planetary mechanism of the loader drive axle are expressed by Equations (5) and (6), respectively:

$$
\beta \mathfrak{z} = \frac{V - W}{W} = \frac{\omega\_6 - \omega\_2}{\omega\_4 - \omega\_2} = -1 \tag{5}
$$

$$
\beta'\_{~3} = \frac{P}{M - P} = \frac{\omega\_6}{\omega\_4} = k
\tag{6}
$$

where the virtual split power ratio β<sup>3</sup> of component 3 of the differential planetary mechanism of the loader drive axle is defined as the ratio between the virtual power passing to component 6 and the virtual power passing to component 4 (Figure 3b). *V* − *W* and *W* are the virtual power value flowing to member 6 and 4 and *P* and *M* − *P* are the power value flowing to member 6 and 4. ω*<sup>i</sup>* is the rotational angular velocity of component *i.* Specifically pointed out, the virtual power was measured by supposing that the observer stood on the carrier of a planetary gear train rotating at an angular velocity ω*i.* The speed ratio *k* was defined as the ratio of the rotational angular velocity of component 6 (ω6) to component 4 (ω4). It can be found from Equations (5) and (6) that *V* = 0, *P* = *<sup>k</sup>* <sup>1</sup>+*kM*.

*Energies* **2019**, *12*, 4540

The virtual power ratio α<sup>4</sup> of component 4 of the differential planetary mechanism of the loader drive axle is defined as the ratio between the virtual power and the power passing through component 4 [9] (Figure 3a,b), and established the relationship between relative angular velocity and angular velocity, as shown in Equation (7).

$$\alpha\_4 = \frac{\mathcal{W}}{M - P} = \frac{\omega\_4 - \omega\_2}{\omega\_4} = \text{sgn}\alpha\_4 \frac{\omega\_3 \mathbf{z}\_3}{\omega\_4 \mathbf{z}\_4} = \text{sgn}\alpha\_4 h i \tag{7}$$

where sgnω<sup>4</sup> = ⎧ ⎪⎪⎪⎨ ⎪⎪⎪⎩ 1 Turn right 0 Straight −1 Turn left , *h* = <sup>ω</sup><sup>3</sup> ω4 , *i* = *<sup>z</sup>*<sup>3</sup> *z*4 It can be obtained from Equation (7) that

$$\mathcal{W} = \text{sgn}\omega\_4 \text{hi}(M - P) = \text{sgn}\omega\_4 \frac{\text{hi}}{1 + k} \mathcal{M} \tag{8}$$

**Figure 3.** Power flows of differential planetary mechanism without considering the gear mesh power loss.

The diagram of the differential planetary mechanism of the loader drive axle shown in Figure 2 can be converted into the diagram shown in Figure 4 when considering the gear mesh power loss. The gear mesh power loss is proportional to the input power, and the proportional ratio is called the mesh power loss coefficient, so *L*1, *L*2, and *L*<sup>3</sup> indicate the mesh power loss passing through the gear pairs *G*1, *G*2, and *G*3, respectively, which are expressed by Equations (9)–(11), respectively:

$$L\_1 = \lambda\_1 M \tag{9}$$

$$L\varphi = \lambda \varphi \mathcal{W} \tag{10}$$

$$L\_3 = \lambda\_3 P \tag{11}$$

where λ*<sup>i</sup>* represents the mesh power loss coefficient.

**Figure 4.** Power flows of differential planetary mechanism considering the gear mesh power loss.

Considering the gear mesh power loss, the split power ratio β <sup>3</sup> of component 3 and the virtual power ratio α<sup>4</sup> of component 4 of the differential planetary mechanism of the loader drive axle are expressed by Equations (12) and (13), respectively (Figure 4a,b):

$$
\beta'\_{\,3} = \frac{P}{M - P - L\_1} = k \tag{12}
$$

$$\alpha\_4 = \frac{W - L\_2}{M - P - L\_1 - L\_2} = \text{sgn}\omega\_4 h i \tag{13}$$

It can be obtained from Equations (12) and (13) that

$$P = \frac{k}{1+k}(1-\lambda\_1)M\tag{14}$$

$$\mathcal{W} = \left\{ \text{sgn}\omega\_4 \frac{hi(1-\lambda\_1)\mathcal{M}}{1+k} / \left[1 + (\text{sgn}\omega\_4hi - 1)\lambda\_2\right] \right\} \tag{15}$$

From Equations (9)–(11), (14), and (15),

$$\frac{\sum L}{M} = \frac{L\_1 + L\_2 + L\_3}{M} = \lambda\_1 + \lambda\_2 \left\{ \text{sgn}\omega\_4 \frac{\text{hi}(1 - \lambda\_1)}{1 + k} / \left[ 1 + (\text{sgn}\omega\_4 \text{hi} - 1)\lambda\_2 \right] \right\} + \lambda\_3 \frac{k(1 - \lambda\_1)}{1 + k} \tag{16}$$

Accordingly, the transmission efficiency η<sup>2</sup> of the differential planetary mechanism of the loader drive axle is expressed as Equation (17):

$$\eta\_2 = 1 - \frac{\sum L}{M} = 1 - \left\{ \lambda\_1 + \lambda\_2 \left\{ \text{sgn}\omega\_4 \frac{\text{hi}(1-\lambda\_1)}{1+k} / [1 + (\text{sgn}\omega\_4 \text{hi} - 1)\lambda\_2] \right\} + \lambda\_3 \frac{k(1-\lambda\_1)}{1+k} \right\} \tag{17}$$

#### 2.2.3. Research on Transmission Efficiency of the Wheel Planetary Reducing Gear

The diagram of the wheel planetary reducing gear of the loader drive axle shown in Figure 2 can be converted into the diagram shown in Figure 5 without considering the gear mesh power loss. The virtual power ratio α<sup>4</sup> of component 4 of the wheel planetary reducing gear of the loader drive axle is defined as the ratio between the virtual power and the power passing through component 4 [9], and expressed by Equation (18) (Figure 5a,b):

$$
\alpha\_4 = \frac{M - V}{M} = \frac{\omega\_4 - \omega\_8}{\omega\_4} = \frac{p}{1 + p} \tag{18}
$$

It can be obtained from Equation (18) that

$$V = \frac{1}{1+p}M\tag{19}$$

**Figure 5.** Power flows of wheel planetary reducing gear without considering the gear mesh power loss.

The diagram of wheel planetary reducing gear of the loader drive axle shown in Figure 2 can be converted into the diagram shown in Figure 6 when considering the gear mesh power loss. *L*<sup>6</sup> and *L*<sup>7</sup> indicate the mesh power loss passing through the gear pairs *G*<sup>6</sup> and *G*7, respectively, which are expressed by Equations (20) and (21), respectively:

$$L\theta = \lambda \theta M\tag{20}$$

$$L\_7 = \lambda\_7 (M - V - L\_6) \tag{21}$$

**Figure 6.** Power flows of wheel planetary reducing gear considering the gear mesh power loss.

Considering the gear mesh power loss, the virtual power ratio α<sup>4</sup> of component 4 of the wheel planetary reducing gear of the loader drive axle is expressed by Equation (22) (Figure 6a,b):

$$
\alpha\_4 = \frac{M - V}{M} = \frac{\omega\_4 - \omega\_8}{\omega\_4} = \frac{p}{1 + p} \tag{22}
$$

From Equation (22),

$$V = \frac{1}{1+p}M\tag{23}$$

Accordingly, the transmission efficiency η<sup>3</sup> of the wheel planetary reducing gear of the loader drive axle is expressed by Equation (24):

$$\eta\_3 = 1 - \frac{\sum L}{M} = 1 - (\lambda\_6 + \lambda\_7 \frac{p}{1+p} - \lambda\_6 \lambda\_7) \tag{24}$$

#### **3. Results and Discussions**

According to Equations (2), (17), and (24), the transmission efficiency of the main reducing gear, the differential planetary mechanism, and the wheel planetary reducing gear of the loader drive axle can be expressed by Equation (25):

$$
\eta = \eta\_1 \eta\_2 \eta\_3 = f(\mu, \phi, \psi, E, D, h, i, \lambda, k, \mathfrak{p}) \tag{25}
$$

To analyze the influence of the mesh friction coefficient μ of the main reducing gear on the transmission efficiency η of the loader drive axle, the curve of the transmission efficiency η of the loader drive axle and the mesh friction coefficient μ of the gear teeth were drawn under the condition that the other parameters shown in Equation (25) were unchanged (Figure 7). According to Figure 7, the transmission efficiency η of the loader drive axle decreased when the mesh friction coefficient μ of the main reducing gear increased under the condition that the other parameters were unchanged. Simultaneously, the mesh friction coefficient μ of the main reducing gear was 0.05, and the transmission efficiency η of the loader drive axle was 0.707.

**Figure 7.** Relationship of the transmission efficiency η of the loader drive axle and the mesh friction coefficient μ.

To determine if the transmission efficiency η of the loader drive axle is affected by the mesh power loss coefficient λ*i*, the curve of the transmission efficiency η of the loader drive axle and the mesh power loss coefficient λ*<sup>i</sup>* was drawn under the condition that the other parameters shown in Equation (25) were unchanged (Figure 8). According to Figure 8, the transmission efficiency η of the loader drive axle decreased as the mesh power loss coefficient λ*<sup>i</sup>* increased. That is, the transmission efficiency η of the loader drive axle decreased as the mesh power loss coefficients λ<sup>1</sup> and λ<sup>6</sup> increased, and it also decreased as the mesh power loss coefficients λ<sup>3</sup> and λ<sup>7</sup> increased. As the mesh power loss coefficients λ<sup>3</sup> and λ<sup>7</sup> increased, the transmission efficiency η of the loader drive axle decreased, but the degree of decrease was slower than λ<sup>1</sup> and λ6; as the mesh power loss coefficient λ<sup>2</sup> increased, it initially decreased slowly and then sharply.

**Figure 8.** Relationship of the transmission efficiency η of the loader drive axle and the mesh power loss coefficient λ*i*.

Here, the speed ratio γ was defined as the ratio of the rotational angular velocity of component 2 (ω2) to component 4 (ω4), as shown in Equation (26):

$$
\gamma = \frac{\alpha\_2}{\alpha\_4} \tag{26}
$$

To determine if the transmission efficiency η of the loader drive axle is affected by the speed ratio γ, the curve of the transmission efficiency η of the loader drive axle and the speed ratio γ were drawn under the condition that the other parameters shown in Equation (25) were unchanged (Figure 9). According to Figure 9, the transmission efficiency η of the loader drive axle increased with the increase of the speed ratio γ when the other parameters were unchanged, but the transmission efficiency η ranged from 0.7055 to 0.7090; that is, the modification of the speed ratio γ had no obvious effect on the transmission efficiency η of the loader drive axle.

**Figure 9.** Relationship of the transmission efficiency η of the loader drive axle and the speed ratio γ.

#### **4. Conclusions**


Based on the above research of this paper, the influence of rolling element bearings and temperatures of the gears on the transmission efficiency of the loader drive axle was carried out, so that the influence of multiple factors on the transmission efficiency of the loader drive axle can be studied more comprehensively.

**Author Contributions:** J.L. performed the formula derivation, data analyses, and wrote the manuscript; T.M. and Q.H. proposed the concept of the study; T.W. and W.C. contributed to writing, reviewing, and editing.

**Funding:** This research was funded by the Research Project of Guangdong Higher Education Institute (19GYB014), the Zhaoqing University Quality Engineering and Teaching Reform Project (zlgc201751 and zlgc201829), and the Guangdong Youth Innovation Talents Project (2018KQNCX290).

**Acknowledgments:** The authors wish to express their gratitude to the School of Mechanical and Automotive Engineering, Zhaoqing University, and express their gratitude for the funding provided by the Guangdong Youth Innovation Talents Project.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
