*Supporting open minds since 2005*

Mass Spectrometry - Future Perceptions and Applications http://dx.doi.org/10.5772/intechopen.78510 Edited by Ganesh Shamrao Kamble

#### Contributors

Vijayabhaskar Veeravalli, Giulio Tarro, Alessandra Rossi, Moreno Paolini, Xianquan Zhan, Tian Zhou, Oscar Núñez, Guillem Campmajó, Nerea Núñez, Najeeb Punnakayathil, Gunawan Witjaksono, Sagir Alva

#### © The Editor(s) and the Author(s) 2019

The rights of the editor(s) and the author(s) have been asserted in accordance with the Copyright, Designs and Patents Act 1988. All rights to the book as a whole are reserved by INTECHOPEN LIMITED. The book as a whole (compilation) cannot be reproduced, distributed or used for commercial or non-commercial purposes without INTECHOPEN LIMITED's written permission. Enquiries concerning the use of the book should be directed to INTECHOPEN LIMITED rights and permissions department (permissions@intechopen.com).

Violations are liable to prosecution under the governing Copyright Law.

Individual chapters of this publication are distributed under the terms of the Creative Commons Attribution 3.0 Unported License which permits commercial use, distribution and reproduction of the individual chapters, provided the original author(s) and source publication are appropriately acknowledged. If so indicated, certain images may not be included under the Creative Commons license. In such cases users will need to obtain permission from the license holder to reproduce the material. More details and guidelines concerning content reuse and adaptation can be found at http://www.intechopen.com/copyright-policy.html.

#### Notice

Statements and opinions expressed in the chapters are these of the individual contributors and not necessarily those of the editors or publisher. No responsibility is accepted for the accuracy of information contained in the published chapters. The publisher assumes no responsibility for any damage or injury to persons or property arising out of the use of any materials, instructions, methods or ideas contained in the book.

First published in London, United Kingdom, 2019 by IntechOpen IntechOpen is the global imprint of INTECHOPEN LIMITED, registered in England and Wales, registration number: 11086078, 7th floor, 10 Lower Thames Street, London, EC3R 6AF, United Kingdom Printed in Croatia

British Library Cataloguing-in-Publication Data A catalogue record for this book is available from the British Library

Additional hard and PDF copies can be obtained from<EMAIL_ADDRESS>
Mass Spectrometry - Future Perceptions and Applications Edited by Ganesh Shamrao Kamble p. cm. Print ISBN 978-1-83962-756-9 Online ISBN 978-1-83962-757-6 eBook (PDF) ISBN 978-1-83962-758-3


**1. Introduction**

Multivariate Cauchy distribution (MCD) belongs to the elliptical symmetric distributions [1] and is a special case of the multivariate *t*-distribution [2] and the multivariate stable distribution [3]. MCD has been recently used in several signal and image processing applications for which non-Gaussian models are needed. To name a few of them, in speckle denoizing, color image denoizing, watermarking, speech enhancement, among others. Sahu et al. in [4] presented a denoizing method for speckle noise removal applied to a retinal optical coherence tomography (OCT) image. The method was based on the wavelet transform where the sub-bands coefficients were modeled using a Cauchy distribution. In [5], a dual tree complex wavelet transform (DTCWT)-based despeckling algorithm was proposed for synthetic aperture radar (SAR) images, where the DTCWT coefficients in each subband were modeled with a multivariate Cauchy distribution. In [6], a new color image denoizing method in the contourlet domain was suggested for reducing noise in images corrupted by Gaussian noise where the contourlet subband coefficients were described by the heavy-tailed MCD. Sadreazami et al. in [7] put forward a novel multiplicative watermarking scheme in the contourlet domain where the watermark detector was based on the bivariate Cauchy distribution and designed to capture the across scale dependencies of the contourlet coefficients. Fontaine et al. in [8] proposed a semi-supervised multichannel speech enhancement system where both speech and noise follow the heavy-tailed multi-variate complex Cauchy distribution.

Kullback–Leibler divergence (KLD), also called relative entropy, is one of the most fundamental and important measures in information theory and statistics [9,10]. KLD was first introduced and studied by Kullback and Leibler [11] and Kullback [12] to measure the

**Citation:** Bouhlel, N.; Rousseau, D. A Generic Formula and Some Special Cases for the Kullback–Leibler Divergence between Central Multivariate Cauchy Distributions. *Entropy* **2022**, *24*, 838. https:// doi.org/10.3390/e24060838

Academic Editors: Karagrigoriou Alexandros and Makrides Andreas

Received: 22 May 2022 Accepted: 15 June 2022 Published: 17 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

divergence between two probability mass functions in the case of discrete random variables and between two univariate or multivariate probability density functions in the case of continuous random variables. In the literature, numerous entropy and divergence measures have been suggested for measuring the similarity between probability distributions, such as Rényi [13] divergence, Sharma and Mittal [14] divergence, Bhattacharyya [15,16] divergence and Hellinger divergence measures [17]. Other general divergence families have been also introduced and studied like the *φ*-divergence family of divergence measures defined simultaneously by Csiszár [18] and Ali and Silvey [19] where the KLD measure is a special case, the Bregman family divergence [20], the R-divergences introduced by Burbea and Rao [21–23], the statistical *f*-divergences [24,25] and recently the new family of a generalized divergence called the (*h*, *φ*)-divergence measures introduced and studied in Menéndez et al. [26]. Readers are referred to [10] for details about these divergence family measures.

KLD has a specific interpretation in coding theory [27] and is therefore the most popular and widely used as well. Since information theoretic divergence and KLD in particular are ubiquitous in information sciences [28,29], it is therefore important to establish closed-form expressions of such divergence [30]. An analytical expression of the KLD between two univariate Cauchy distributions was presented in [31,32]. To date, the KLD of MCDs has no known explicit form, and it is in practice either estimated using expensive Monte Carlo stochastic integration or approximated. Monte Carlo sampling can efficiently estimate the KLD provided that a large number of independent and identically distributed samples is provided. Nevertheless, Monte Carlo integration is a too slow process to be useful in many applications. The main contribution of this paper is to derive a closed-form expression for the KLD between two central MCDs in a general case to benchmark future approaches while avoiding approximation using expensive Monte Carlo (MC) estimation techniques. The paper is organized as follows. Section 2 introduces the MCD and the KLD. Section 3 gives some definitions and propositions related to a multiple power series used to compute the closed-form expression of the KLD between two central MCDs. In Sections 4 and 5, expressions of some expectations related to the KLD are developed by exploiting the propositions presented in the previous section. Section 6 demonstrates some final results on the KLD computed for the central MCD. Section 7 presents some particular results such as the KLD for the univariate and the bivariate Cauchy distribution. Section 8 presents the implementation procedure of the KLD and a comparison with Monte Carlo sampling method. A summary and some conclusions are provided in the final section.


*Article* **Study on Dynamic Characteristics of a Rotating Sandwich Porous Pre-Twist Blade with a Setting Angle Reinforced by Graphene Nanoplatelets**

**Jiapei Peng 1 , Lefa Zhao 2,\* and Tianyu Zhao 1,\***

**Abstract:** Lightweight blades with high strength are urgently needed in practical rotor engineering. Sandwich structures with porous core and reinforced surfaces are commonly applied to achieve these mechanical performances. Moreover, blades with large aspect ratios are established by the elastic plate models in theory. This paper studies the vibration of a rotating sandwich pre-twist plate with a setting angle reinforced by graphene nanoplatelets (GPLs). Its core is made of foam metal, and GPLs are added into the surface layers. Supposing that nanofillers are perfectly connected with matrix material, the effective mechanical parameters of the surface layers are calculated by the mixing law and the Halpin–Tsai model, while those of the core layers are determined by the open-cell scheme. The governing equation of the rotating plate is derived by employing the Hamilton principle. By comparing with the finite element method obtained by ANSYS, the present model and vibration analysis are verified. The material and structural parameters of the blade, including graphene nanoplatelet (GPL) weight faction, GPL distribution pattern, porosity coefficient, porosity distribution pattern, length-to-thickness ratio, length-to-width ratio, setting angle and pre-twist angle of the plate are discussed in detail. The finds provide important inspiration in the designing of a rotating sandwich blade.

**Keywords:** sandwich plate; graphene nanoplates; vibration; rotating; porosity

**MSC:** 74K20

#### **1. Introduction**

Blade structures are commonly applied in several rotating machineries, such as aero engines and wind turbines. As the mechanical performance of the rotor structure is required to improve, the blade structure is designed to be lighter and thinner. Blades are extremely susceptible to vibration damage, which accounts for a very large proportion of the failure of the blade. So, it is very meaningful to study vibration behaviors of blades. In previous studies, blades are often simplified as rotating plates or rotating beams, which have been studied comprehensively [1–10]. Based on Chebyshev polynomials, Gen et al. [11] studied the nonlinear dynamic behavior of blades with variable rotating speed. By using the Galerkin method, Avramov et al. [12] investigated the bending-torsional nonlinear vibration behavior of a rotating beam with asymmetric cross-section. By adopting the Ritz method, Mcgee [13] calculated natural frequencies of cantilever parallelepipeds that were skewed and twisted in the meantime. Yao et al. [14] studied the nonlinear dynamic behavior of high-speed rotating plates. Based on the Hamilton principle, Xu et al. [15] studied the nonlinear vibration of a rotating cantilever beam in the magnetic field. Wang et al. [16] studied the nonlinear stability of rotating blades. Hashemia et al. [17] presented the finite element formulas for vibration analysis of rotating thick plates. Li et al. [18] studied the

**Citation:** Peng, J.; Zhao, L.; Zhao, T. Study on Dynamic Characteristics of a Rotating Sandwich Porous Pre-Twist Blade with a Setting Angle Reinforced by Graphene Nanoplatelets. *Mathematics* **2022**, *10*, 2814. https://doi.org/10.3390/ math10152814

Academic Editors: Xiang Li, Shuo Zhang and Wei Zhang

Received: 26 April 2022 Accepted: 29 July 2022 Published: 8 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

free vibration of variable rotating cantilever rectangular functionally graded (FG) plates moving in a wide range. By using the Donner shell theory, Shakour [19] studied vibrations of a rotating FG conical shell in different thermal environment. Based on the Hamilton principle, Qin et al. [20] investigated the bent-bending coupled vibration of rotating composite thin-walled beams under aerodynamic and humid conditions. Yutaek et al. [21] proposed a dynamic model of rotating pre-torsional blades with variable cross-section. Yang et al. [22] promoted the FG plate theory proposed by Mian and Spencer from two aspects. Arumugam et al. [23] studied vibration behaviors of a planar layered laminate with rotational effects. Tuzzi et al. [24] conducted a study of coupling vibrations between shaft-bending and disk-zero nodal diameter mode vibration in a flexible shaft–disk assembly. Based on classical plate theory and Hamilton principle, Sun et al. [25] investigated the vibration characteristics of rotating blades with a staggered angle.

Due to the influence of aerodynamic loads on the blade during rotation, higher requirements are put forward for the structural strength of the blade. The traditional blade material has not met the needs. GPLs have been very popular as a reinforced material in the past few years and widely used in various engineering fields, such as physics, electrical engineering, materials science and nanoengineering applications [26–30]. Its structural stiffness and lightweight meet the requirements, which has become the focus of scholars' research. Zhao et al. [31] studied the free vibration of GPL reinforced rotating FG pre-twisted blade-shaft assembly. Li et al. [32] carried out nonlinear vibration and dynamic buckling analysis of laminated FG porous plates reinforced by GPLs. Yang et al. [33] studied the buckling and free vibration behavior of FG porous nanocomposite plates reinforced with GPLs. Based on the Timoshenko beam theory, Chen et al. [34] analyzed the elastic buckling and static bending of shear deformable FG porous beams. Wu et al. [35] studied the parameter instability of FG-GPL reinforced nanocomposite plate under the action of periodic uniaxial plane internal force in uniform temperature field. Lei et al. [36] used element-free -Ritz method to conduct free vibration analysis on single-walled carbon nanotube reinforced FG nanocomposite plate. Yang et al. [37] studied the buckling and post-buckling behavior of FG multilayer nanocomposite beams reinforced by low content GPLs on elastic foundation. Wang et al. [38] studied the eigenvalue buckling of multilayer FG cylindrical shells reinforced with GPLs by using finite element method.

To sum up, GPL reinforced porous rotor structures are suitable for high-speed rotating blades to achieve high strength or light weight. Moreover, the sandwich blade structure with GPL reinforced surfaces and porous core could be designed to avoid air turbulence. This paper investigated the theoretical modeling and free vibration of a pre-twist sandwich blade. To the authors0 knowledge, almost no study focused on the vibration analysis of a rotating sandwich blade with GPL reinforced porous core. Based on the Kirchhoff plate theory, the rotating pre-twist blade with a setting angle is established in this paper. The equations of motion are derived by applying the Hamilton principle and then solved by the using the hypothetical modal method. The effects of material and structural parameters of the sandwich blade on its free vibration are examined in detail. The findings shed a bright light on the design of sandwich blades with GPL reinforced porous core to achieve better mechanical performance.

#### **2. Theoretical Model**

Figure 1 shows a rotating pre-twist blade with a setting angle, which attached to a rigid body. To describe the motion and displacements of the blade, three coordinate systems are carried out. The *O-XYZ* coordinate system is established on the rigid body, while the reference coordinate system *o-xyz* is fixed on the blade root. Moreover, another reference coordinate system *o* 0 -*uvw* is fixed on the blade at a distance *x* away from the point o. The radius of the rigid body is *Rd*, and it rotates around the *Y* axis at a constant speed Ω. The length, width and thickness of the blade are *a*, *b* and *h*, respectively. Moreover, *h*c and *h*f are the thicknesses of the porous core and the GPL reinforced face layer, respectively. It can be known that *h* = *hc* + 2 *hf* . Then, define *m* = *h*c/*h*f , which is the ratio of the core to the surface

layer. The setting angle and pre-twist angle of the blade are *θ* and *ϕ*, respectively. The twist angle at an arbitrary point on the blade is *β*(*x*) = *θ* + *kx*, where *k* = (*ϕ* − *θ*)/*a*. The twist angle at an arbitrary point on the blade is β() + *x* =θ *kx* , where *k a* = − ( )/ ϕ θ.

*h*f are the thicknesses of the porous core and the GPL reinforced face layer, respectively. It can be known that *h* = *hc* + 2 *hf*. Then, define *m* = *h*c/*h*f, which is the ratio of the core to the surface layer. The setting angle and pre-twist angle of the blade are *θ* and *φ*, respectively.

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 3 of 18

**Figure 1.** A rotating pre-twist blade with a setting angle attached to a rigid body. **Figure 1.** A rotating pre-twist blade with a setting angle attached to a rigid body.

#### **3. Material Properties 3. Material Properties**

O*P*.

given by

shown as

where η*T* and η

and *T* ξ

 and *L* ξ

As shown in Figure 2, three different porosity distribution patterns are considered in this paper. Plus, it can be seen clearly that the material parameters vary along the blade thickness direction. Pattern X means that more pores are distributed around the edges; Pattern O tells more pores are distributed in the middle surface; Pattern U gives the uniform porosity distribution. Based on the open-cell scheme [39], the material properties are As shown in Figure 2, three different porosity distribution patterns are considered in this paper. Plus, it can be seen clearly that the material parameters vary along the blade thickness direction. Pattern X means that more pores are distributed around the edges; Pattern O tells more pores are distributed in the middle surface; Pattern U gives the uniform porosity distribution. Based on the open-cell scheme [39], the material properties are

$$\begin{aligned} \text{Pattern} \, \mathsf{X}: \begin{cases} \begin{aligned} &E(Z) = E\_1[1 - e\_0 \cos(\pi Z/h)] \\ \nu(Z) = \nu\_1[1 - e\_0 \cos(\pi Z/h)] \\ \rho(Z) = \rho\_1[1 - e\_m \cos(\pi Z/h)] \end{aligned} \end{aligned} \\ \text{Pattern} \, \mathsf{U}: \begin{cases} &E(Z) = E\_1\alpha \\ \nu(Z) = \nu\_1\alpha \\ \rho(Z) = \rho\_1\alpha' \end{cases} \end{aligned} \tag{1}$$

$$\begin{aligned} \text{Pattern} \, \mathsf{O}: \begin{cases} &E(Z) = E\_1[1 - e\_0^\*(1 - \cos(\pi Z/h))] \\ \nu(Z) = \nu\_1[1 - e\_0^\*(1 - \cos(\pi Z/h))] \\ \rho(Z) = \rho\_1[1 - e\_m^\*(1 - \cos(\pi Z/h))] \end{cases} \end{aligned} \tag{2}$$

Pattern O: π νν π ρρ π ∗ ∗ ∗ = −− = −− = −− ( ) ( ( ) ( ( 1 0 1 0 1 1 1 cos( / )) 1 1 cos( / )) ( ) 1 1 cos( / )) *m EZ E e Z h Z e Zh Z e Zh* where *E*(*Z*), *ρ*(*Z*), *ν*(*Z*) are the effective elasticity modulus, mass density, Poisson's ratio of the core layer, respectively.*E*1, *ρ*1, *ν*1 are the corresponding material parameters of core layer without pores. *e*0, *α*,*e* ∗ 0 are the porosity coefficients of the three distribution patterns. *em*, *α* 0 ,*e* ∗ *m* are the mass density coefficients of the three distribution patterns. *Mathematics* **2022**, *10*, x FOR PEER REVIEW 4 of 18

**Figure 2.** Porosity distribution patterns of the core layer: (**a**) Pattern X*P*; (**b**) Pattern U*P*; (**c**) Pattern **Figure 2.** Porosity distribution patterns of the core layer: (**a**) Pattern X*P*; (**b**) Pattern U*P*; (**c**) Pattern O*P*.

The relationship between porosity coefficients and mass density coefficients can be

0

− − = − −

1 1 cos( / ) 1 1 cos( / )

α

0.1 0.9361 0.1738 0.2 0.8716 0.3442 0.3 0.8064 0.5103 0.4 0.7404 0.6708

GPLs are added to the surface layers. Due to high precision, the Halpin–Tsai model [40,41] is used to calculate the surface layers' effective elasticity modulus. The formula is

> + + = − − T 3 5 1 1 +

> > ξ

*L*

− + − +

*E E E E E E E E*

*GPL M GPL M GPL M GL M P T*

ξ

81 81

*V V E E*

ξη

η

η

*L*

=

=

η

*T*

ξ

=2 =2 *L GPL GPL T GPL GPL*

where *EM* and *EGPL* are the elasticity modulus of the matrix and GPLs. *LGPL* , *wGPL* and

*L t*

Based on the mixing rule, the mass density and Poisson's ratio are

ξ

*e Zh e Zh*

 π

0

 π

\*

 ξη

*L L GPL T GPL*

*L GPL T GPL*

 η

*V V* (3)

*w t* (5)

0 e

(2)

(4)

 α

*m*

*m*

*S*

*L* are given by

can be expressed as

*GPL t* are the length, width and thickness of GPLs.

*M*

α

 ′ =

**Table 1.** Relationship of porosity coefficients.

e0 π

− =−

π

\* \*

1 cos( / ) 1 cos( / )

*e Zh e Zh*

The relationship between porosity coefficients and mass density coefficients can be given by

$$\begin{cases} 1 - e\_{\mathfrak{m}} \cos(\pi \mathbf{Z}/h) = \sqrt{1 - e\_0 \cos(\pi \mathbf{Z}/h)} \\ \mathbf{a}' = \sqrt{\mathbf{a}} \\ 1 - e\_{\mathfrak{m}}^\* [1 - \cos(\pi \mathbf{Z}/h)] = \sqrt{1 - e\_0^\* [1 - \cos(\pi \mathbf{Z}/h)]} \end{cases} \tag{2}$$

Table 1 shows the relationship of porosity coefficients of the three patterns.

**Table 1.** Relationship of porosity coefficients.

GPLs are added to the surface layers. Due to high precision, the Halpin–Tsai model [40,41] is used to calculate the surface layers' effective elasticity modulus. The formula is shown as

$$E\_S = E\_M \left[ \frac{3}{8} \left( \frac{1 + \tilde{\xi}\_L \eta\_L V\_{GPL}}{1 - \eta\_L V\_{GPL}} \right) + \frac{5}{8} \left( \frac{1 + \tilde{\xi}\_T \eta\_T V\_{GPL}}{1 - \eta\_T V\_{GPL}} \right) \right] \tag{3}$$

where *ηT* and *ηL* are given by

$$\begin{cases} \begin{array}{c} \eta\_L = \frac{E\_{GPL} - E\_M}{E\_{GPL} + \xi\_L E\_M} \\ \eta\_T = \frac{E\_{GPL} - E\_M}{E\_{GPL} + \xi\_T E\_M} \end{array} \tag{4} \end{cases} \tag{4}$$

and *ξT* and *ξL* can be expressed as

$$\begin{cases} \ ^\sharp \xi\_L = 2L\_{\rm GPL}/t\_{\rm GPL} \\ \ ^\sharp \xi\_T = 2w\_{\rm GPL}/t\_{\rm GPL} \end{cases} \tag{5}$$

where *EM* and *EGPL* are the elasticity modulus of the matrix and GPLs. *LGPL*, *wGPL* and *tGPL* are the length, width and thickness of GPLs.

Based on the mixing rule, the mass density and Poisson's ratio are

$$\begin{cases} \rho\_{\rm s} = V\_{\rm GPL} \rho\_{\rm GPL} + (1 - V\_{\rm GPL}) \rho\_M\\ \upsilon\_{\rm s} = V\_{\rm GPL} \upsilon\_{\rm GPL} + (1 - V\_{\rm GPL}) \upsilon\_M \end{cases} \tag{6}$$

where *υGPL* and *υm* are the Poisson's ratio of matrix and GPLs, respectively. *ρGPL* and *ρm* are the mass density of matrix material GPLs, respectively. *VGPL* is the volume fraction of GPLs in the plate, expressed as

$$V\_{\rm GPL} = \frac{W\_{\rm GPL}}{W\_{\rm GPL} + \rho\_{\rm GPL}(1 - W\_{\rm GPL})/\rho\_M} \tag{7}$$

where, *WGPL* is the weight fraction of GPLs, obtained as

$$\mathcal{W}\_{\text{GPL}} = \begin{cases} \frac{4Z^2}{h\_f^2} \lambda\_1 \mathcal{W}\_0 & \text{Pattern} \mathcal{X}\_p\\ \lambda\_2 \mathcal{W}\_0 & \text{Pattern} \mathcal{U}\_p\\ \left(1 - \frac{4Z^2}{h\_f^2}\right) \lambda\_3 \mathcal{W}\_0 & \text{Pattern} \mathcal{O}\_p \end{cases} \tag{8}$$

where *W*0 is the GPL characteristic value; *λ*1, *λ*2, *λ*3 are the GPL weight fraction indices determined by the GPL weight distribution patterns, as listed in Table 2.

**Table 2.** Values of GPL weight fraction indices. **Table 2.** Values of GPL weight fraction indices. λ

where

ρ

υ*GPL* and

υ

tion of GPLs in the plate, expressed as

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 5 of 18

ρρ

 = +− = +−

s (1 )

*V V*

*m* are the mass density of matrix material GPLs, respectively. *VGPL* is the volume frac-

ρ

λ

*Z W h*

2 2 1 0

 −

2 2 3 0

*f*

*h*

2 0

*f*

λ

determined by the GPL weight distribution patterns, as listed in Table 2.

 =

 ( ) = + −1 *GPL*

λ

*Z W*

4 1 Patte

λλλ

*W*

*GPL GPL GPL M*

X 4 Pattern

*s GPL GPL GPL M*

*m* are the Poisson's ratio of matrix and GPLs, respectively.

υυ

*GPL*

*V*

*L*

*W W*

*GP*

where *W*0 is the GPL characteristic value; 123

where, *WGPL* is the weight fraction of GPLs, obtained as

 ρ

> ρ

Pattern

*W W* (7)

U

*p*

*p*

, , are the GPL weight fraction indices

*p*

rnO

*V V* (6)

ρ*GPL* and

(8)

 υ

(1 ) *GPL GPL GPL M*

The corresponding GPL distribution patterns are as shown in Figure 3. The corresponding GPL distribution patterns are as shown in Figure 3.

**Figure 3.** GPL distribution patterns in the surface layer of the plate: (**a**) Pattern XG; (**b**) Pattern UG; (**c**) Pattern OG. **Figure 3.** GPL distribution patterns in the surface layer of the plate: (**a**) Pattern XG; (**b**) Pattern UG; (**c**) Pattern OG.

#### **4. Equations of Motion 4. Equations of Motion**

To describe the motion of the rotating plate, three coordinate systems are established. The coordinates of an arbitrary point on the plate in the local coordinate system are =+ +++ ( ) () *R x R i yj z w k d* . To describe the motion of the rotating plate, three coordinate systems are established. The coordinates of an arbitrary point on the plate in the local coordinate system are *R* = (*x* + *Rd* )*i* + *yj* + (*z* + *w*)*k*.

Then, the coordinate of the same point on the plate in the fixed coordinate system is Then, the coordinate of the same point on the plate in the fixed coordinate system is

$$R = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos\beta(\mathbf{x}) & -\sin\beta(\mathbf{x}) \\ 0 & \sin\beta(\mathbf{x}) & \cos\beta(\mathbf{x}) \end{pmatrix} \begin{pmatrix} \mathbf{x} + R\_d \\ y \\ z + w \end{pmatrix} \tag{9}$$

The speed of the plate is *V* = . *R* + (Ω*i* × *R*), where Ω is the rotating speed. Then, the kinetic energy of the rotating blade can be obtained as

$$T\_P = \frac{1}{2}\rho h \int\_0^b \int\_0^a \left[ w^2 \Omega^2 + \left(\frac{\partial w}{\partial t}\right)^2 + (R+x)^2 - 2w(\cos\beta(x))\frac{\partial w}{\partial t}(R+x)\Omega + y^2 \right] dx dy \tag{10}$$

Considering the gyroscopic effect, the rotating strain potential energy is

$$\mathcal{U}\_{1} = \frac{\mathbb{E}h^{3}}{24(1-\nu^{2})} \int\_{0}^{b} \int\_{0}^{a} \left\{ \frac{1}{\mathcal{B}^{4}} \left( \frac{\partial^{2}w}{\partial x^{2}} + \frac{\partial^{2}w}{\partial y^{2}} \right)^{2} - \frac{2(1-\nu)}{\mathcal{B}^{2}} \left[ \frac{\partial^{2}w}{\partial x^{2}} \frac{\partial^{2}w}{\partial y^{2}} - \left( \frac{\partial^{2}w}{\partial x \partial y} \right)^{2} \right] \right\} d\mathbf{x} dy \tag{11}$$

where *fc*1 and *fc*2 are the centrifugal forces in two directions of the blade, expressed as

$$f\_{c1} = \rho h \Omega^2 \int\_{\chi}^{a} (R + x) dx \tag{12}$$

$$f\_{c2} = \rho h \Omega^2 \int\_y^{\frac{b}{2}} y \left(\sin\left(\beta(x)\right)^2\right) dy\tag{13}$$

The total potential energy can be given by

$$\begin{split} \mathcal{U}\_{P} &= \frac{Eh^{3}}{24(1-v^{2})} \int\_{0}^{b} \int\_{0}^{a} \left\{ \frac{1}{\mathcal{S}^{4}} \left( \frac{\partial^{2}w}{\partial x^{2}} + \frac{\partial^{2}w}{\partial y^{2}} \right)^{2} - \frac{2(1-v)}{\mathcal{B}^{2}} \left[ \frac{\partial^{2}w}{\partial x^{2}} \frac{\partial^{2}w}{\partial y^{2}} - \left( \frac{\partial^{2}w}{\partial x \partial y} \right)^{2} \right] \right\} \mathrm{d}x \mathrm{d}y \\ &+ \frac{1}{2} \rho h \Omega^{2} \int\_{0}^{b} \int\_{0}^{a} \left\{ \left[ R(a-x) + \frac{1}{2} (a^{2} - x^{2}) \right] \left( \frac{\partial w}{\partial x} \right)^{2} + (b^{2} - y^{2}) (\sin(\beta(x)))^{2} \right\} \mathrm{d}x \mathrm{d}y \end{split} \tag{14}$$

Applying the Hamiltonian variational principle *δ* R *t*1 *t*0 (*TP* − *UP*)d*t* = 0, the equations of motion can be derived as

$$\begin{aligned} \int\_{0}^{b} \int\_{0}^{a} \int\_{0}^{b} & \left[ \begin{aligned} &\rho h \Big(\omega^{2} + \Omega^{2}\Big) \mathcal{W} \delta \mathcal{W} \\ & - \frac{1}{2} \rho h \mathcal{Q}^{2} \Bigg] \left[ \begin{aligned} &\left[ R(a-x) + \frac{1}{2} (a^{2} - x^{2}) \right] \frac{\partial \mathcal{W}}{\partial x} \delta \left( \frac{\partial \mathcal{W}}{\partial x} \right) \\ & + (\frac{b^{2}}{4} - y^{2}) (\sin \beta(x))^{2} \frac{\partial \mathcal{W}}{\partial y} \delta \left( \frac{\partial \mathcal{W}}{\partial y} \right) \end{aligned} \right] \\ & - \frac{D}{B^{2}} \Bigg[ & + \frac{\partial^{2} \mathcal{W}}{\partial x^{2}} \delta \left( \frac{\partial^{2} \mathcal{W}}{\partial x^{2}} \right) + \frac{\partial^{2} \mathcal{W}}{\partial y^{2}} \delta \left( \frac{\partial^{2} \mathcal{W}}{\partial x^{2}} \right) \\ & + \frac{D(1-v)}{B^{2}} \Bigg[ & + \frac{D(1-v)}{B^{2}} \Bigg[ \begin{aligned} &\frac{\partial^{2} \mathcal{W}}{\partial x^{2}} \delta \left( \frac{\partial^{2} \mathcal{W}}{\partial x^{2}} \right) + \frac{\partial^{2} \mathcal{W}}{\partial y^{2}} \delta \left( \frac{\partial^{2} \mathcal{W}}{\partial x^{2}} \right) \\ & - 2 \frac{\partial^{2} \mathcal{W}}{\partial x \partial y} \delta \left( \frac{\partial^{2} \mathcal{W}}{\partial x \partial y} \right) \end{aligned} \end{aligned} \right] \tag{15}$$

According to the approximation method of combined series of beam function, the vibration mode expression, satisfying all displacement boundary conditions, is

$$\mathcal{W}(\mathbf{x}, \mathbf{y}) = \sum\_{m=1}^{M} \sum\_{n=1}^{N} A\_{mn} \phi\_m(\mathbf{x}) \varphi\_n(\mathbf{y}) \tag{16}$$

The fixed-free beam function *φm*(*x*) and the free-free beam function *ϕn*(*y*) are given by

$$\begin{cases} \begin{aligned} \phi\_m(\mathbf{x}) &= \cosh(a\_m \mathbf{x}) - \cos(a\_m \mathbf{x}) - c\_m[\sinh(a\_m \mathbf{x}) - \sin(a\_m \mathbf{x})] \\ \phi\_n(y) &= \cosh(\beta\_n y) + \cos(\beta\_n y) - d\_n[\sinh(\beta\_n y) + \sin(\beta\_n y)] \\ c\_m &= \frac{\cos(a\_m a) + \cosh(a\_m a)}{\sin(a\_m a) + \sinh(a\_m a)}, d\_n = \frac{\cos(\beta\_n b) - \cosh(\beta\_n b)}{\sin(\beta\_n b) - \sinh(\beta\_n b)} \\ \cosh(a\_m a)\cos(a\_m a) + 1 &= 0, \cosh(\beta\_n b)\cos(\beta\_n b) - 1 = 0 \end{aligned} \tag{17}$$

By substituting the Equation (16) into Equation (15) and making the variational coefficient *δAmn* zero, the frequency equation can be obtained as

$$\begin{bmatrix} \rho h(\omega^{2} + \Omega^{2}) \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n}(x) \rho\_{n}(y) \rho\_{l}(x) \rho\_{j}(y) \\\\ - \frac{D}{\mathcal{F}} \left[ \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n}(x) \rho\_{n}(y) \rho\_{n}^{\*}(x) \rho\_{j}(y) + \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n}(x) \rho\_{n}^{\*}(y) \rho\_{l}(x) \rho\_{j}^{\*}(y) \right.\\\\ \left. \left. + \frac{M}{\mathcal{F}} \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n} \sigma^{\*}(x) \rho\_{n}(y) \rho\_{l}(x) \rho\_{j}^{\*}(y) + \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n}(x) \rho\_{n}^{\*}(y) \rho\_{l}^{\*}(x) \rho\_{j}(y) \right] \right.\\\\ \left. + \frac{D(1-\varepsilon)}{\mathcal{F}} \left[ \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n} \sigma^{\*}(x) \rho\_{n}(y) \rho\_{l}^{\*}(y) \right.\\\\ \left. + \frac{M}{\mathcal{F}} \sum\_{n=1}^{M} A\_{nm} \mathfrak{p}\_{n} (x) \rho\_{n}^{\*}(y) \rho\_{l}^{\*}(x) \rho\_{l}(y) \right.\\\\ \left. - 2 \sum\_{n=1}^{M} \sum\_{n=1}^{N} A\_{nm} \mathfrak{p}\_{n} (x) \rho\_{n}^{\*$$

where, *ω* is the natural frequency of the plate.

#### **5. Results and Discussions**

Before discussing the parameters, we need to conduct convergence analysis to verify correctness of theoretical solution. The structural and material parameters are shown as Table 3.

**Table 3.** Structural and material parameters.

Table 4 shows variations of the first four frequencies with different mode numbers. It is seen that the mode number tends to be convergent at *M* = *N* = 15, which will be adopted in the following discussion. *Mathematics* **2022**, *10*, x FOR PEER REVIEW 8 of 18

Table 4 shows variations of the first four frequencies with different mode numbers.

**Table 4.** Frequencies with different mode numbers. It is seen that the mode number tends to be convergent at *M* = *N* = 15, which will be

The comparison between theoretical results and finite element results obtained from ANSYS of the first four-order frequencies and corresponding vibration modes is depicted in Figure 3 and Table 4, the Solid 186 element in ANSYS is applied. Plus, the numbers of elements and nodes are 11,616 and 23,735, respectively. The convergence analysis for the element number is conducted as given in Figure 4 and Table 5. It can be clearly seen from Figure 5 and Table 6 that the error is less than 2%, and the vibration modes are in good agreement, which indicates that the present model and analysis are correct. Figure 5 and Table 6 that the error is less than 2%, and the vibration modes are in good agreement, which indicates that the present model and analysis are correct. The idea of controlling variable method is adopted for the argument discussion in this article. The parameters are demonstrated in the Table 7. Moreover, the GPL distribution and pore distribution are Pattern X, and the porosity coefficient *e*0 = 0.1, if not specified.

**Figure 4.** Finite element grid diagram. (**a**) *Ne* = 2868; (**b**) *Ne* = 4092; (**c**) *Ne* = 6523; (**d**) *Ne* = 11,616. **Figure 4.** Finite element grid diagram. (**a**) *Ne* = 2868; (**b**) *Ne* = 4092; (**c**) *Ne* = 6523; (**d**) *Ne* = 11,616.

**Frequency (Hz)** *M* **=** *N* **= 5** *M* **=** *N* **= 10** *M* **=** *N* **= 15**

Second 725.4 724.2 723.6 Third 2709.7 2707.3 2704.7 Fourth 4036.4 4033.2 4025.1

**Frequency (Hz)** *N***e = 2868** *Ne* **= 4092** *Ne =* **6523** *Ne* **= 11,616**  First 116.04 116.07 116.06 116.06 Second 720.91 720.61 720.39 720.29 Third 2678.4 2674.3 2670.8 2668.8

**Table 5.** Frequencies with different finite element numbers.

**Table 4.** Frequencies with different mode numbers.

**Table 5.** Frequencies with different finite element numbers. **Table 6.** Comparison between theoretical and finite element results. **Frequency (Hz) Theoretical Result ANSYS Error**

Fourth 4001.0 3994.2 3989.3 3986.7

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 9 of 18

**Figure 5.** Comparison of the first-four vibration modes between theoretical and finite element results: (**a**,**b**) are the first vibration mode; (**c**,**d**) are the second vibration mode; (**e**,**f**) are the third vibration mode; (**g**,**h**) are the fourth vibration mode. **Figure 5.** Comparison of the first-four vibration modes between theoretical and finite element results: (**a**,**b**) are the first vibration mode; (**c**,**d**) are the second vibration mode; (**e**,**f**) are the third vibration mode; (**g**,**h**) are the fourth vibration mode.

**Table 6.** Comparison between theoretical and finite element results.

The idea of controlling variable method is adopted for the argument discussion in this article. The parameters are demonstrated in the Table 7. Moreover, the GPL distribution and pore distribution are Pattern X, and the porosity coefficient *e*0 = 0.1, if not specified.

**Table 7.** Values of structural and material parameters. *Mathematics* **2022**, *10*, x FOR PEER REVIEW 10 of 18

Figure 6 shows variations of the first four frequencies with the rotating plate in the case of different GPL weight fraction. The frequencies of the plate gradually increase as the rotating speed rises. Moreover, it can be clearly seen that adding more GPLs into the matrix leads to the higher frequencies. Therefore, we can add more GPLs in a small amounts to the surface layers of the plate to improve its mechanical performance. Figure 6 shows variations of the first four frequencies with the rotating plate in the case of different GPL weight fraction. The frequencies of the plate gradually increase as the rotating speed rises. Moreover, it can be clearly seen that adding more GPLs into the matrix leads to the higher frequencies. Therefore, we can add more GPLs in a small amounts to the surface layers of the plate to improve its mechanical performance.

**Figure 6.** Variations of frequencies with rotating speed under different GPL weight fractions. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

Figure 7 depicts variations of the first four frequencies with the rotating plate in the case of different porosity coefficients. Obviously, the frequencies of the plate are decreased with larger porosity coefficients, which means that the core layer has more pores. It can significantly reduce the

Variations of the first four frequencies with the rotating plate in the case of different

**Figure 6.** Variations of frequencies with rotating speed under different GPL weight fractions. (**a**)

Figure 7 depicts variations of the first four frequencies with the rotating plate in the case of different porosity coefficients. Obviously, the frequencies of the plate are decreased with larger porosity coefficients, which means that the core layer has more pores. It can significantly reduce the weight, and the impact on the vibration frequency of the plate is relatively small. So, it is necessary to choose the porous core to achieve weight

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 11 of 18

First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

weight, and the impact on the vibration frequency of the plate is relatively small. So, it is necessary to choose the porous core to achieve weight reduction. reduction.

**Figure 7.** Variations of frequencies with rotating speed under different porosity coefficients*.* (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 7.** Variations of frequencies with rotating speed under different porosity coefficients. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

ratios of surface-thickness and core-thickness are exhibited in Figure 8. As the ratios continue to increase, frequencies of the plate increase significantly. It tells those thicker surfaces give a hand in enhancing the structural stiffness. Variations of the first four frequencies with the rotating plate in the case of different ratios of surface-thickness and core-thickness are exhibited in Figure 8. As the ratios continue to increase, frequencies of the plate increase significantly. It tells those thicker surfaces give a hand in enhancing the structural stiffness. *Mathematics* **2022**, *10*, x FOR PEER REVIEW 12 of 18

**Figure 8.** Variations of frequencies with rotating speed under different ratios of core layer thickness and surface layer thickness. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. Figure 9 displays variations of the first four frequencies with the rotating plate in the **Figure 8.** Variations of frequencies with rotating speed under different ratios of core layer thickness and surface layer thickness. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

case of different GPL distribution patterns. The plate with GPL distribution pattern X has the highest natural frequencies. To obtain better mechanical performance, adding more GPLs around the surfaces is a good choice. Figure 9 displays variations of the first four frequencies with the rotating plate in the case of different GPL distribution patterns. The plate with GPL distribution pattern X has the highest natural

(**a**) (**b**)

frequency.

(**a**) (**b**)

(**c**) (**d**) **Figure 8.** Variations of frequencies with rotating speed under different ratios of core layer thickness and surface layer thickness. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth

Figure 9 displays variations of the first four frequencies with the rotating plate in the case of different GPL distribution patterns. The plate with GPL distribution pattern X has

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 12 of 18

frequencies. To obtain better mechanical performance, adding more GPLs around the surfaces is a good choice. the highest natural frequencies. To obtain better mechanical performance, adding more GPLs around the surfaces is a good choice.

**Figure 9.** Variations of frequencies with rotating speed under different GPL distribution patterns. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 9.** Variations of frequencies with rotating speed under different GPL distribution patterns. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 9.** Variations of frequencies with rotating speed under different GPL distribution patterns. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

case of different porosity distribution patterns. It is obvious that porosity distribution pattern X provides the greater frequencies than the other two porosity distribution patterns. One can tell that setting more pores around the edges and arranging larger pores in the middle plane are conducive to achieve better structural stiffness. Figure 10 shows variations of the first four frequencies with the rotating plate in the case of different porosity distribution patterns. It is obvious that porosity distribution pattern X provides the greater frequencies than the other two porosity distribution patterns. One can tell that setting more pores around the edges and arranging larger pores in the middle plane are conducive to achieve better structural stiffness. Figure 10 shows variations of the first four frequencies with the rotating plate in the case of different porosity distribution patterns. It is obvious that porosity distribution pattern X provides the greater frequencies than the other two porosity distribution patterns. One can tell that setting more pores around the edges and arranging larger pores in the middle plane are conducive to achieve better structural stiffness.

Figure 10 shows variations of the first four frequencies with the rotating plate in the

(**c**) (**d**)

terns. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. Variations of the first four frequencies with the rotating plate in the case of length-to-**Figure 10.** Variations of frequencies with rotating speed under different porosity distribution patterns. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 10.** Variations of frequencies with rotating speed under different porosity distribution patterns. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

thickness ratio of the plate are given in Figure 11, where the blade length is constant. Lower length-to-thickness ratio of the plate shows greater frequencies. It indicates that

Variations of the first four frequencies with the rotating plate in the case of length-tothickness ratio of the plate are given in Figure 11, where the blade length is constant. Lower length-to-thickness ratio of the plate shows greater frequencies. It indicates that mance.

mance.

Variations of the first four frequencies with the rotating plate in the case of length-to-thickness ratio of the plate are given in Figure 11, where the blade length is constant. Lower length-to-thickness ratio of the plate shows greater frequencies. It indicates that thinner blades should be designed in rotating machines for better mechanical performance. additional requirements. Figure 14 shows variations of the first four frequencies with the rotating plate in the case of different setting angles. Visibly, raising the setting angle will increase the frequencies of the plate. We can know that the blade with larger setting angle contributes to improve the structural performance. case of different length-to-width ratio of the plate, where the blade width is constant. One can see that frequencies decrease dramatically with increase of the length-to-width ratio of the plate. It can be told that the blade length should be closer to the blade width to enhance the structural performance of the plate. Figure 13 depicts variations of the first four frequencies with the rotating plate in the

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 14 of 18

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 14 of 18

enhance the structural performance of the plate.

thinner blades should be designed in rotating machines for better mechanical perfor-

Figure 12 plots variations of the first four frequencies with the rotating plate in the case of different length-to-width ratio of the plate, where the blade width is constant. One can see that frequencies decrease dramatically with increase of the length-to-width ratio of the plate. It can be told that the blade length should be closer to the blade width to

Figure 13 depicts variations of the first four frequencies with the rotating plate in the case of different pre-twist angles. Increasing the pre-twist angle of the plate tends to reduce the frequencies. It should be noted that the blade should be flatter in the absence of

Figure 12 plots variations of the first four frequencies with the rotating plate in the

thinner blades should be designed in rotating machines for better mechanical perfor-

case of different pre-twist angles. Increasing the pre-twist angle of the plate tends to re-

**Figure 11.** Variations of frequencies with rotating speed under different length-to-thickness ratio of the plate. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 11.** Variations of frequencies with rotating speed under different length-to-thickness ratio of the plate. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

Figure 12 plots variations of the first four frequencies with the rotating plate in the case of different length-to-width ratio of the plate, where the blade width is constant. One can see that frequencies decrease dramatically with increase of the length-to-width ratio of the plate. It can be told that the blade length should be closer to the blade width to enhance the structural performance of the plate. (**c**) (**d**) **Figure 11.** Variations of frequencies with rotating speed under different length-to-thickness ratio of the plate. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

**Figure 12.** Variations of frequencies with rotating speed under different length-to-width ratio of the plate. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 12.** Variations of frequencies with rotating speed under different length-to-width ratio of the plate. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

(**a**) (**b**)

(**c**) (**d**) **Figure 13.** Variations of frequencies with rotating speed under different pre-twist angles. (**a**) First

(**a**) (**b**)

frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

Figure 13 depicts variations of the first four frequencies with the rotating plate in the case of different pre-twist angles. Increasing the pre-twist angle of the plate tends to reduce the frequencies. It should be noted that the blade should be flatter in the absence of additional requirements. (**c**) (**d**) **Figure 12.** Variations of frequencies with rotating speed under different length-to-width ratio of the plate. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 15 of 18

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 15 of 18

**Figure 13.** Variations of frequencies with rotating speed under different pre-twist angles. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 13.** Variations of frequencies with rotating speed under different pre-twist angles. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

Figure 14 shows variations of the first four frequencies with the rotating plate in the case of different setting angles. Visibly, raising the setting angle will increase the frequencies of the plate. We can know that the blade with larger setting angle contributes to improve the structural performance. (**c**) (**d**) **Figure 13.** Variations of frequencies with rotating speed under different pre-twist angles. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

**Figure 14.** Variations of frequencies with rotating speed under different setting angles. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency. **Figure 14.** Variations of frequencies with rotating speed under different setting angles. (**a**) First frequency; (**b**) second frequency; (**c**) third frequency; (**d**) fourth frequency.

#### **6. Conclusions**  This paper studies free vibration analysis of a rotating sandwich porous blade rein-**6. Conclusions**

been drawn as follows:

structure.

the manuscript.

erately thick shallow panels. *Mech. Adv. Mater. Struct.* **2022**, *29*, 594–612.

*Phys.* **2022**, *12*, 97. https://doi.org/10.1007/s13324-022-00697-3.

twist blade-disk assembly with a setting angle. *Appl. Math. Model.* **2021**, *93*, 578–596.

reinforced disk-shaft assembly with whirl motion. *Int. J. Mech. Sci.* **2021**, *197*, 106335.

eccentric mass. *Mech. Adv. Mater. Struct.* **2021**. https://doi.org/10.1080/15376494.2021.1904525.

shell-plate structure with whirl motion. *Appl. Math. Model.* **2022**, *110*, 618–632.

**References**

choice to obtain better structural stiffness.

1. Mohammad, R.P.; Amir, R.M. Hygro-thermo-elstic nonlinear analysis of functionally gaded porous composite thin and mod-

2. Zhao, T.Y.; Ma, Y.; Zhang, H.Y.; Pan, H.G.; Cai, Y. Free vibration analysis of a rotating graphene nanoplatelet reinforced pre-

3. Zhao, T.Y.; Yan, K.; Li, H.W.; Wang, X. Study on theoretical modeling and vibration performance of an assembled cylindrical

4. Zhao, T.Y.; Cui, Y.S.; Pan, H.G.; Yuan, H.Q.; Yang, J. Free vibration analysis of a functionally graded graphene nanoplatelet

5. Zhao, T.Y.; Cui, Y.S.; Wang, Y.Q.; Pan, H.G. Vibration characteristics of graphene nanoplatelet reinforced disk-shaft rotor with

6. Zhao, T.Y.; Li, K.; Ma, H. Study on dynamic characteristics of a rotating cylindrical shell with uncertain parameters. *Anal. Math.*

forced with GPLs. The effective material parameters of the sandwich plate are obtained by the mixing law and the Halpin–Tsai model. In addition, the frequencies can be calculated based on the Hamilton principle and the method of hypothesis modes. Plus, the present results are verified by the finite element method. Some useful conclusions have This paper studies free vibration analysis of a rotating sandwich porous blade reinforced with GPLs. The effective material parameters of the sandwich plate are obtained by the mixing law and the Halpin–Tsai model. In addition, the frequencies can be calculated based on the Hamilton principle

(1) Better mechanical performance can be achieved by adding more GPLs into the surface layer in a small content and making more GPLs on the edges of the plate. (2) Setting less pores in the core layer and arranging more smaller pores around edges

(4) Shorter and thinner blade should be designed to obtain better mechanical properties. (5) Appropriately increasing setting angle and reducing pre-twist angle can is a good

**Author Contributions:** Conceptualization, J.P. and T.Z.; Funding acquisition, L.Z. and T.Z.; Investigation, J.P., L.Z. and T.Z.; Methodology, J.P. and L.Z.; Supervision, T.Z.; Writing—original draft, J.P.; Writing—review & editing, T.Z. All authors have read and agreed to the published version of

**Funding**: This research was funded by the major project of aeroengine and gas turbine of China (HT-J2019-IV-0016-0084), the National Natural Science Foundation of China (No. 51805076, No. U1708255 and No. 52005086), the National Science and Technology Major Project of China (J2019-I-0008-0008) and the Fundamental Research Funds for the Central Universities of China (N2105013).

**Conflicts of Interest:** The authors declare no conflict of interest in preparing this article.

and the method of hypothesis modes. Plus, the present results are verified by the finite element method. Some useful conclusions have been drawn as follows:

**Author Contributions:** Conceptualization, J.P. and T.Z.; Funding acquisition, L.Z. and T.Z.; Investigation, J.P., L.Z. and T.Z.; Methodology, J.P. and L.Z.; Supervision, T.Z.; Writing—original draft, J.P.; Writing—review & editing, T.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Major Project of Aeroengine and Gas Turbine of China (HT-J2019-IV-0016-0084), the National Natural Science Foundation of China (No. 51805076, No. U1708255 and No. 52005086), the National Science and Technology Major Project of China (J2019-I-0008-0008) and the Fundamental Research Funds for the Central Universities of China (N2105013).

**Conflicts of Interest:** The authors declare no conflict of interest in preparing this article.

#### **References**


**5. Main methods of disperse dyeing**

Disperse dyes are the most important class of dye used in dyeing polyester fibres and provide a wide range of hues with good build-up and fastness properties adequate for most purposes. The rate of dyeing may be raised to the level of commercial acceptability, either by raising the working temperature to the region of 130C, or by dyeing at the boil in the presence of an accelerating agent or 'carrier' (Waters, 1950). Disperse dyes can also be transferred to polyester fibres under dry conditions by impregnating the goods with a suitable dispersion, drying and then baking at temperatures in the range 190-220C (Ingamells, 1993).


**3.8 RQ6: Who are the perpetrators of human trafficking, migrant smuggling, and human enslavement?**

Human trafficking, migrant smuggling, and human enslavement are perpetrated by unregistered international organized networks of individuals, mafias, and organizations that generate substantial illicit proceeds, often laundered into licit markets, with a too-low risk of prosecution and confiscation. These illegal and criminal activities are lucrative crimes that take advantage of people who are vulnerable, desperate, or simply seeking a better life.

The perpetrators include:


**1. Introduction**

The food industry is one of the major manufacturing sectors, representing 15% of the sales in the European Union and more than 1,120,000 million euros [1]. However, the rapid growth of global population and its corresponding consequences, such as the increasing demand of food production, caused an intensive use of energy resources [2]. In this sense, a considerable use of energy, estimated at 30% of the total energy consumption, is attributed to agri-food systems [1]. Given the strong dependency on fossil fuels of our current energy system, such energy consumption is responsible for 20–30% of total anthropogenic greenhouse gases (GHG) [3]. This fact results in the depletion of fossil resources and GHG emissions being of major concern. Energy is required at every stage in food production, including the cultivation and harvesting of crops, animal husbandry, transportation and distribution, and food processing for consumption, with the agricultural stage as the most critical [4]. This causes a low efficiency for food systems, with around 10 kcal of fossil fuel energy to produce 1 kcal

of food [1]. Some foods are very efficient in their use of resources to produce a nutritious meal, such as fruits and vegetables, while others are very inefficient, such as animal-based products. In general, animal foods require eight times more energy per calorie than plant-based commodities [5].

Concerns regarding inefficiency are wider when food losses and waste (FLW) are taken into account. According to the Organization for Economic Co-operation and Development (OECD) [6], more than one-third of the food produced worldwide is lost along the supply chain, involving around 38% of the energy consumed in its production. In fact, FLW was identified as a major global concern putting at risk environmental, economic, and health security [7]. This food produced is, therefore, lost or wasted from initial agricultural production down to final household consumption. Since the embodied amount of energy builds up along the chain, the later the waste occurs, the higher the energy waste and the related GHG emissions will be. In addition, due to the expected global population increase, a 60% rise in food production was forecasted for the year 2050, entailing a 50% rise in global energy consumption [8]. Therefore, the efficient use of energy is a necessary step toward reducing environmental hazards, preventing destruction of natural resources, and ensuring food sustainability [9].

With the aim of attaining sustainability, many analytical tools are applied individually or combined to determine an appropriate productive performance. Life-cycle assessment (LCA) is a powerful tool that assesses the environmental impacts of products, processes, and services; it gained in acceptance since it first appeared in the 1990s and is today well established [10]. Many authors used LCA to assess the energy efficiency of their processes. For instance, Marique and Rossi [11], Ingrao et al. [12], and Berg and Fuglseth [13] applied LCA to energy uses in the building sector, while Li and Feng [14] used it to improve the energy recovery from sewage sludge comparing different pathways. Sundaram et al. [15] assessed the energy efficiency of the production of gasoline from biogas and pyrolysis oil using an LCA approach. Regarding the agri-food sector, Carrasquer et al. [16] created a new indicator to estimate the water and energy efficiency in agro-industries based on life-cycle thinking methodologies. Skunca et al. [17] performed the LCA of the chicken meat chain to estimate the cumulative energy demand (CED), while Pires-Gaspar et al. [18] carried out an energy-LCA of Portuguese peach production using energy efficiency indicators.

On the other hand, additional methods to compute which resource/technology combinations provide the highest amounts of net energy to society are reported as energy return ratios (ERRs) [19]. The most well-known indicator is the energy return on investment (EROI). In fact, its use is widespread within the energy sector to determine the energy that is returned from an energy-collecting process as compared to the embodied energy to provide this energy [20]. However, in recent years, some authors highlighted the benefits of calculating EROI ratios to monitor energy return in food systems [21]. For instance, Laso et al. [22,23], Vázquez-Rowe et al. [20], Tyedmers [24], and Ramos et al. [25] applied the EROI ratio to fishery systems. Cancino-Espinoza et al. [26] computed the EROI of organic quinoa, and Pérez Neira et al. [27] estimated the EROI of tomato production.

In contrast to traditional EERs, the use of the data envelopment analysis (DEA) enables the identification of energy-saving targets for energy-inefficient processes [28]. DEA is a non-parametric multi-input/output linear approach for the calculation of the relative efficiency of a set of comparable decision-making units (DMUs) [29]. In recent years, the application of DEA increased in the fields of environmental and energy. For instance, DEA was applied in the building sector [30–32], the power industry [33,34], the food production sector [1,9,35], and for agricultural production [36–38]. Moreover, it was also used with economic and social variables [3]. The combination of LCA and DEA for the identification and quantification of potential environmental consequences of operational inefficiencies experienced a pronounced increase, and enriched the interpretation of results [39,40]. Several previous studies used the joint LCA/DEA method in different agricultural production systems: fisheries [41,42], grapes for winemaking [39], wheat [43], rice [44], and dairy farms [45,46].

In this contribution, a hybrid approach which combines EROI, LCA, and DEA is proposed to assess the energy and environmental efficiency of agri-food systems. The capabilities of this combined approach are illustrated through its application to the Spanish agri-food system in order to explore the operational inefficiencies and improvement actions.

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

Figure 1 depicts the methodology described in this analysis, which attempted to assess the energy and environmental efficiency of Spanish agri-food system using a life-cycle approach. Firstly, an LCA was conducted to estimate the primary energy demand (PED) and related global warming potential (GWP). Secondly, the EROI was estimated to describe the energy efficiency of food categories under study. Finally, DEA was synergistically combined with LCA results to assess the efficiency of the system.

**Figure 1.** Schematic representation of the methodology applied.

#### *2.1. Life-Cycle Assessment (LCA)*

LCA is a powerful tool for performing the potential environmental assessment of impacts and resource consumption throughout a product's life cycle [47]. In this regard, LCA became one of the most relevant methodologies to help organizations perform their activities in the most environmental friendly way along the whole value chain.

In the current study, LCA was conducted following the recommendations of International Organization for Standardization (ISO) 14040 and 14044 international standards [48,49], and was divided into four stages.

#### (i). Goal and Scope Definition

This stage consisted of a detailed definition of the intended application of the study in terms of the system boundaries and functional unit (FU) [50,51]. Moreover, allocation procedures, cut-off rules, and assumptions were also defined in this step.

#### (ii). Life-Cycle Inventory (LCI)

In this stage, all relevant energy, water, and material consumptions, and emissions, effluents, and residues of the process in a specific temporality and geography were collected [52]. Therefore, in this stage, it was necessary to determine the FLW flows. Moreover, in this study, the LCI also collected data about the nutritional content of the different food categories.

#### (iii). Life-Cycle Impact Assessment (LCIA)

The LCIA step transformed the inputs and outputs of the LCI into environmental impact over all of the stages involved in the supply chain. In this study, the impact assessment method selected followed the International Reference Life Cycle Data Product Environmental Footprint (ILCD/PEF) recommendations v1.09 for determining the global warming potential (GWP). The consumption of primary energy resources (net calorific value, PED) was determined according to PE International's [53] life-cycle inventory. Once the PED and GWP embodied along the agri-food supply chain (FSC) were determined, the embodied energy and GHG emissions in FLW could be determined using FLW rates from Gustvasson et al. [54], as indicated Equations (1) and (2).

$$\text{FLW}\_{\text{EE},i} = \sum\_{\mathbf{j}=1}^{\mathbf{j}} \left( \text{PED}\_{\mathbf{i},\mathbf{j}} \right) \alpha\_{\mathbf{i},\mathbf{j}} \tag{1}$$

$$\text{FLW}\_{\text{GWP},i} = \sum\_{j=1}^{j} \left( \text{GWP}\_{i,j} \right) \alpha\_{i,j} \tag{2}$$

where PEDi,j is the PED in stage j for food category i, αi,j is the percentage of FLW in stage j for food category i, and GWPi,j is the GWP of stage j for food category i (see the example in Section S4 of the Supplementary Materials).

#### (iv). Interpretation of the Results

#### *2.2. Energy Return on Investment (EROI)*

The concept of EROI is part of the field of net energy analysis (NEA), and is one way of measuring and comparing the net energy availability from different energy sources and processes. In general, EROI can be defined as "the ratio between the energy returned from an energy-gathering activity compared to the embodied energy in that process" [55]. Although this concept was used initially to develop an energy-focused approach to the economy [56], the concept was adapted to calculate ratios between food energy output and food production energy inputs. In this sense, the most commonly used EROI perspective is the human-edible food energy return on industrially energy investment, as indicated in Equation (3). This provides an anthropocentric perspective on the non-renewable resource efficiencies of competing food production technologies [5]. Therefore, the higher the EROI of a food system is, the more "valuable" it is in terms of producing (nutritionally) useful energy output.

$$\text{EROI} = \frac{\sum\_{\text{j=1}}^{\text{J}} \text{Food}\_{\text{i,j}} \times \text{NC}\_{\text{i,j}}}{\sum\_{\text{j=1}}^{\text{J}} \text{PED}\_{\text{i,j}}},\tag{3}$$

where Foodi,j represents the amount of the food category i available in stage j after withdrawing the FLW, NCi,j is the nutritional content of food category i in supply chain j expressed in kcal, and PEDi,j is the PED in stage j for food category i calculated using LCA.

#### *2.3. Data Envelopment Analysis (DEA)*

DEA is a non-parametric method widely used to compare the inputs and outputs of a set of homogenous DMUs. This method focuses on evaluating the performance of DMUs based on the evaluation of relative efficiency of comparable DMUs by estimating an empirical efficient boundary [9].

Charnes, Cooper, and Rhodes (CCR) introduced the DEA method for the first time. The original CCR model was applicable to the assumption of constant returns to scale (CRS) [57]. In this sense, the CCR model considers that the efficiency frontier is a straight line intersecting the point of origin and the best performer(s), as shown in Figure 2a. The best performer is determined by the highest ratio of output to input; thus, as DMU2 has this condition, it is considered as the reference DMU for all other units [38]. The remaining four DMUs in Figure 2a are inefficient since they use more amounts of input than DMU P2 to produce one unit of output.

Banker et al. [58] modified the CCR model by introducing the so-called "convexity constraint", which changed the efficiency frontier from being a straight line to a convex hull (see Figure 2b). The new model, referred to as BCC, was built based on variable returns to scale (VRS). This model presents two advantages as compared to the CCR model: (i) more units can potentially be considered efficient, and (ii) inefficient units can be compared to more appropriate peers [38]. As shown in Figure 2b, using the BCC model, a higher number of units are efficient (DMU1, DMU2 and DMU3).

**Figure 2.** Graphical representation of the efficiency frontier of the (**a**) Charnes, Cooper, and Rhodes (CCR) model, and (**b**) the CCR model modified by Banker et al. (BCC model). Adapted from Hosseinzadeh-Bandbafha et al. [38].

The slacks-based measure (SBM) model is a typical extension of DEA. Although CCR and BCC approaches calculate the efficiencies of the DMUs based on the proportional decrease/increase of inputs/outputs, the SBM model uses the input excess and output shortfall of each DMU to measure its efficiency [59]. Moreover, it is possible to create and estimate models that provide input-oriented or output-oriented approaches for both CRS and VRS envelopments. An input-oriented model attempts to reduce input variables while remaining within the envelopment space. On the other hand, an output-oriented model increases output variables while remaining within the envelopment space [37]. In this study, an input-oriented approach was selected because multiple inputs were used, while there was only one output.

The DEA linear programming (LP) model (M1) deployed to generate CCR efficiency factors of DMUs is as follows [54]:

M1 (CCR model input-oriented to be solved for each DMU k0):

$$\begin{aligned} \mathsf{min}\mathsf{e}\_{\mathsf{CCR}}(\mathsf{k}\_{0}) &= \sum\_{\mathbf{j}=1}^{\mathsf{n}} \mathsf{u}\_{\mathsf{j}} \mathsf{y}\_{\mathsf{jk}\_{0}};\\ \mathsf{s.t.}\\ \sum\_{\mathbf{i}=1}^{\mathsf{m}} \mathsf{v}\_{\mathsf{i}} \mathsf{x}\_{\mathsf{ik}\_{0}} &= 1;\\ \sum\_{\mathbf{j}=1}^{\mathsf{n}} \mathsf{u}\_{\mathsf{j}} \mathsf{y}\_{\mathsf{jk}\_{0}} - \sum\_{\mathbf{i}=1}^{\mathsf{m}} \mathsf{v}\_{\mathsf{i}} \mathsf{x}\_{\mathsf{ik}\_{0}} &\leq 0, \quad \mathsf{k} = 1, \dots, \mathsf{K}; \ \mathsf{j} = 1, \dots, \mathsf{n}; \ \mathsf{i} = 1, \dots, \mathsf{m};\\ \mathsf{u}\_{\mathsf{j}} &\geq 0, \\ \mathsf{v}\_{\mathsf{i}} &\geq 0, \end{aligned}$$

where uj is the weight for output j, vi is the weight for input i, m is the number of inputs, n is the number of outputs, k is the number of DMUs, yjk is the amount of output j of DMU k, and xjk is the amount of input i of DMU k.

On the other hand, the DEA LP model deployed to generate BCC efficiency factors of DMUs is as follows [58]:

M2 (BCC model input-oriented to be solved for each DMU k0):

$$\begin{aligned} \min & \theta\_{\text{BCC}}(\mathbf{k}\_0) = \sum\_{\mathbf{j}=1}^n \mathbf{u}\_{\mathbf{j}} \mathbf{y}\_{\mathbf{jk}\_0} - \mathbf{u}(\mathbf{k}\_0); \\ \text{s.t.} & \\ \sum\_{\mathbf{i}=1}^m \mathbf{v}\_{\mathbf{i}} \mathbf{x}\_{\mathbf{ik}\_0} &= 1; \\ \sum\_{\mathbf{j}=1}^n \mathbf{u}\_{\mathbf{j}} \mathbf{y}\_{\mathbf{jk}\_0} - \sum\_{\mathbf{i}=1}^m \mathbf{v}\_{\mathbf{i}} \mathbf{x}\_{\mathbf{ik}\_0} - \mathbf{u}\_{(\mathbf{k}\_0)} &\le 0, \\ \mathbf{k} &= 1, \dots, \mathbf{K}; \ \mathbf{j} = 1, \dots, \mathbf{n}; \ \mathbf{i} = 1, \dots, \mathbf{m}; \\ \mathbf{u}\_{\mathbf{j}} &\ge 0, \\ \mathbf{v}\_{\mathbf{i}} &\ge 0, \end{aligned}$$

where uj is the weight for output j, vi is the weight for input i, m is the number of inputs, n is the number of outputs, k is the number of DMUs, yjk is the amount of output j of DMU k, and xjk is the amount of input i of DMU k.

To combine the LCA and DEA methods, a five-step procedure was followed [60,61]: (i) individual LCI for each of the DMUs; (ii) LCIA for every DMU included in the inventory (in this case, using the PED and GWP indicators, as explained in Section 2.1); (iii) determination of the operational efficiency for each DMU; (iv) LCIA of the target DMUs; (v) quantification of the environmental consequences of operational inefficiencies [42].

#### **3. Case Study: The Spanish Food Basket**

#### *3.1. Goal and Scope*

The Spanish Agency of Food Security and Nutrition (AECOSAN), along with the Spanish Society of Community Nutrition (SENC), is the institution that sets the nutritional recommendations for the Spanish population. These recommendations take into account the main diet-related public health challenges, as well as cultural habits of the Spanish population, to promote healthy eating and regular physical activity [62].

In this study, a basket of products was selected based on the consumption data reported by the Ministry of Agriculture, Fisheries, and Food (MAPAMA) [63]. These food commodities were classified according to eleven categories following the Food and Agriculture Organization of the United Nations (FAOSTAT) classification: eggs, meat and animal fat, fish and seafood, dairy, cereals, sweets, pulses, vegetable oils, vegetables, fruits, and roots and tubers. The food commodities considered in each category are collected in Table 1.

The goal and scope of this study was to assess the energetic and environmental efficiency of the Spanish food basket by means of the combination of LCA and DEA approaches. An additional goal of this analysis was to determine whether the use of DEA methodology is a viable tool for assessing the energy efficiency of food systems with FLW. The results are expected to provide an overview regarding the energy efficiency of the different food categories under study with the aim of identifying the most inefficient and the stages where improvement measures should be applied.

**Table 1.** Food commodities included in the study.

#### *3.2. Function, Functional Unit, and System Boundaries*

The methodology proposed included three main parts: LCA, EROI, and the combination of LCA and DEA. Therefore, the function of the study was the estimation of the environmental impacts of the Spanish food basket and the determination of its energy efficiency. To quantify this function, it was necessary to define an FU to which inputs and outputs would be referred. In this case, the FU was described as the food basket with the representative food products consumed by a Spanish citizen in a year, covering the daily energy requirement of 2000 kcal of an adult.

The system boundaries comprise the entire supply chain of a food system following Garcia-Herrero et al. [4] (see Figure 3), i.e., agricultural production, postharvest and storage, industrial processing, distribution (i.e., retail/wholesale), and consumption. The consumption stage was divided into household consumption and related extradomestic consumption.

#### *3.3. Life-Cycle Inventory (LCI)*

An extensive LCI was built up using data from the literature (for more details, see Batlle-Bayer et al. [64]) and the PE GaBi database [53]. The food balance sheet and the unavoidable FLW percentages of each food category were collected from García-Herrero et al. [3], where production values were mainly sourced from Eurostat [65–68]. Imported food products were taken into account based on trade statistics compiled by the Spanish Tax Agency (AEAT) [69]. The food balance sheet for Spain in 2015 and the FLW percentages are available in Sections S1 and S2 in the Supplementary Materials (SM).

For transportation, 400-km and 100-km distances to wholesale and retailers were assumed, respectively. Electricity consumption for retail storage was considered, assuming two days of storage for products requiring cooling conditions, and 15 days under freezing conditions. Transport to home was estimated based on Milà i Canals et al. [70]. Regarding home storage, data from the LCA Food

Database was used [71]. For cooking, energy factors of Foster et al. [72] were taken (for more details, see Batlle-Bayer et al. [64]).

On the other hand, nutritional data for the EROI estimation were obtained from the food composition tables of the Institute for Education in Nutrition and Dietetics from Spain (CESNID) [73]. Such tables are registered in the FAO's International Network of Food Data Systems [74]. The nutritional content for each food category under study is collected in Section S3 in the Supplementary Materials.

**Figure 3.** System boundaries for the Spanish food basket using a life-cycle assessment (LCA)/data envelopment analysis (DEA) approach.

#### *3.4. Selection of the DEA Model*

The selected DEA model was the slacks-based measure of efficiency (SBM). The choice of model was based on its flexibility concerning the computation of the DMUs irrespective of the units of measure used for the different inputs/outputs [42]. Moreover, for the sake of further discussion, the CRS and VRS approaches were computed in parallel in the current study under an input-oriented approach.

#### *3.5. Inputs and Output Selection for the DEA Matrix*

One DEA matrix composed of 11 DMUs was computed in this study. The number of DMUs under study must satisfy Equation (4). Each DMU included two different inputs (see Figure 3). As input 1, the primary energy invested in food production was computed in MJ/cap/day, whereas the GWP was considered as input 2. As mentioned previously, an input-oriented approach was based on the assumption that inputs have to be minimized and outputs have to be maximized. However, in some situations, undesirable (bad) inputs and outputs may be presented in the production process [75]. In this case, the GWP was a bad output that was computed as an input to be minimized. Finally, the nutritional energy supply was computed in kcal/cap/day as output.

$$m \ge \max\left\{m \times s, 3 \times (m+s)\right\},\tag{4}$$

where *m* is the number of inputs used in the DEA study, and *s* is the number of outputs involved.

#### **4. Results and Discussion**

Following the methodology described in Section 2, the LCA of the 11 food categories selected was performed, and the PED and GWP were calculated. These values were obtained from a previous study developed by Batlle-Bayer et al. [64]. As shown in Table 2, the total energy invested in food production daily amounts to around 19,500 kcal for an average Spanish citizen. This involves an emission of close to 4 kg of CO2 equivalent per day, with meat production as the most responsible: 28% for PED and 42% in terms of CO2 equivalent emissions. Dairy, and fish and seafood are the second main categories contributing to GWP, each representing 12% of total emissions. Results suggest a high correlation between the investment of energy along the FSC and the generation of GHG emissions, owing to the high dependency of the energy matrix on fossil fuels. Conversely, this is not directly related to the nutritional energy provided to consumers. As displayed in Table 1, the major contributors to the Spanish average diet are vegetable oils and cereals, responsible for nearly half of the energy supply. However, the contribution of meat to nutritional energy supply is relegated to less than 12%, with eggs as the food category providing the lowest amount of nutritional energy to consumers (2%). These results are in agreement with data reported by the European Commission [76], which confirm that livestock products, such as meat, fish and seafood, and dairy products, incorporate a substantial amount of energy; in addition, the consumption of dairy and meat products has a major role in GHG emissions.

**Table 2.** Energy invested from cradle to plate for the food categories under study, related CO2 eq. emissions, and nutritional energy provided to consumer. PED—primary energy demand; GWP—global warming potential; FLW—food losses and waste.

When considering FLW, it was observed that embodied energy waste duplicates the daily energy supplied to consumers (4685 kcal/cap/d). As shown in Figure 4a, the largest contributor to this fact is the meat category, accounting for 25%. This is directly related to the unnecessary emission of CO2 eq., which also contributes 40% to the total GWP (Figure 4b). It is followed by cereals and fish and seafood categories, for which a 22% and 18% embodied energy waste was estimated, respectively, whereas a 15% and 12% of embodied GWP waste was calculated.

**Figure 4.** Contribution of each food category in the average Spanish diet to the (**a**) embodied energy, and (**b**) embodied global warming potential (GWP).

These findings suggest the need for estimating the efficiency of the Spanish agri-food system. As a first approach, the EROI indicator was used to estimate the energetic efficiency. As displayed in Table 3, the maximum EROI was observed for sweets and vegetable oils, which means that such food categories are the most efficient from an energy perspective. These results also suggest the need of studying more in-depth nutrients that should be either encouraged or limited, in order to provide a nutritional efficiency perspective. Results also show that animal-based products such as eggs, meat, and fish present the lowest EROI values (3.1–4.0%), which agree with the study of Pelletier et al. [5], where an eightfold difference between animal- and plant-based products was determined. An exception was observed in our results for vegetables. This is mainly due to the inclusion of processed commodities in this category, particularly tomato sauce, which exhibits an EROI value of 1.2%.

**Table 3.** Energy return on investment (EROI) values for food categories in the Spanish diet.

However, these efficiency results do not consider the environmental impacts of food production along the supply chain. For this reason, DEA was combined with LCA, allowing a wider scope of the efficiency of the Spanish agri-food system to be provided. To perform an LCA/DEA analysis, it is firstly necessary to elaborate the DEA matrix. Table 4 represents the DEA matrix, composed of 11 DMUs, with each one represented by two inputs and one output, as explained in Section 3.4.

DEA Frontier was the software used for the computation of the DEA matrix [77]. Results using CRS and VRS perspectives were compared and are shown in Figure 5, which includes the efficiency score for each DMU. A food category was considered inefficient with an efficiency score φ < 1, whereas ϕ = 1 represented an efficient food category. According to that, when the CRS perspectives were considered, one food category was efficient (sweets), while a 0.97 efficiency was observed for vegetable oils (see Figure 5a). The CRS approach presented a wide range of efficiencies, between 0.06 and 0.97. Egg production was the least efficient food category, followed by fish, meat, and vegetables, with similar efficiencies ranging from 0.06 to 0.08. On the other hand, when the VRS model was applied, a total of two food categories were considered efficient (sweets and vegetable fats; see Figure 5b). Despite the fact that the VRS approach presented an additional efficient DMU, the range of efficiencies was lower than when using a CRS perspective, ranging from 0.06 to 0.56. Finally, the average efficiencies (including efficient DMUs) were similar: 0.33 ± 0.34 for CRS and 0.40 ± 0.34 for VRS. Some authors state that using a VRS approach rather than a CRS formulation leads to a higher number of efficient DMUs since the constraint set for CRS is less restrictive and, consequently, lower efficiency scores are possible [78]. In addition, the convexity constraint added in the VRS simply guarantees that each DMU is only compared to others of similar size [43]. Therefore, considering that the unit of reference (i.e., DMU) is, in each case, a different food category and that the food categories inventoried in this study present different characteristics, we considered that the use of the VRS approach should prevail.

**Figure 5.** Efficiency (φ) score of the selected food categories (i.e., decision-making units (DMUs)) using the (**a**) constant return to scale (CRS) model, and (**b**) variable return to scale (VRS) model.

The previous findings agree with the pattern observed for the EROI assessment. While this can be taken as a verification of both approaches, DEA also provides a further exploration of the operational inefficiencies and the improvement actions of the system under study. Hence, the DEA model was also used to formulate new virtual and efficient values for the inputs of the inefficient DMUs, by projecting the inefficient scores on the efficient targets established. Both PED and GWP were subjected to minimization while maintaining the same nutritional supply.

Figures 6 and 7 compare the current environmental impacts of the different DMUs with the target PED and GWP calculated with the LCA/DEA methodology. In other words, it presents the environmental savings that would be reached in these DMUs if they were to operate under the efficient conditions projected in the computed matrix. In this regard, when the CRS perspective was used, primary energy invested savings ranged from 2.83% to 96.7%, while, when the VRS approach was used, reductions ranged from 52.9% to 95.3%. The highest reductions were obtained by the most inefficient food categories (i.e., vegetables, meat, fish, and eggs). On the other hand, GWP savings followed the same trend, ranging from 19% to 96.4% for CRS and from 43.6% to 95.8% for VRS. Similarly, the highest reductions were also reached by vegetables, meat, fish, and eggs.

**Figure 6.** Primary energy invested (PED) for original DMUs (blue bar) and virtual targets for the CRS model (orange bar) and VRS model (green bar). The lines represent the overall environmental improvements for the CRS model (orange line) and VRS model (green line).

**Figure 7.** Global warming potential (GWP) for original DMUs (blue bar) and virtual targets for the CRS model (orange bar) and VRS model (green bar). The lines represent the overall environmental improvements for the CRS model (orange line) and VRS model (green line).

The results suggest that the planification of citizens' diet may affect the quantity of embodied energy waste and related environmental impacts. For this reason, specific strategies should be addressed to the categories revealing larger inefficiency scores, such as meat, eggs, fish, and vegetables, which all have better nutritional reputation than sweets or processed cereals. Figure 8 displays the contribution of each stage of the FSC to the PED for eggs, meat, fish and seafood, and vegetable categories. As stated by Batlle-Bayer et al. [64], in terms of the life-cycle processes, the primary production phases were the major contributors to the PED for eggs (85%), meat (65%), and fish and seafood (59%), while the household stage was the main contributor for vegetables (52%). These results are related to the production of feed in the case of eggs and meat production [60]. According to Ghasempour and Ahmadi [79], the production of feed (corn, soybean, and wheat), as well as the poultry equipment, used in the production of eggs is associated wuth the most energy consumption, whereas Skunca et al. [17] suggested the use of grain legumes as a protein source in feed instead of soybean-based ingredients. This fact could decrease the resource consumption since the cultivation of grain legumes does not require mineral fertilizer application. Moreover, in the case of beef, the use of electricity and diesel during the suckling cow–calf stage was also significant in the primary production phase [64].

Other authors, such as Avadí et al. [80,81], Laso et al. [22], and Vázquez-Rowe et al. [82] assessed the life cycle of processed fish products identifying the fishery stage as one of the most energy-intensive due to the use of fossil fuels. Batlle-Bayer et al. [64] analyzed eight fish and seafood species (i.e., mussels, shrimps and prawns, Atlantic mackerel, European hake, sardines, salmon, tuna, and octopus), identifying the use of diesel as the main contributor (around 98%) in the fishery stage, and the use of electricity (around 97%) in the culture stage (mussels farming).

On the other hand, for vegetables, several foods were considered (tomato, lettuce, and veggies). The fact that household consumption was the stage with the highest contribution to PED was due to the use of energy when cooking veggies (i.e., boiling and frying) and cold storage for preserving their shelf life. However, the agricultural production stage for vegetables was also important, representing 26% of total PED, due to the cultivation of lettuce, in particular, the use of diesel for machinery. In fact, although the amount of PED required for meat, fish and seafood, and vegetable consumption is similar (Table 4), it was distributed in different stages. Finally, the processing stage represented 19% of the total impact due to the inclusion of processed products, such as tomato processed products, which use fuel oil and electricity in their processing.

**Figure 8.** Distribution of the primary energy demand (PED) between the different stages of the food supply chain for the most inefficient food categories: eggs, meat, fish and seafood, and vegetables.