The following article has been accepted by Applied Physics Letters. After it
is published, it will be found at: https://aip.scitation.org/journal/apl

# Elastic three-dimensional metaframe for selective wave filtering and
polarization control

J. M. De Ponti<EMAIL_ADDRESS>Dept. of Civil and Environmental
Engineering, Politecnico di Milano, Piazza Leonardo da Vinci 32 20133 Milano
Italy  E. Riva  Dept. of Mechanical Engineering, Politecnico di Milano, Via
Giuseppe La Masa 1 20156 Milano Italy  F. Braghin  Dept. of Mechanical
Engineering, Politecnico di Milano, Via Giuseppe La Masa 1 20156 Milano Italy
R. Ardito  Dept. of Civil and Environmental Engineering, Politecnico di
Milano, Piazza Leonardo da Vinci 32 20133 Milano Italy

###### Abstract

We experimentally achieve selective wave filtering and polarization control in
a three-dimensional elastic frame embedding local resonators. By connecting
multi-resonating elements to a frame structure, a complete low-frequency,
subwavelength bandgap with strong selective filtering properties is obtained.
Theory and experiments demonstrate the metaframe capability to selectively
stop transverse waves while allowing longitudinal wave propagation, as in
’fluid-like’ elasticity. This peculiar behaviour, together with the complete
bandgap structure, may open opportunities in the context of wave control,
envisaging concurrent applications for three-dimensional filters and elastic
wave polarizers.

††preprint: APS/123-QED

A dominant line of research in the wave physics focuses on the investigation
of materials and structures to control the propagation of waves. Periodic
media have been demonstrated particularly influential to manipulate waves
exploiting Bragg scattering bandgaps, creating the fields of photonic crystals
in optics [1] and their acoustic counterparts in phononic structures [2, 3,
4]. In the context of elastic waves, notable effort has been devoted to the
definition of low frequency and subwavelength devices to comply with ambient
and most widespread vibration spectra. To push the operational regime towards
lower frequencies, the exploitation of local resonance has received
considerable attention [5, 37, 7, 8], leading to the definition of
metamaterials [9, 10]. Potential applications of elastic metamaterials could
be implemented at any lengthscale, for seismic wave protection, [11, 12, 13,
14, 15, 16], vibration absorption [17, 18, 19, 20], nondestructive evaluation
[21] and several nuances for wave enhancement and manipulation [22, 23, 24]
recently inspired by topological insulators[25, 26, 27]. Metamaterials are
often employed in combination with smart materials to obtain multifunctional
devices for enhanced sensing or energy harvesting applications [28, 29, 30].
The effect of local resonance on the dispersion properties in continuous
elastic substrates has been widely investigated in the past through analytical
and approximated formulations in beams [31, 32], plates [33, 35, 34] and
halfspaces [35, 36], amongst others. Motivated by the idea to control elastic
waves in a complete three-dimensional setting, we propose a metaframe able to
provide selective wave filtering and polarization control. Moreover, we show
that simple lumped models can be suitably employed to predict the behaviour of
such structures, providing a rapid analytical design tool. We show that the
attenuation performance of locally resonant structures strongly depends on the
relative stiffness ratio (frame vs. resonator), which should be considered an
additional parameter to the conventional total mass ratio [32]. Even if a
complete bandgap is obtained, the imaginary dispersion curves are remarkably
different for longitudinal and transverse waves; this behavior is linked to
the stiffness ratio between the frame and the resonator, which affects the
attenuation properties of the structure. More specifically, we demonstrate
selective attenuation of transverse waves, while allowing longitudinal wave
propagation, as in ’fluid-like’ elasticity [37, 38]. The possibility to create
polarization bandgaps by a selective suppression of each vibration mode has
been demonstrated in stiff-in-soft structures with a single resonator per
cell, obtaining different bandgaps for flexural, longitudinal or torsional
modes [37], or leveraging multiple resonances [38]. Modes can be also isolated
on the basis of their polarization by creating distinct topological interfaces
[39] supporting heterogeneous helical-valley edge waves. Differently from
these approaches, we show a three-dimensional metaframe with the coexistence
of attenuation and trasmission of transverse and longitudinal waves, both
inside a common bandgap.

Figure 1: Elastic three-dimensional metaframe. The structure is obtained by
connecting a set of resonators to the nodes of an elastic frame. In the bottom
part, a detailed view of the unit cell made of a frame with 8 resonators
attached is shown, together with the main geometrical dimensions and the high
symmetry points used to compute the dispersion curves. Figure 2: (a) Numerical
dispersion curves of the metaframe obtained through FEM, superimposed to the
dispersion of a plain elastic frame without the resonators. A low frequency
bandgap with gap-mid gap ratio of $35.2\%$ is observed in the range
$58.04-83.28$ $Hz$ for both longitudinal and transverse waves. Insets show the
vibration modes of the plain frame and the metaframe at the $\Gamma$ point of
the bandgap edges. The bandgap opens with an in-phase response of the
resonators, while the closure is related to out-of-phase oscillations. (b)
Wave propagation analyses on an infinitely long elastic frame embedding the
metaframe. Selective wave filtering is observed: for the same input,
longitudinal waves propagate, while transverse waves are attenuated.

Although this mechanism may resemble an elastic three-dimensional analogue of
optical waveguide polarizers in photonics [40, 41], the concept of wavelength
selective filtering is caused here by a complete bandgap with an imaginary
component strongly dependent on the wave polarization. This can be regarded as
an additional feature coming from the higher complexity of wave polarization
in elasticity, which allows creating complete bandgaps with strong selective
filtering properties. This mechanism can be employed for the creation of
three-dimensional single phase selective filters, or elastic wave polarizers.

On the basis of previous research results [42], we consider the metaframe
depicted in Fig. 1, obtained connecting local resonators to an elastic frame.
The finite structure is defined by periodically repeating a unit cell
containing 8 resonators made of a thin elastic beam equipped with a sphere on
the tip, and connected to a frame node. Such arrangement allows to increase
the resonating mass, while keeping the stiffness of the resonators low.

The unit cell, of size $a=50$ mm, includes a frame with beams of square cross
section of dimension $w_{frame}=0.08a$, while the resonators are characterized
by a cylindrical elastic connection of diameter $w_{res}=0.04a$, length
$l_{res}=0.14a$ and a sphere of radius $r_{sphere}=0.21a$. The entire
structure is made of Nylon $PA12$, with Young’s modulus $E_{PA12}=1586$ MPa,
Poisson’s ratio $\nu=0.4$ and density $\rho=1000$ kg/m3.
The numerical dispersion curves of the metaframe, obtained through Abaqus®,
are displayed in Fig. 2 (a), along with the graphical representation of a
number of Bloch modes. A complete low- frequency bandgap in the range
$58.04-83.28$ Hz, with gap-mid gap ratio [42] of $35.2\%$ is observed. In
terms of normalized frequency, i.e. in units of the resonance frequency, we
obtain a bandgap in the range $0.97-1.39$. In correspondence of the bandgap
edges, the Bloch modes exhibit a localized flexural motion of the resonators,
with in-phase and out-of-phase response for the mode below and above,
respectively. At this step, we also observe that the number of propagating
wave modes that collapse in a bandgap can be dynamically affected in a
different manner from the resonator motion, due to the associated equivalent
stiffness, which is very different for longitudinal and transverse waves. Fig.
2 (b) shows the wave propagation along an infinitely long elastic frame with
the insertion of the metaframe composed of 9 unit cells. Spurious edge
reflections are avoided by imposing absorbing boundary conditions at the end
of the frame along the wave propagation direction. They are implemented in
Abaqus® using the ALID (Absorbing Layers using Increasing Damping) method,
adopting a cubic law function for mass proportional Rayleigh damping [43]. By
inspecting the wave amplitude after the metaframe it can be noticed that, for
a constant input excitation, longitudinal waves propagate, while transverse
waves are attenuated. This phenomenon is observed independently on the
metaframe dimension, as can be noticed considering larger samples (see
Supplementary Material). While the influence of different resonator modes on
the attenuation performances has already been investigated [33, 35], comparing
for e.g. flexural and longitudinal resonances, we focus our attention on a
single mode of the resonator, changing the stiffness of the frame only. This
is possible thanks to the three-dimensional symmetry of the cell, which
guarantees a flexural polarization of the resonators in accord with both
longitudinal and transverse waves.

To investigate on this aspect, we adopt a simplified lumped model to predict
the bandgap characteristics, and the associated level of attenuation. We
consider a periodic spring-mass chain, made of unit cells of stiffness $k$ and
mass $M$, with attached a resonator of stiffness $k_{R}$ and mass $m_{R}$, as
shown in the inset of Fig. 3. The quantity $k$ represents the stiffness of the
plain frame, while $m$ is the overall non-resonating mass. We adopt $2$
degrees of freedom to describe the translation of the main mass $M$ and of the
resonator $m_{R}$.

Figure 3: Analytical dispersion curve for the metaframe, superimposed to the
numerical dispersion (scattered points) from FEM. Even if the bandgap exists
for both longitudinal and transverse waves, the imaginary part of the
dispersion curves is remarkably different, reporting almost negligible
attenuation for longitudinal waves.

By defining $\Omega^{2}=k/M$, $\omega_{R}^{2}=k_{R}/m_{R}$, and introducing
the non-dimensional quantity $\mu=\kappa d$, being $d$ the unit cell size and
$\kappa$ the wavenumber, the dispersion relation is:

$\cos(\mu)=1-\frac{\omega^{2}}{2\Omega^{2}}\left[1+\frac{m_{R}}{M}\frac{\omega_{R}^{2}}{\omega_{R}^{2}-\omega^{2}}\right]$
(1)

Evanescent waves are present if $\cos(\mu)\geq 1$ and $\cos(\mu)\leq-1$. By
imposing these conditions we define, after some algebraic manipulations (see
Supplementary Material), the limits of the bandgap as:

4Ω2+ωR2(1+mRM)-[4Ω2+ωR2(1+mRM)]2-16ωR2Ω22≤ω2≤ωR2(1+mRM)ω2≥4Ω2+ωR2(1+mRM)+[4Ω2+ωR2(1+mRM)]2-16ωR2Ω22

(2)

The equations reveal that the bandgap does not open at the resonance frequency
$\omega_{R}$, but at a frequency which depends also on the inertial and
elastic properties of the main structure.

The longitudinal (axial) and transverse (bending) stiffness of the frame can
be easily computed by means of well-known techniques of structural mechanics
(see Supplementary Material). Moreover, we model the resonator as a Timoshenko
elastic beam with a lumped rigid mass, assuming as Lagrangian coordinates the
translation of the tip of beam and the rotation of the mass. By doing so, we
also neglect the mass of the beam ($0.46\%$ of the one of the sphere) and
assume an infinite longitudinal (axial) stiffness. By solving the eigenvalue
problem for the $2$ degrees of freedom model of the resonator, we obtain that
$82.4\%$ of the mass of the sphere participates to the motion, at the
resonance frequency $f_{R}=60.7$ Hz. Having computed the mass and stiffness
parameters, we get from Eq. 1 the analytical dispersion curve reported in
Fig.3. A good agreement between the lumped and the three-dimensional finite
element model is observed. The simplified model provides accurate predictions
in the frequency range of interest, with an error of approximately $4\%$ for
the bandgap opening and closing frequencies. The plot shows clearly that the
bandgap opening frequency for the transverse excitation is significantly
smaller than the case of longitudinal wave. As a consequence, there is a
frequency range in which transverse waves are evanescent and longitudinal
waves are amplified; the relative width of such a range, with respect to the
mid range frequency, is equal to $13\%$. More importantly, the lumped model is
used to easily evaluate the imaginary part of the dispersion relation, which
quantifies the actual level of attenuation provided within the bandgap [3]. It
can be noticed that transverse waves are strongly attenuated, while the
longitudinal ones are slightly affected by the bandgap, except for a narrow
region centered at resonance frequency of the resonator. This feature allows
the metaframe to selectively attenuate transverse waves and enable
longitudinal wave propagation. That is, for a mixed longitudinal and
transverse polarization of the input at the bandgap frequencies, one can make
a filter for transverse wave components, while keeping unaltered the
longitudinally polarized waves. The different attenuation levels are
attributed to the stronger mismatch between the axial stiffness of the frame
and the flexural stiffness of the resonators. In terms of waves, a higher
wavelength is obtained, thus increasing the subwavelength requirements of the
system.

The theoretical results are hereafter validated experimentally. Prototypes are
built using Selective Laser Sintering (SLS), an additive manufacturing process
that has no need of support material (a critical issue for the realization of
suspended parts, like the resonators in this case) while guaranteeing good
precision. The finite structure is composed of $9$ unit cells, resulting in a
cube of dimension $150$ mm. The experimental setup is shown in Fig. 4.

Figure 4: Experimental setup. The excitation is provided through an
electrodynamic shaker. The metaframe is glued to an aluminium thick plate
which is rigidly connected to the shaker. Longitudinal excitation is provided
by vertically placing the plate onto the shaker. For transverse excitation,
the shaker is rotated of $90^{\circ}$ and the prototype is placed on a
$L$-shaped fixture. Input and output accelerations are measured through
accelerometers glued on the rigid plate and on the metaframe, respectively.

The metaframe is glued to an aluminium thick plate which is connected to a LDS
v406 electrodynamic shaker to provide excitation. The aluminium plate is
sufficiently thick to avoid relevant deformations for the input frequency,
i.e. it is assumed as a perfectly rigid connection between the metaframe and
the shaker. Two configurations are tested, imposing a longitudinal and a
transverse polarization of the wave respectively. Both cases consider the same
direction of the wave propagation, $b_{1}$ in Fig.1, but different
polarization, i.e. $b_{1}$ for longitudinal waves and $b_{3}$ for transverse
waves. For the longitudinal wave, the plate is directly mounted onto the
shaker which provides a vertical excitation. In contrast, the transverse wave
is obtained by rotating the shaker of $90^{\circ}$ and connecting the plate to
a $L$-shaped structure. The experimental acceleration spectra are measured
using PCB Piezotronics 352C33 accelerometers with a sensitivity of $10$ mV/$g$
and resonant frequency higher than $70$ kHz. Since the accelerometers are able
to measure the acceleration along one direction only, we rotated them of
$90^{\circ}$ to measure transverse waves. The input is measured through an
accelerometer glued to the thick plate, while the output is measured on the
frame structure at the opposite side, as indicated by red circles and arrows
in Fig. 4. The excitation signal is provided through a KEYSIGHT 33500B
waveform generator, which synchronously starts with the measurement system and
consists in a linear swept-frequency cosine signal modulated by a Hann window.
A frequency range between $50$ Hz to $150$ Hz is considered, filtering out the
other frequency components. Data are then post-processed in order to estimate
the transmission spectra, defined as $\tau\>$[dB] $=20$
log${}_{10}(\ddot{w}_{OUT}/\ddot{w}_{IN})$.

Figure 5: Numerical FEM viscoelastic (continuous lines) and experimental
(dashed-dot lines) transmission spectra for longitudinal and transverse waves.
Longitudinal waves (red lines) are slightly attenuated inside the expected
bandgap region; maximum attenuation is detected at the resonance frequency of
the resonators, marked by the continuous black vertical line. A stronger
attenuation is measured for transverse waves (blue lines), in agreement with
previous dispersion curves predictions.

Fig. 5 shows the numerical and experimental transmission diagrams for
longitudinal and transverse waves. To investigate the results, we carry out
FEM viscoelastic transmission analyses [44, 45], using the Standard Linear
Solid – Maxwell type – model with $\tau_{Maxwell}=1.5e-4$ s and a Young’s
modulus $E_{Maxwell}=315\>$MPa. The experimental data are in satisfactory
agreement with the computations, corroborating the different attenuation
performance of the structure for longitudinal and transverse waves. The
different level of attenuation, and the slight shift in the position of the
transmission peaks is attributed to the relative simplicity of the
viscoelastic model with respect to the actual material behaviour, as well as
experimental misalignments that prevent the possibility to excite longitudinal
and transverse waves separately.

In conclusions, we have experimentally demonstrated potential advantages in
using a three-dimensional metaframe for the creation of a complete, low-
frequency, subwavelength bandgap. The strong mismatch between the longitudinal
and transverse stiffness of the frame allows to achieve a different control of
the waves depending on the polarization direction. Strong attenuation of
transverse waves is achieved, while longitudinal wave propagation is mostly
unaffected by the bandgap, except from a sharp region close to the resonant
frequency of the resonators. This behavior can be suitably employed for
applications involving selective wave filtering and polarization control.

## Supplementary Material

See the supplementary material for a detailed description of the analytical
and numerical models.

## Acknowledgements

The support of the H2020 FET-proactive project MetaVEH under grant agreement
No. 952039 is acknowledged. We also gratefully acknowledge the Italian
Ministry of Education, University and Research for the support provided
through the Project “Department of Excellence LIS4.0—Lightweight and Smart
Structures for Industry 4.0.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * [1] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, R. D. Meade, Photonic Crystals, Molding the Flow of Light, 2nd ed., Princeton University Press, Princeton, NJ (2008).
  * [2] M. S. Kushwaha, P. Halevi, L. Dobrzynski, B. Djafari-Rouhani, Acoustic band structure of periodic elastic composites, Phys. Rev. Lett. 71, 2022 (1993).
  * [3] V. Laude, Phononic Crystals, Artificial Crystals for Sonic, Acoustic, and Elastic Waves, De Gruyter Studies in Mathematical Physics, 26 (2015).
  * [4] M. I. Hussein , M. J. Leamy, M. Ruzzene, Dynamics of Phononic Materials and Structures: Historical Origins, Recent Progress, and Future Outlook, Appl. Mech. Rev., 66(4): 040802 (2014).
  * [5] Z. Liu, X. Zhang, Y. Mao, Y. Y. Zhu, Z. Yang, C. T. Chan, P. Sheng, Locally Resonant Sonic Materials, Science, Vol. 289, Issue 5485, pp. 1734-1736 (2000).
  * [6] G. Ma, P. Sheng, Acoustic metamaterials: From local resonances to broad horizons, Sci. Adv., Vol. 2, no. 2, e1501595 (2016).
  * [7] A. E. Miroshnichenko, S. Flach, Y. S. Kivshar, Fano resonances in nanoscale structures, Rev. Mod. Phys., 82, 2257–98 (2010).
  * [8] F. Lemoult, M. Fink, G. Lerosey, Acoustic Resonators for Far-Field Control of Sound on a Subwavelength Scale, Phys. Rev. Lett., 107, 064301 (2011).
  * [9] R. V. Craster, S. Guenneau, Acoustic Metamaterials Negative Refraction, Imaging, Lensing and Cloaking, Springer Series in Materials Science, ISBN 978-94-007-4812-5 (2013).
  * [10] R. V. Craster, S. Guenneau, World Scientific Handbook of Metamaterials and Plasmonics: Volume 2: Elastic, Acoustic and Seismic Metamaterials (Singapore: World Scientific) (2017).
  * [11] S. Brûlé, E. H. Javelaud, S. Enoch, S. Guenneau, Experiments on seismic metamaterials: molding surface waves, Phys. Rev. Lett., 112, 133901 (2014).
  * [12] Y. Achaoui, T. Antonakakis, S. Brule, R. V. Craster, S. Enoch, S. Guenneau, Clamped seismic metamaterials: Ultra-low frequency stop bands, New J. Phys., 19, 063022 (2017).
  * [13] G. Finocchio, O. Casablanca, G. Ricciardi, U. Alibrandi, F. Garescì, M. Chiappini, B. Azzerboni, Seismic metamaterials based on isochronous mechanical oscillators Appl. Phys. Lett., 104, 191903 (2014).
  * [14] A. Colombi, P. Roux, S. Guenneau, P. Gueguen, R. V. Craster, Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances, Sci. Rep., 6, 19238 (2016).
  * [15] M. Miniaci, A. Krushynska, F. Bosia, N. M. Pugno, Large scale mechanical metamaterials as seismic shields, New J. Phys., 18, 083041 (2016).
  * [16] A. Colombi, D. Colquitt, P. Roux, S. Guenneau, R. V. Craster, A seismic metamaterial: The resonant metawedge, Sci. Rep., 6, 27717 (2016).
  * [17] A. Colombi, V. Ageeva, R. J. Smith, A. Clare, R. Patel, M. Clark, D. Colquitt, P. Roux, S. Guenneau, R. V. Craster, Enhanced sensing and conversion of ultrasonic rayleigh waves by elastic metasurfaces. Sci. Rep. 7, 6750 (2017).
  * [18] K. H. Matlack, A. Bauhofer, S. Krödel, A. Palermo, C. Daraio, Composite 3D-printed metastructures for low-frequency and broadband vibration absorption, Proc. Natl Acad. Sci., 113, 8386–90 (2016).
  * [19] N. Kaina, M. Fink, G. Lerosey, Composite media mixing Bragg and local resonances for highly attenuating and broad bandgaps, Sci. Rep., 3, 3240 (2013).
  * [20] M. Moscatelli, R. Ardito, L. Driemeier, C. Comi, Band-gap structure in two- and three-dimensional cellular locally resonant materials, J. Sound Vib., 454, pp. 73-84 (2019).
  * [21] M. Molerón, C. Daraio, Acoustic metamaterial for subwavelength edge detection, Nat. Commun. 6, 8037 (2015).
  * [22] V. Romero-García, R. Picó, A. Cebrecos, V. J. Sánchez-Morcillo, K. Staliunas, Enhancement of sound in chirped sonic crystals, Appl. Phys. Lett. 102, 091906 (2013).
  * [23] Y. Chen, H. Liu, M. Reilly, H. Bae, M. Yu, Enhanced acoustic sensing through wave compression and pressure amplification in anisotropic metamaterials, Nat. Commun. 5, 5247 (2014).
  * [24] Y. Liu, Z. Liang, F. Liu, O. Diba, A. Lamb, J. Li, Source illusion devices for flexural lamb waves using elastic metasurfaces, Phys. Rev. Lett. 119, 034301 (2017).
  * [25] S. H. Mousavi, A. B. Khanikaev, Z. Wang, Topologically protected elastic waves in phononic metamaterials, Nat. Commun., Vol. 6, 8682 (2015).
  * [26] R. Süsstrunk, S. D. Huber, Observation of phononic helical edge states in a mechanical topological insulator, Science, Vol. 349, Issue 6243, pp. 47-50 (2015).
  * [27] G. J. Chaplain, J. M. De Ponti, G. Aguzzi, A. Colombi, R. V. Craster, Topological Rainbow Trapping for Elastic Energy Harvesting in Graded Su-Schrieffer-Heeger Systems, Phys. Rev. Applied, 14, 054035 (2020).
  * [28] C. Sugino, A. Erturk, Analysis of multifunctional piezoelectric metastructures for low-frequency bandgap formation and energy harvesting, J. Phys. D: Appl. Phys., 51, 215103 (2018).
  * [29] J. M. De Ponti, A. Colombi, R. Ardito, F. Braghin, A. Corigliano, R. V. Craster, Graded elastic metasurface for enhanced energy harvesting, New J. Phys., 22, 013013 (2020).
  * [30] J. M. De Ponti, A. Colombi, E. Riva, R. Ardito, F. Braghin, A. Corigliano, R. V. Craster, Experimental investigation of amplification, via a mechanical delay-line, in a rainbow-based metamaterial for energy harvesting, Appl. Phys. Lett. 117, 143902 (2020).
  * [31] E. Baravelli, M. Ruzzene, Internally resonating lattices for bandgap generation and low-frequency vibration control, J. Sound Vib., Vol. 332, Issue 25, pp. 6562-6579 (2013).
  * [32] C. Sugino, S. Leadenham, M. Ruzzene, A. Erturk, On the mechanism of bandgap formation in locally resonant finite elastic metamaterials, J. Appl. Phys., 120, 134501 (2016).
  * [33] E. G. Williams, P. Roux, M. Rupin, W. A. Kuperman, Theory of multiresonant metamaterials for $A_{0}$ lamb waves, Phys. Rev. B, 91, 104307 (2015).
  * [34] J. Y. Yoritomo, R. L. Weaver, On band gap predictions for multiresonant metamaterials on plates, J. Acoust. Soc. Am. 139, 1282 (2016)
  * [35] D. J. Colquitt, A. Colombi, R. V. Craster, P. Roux, S. R. L. Guenneau, Seismic metasurfaces: Sub-wavelength resonators and Rayleigh wave interaction, J. Mech. Phys. Solids, Vol. 99, pp. 379-393 (2017).
  * [36] C. Comi, M. Moscatelli, J. J. Marigo, Two scale homogenization in ternary locally resonant metamaterials, Mat. Phys. Mech., 44 (1), pp. 8-18 (2020).
  * [37] G. Ma, C. Fu, G. Wang, P. del Hougne, J. Christensen, Y. Lai, P. Sheng , Polarization bandgaps and fluid-like elasticity in fully solid elastic metamaterials, Nat. Commun. 7, 13536 (2016).
  * [38] Y. Lai, Y. Wu, P. Sheng, Z. -Q. Zhang, Hybrid elastic solids, Nat. Mat. 10, pp. 620–624 (2011).
  * [39] M. Miniaci, R. K. Pal, R. Manna, M. Ruzzene, Valley-based splitting of topologically protected helical waves in elastic plates, Phys. Rev. B 100, 024304 (2019).
  * [40] R. K. Sinha, Y. Kalra, Design of optical waveguide polarizer using photonic band gap, Optics Express, 14, 22, pp. 10790-10794 (2006).
  * [41] J. Zimmermann, M. Kamp, A. Forchel, R. März, Photonic crystal waveguide directional couplers as wavelength selective optical filters, Optics Communications, 230, 4–6, pp. 387-392 (2004).
  * [42] L. D’Alessandro, E. Belloni, G. D’Alo, L. Daniel, R. Ardito, A. Corigliano, F. Braghin, Modelling and experimental verification of a single phase three-dimensional lightweight locally resonant elastic metamaterial with complete low frequency bandgap, 11th International Congress on Engineered Material Platforms for Novel Wave Phenomena, Metamaterials 2017, art. no. 8107842, pp. 70-72 (2017).
  * [43] P. Rajagopal, M. Drozdz, E. A. Skelton, M. J. S. Lowe, R. V. Craster, On the use of absorbing layers to simulate the propagation of elastic waves in unbounded isotropic media using commercially available Finite Element packages, NDT $\&$ E International 51, pp 30-40 (2012).
  * [44] A. Krushynska, V. G. Kouznetsova, M. G. D. Geers, Visco-elastic effects on wave dispersion in three-phase acoustic metamaterials, J. Mech. Phys. Solids, 96, pp. 29-47 (2016).
  * [45] L. D’Alessandro, A. Krushynska, R. Ardito, N. M. Pugno, A. Corigliano, A design strategy to match the band gap of periodic and aperiodic metamaterials, Sci. Rep., 10 (1), art. no. 16403 (2020).


# Instantaneous Frequency Estimation in Unbalanced Systems Using Affine
Differential Geometry

Ali Alshawabkeh, Georgios Tzounas, , and Federico Milano A. Alshawabkeh, G.
Tzounas and F. Milano are with School of Electrical and Electronic
Engineering, University College Dublin, Dublin, D04V1W8, Ireland. E-mails:
<EMAIL_ADDRESS>{georgios.tzounas<EMAIL_ADDRESS>work is supported by the Sustainable Energy Authority of Ireland (SEAI) by
funding A. Alshawabkeh and F. Milano under project FRESLIPS, Grant No.
RDD/00681.

###### Abstract

The paper discusses the relationships between electrical quantities, namely
voltages and frequency, and affine differential geometry ones, namely affine
arc length and curvature. Moreover, it establishes a link between frequency
and time derivatives of voltage, through the utilization of affine
differential geometry invariants. Based on this link, a new instantaneous
frequency estimation formula is proposed, which is particularly suited for
unbalanced systems. An application of the proposed formula to single-phase
systems is also provided. Several numerical examples based on balanced,
unbalanced, as well as single-phase systems illustrate the findings of the
paper.

###### Index Terms:

Frequency estimation, affine differential geometry, instantaneous frequency,
unbalanced systems, curvature, phase-locked loop (PLL).

##  I Introduction

###  I-A Motivation

On-line frequency estimation in electric power systems is known to be a
complex problem, especially in conditions characterized by the presence of
harmonics and imbalances of amplitude or phase angle within the measured
signals. Aiming in particular to tackle the complexities posed by unbalanced
systems, this paper introduces a technique for on-line frequency estimation,
based on the principles of affine differential geometry.

###  I-B Literature Review

The problem of frequency estimation has been studied for many years and
several solution approaches have been reported in the literature, e.g. see [1,
2, 3, 4, 5, 6, 7, 8, 9, 10]. These approaches rely on a variety of methods,
including phase-locked loops (PLLs) [1, 2], discrete Fourier transform [3, 4],
Kalman filtering [5, 6], least-squares [7, 8], adaptive notch filters [9, 10],
etc.

For grid synchronization and control applications in particular, PLLs are a
popular solution due their performance characteristics and straightforward
structure and practical implementation. Three-phase PLLs, for example, are
widely utilized to provide on-line phase and frequency estimations in grid-
connected power converters. With this regard, a conventional PLL configuration
in three-phase system applications is the synchronous reference frame (SRF)
PLL, which relies on transforming input voltages to the $\rm dq$ synchronous
reference frame and on regulating the frame’s angular position so that either
the $\rm d$\- or $\rm q$-axis component is zero. The analogous of SRF-PLL for
single-phase systems is the quadrature signal generation (QSG)-based PLL.
Given a single-phase voltage signal, a QSG-PLL defines a second dimension
through a fictitious quadrature signal, required to enable the application of
the Park transform (and thus the formulation of $\rm dq$-axis voltage
components). The simplest approach to do that is using a transport delay of
$T/4$, where $T$ is the period of the fundamental frequency, e.g. see [11].

Other approaches are based on the inverse Park transform [12], the Hilbert
transform [13], and on second-order generalized integrators [14]. Although
they provide robust frequency estimations under balanced conditions, they are
also known to perform poorly for unbalanced systems, wherein they often result
in estimations with sinusoidal ripple errors [15, 16, 17, 18, 19]. Reducing
the bandwidth helps mitigate this issue and refine accuracy, but also
compromises dynamic performance [20]. In this regard, careful tuning of PLL
parameters is essential to achieve good trade-offs between dynamic performance
and estimation accuracy. Efforts to improve the performance of PLLs under
unbalanced conditions include, among other studies, [1, 21, 22, 17].

PLLs belong to the broad family of time-domain methods. In this paper we also
focus on time-domain methods but approach the problem of frequency estimation
from an unconventional perspective, that is based on the theory of
differential geometry. The starting assumption is that any voltage vector can
be perceived as the velocity of a point on a space curve and, as such, be
analyzed using differential geometrical invariants. In our recent work on the
topic, we described the definition of these curves in a Euclidean space and,
by applying the Frenet–Serret formulas, we derived a correspondence between
curvature and instantaneous electrical frequency [23, 24, 25, 26]. Despite
providing accurate frequency estimations for balanced systems, the curvature
obtained in these works is time-varying in stationary unbalanced operating
conditions, a result that clearly does not align well with the notion of
angular frequency of stationary ac signals.

In this paper, we aim at solving this issue through an alternative theory of
differential geometry of curves, namely through affine differential geometry.
This theory has found applications in various areas, such as control of
mechanical systems [27], computer vision [28], and motion identification [29].
But there has been, to the best of our knowledge, no application to power
system analysis or frequency estimation.

###  I-C Contributions

The specific contributions of the paper are as follows.

  * •

A derivation of the expressions for the affine arc length and affine curvature
in terms of the voltage of an ac system.

  * •

An approximated yet accurate formula of the instantaneous angular frequency of
a three-phase voltage as a function of affine geometrical invariants.

  * •

A demonstration of the effectiveness of the proposed formula to serve as an
instantaneous frequency estimation technique for unbalanced three-phase
systems, as well as for single-phase systems.

The last two points are fully supported through a variety of examples, which
are provided in the case study section. The examples show, in particular, that
the proposed expression yields a more precise estimation of the instantaneous
frequency in unbalanced systems, compared to PLLs and the Frenet-frame based
method from [23].

###  I-D Paper Organization

The remainder of the paper is organized as follows. Section II provides an
overview of basic concepts from affine geometry. These concepts are essential
for the derivation of the theoretical results of the paper presented in
Section III. Section IV tests the proposed approach through analytical
examples, as well as through a case study based on a fully-fledged EMT model
of the IEEE 39-bus system. Finally, Section V draws relevant conclusions.

##  II Outlines of Affine Differential Geometry

This section provides a brief overview of affine quantities that are relevant
for the derivation of the instantaneous frequency formula for ac systems
presented in Section III. The interested reader can find a comprehensive
presentation of the theory of affine differential geometry in [30].

Affine geometry can be defined as a Euclidean geometry without measuring
distances or angles. In other words, it is a Euclidean geometry whose metric
structure has been removed [29]. For the affine plane, a pair of non-collinear
vectors determines a parallelogram whose area is given by the determinant of
the two-by-two matrix formed by these vectors.

Let us consider a smooth parametric curve in the plane:

$\boldsymbol{x}(t)=x_{1}(t)\,\boldsymbol{e}_{\rm
1}+x_{2}(t)\,\boldsymbol{e}_{\rm 2}\,,$ (1)

where $x_{1}(t),x_{2}(t):\mathbb{R}\mapsto\mathbb{R}$ are smooth and
$\boldsymbol{e}_{\rm 1}$ and $\boldsymbol{e}_{\rm 2}$ form an orthogonal basis
of the plane. Let us also assume that the curve $\boldsymbol{x}$ does not have
inflection points, that is, the magnitude of the bracket operator

$[\dot{\boldsymbol{x}}(t),\ddot{\boldsymbol{x}}(t)]\neq 0,\quad\forall t\,,$
(2)

never vanishes. In (2), $\dot{\boldsymbol{x}}=d\boldsymbol{x}/dt$ and
$\ddot{\boldsymbol{x}}=d^{2}\boldsymbol{x}/dt^{2}$, and the bracket operator
of two vectors, say $[\boldsymbol{a},\boldsymbol{b}]$, with
$\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{2}$, is defined as:

$[\boldsymbol{a},\boldsymbol{b}]={\rm det}\left|\begin{matrix}a_{1}&b_{1}\\\
a_{2}&b_{2}\end{matrix}\right|=a_{1}b_{2}-b_{1}a_{2}\,,$ (3)

The affine arc length or equi-affine arc length, indicated with $\sigma$, is
defined as:

$\sigma(t)=\int_{t_{0}}^{t}[\dot{\boldsymbol{x}}(t),\ddot{\boldsymbol{x}}(t)]^{1/3}dt\,,$
(4)

or, equivalently:

$\dot{\sigma}(t)=\frac{d\sigma(t)}{dt}=[\dot{\boldsymbol{x}}(t),\ddot{\boldsymbol{x}}(t)]^{1/3}\,.$
(5)

A curve $\boldsymbol{x}$ is said to be parameterized with $\sigma$ if, for all
$\sigma$, it satisfies the condition:

$[\boldsymbol{x}^{\prime}(\sigma),\boldsymbol{x}^{\prime\prime}(\sigma)]=1\,,$
(6)

where $\boldsymbol{x}^{\prime}=d\boldsymbol{x}/d\sigma$ is the affine tangent
and $\boldsymbol{x}^{\prime\prime}=d^{2}\boldsymbol{x}/d\sigma^{2}$ is the
affine normal. Applying the chain rule, $\boldsymbol{x}^{\prime}$ becomes:

$\boldsymbol{x}^{\prime}(\sigma(t))=\frac{d\boldsymbol{x}}{d\sigma}=\frac{d\boldsymbol{x}}{dt}\frac{dt}{d\sigma}=\frac{\dot{\boldsymbol{x}}(t)}{[\dot{\boldsymbol{x}}(t),\ddot{\boldsymbol{x}}(t)]^{1/3}}\,,$
(7)

and, differentiating (6) with respect to $\sigma$, one obtains:

$[\boldsymbol{x}^{\prime}(\sigma),\boldsymbol{x}^{\prime\prime\prime}(\sigma)]=0\,.$
(8)

This result implies that $\boldsymbol{x}^{\prime}$ and
$\boldsymbol{x}^{\prime\prime\prime}$ are linearly independent, leading to the
relationship:

$\boldsymbol{x}^{\prime\prime\prime}(\sigma)=-\kappa_{a}(\sigma)\,\boldsymbol{x}^{\prime}(\sigma)\,,$
(9)

where $\kappa_{a}$ is the affine curvature or equi-affine curvature of the
curve $\boldsymbol{x}$ and is defined as:

$\kappa_{a}(\sigma)=[\boldsymbol{x}^{\prime\prime}(\sigma),\boldsymbol{x}^{\prime\prime\prime}(\sigma)]\,.$
(10)

The affine curvature is represented by the area of the parallelogram formed by
the vectors $\boldsymbol{x}^{\prime\prime}$ and
$\boldsymbol{x}^{\prime\prime\prime}$.

It is relevant for the following discussion on the estimation of the
instantaneous frequency to note that for non-singular conic sections,
$\kappa_{a}$ is constant, as follows [30, 31]:

  * •

for $\kappa_{a}=0$, the curve is a parabola;

  * •

for $\kappa_{a}>0$, the curve is an ellipse;

  * •

for $\kappa_{a}<0$, the curve is an hyperbola.

In the next section, we consider the specific case of the ellipse, that is
$\kappa_{a}>0$.

##  III Voltage in the Affine Plane

We adopt the assumption made in [25], that is, the magnetic flux
$\boldsymbol{\varphi}$ is the position of a point on a space curve in
generalized coordinates and, from Faraday’s law, the speed of such a point is
the voltage, as follows:

$\boldsymbol{\varphi}(t)\equiv-\boldsymbol{x}(t)\,\quad\Rightarrow\quad\boldsymbol{v}(t)=-\dot{\boldsymbol{\varphi}}(t)\equiv\dot{\boldsymbol{x}}(t)\,.$
(11)

Reference [23] shows that one can express electrical quantities such as
voltage and current in terms of the coordinates of the Frenet frame and of
geometric invariants such as arc length, curvature and torsion [23]. In the
same vein, but using the definitions of coordinates, arc length and curvature
given by affine differential geometry, this section derives a new formula for
the instantaneous frequency of electrical quantities. In the remainder of this
section, we discuss exclusively voltages, but the same procedure and results
can be obtained using currents. We consider two scenarios, namely unbalanced
three-phase systems; and single-phase systems.

###  III-A Three-Phase Unbalanced Voltages

Let’s assume that the phases $\rm abc$ of a three-phase voltage
$\boldsymbol{v}(t)$ constitute a set of orthogonal coordinates:

$\boldsymbol{v}(t)=v_{\rm a}(t)\,\boldsymbol{e}_{\rm a}+v_{\rm
b}(t)\,\boldsymbol{e}_{\rm b}+v_{\rm c}(t)\,\boldsymbol{e}_{\rm c}\,.$ (12)

In order to employ the theory described in the previous section, which applies
to curves in two dimensions, we first need to transform the voltage
$\boldsymbol{v}(t)$ into the following shape:

$\boldsymbol{v}(t)=v_{1}(t)\,\boldsymbol{e}_{\rm
1}+v_{2}(t)\,\boldsymbol{e}_{\rm 2}\,.$ (13)

This is conveniently achieved by applying the Clarke transform to (12) and
taking the $\alpha$ and $\beta$ components, as follows:

$\begin{bmatrix}v_{\alpha}(t)\\\
v_{\beta}(t)\end{bmatrix}=\frac{2}{3}\begin{bmatrix}1&-\frac{1}{2}&-\frac{1}{2}\\\
0&\frac{\sqrt{3}}{2}&-\frac{\sqrt{3}}{2}\\\ \end{bmatrix}\begin{bmatrix}v_{\rm
a}(t)\\\ v_{\rm b}(t)\\\ v_{\rm c}(t)\end{bmatrix}.$ (14)

Thus, the components of the voltage in (13) are:

$v_{1}(t)=v_{\alpha}(t)\,,\qquad v_{2}(t)=v_{\beta}(t)\,.$ (15)

####  III-A1 Stationary Sinusoidal Voltages

We discuss in this section a “base case” scenario for which the theory of
affine differential geometry allows obtaining the exact value of the frequency
of the voltage. This is the case of an unbalanced stationary sinusoidal (i.e.,
without harmonics) three-phase voltage. As discussed above, since we apply the
Clarke transform, the components of the voltage vector in (13) are:

$v_{1}(t)=V_{1}\cos\theta(t)\,,\qquad v_{2}(t)=V_{2}\sin\theta(t)\,,$ (16)

where $V_{1}$ and $V_{2}$ are constant with $V_{1}\neq V_{2}$, and:

$\theta(t)=\omega_{o}t+\theta_{o}\,,$ (17)

where $\omega_{o}$ is the fundamental synchronous reference frequency of the
system and $\theta_{o}$ is constant and depends on the chosen phase angle
reference of the system.

With the equivalence given in (11), the time derivative of the affine arc
length $\dot{\sigma}$ in (5) can be written as:

$\dot{\sigma}=[\boldsymbol{v}(t),\dot{\boldsymbol{v}}(t)]^{1/3}=(\omega_{o}V_{1}V_{2})^{1/3}\,.$
(18)

Note that while $\boldsymbol{v}$ and $\dot{\boldsymbol{v}}$ depend on time,
$\dot{\sigma}$ does not. Then, imposing that the components of the voltage are
as those given in (16), one obtains:

$\boldsymbol{x}^{\prime}(t)=\frac{\boldsymbol{v}(t)}{\dot{\sigma}}\,,\quad\boldsymbol{x}^{\prime\prime}(t)=\frac{\dot{\boldsymbol{v}}(t)}{\dot{\sigma}^{2}}\,,\quad\boldsymbol{x}^{\prime\prime\prime}(t)=\frac{\ddot{\boldsymbol{v}}(t)}{\dot{\sigma}^{3}}\,,$
(19)

where

$\displaystyle\boldsymbol{v}(t)$
$\displaystyle=V_{1}\cos\theta(t)\,\boldsymbol{e}_{\rm
1}+V_{2}\sin\theta(t)\,\boldsymbol{e}_{\rm 2}\,,$ (20)
$\displaystyle\dot{\boldsymbol{v}}(t)$
$\displaystyle=-\omega_{o}V_{1}\sin\theta(t)\,\boldsymbol{e}_{\rm
1}+\omega_{o}V_{2}\cos\theta(t)\,\boldsymbol{e}_{\rm 2}\,,$
$\displaystyle\ddot{\boldsymbol{v}}(t)$
$\displaystyle=-\omega^{2}_{o}V_{1}\cos\theta(t)\,\boldsymbol{e}_{\rm
1}-\omega^{2}_{o}V_{2}\sin\theta(t)\,\boldsymbol{e}_{\rm 2}\,.$

Then, using (10), (18) and (19), the expression of the affine curvature
$\kappa_{a}$ becomes:

$\kappa_{a}=\frac{1}{\dot{\sigma}^{5}}[\dot{\boldsymbol{v}}(t),\ddot{\boldsymbol{v}}(t)]=\frac{\omega_{o}^{3}V_{1}V_{2}}{\dot{\sigma}^{5}}\,,$
(21)

where $\kappa_{a}$ is constant, which is as expected since (16) describes an
ellipse in the plane $(v_{1},v_{2})$. Merging (18) and (21) we obtain:

$\omega_{o}=\sqrt{\kappa_{a}}\,\dot{\sigma}=\sqrt{\frac{[\dot{\boldsymbol{v}}(t),\ddot{\boldsymbol{v}}(t)]}{[\boldsymbol{v}(t),\dot{\boldsymbol{v}}(t)]}}\,.$
(22)

Equation (22) indicates that, in order to calculate the angular frequency of
the voltage in unbalanced conditions, it suffices to measure $\boldsymbol{v}$
and estimate its first and second time derivatives.

####  III-A2 Transient Voltages

Section III-A1 considers an ideal scenario for which the magnitude and the
angular frequency of a three-phase voltage are constant. This scenario leads
to a compact and elegant analytical result. However, such scenario is hardly
found in practice, wherein the presence of noise, harmonics, and transient
conditions prevents obtaining a general explicit expression for the
instantaneous frequency.

Under certain conditions, however, it is still possible to utilize the results
of Section III-A1 for a voltage of time-varying angular frequency and/or
magnitudes of $v_{1}$ and $v_{2}$. Consider a time-varying voltage vector:

$\boldsymbol{v}(t)=V_{1}(t)\cos\vartheta(t)\,\boldsymbol{e}_{\rm
1}+V_{2}(t)\sin\vartheta(t)\,\boldsymbol{e}_{\rm 2}\,,$ (23)

where $\vartheta(t)=\omega_{o}t+\phi(t)$. The conditions so that (22) holds
for a voltage $\boldsymbol{v}(t)$ in the form of (23) are:

$\displaystyle\frac{d^{h}}{dt^{h}}\phi(t)$
$\displaystyle\ll\omega^{h}_{o}\,,\quad h=1,2\,,$ (24)
$\displaystyle\frac{d^{h}}{dt^{h}}\frac{V_{i}(t)}{\langle V_{i}\rangle}$
$\displaystyle\ll\omega^{h}_{o}\,,\quad i,h=1,2\,,$ (25)

where $\langle\cdot\rangle$ denotes the average value. Condition (24) for
$h=1$ indicates that the variation of the instantaneous frequency of the
voltage is close to the synchronous reference angular frequency of the system;
for $h=2$, (24) imposes a boundary to the rate of change of frequency (RoCoF);
and conditions (25) impose that the variations of the radial frequency (see
definition in [25]) are small compared to the fundamental frequency of the
grid. All these assumptions are generally well satisfied in power systems.

Conditions (24) and (25) are sufficient for (22) to hold at least as a first
order approximation. In fact, the first time derivative of the voltage vector
in (23) is:

$\dot{\boldsymbol{v}}=(\dot{V}_{1}\cos\vartheta-\dot{\vartheta}V_{1}\sin\vartheta)\,\boldsymbol{e}_{\rm
1}+(\dot{V}_{2}\sin\vartheta+\dot{\vartheta}V_{2}\cos\vartheta)\,\boldsymbol{e}_{\rm
2}\,,$

and the second time derivative is:

$\displaystyle\ddot{\boldsymbol{v}}=$
$\displaystyle-\big{(}\ddot{\vartheta}V_{1}\sin\vartheta+\dot{\vartheta}^{2}V_{1}\cos\vartheta+\dot{\vartheta}\dot{V}_{1}\sin\vartheta-\ddot{V}_{1}\cos\vartheta\big{)}\,\boldsymbol{e}_{\rm
1}$
$\displaystyle+\big{(}\ddot{\vartheta}V_{2}\cos\vartheta+\dot{\vartheta}^{2}V_{2}\sin\vartheta+\dot{\vartheta}\dot{V}_{2}\cos\vartheta-\ddot{V}_{2}\sin\vartheta\big{)}\,\boldsymbol{e}_{\rm
2}\,,$

where the dependency on time has been omitted for economy of notation.

It is straightforward to show that by applying (24) and (25), the voltage
derivatives can be approximated with the second and third equations of (20)
and, hence, the instantaneous frequency can be approximated using (22). In
summary, (24) and (25) lead to the following approximated expression of the
instantaneous frequency of a time-varying unbalanced voltage:

$\dot{\vartheta}(t)\approx\boxed{\omega_{a}(t)=\sqrt{\frac{[\dot{\boldsymbol{v}}(t),\ddot{\boldsymbol{v}}(t)]}{[\boldsymbol{v}(t),\dot{\boldsymbol{v}}(t)]}}}$
(26)

The expression of $\omega_{a}$ in (26) is the main result of this work. We
test the accuracy of (26) through a variety of examples and a case study in
Section IV.

###  III-B Single-Phase Voltages

In this section we consider a single-phase voltage with instantaneous value
$v(t)$. To apply the theory described in Section II, we first need to
transform $v(t)$ into the shape of (13). To this aim, we construct the second
dimension by employing the voltage derivative. That is:

$v_{1}(t)=v(t)\,,\qquad v_{2}(t)=\dot{v}(t)\,.$ (27)

Since the time derivative of sinusoidal signals gives a $90^{\circ}$ rotation,
using the time derivative is equivalent to defining a quadrature axis (see
also discussion on quadrature signal generation in Section I).

####  III-B1 Stationary Sinusoidal Voltages

The result obtained in the previous section can be easily extended to a
stationary sinusoidal single-phase voltage using (27). Let the voltage be:

$v(t)=V\cos\theta(t)\,,$ (28)

where $V$ is constant and $\theta$ is defined in (17). Then, from (27), the
components of the voltage vector are:

$v_{1}(t)=V\cos\theta(t)\,,\qquad v_{2}(t)=-\omega_{o}V\sin\theta(t)\,.$ (29)

Substituting $V_{1}=V$ and $V_{2}=\omega_{o}V$ in (18) and (21), one obtains:

$\dot{\sigma}=(\omega_{o}V)^{2/3}\,,\quad\kappa_{a}=\frac{\omega_{o}^{4}V^{2}}{\dot{\sigma}^{5}}\,.$
(30)

Apart from the fact that calculation of $\ddot{\boldsymbol{v}}(t)$ in this
case requires the additional step of computing the third derivative of $v(t)$,
equation (22) holds and allows estimating the angular frequency also for a
single-phase voltage.

####  III-B2 Transient Voltages

Consider a time-varying single-phase voltage:

$v(t)=V(t)\cos\vartheta(t)\,,$ (31)

where $\vartheta(t)=\omega_{o}t+\phi(t)$. The voltage vector is defined as:

$\displaystyle\boldsymbol{v}(t)=$ $\displaystyle
V(t)\cos\vartheta(t)\,\boldsymbol{e}_{\rm 1}$ (32) $\displaystyle+$
$\displaystyle[\dot{V}(t)\cos\vartheta(t)-V(t)\dot{\vartheta}(t)\sin\vartheta(t)]\,\boldsymbol{e}_{\rm
2}\,.$

If one assumes:

$\displaystyle\frac{d^{h}}{dt^{h}}\phi(t)$
$\displaystyle\ll\omega^{h}_{o}\,,\quad h=1,2,3\,,$ (33)
$\displaystyle\frac{d^{h}}{dt^{h}}\frac{V(t)}{\langle V_{i}\rangle}$
$\displaystyle\ll\omega^{h}_{o}\,,\quad h=1,2,3\,,$ (34)

then (26) is also a good approximation of the instantaneous frequency of the
time-varying single-phase voltage in (31).

##  IV Case Studies

The examples discussed in this section are aimed at illustrating the accuracy
of (26) for both unbalanced three- and single-phase voltages in various non-
sinusoidal and non-stationary conditions. The proposed approach is compared
with a conventional technique, namely the SRF-PLL, as well as with the Frenet
frame-based method proposed in [23]. The first section of the study focuses on
three-phase systems under various conditions. Both balanced and unbalanced
cases are discussed. The IEEE 39-bus system is also utilized to illustrate the
accuracy of the proposed method under transient unbalanced conditions.
Finally, an example discussing the frequency estimation for a non-stationary
single-phase voltage is presented.

In all figures shown in this section, $\omega_{a}$ represents the estimated
angular frequency derived from the proposed approach, while $\omega_{\kappa}$
represents the estimated angular frequency obtained using the Frenet frame-
based method from [23], as follows:

$\omega_{\kappa}(t)=\frac{[\boldsymbol{v}(t),\dot{\boldsymbol{v}}(t)]}{|\boldsymbol{v}(t)|^{2}}\,.$
(35)

The three-phase voltage trajectories are given in per unit (pu) with respect
to a base of 12 kV, and $\omega_{a}$ and $\omega_{\kappa}$ are in pu with
respect to $\omega_{o}=100\pi$ rad/s. Finally, note that in all examples, (26)
is calculated using a sampling of the voltage signal, transforming it through
the Clarke transform and then evaluating numerically the time derivatives of
the $\alpha$ and $\beta$ components.

###  IV-A Three-Phase Voltage

Let us consider the three-phase voltage vector given in (12), which we repeat
here for convenience:

$\boldsymbol{v}(t)=v_{\rm a}(t)\,\boldsymbol{e}_{\rm a}+v_{\rm
b}(t)\,\boldsymbol{e}_{\rm b}+v_{\rm c}(t)\,\boldsymbol{e}_{\rm c}\,,$ (36)

with components:

$\displaystyle v_{\rm a}(t)$ $\displaystyle=V_{\rm
a}\sin(\omega_{o}t+\phi_{\rm a}(t))\,,$ (37) $\displaystyle v_{\rm b}(t)$
$\displaystyle=V_{\rm b}\sin(\omega_{o}t+\phi_{\rm b}(t)-\zeta_{\rm b})\,,$
$\displaystyle v_{\rm c}(t)$ $\displaystyle=V_{\rm
c}\sin(\omega_{o}t+\phi_{\rm c}(t)+\zeta_{\rm c})\,.$

Recall that the proposed approach operates in two dimensions, and that we use
Clarke transform to convert (12) to the two-dimensional $(\alpha,\beta)$
plane.

####  IV-A1 Balanced Three-Phase Voltage

We discuss two examples: the first example involves a stationary condition,
whereas the second example considers a signal with time-varying voltage
magnitude. In both examples, the angular frequency is constant and equal to
$\omega_{o}$. The parameters used are:

  * •

E1: $V_{i}=12$ kV, with $\omega_{o}=100\pi$ rad/s, $\phi_{i}=0$ and
$\zeta_{\rm b}=\zeta_{\rm c}=\frac{2\pi}{3}$ rad.

  * •

E2: $V_{i}=12+3\sin(\pi t)$ kV, with $\omega_{o}=100\pi$ rad/s, $\phi_{i}=0$
and $\zeta_{\rm b}=\zeta_{\rm c}=\frac{2\pi}{3}$ rad.

Figure 1 shows the phase voltages, the geometric frequencies $\omega_{a}$ and
$\omega_{\kappa}$, as well as the instantaneous frequency
$\omega_{\scriptscriptstyle{\rm PLL}}$ obtained with a conventional SRF-PLL
for E1 and E2. As expected, since the voltage is balanced and the curve in the
plane $(\alpha,\beta)$ is a circle, there is a perfect match between the
estimations obtained with the two differential geometry-based methods. Note,
however, that the two approaches return the right results for two different
reasons: the Frenet frame-based formula returns a constant $\omega_{\kappa}$
because the circle has a constant curvature; whereas the proposed affine
differential geometry approach returns a constant $\omega_{a}$ because the
circle is a special case of an ellipse. We note that the conventional PLL also
works well in this balanced-voltage case.

(a) E1: Voltage components

(b) E1: Geometric frequency

(c) E2: Voltage components

(d) E2: Geometric frequency

Figure 1: Balanced three-phase voltages and estimated frequency.

####  IV-A2 Unbalanced Three-Phase Voltage

Three examples of unbalanced voltages with constant angular frequency
$\omega_{o}$ are considered in this section. The first example involves
unequal constant voltage magnitudes; the second example examines a system with
unequal and time-varying voltage magnitudes; and the third example examines a
system with unequal phase displacements. The following three cases are
considered:

  * •

E3: $V_{\rm a}=V_{\rm c}=12$ kV, $V_{\rm b}=8$ kV, with $\omega_{o}=100\pi$
rad/s, $\phi_{i}=0$ and $\zeta_{\rm b}=\zeta_{\rm c}=\frac{2\pi}{3}$ rad.

  * •

E4: $V_{\rm a}=V_{\rm c}=12+3\sin(\pi t)$ kV, $V_{\rm b}=8+2\sin(2\pi t)$ kV,
with $\omega_{o}=100\pi$ rad/s, $\phi_{i}=0$ and $\zeta_{\rm b}=\zeta_{\rm
c}=\frac{2\pi}{3}$ rad.

  * •

E5: $V_{\rm a}=V_{\rm b}=V_{\rm c}=12$ kV, with $\omega_{o}=100\pi$ rad/s,
$\phi_{i}=0$ and $\zeta_{\rm b}=\frac{-2\pi}{3},\zeta_{\rm
c}=\frac{1.5\pi}{3}$ rad.

Figure 2 shows the voltage components and estimated geometric and PLL
frequencies for examples E3-E5. In all these examples, the curves in the
$(\alpha,\beta)$ plane are ellipses. This means that the curvature obtained
using the Frenet frame is time-varying and periodic, thus leading to a time-
varying and periodic $\omega_{\kappa}$. Moreover, the conventional PLL also
outputs a time-varying frequency in the form of a significant ripple around
the frequency $\omega_{o}$. On the other hand, the proposed affine geometry
formula returns a constant $\omega_{a}$ equal to $\omega_{o}$ (in pu), which
is consistent with the expected in this case result.

(a) E3: Voltage components

(b) E3: Geometric frequency

(c) E4: Voltage components

(d) E4: Geometric frequency

(e) E5: Voltage components

(f) E5: Geometric frequency

Figure 2: Unbalanced three-phase voltages and estimated frequencies.

###  IV-B Three-Phase Voltage with Time-Varying Frequency

Two examples of three-phase system with varying angular frequency are
considered in this section. The first example considers a voltage with angular
frequency that varies periodically around its average value. This example is
used to resemble the transient behavior of the voltage following a contingency
in a power system, where voltage phase angle oscillations arising due to
electro-mechanical swings of synchronous machines are poorly damped and thus
sustain for a relatively long time. The second example is an extreme and
uncommon situation in power systems, where the components of the three-phase
voltage are time-varying and have unequal angular frequencies. The following
parameters are used:

  * •

E6: $V_{i}=12$ kV, with $\omega_{o}=100\pi$ rad/s and $\zeta_{\rm
b}=\zeta_{\rm c}=\frac{2\pi}{3}$ rad and $\phi_{i}(t)=\pi\sin(0.4\pi t)$ rad.

  * •

E7: $V_{i}=12$ kV, with $\omega_{o}=100\pi$ rad/s and $\zeta_{\rm
b}=\zeta_{rmc}=\frac{2\pi}{3}$ rad and $\phi_{\rm a}(t)=\phi_{\rm
b}(t)=\pi\sin(0.4\pi t),\phi_{\rm c}(t)=1.1\pi\sin(0.4\pi t)$ rad.

Figure 3(a) shows the estimated frequency with the proposed formula, the
frequency estimated with the PLL, and the geometrical frequency
$\omega_{\kappa}$, for E6. Despite the approximations imposed by assuming (24)
and (25), we note that (26) is able to precisely track the exact instantaneous
frequency. In this example, also $\omega_{\kappa}$ and
$\omega_{\scriptscriptstyle{\rm PLL}}$ track well $\rm IF$. On the other hand,
for E7, while $\omega_{a}$ and $\omega_{\scriptscriptstyle{\rm PLL}}$ still
track well the exact frequency, $\omega_{\kappa}$ shows significant
fluctuations. This behavior is illustrated in Fig. 3(b).

(a) E6

(b) E7

Figure 3: Estimated angular frequency.

###  IV-C Stationary Three-Phase Voltage with Harmonics

As a last example on three-phase voltage, we discuss the effect of harmonics
on the estimation of the frequency based on (26). A fundamental condition for
the affine differential geometry approach to work properly is that (2) is
satisfied at all times. Harmonics, however, introduce inflection points, that
is, points for which $[\dot{\boldsymbol{x}},\ddot{\boldsymbol{x}}]\leq 0$.
Moreover, since the term $[\dot{\boldsymbol{x}},\ddot{\boldsymbol{x}}]$
appears in the denominator of (26), this leads to numerical issues. Figure
4(b) shows the performance of (26) as well as of the PLL and Frenet-based
estimated frequencies for a stationary balanced three-phase voltage. For the
fundamental frequency, $V_{i}=12$ kV, with $\omega_{o}=100\pi$ rad/s and
$\zeta_{\rm b}=\zeta_{\rm c}=\frac{2\pi}{3}$ rad and $\phi_{i}(t)=0$ are
assumed. Then, 7-th and 11-th harmonics are added with magnitudes $0.02V_{i}$
and $0.01V_{i}$, respectively. To overcome numerical issues, we have set
$\omega_{a}=0$ if $[\dot{\boldsymbol{x}},\ddot{\boldsymbol{x}}]\leq 0$. As
expected, in this scenario, $\omega_{a}$ shows a poor performance. The best
estimation is obtained with the PLL.

(a) Voltage components

(b) Estimated angular frequency

Figure 4: E8: Three-phase voltage with harmonics and estimated angular
frequency.

###  IV-D IEEE 39-Bus System in Unbalanced Conditions

In this section, the accuracy of (26) is tested using a fully-fledged EMT
model of the IEEE 39-bus system. The system setup is the same as the
unbalanced scenario described in [24], where it was utilized to show the
performance of the Frenet frame-based frequency estimation. In particular, the
power consumption of all 19 loads of the system is unbalanced, with imbalances
ranging from 5 to 10% on one of the phases. A three-phase fault is simulated
at bus 4. The fault occurs at $t=0.2$ s and is cleared at $t=0.3$ s. The
behavior of the three-phase voltage at bus 26 following the contingency is
illustrated in Fig. 5.

Figure 5: Voltage at bus 26 of IEEE 39-bus system following a three-phase
fault at bus 4.

Figure 6 shows the results of the frequency estimation for the unbalanced
voltage at bus 26. The $\omega_{a}$ is more accurate than
$\omega_{\scriptscriptstyle{\rm PLL}}$. Note that (26) has been evaluated
using numerical differentiation of the time series of the three-phase voltage
components that have a constant time sampling rate of $10$ ms. The numerical
derivatives are filtered using a second order Butterworth digital filter and a
IIR filter. Note that, in this scenario, $\omega_{\kappa}$ shows a bigger
ripple than $\omega_{\scriptscriptstyle{\rm PLL}}$ and, thus $\omega_{\kappa}$
is omitted in Fig. 6 for clarity. The interested reader can find a
comprehensive comparison between PLL frequency estimations and
$\omega_{\kappa}$ in [24].

Figure 6: Estimated frequency for IEEE 39-bus system for unbalanced time-
varying voltage.

###  IV-E Single-Phase Voltage

This last example illustrates the performance of the proposed formula (26)
when applied to a single-phase voltage with time-varying angular frequency
$\omega t+\phi(t)$ and constant amplitude $V$:

$v(t)=V\sin(\omega_{o}t+\phi(t))\,.$ (38)

The parameters considered for this example are: $V=12$ kV, $\omega_{o}=100\pi$
rad/s and $\phi(t)=0.05\omega_{o}e^{-t}(1-\cos(\pi t))$ rad.

As discussed in Section III-B, we construct the second dimension by using the
derivative of the original signal, as in (27). Figure 7 illustrates the
accuracy of (26) in matching the actual analytical value of the instantaneous
frequency, that is, $\rm IF=\omega_{o}+\dot{\phi}$.

Figure 7: Estimated frequency, analytical instantaneous frequency ($\rm IF$)
and frequency estimated using a phase shift for the single-phase voltage.

Figure 7 also shows the frequency estimated using a conventional PLL where the
quadrature signal is obtained using $v(t-\tau)$, where the transport delay is
$\tau=0.25\,T=0.5\,\pi/\omega_{o}$. Despite the approximations resulting from
the assumptions (33) and (34), also in this case the proposed approach shows
very good accuracy, whereas the PLL shows some ripples due to the fact that
the quadrature signal is not exact because the frequency is time-varying.

##  V Conclusions

This paper presents an approach based on affine differential geometry to
estimate the angular frequency of unbalanced three-phase voltages as well as
of single-phase voltages. The main contribution of this work is the
approximated formula (26), that estimates the angular frequency of the voltage
through measurements of the voltage, and calculation of the first and second
derivatives of the voltage itself.

Approximations based on the nature of typical power system transients are
assumed in order to achieve a compact explicit expression of the proposed
angular frequency estimation formula. Then, a variety of examples are provided
to demonstrate the adequateness of such approximations and the performance of
the proposed formula. When compared to PLLs as well as to the Frenet frame-
based estimation from [23], the proposed formula proves to be accurate and
robust in balanced and unbalanced conditions, as well as for voltages of time-
varying magnitude and frequency.

Future work will focus on testing the proposed formula with real-world
measurements and on extending its formulation to process signals with harmonic
contents as well as with multi-phase systems with a number of phases higher
than three.

## References

  * [1] B. Liu, F. Zhuo, Y. Zhu, H. Yi, and F. Wang, “A three-phase PLL algorithm based on signal reforming under distorted grid conditions,” _IEEE Transactions on Power Electronics_ , vol. 30, no. 9, pp. 5272–5283, 2014.
  * [2] M. Eskandari and A. V. Savkin, “Robust PLL synchronization unit for grid-feeding converters in micro/weak grids,” _IEEE Transactions on Industrial Informatics_ , vol. 19, no. 4, pp. 5400–5411, 2022.
  * [3] P. Romano, M. Paolone, J. Arnold, and R. Piacentini, “An interpolated-DFT synchrophasor estimation algorithm and its implementation in an FPGA-based PMU prototype,” in _IEEE PES General Meeting_. IEEE, 2013, pp. 1–6.
  * [4] J. Song, A. Mingotti, J. Zhang, L. Peretto, and H. Wen, “Fast iterative-interpolated dft phasor estimator considering out-of-band interference,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 71, pp. 1–14, 2022.
  * [5] S. Reza, M. Ciobotaru, and V. G. Agelidis, “Accurate estimation of single-phase grid voltage fundamental amplitude and frequency by using a frequency adaptive linear Kalman filter,” _IEEE Journal of Emerging and Selected Topics in Power Electronics_ , vol. 4, no. 4, pp. 1226–1235, 2016\.
  * [6] X. Nie, “Detection of grid voltage fundamental and harmonic components using Kalman filter based on dynamic tracking model,” _IEEE Transactions on Industrial Electronics_ , vol. 67, no. 2, pp. 1191–1200, 2019.
  * [7] S. Giarnetti, F. Leccese, and M. Caciotta, “Non recursive multi-harmonic least squares fitting for grid frequency estimation,” _Measurement_ , vol. 66, pp. 229–237, 2015.
  * [8] A. Pradhan, A. Routray, and A. Basak, “Power system frequency estimation using least mean square technique,” _IEEE Transactions on Power Delivery_ , vol. 20, no. 3, pp. 1812–1816, 2005.
  * [9] F. Wilches-Bernal, J. Wold, and W. H. Balliet, “A method for correcting frequency estimates for synthetic inertia control,” _IEEE Access_ , vol. 8, pp. 229 141–229 151, 2020.
  * [10] M. Mojiri, M. Karimi-Ghartemani, and A. Bakhshai, “Estimation of power system frequency using an adaptive notch filter,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 56, no. 6, pp. 2470–2477, 2007.
  * [11] L. Hadjidemetriou, Y. Yang, E. Kyriakides, and F. Blaabjerg, “A synchronization scheme for single-phase grid-tied inverters under harmonic distortion and grid disturbances,” _IEEE Transactions on Power Electronics_ , vol. 32, no. 4, pp. 2784–2793, 2016.
  * [12] R. M. Santos Filho, P. F. Seixas, P. C. Cortizo, L. A. Torres, and A. F. Souza, “Comparison of three single-phase PLL algorithms for ups applications,” _IEEE Transactions on Industrial Electronics_ , vol. 55, no. 8, pp. 2923–2932, 2008.
  * [13] P. Hao, W. Zanji, and C. Jianye, “A measuring method of the single-phase ac frequency, phase, and reactive power based on the hilbert filtering,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 56, no. 3, pp. 918–923, 2007.
  * [14] A. Sahoo, J. Ravishankar, and C. Jones, “Phase-locked loop independent second-order generalized integrator for single-phase grid synchronization,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 70, pp. 1–9, 2021.
  * [15] M. S. Reza, F. Sadeque, M. M. Hossain, A. M. Ghias, and V. G. Agelidis, “Three-phase PLL for grid-connected power converters under both amplitude and phase unbalanced conditions,” _IEEE Transactions on Industrial Electronics_ , vol. 66, no. 11, pp. 8881–8891, 2019.
  * [16] H. Karimi, M. Karimi-Ghartemani, and M. R. Iravani, “Estimation of frequency and its rate of change for applications in power systems,” _IEEE Transactions on Power Delivery_ , vol. 19, no. 2, pp. 472–480, 2004.
  * [17] M. E. Meral, “Improved phase-locked loop for robust and fast tracking of three phases under unbalanced electric grid conditions,” _IET Generation, Transmission & Distribution_, vol. 6, no. 2, pp. 152–160, 2012.
  * [18] A. Kulkarni and V. John, “Design of synchronous reference frame phase-locked loop with the presence of dc offsets in the input voltage,” _IET Power Electronics_ , vol. 8, no. 12, pp. 2435–2443, 2015.
  * [19] S. Golestan, M. Monfared, and F. D. Freijedo, “Design-oriented study of advanced synchronous reference frame phase-locked loops,” _IEEE Transactions on Power Electronics_ , vol. 28, no. 2, pp. 765–778, 2012.
  * [20] G. Escobar, C. N. M. Ho, S. Pettersson, M. J. Lopez-Sanchez, and A. A. Valdez-Fernandez, “Cascade three-phase PLL for unbalance and harmonic distortion operation (CSRF-PLL),” in _IECON_. IEEE, 2014, pp. 5489–5493.
  * [21] P. Rodríguez, J. Pou, J. Bergas, J. I. Candela, R. P. Burgos, and D. Boroyevich, “Decoupled double synchronous reference frame PLL for power converters control,” _IEEE Transactions on Power Electronics_ , vol. 22, no. 2, pp. 584–592, 2007.
  * [22] P. Rodríguez, A. Luna, I. Candela, R. Mujal, R. Teodorescu, and F. Blaabjerg, “Multiresonant frequency-locked loop for grid synchronization of power converters under distorted grid conditions,” _IEEE Transactions on Industrial Electronics_ , vol. 58, no. 1, pp. 127–138, 2010.
  * [23] F. Milano, G. Tzounas, I. Dassios, and T. Kerci, “Applications of the Frenet frame to electric circuits,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 69, no. 4, pp. 1668–1680, 2021.
  * [24] F. Milano, G. Tzounas, I. Dassios, M. A. A. Murad, and T. Kërçi, “Using differential geometry to revisit the paradoxes of the instantaneous frequency,” _IEEE Open Access Journal of Power and Energy_ , 2022.
  * [25] F. Milano, “A geometrical interpretation of frequency,” _IEEE Transactions on Power Systems_ , vol. 37, no. 1, pp. 816–819, 2021.
  * [26] F. Milano, “The Frenet frame as a generalization of the Park transform,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 70, no. 2, pp. 966–976, 2022.
  * [27] A. D. Lewis, “The bountiful intersection of differential geometry, mechanics, and control theory,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 1, pp. 135–158, 2018.
  * [28] M. Craizer, T. Lewiner, and J.-M. Morvan, “Parabolic polygons and discrete affine geometry,” in _Brazilian Symposium on Computer Graphics and Image Processing_. IEEE, 2006, pp. 19–26.
  * [29] T. Flash and A. A. Handzel, “Affine differential geometry analysis of human arm movements,” _Biological Cybernetics_ , vol. 96, pp. 577–601, 2007.
  * [30] K. Nomizu and T. Sasaki, _Affine differential geometry: geometry of affine immersions_. Cambridge University Press, 1994.
  * [31] E. Calabi, P. J. Olver, and A. Tannenbaum, “Affine geometry, curve flows, and invariant numerical approximations,” _Advances in Mathematics_ , vol. 124, no. 1, pp. 154–196, 1996.

Ali Alshawabkeh received the B.Sc. in Energy Engineering and Electric Power
from German Jordanian University in 2017 and the M.Sc. in Power and Data
Engineering from Offenburg University, Germany, in 2021. In 2022, he joined
the PhD program in Electrical Power Engineering of University College Dublin,
Ireland. He is funded by Sustainable Energy Authority of Ireland (SEAI) under
project FRESLIPS, Grant No. RDD/00681.
---
|  Georgios Tzounas (M’21) received the Diploma (M.E.) in Electrical and
Computer Engineering from the National Technical Univ. of Athens, Greece, in
2017, and the Ph.D. from University College Dublin (UCD), Ireland, in 2021. In
Jan.-Apr. 2020, he was a visiting researcher at Northeastern Univ., Boston,
MA. From Oct. 2020 to Apr. 2023, he was a postdoctoral researcher with UCD
(2020-2022) and ETH Zürich (2022-2023). Since Apr. 2023, he has been an
Assistant Professor with the School of Electrical and Electronic Engineering
at UCD. His research interests include stability analysis and control of power
systems.
---|---
|  Federico Milano (F’16) received from the University of Genoa, Italy, the ME
and Ph.D. in Electrical Engineering in 1999 and 2003, respectively. In 2013,
he joined University College Dublin, Ireland, where he is a full professor. He
is an IEEE Distinguished Lecturer, a senior editor of the IEEE Transactions on
Power Systems, an IET Fellow, and an editor in chief of the IET Generation,
Transmission & Distribution. He is chair of the IEEE Power System Stability
Controls Subcommittee; chair of the Technical Program Committee of PSCC 2024;
and a member of the Cigre Irish National Committee. His research interests
include power system modeling, control and stability analysis.
---|---