**Modelling, Simulation and Data Analysis in Acoustical Problems**

Special Issue Editors

**Claudio Guarnaccia Lamberto Tronchin Massimo Viscardi**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editors* Claudio Guarnaccia University of Salerno Italy

Lamberto Tronchin University of Bologna Italy

Massimo Viscardi University of Naples Federico II Italy

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/MSDAAP).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03928-284-5 (Pbk) ISBN 978-3-03928-285-2 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

### **Contents**


#### **Fang Wang, Yong Chen and Jianwei Wan**



**Tsukasa Ito, Toshinori Kubota, Takatoshi Furukawa, Hirooki Matsui, Kazunori Futai, Melinda Hull and Seiji Kakehata**

The Role of Powered Surgical Instruments in Ear Surgery: An Acoustical Blessing or a Curse? Reprinted from: *Applied Sciences* **2019**, *9*, 765, doi:10.3390/app9040765 ................ **555**

### **About the Special Issue Editors**

**Claudio Guarnaccia** is an associate professor of Applied Physics at the Department of Civil Engineering, University of Salerno. He has an M.Sc (2004) and Ph.D. (2008) in Physics from the Physics Department "E.R.Caianiello" of the University of Salerno. He has carried out post-doc research studies both in the Physics Department and in the Industrial Engineering Department of the University of Salerno. His research activity is focused on the application of physics to engineering and biomedical problems. The main focus of his research is acoustics. He works with several national and international universities and research institutes. He is a member of several national and international organizations, in the field of Acoustics and Applied Physics. He is an editor-in-chief, associate editor, member of the editorial board, and reviewer for several international journals, indexed in the most important international databases. He is the author of more than 100 papers published in international journals and the proceedings of international conferences with referees. He was invited to several conferences as a "plenary lecturer", to present the results of his research. Since 2005, he is working in teaching, and since 2007 he has been a teaching assistant in all the physics courses at the Engineering Faculty of Salerno University. He is a member of the Council of the Department and he is a delegate of the Director for student careers and mentorship/tutorship.

**Lamberto Tronchin** is Associate Professor in Environmental Physics from the University of Bologna and is recognised internationally as a leading authority on the subject of sound and acoustics. A pianist himself, with a diploma in piano from the Conservatory of Reggio Emilia, Dr Tronchin's principal area of research has been musical acoustics and room acoustics. He is the author of more than 200 papers and was Chair of the Musical Acoustics Group of the Italian Association of Acoustics from 2000 to 2008. Dr. Tronchin is the President of the Italian Section of AESDr. Tronchin is Associate Editor of the Journal of the Audio Engineering Society and a member of the Scientific Committee of the CIARM, the Inter- University Centre of Acoustics and Musical research. He has chaired sessions of architectural and musical acoustics during several international symposiums, been a referee for a number of International journals and is the Chair of the Organising and Scientific Committees of IACMA (International Advanced Course on Musical Acoustics). He was a visiting researcher at the University of Kobe in Japan, a visiting professor at the University of Graz in Austria and a special honoured International Guest at the International Workshop, 'Analysis, Synthesis and Perception of Music Signals', at Jadavpur University of Kolkata, India, in 2005. He has chaired the International Advanced Course on Musical Acoustics (IACMA), organised by the European Association of Acoustics, which was held in Bologna, in 2005. In 2008 and 2009, he gave plenary lectures at international congresses on acoustics in Cambridge (UK), Paris, Vancouver, Prague, Bucharest, Santander, Kos, Malta, Rodi, as well as at WIELS in Bruxelles. Dr Tronchin holds a Masters Degree in Building Engineering and a PhD in Applied Physics (Architectural Acoustics) from the University of Bologna. He has completed advanced courses on the Mechanics of Musical Instruments at CISM, Udine, Italy, and on Noise and Vibration at the University of Southampton in the UK where he has also worked as a visiting researcher. He was granted a post-doctoral scholarship in Room Acoustics. He designed theatres and other buildings, as an acoustic consultant, in collaboration with several architects, among them Richard Meier and Paolo Portoghesi. He is the inventor of an international patent with Alma Mater Studiorum: "Method for artificially reproducing an output signal of a non-linear time invariant system".

**Massimo Viscardi** graduated from the University of Naples "Federico II" with a degree in Aeronautical Engineering in 1994, and a Ph.D. in Aerospace Engineering in 1998. From 2006, he has been a confirmed researcher at the University of Naples "Federico II", and is Assistant Professor of Structural Testing and Assistant Professor of Experimental Vibroacoustic at the same university. His main research fields of interest are: acoustics and vibrations, structural testing, sensing and actuation systems, composite materials, natural fibers applications, NDT/Health monitoring activities and noise and vibration active control. He has been involved in many research activities within national (L297, PIA, Regional research program, PRIN) and international programs, like (EU Funded ASANCA, ASANCA II, MESA, FACE, DAMOCLES, PIROS, AERONEWS, ITEM B, SPAIN). From April 2002 to the end of 2006, he was a Member of the Scientific Committee on Ministry of Industry for the Research and Development program . He has developed specific competencies and has promoted industrial applications for innovative materials and relative testing facilities. He has been scientific coordinator of three international projects and more than 15 research and innovation projects within the national research framework. He is the author of more than 100 scientific works and more than 200 thesis degrees. He held more than 30 training courses in the reference research field at post-graduate and master's level.

### **Preface to "Modelling, Simulation and Data Analysis in Acoustical Problems"**

Modelling and simulation in acoustics is currently gaining in importance. In fact, with the development and improvement of innovative computational techniques and with the growing need for predictive models, an impressive boost has been observed in this domain.

The design of a model requires a proper conversion of reality to functions and parameters. Once the model has been designed, an adequate simulation must be run in terms of modelling and computational parameters. Keeping in mind the limitations and approximations of any model, data analysis, both online and offline, is the last step of this process and can be extremely important to extract the required output from the process. These basic and general concepts can be applied in many acoustical problems. In acoustics, there is a large demand for modelling and simulation in several research and application areas, such as noise control, indoor acoustics, and industrial applications.

These factors led us to propose a special issue about "Modelling, Simulation and Data Analysis in Acoustical Problems", https://www.mdpi.com/journal/applsci/special issues/MSDAAP, as we believe in the importance of these topics in modern acoustics' studies. In total, 81 papers were submitted and 33 of them were published, with an acceptance rate of 37.5%. Among the 33 papers published, two of them were classified as review papers, while the rest are classified as research papers.

According to the number of papers submitted, it can be affirmed that this is a trending topic in the scientific and academic community and this special issue will try to provide a future reference for the research that will be developed in the coming years.

> **Claudio Guarnaccia, Lamberto Tronchin, Massimo Viscardi** *Special Issue Editors*

### *Editorial* **Special Issue on Modelling, Simulation and Data Analysis in Acoustical Problems**

#### **Claudio Guarnaccia 1,\* , Lamberto Tronchin <sup>2</sup> and Massimo Viscardi <sup>3</sup>**


Received: 19 November 2019; Accepted: 23 November 2019; Published: 3 December 2019

#### **1. Introduction**

Modelling and simulation in acoustics is gathering more and more importance nowadays. In fact, with the development and improvement of innovative computational techniques and with the growing need of predictive models, an impressive boost has been observed in this domain. The design of a model needs a proper conversion of reality to functions and parameters. On the other hand, once the model has been designed, an adequate simulation must be run, in terms of modelling and computational parameters. Keeping in mind the limitation and the approximations of any model, the data analysis, both online and offline, is the last step of this process and can be extremely important to extract the required output from the process.

These basic and general concepts can be applied in many acoustical problems. In acoustics, in fact, there is a large demand for modelling and simulation, in several research and application areas, such as noise control, indoor acoustics, industrial applications, etc.

These motivations led us to the proposal of a Special Issue about "Modelling, Simulation and Data Analysis in Acoustical Problems", since we definitely believe in the importance of these topics in modern acoustics studies. In total, 81 papers were submitted and 33 were published, with an acceptance rate of 37.5%. Among the 33 papers published, two of them were classified as review papers, while the rest were classified as research papers. According to the number of papers submitted, it can be affirmed that this is a trendy topic in the scientific and academic community and this special issue will try to be a future reference for research to be developed in the next years.

#### **2. Modelling, Simulation and Data Analysis in Acoustical Problems**

As stated in the introduction, the need of models in acoustical problems is very large and can interest many subareas of acoustics. Withstanding this variety of possible applications and considering the interdisciplinary features of the Special Issue, several topics were studied in the papers submitted to the issue.

The noise control topic is studied by Wang et al. [1] and Zhang et al. [2], regarding respectively noise barrier insertion loss and car silencers.

Medical applications are presented in [3–6]. Hearing issues are studied by Cucis et al., comparing normal-hearing subjects and cochlear implant users, and Ito et al., regarding the effects of surgical instruments in ear surgery. High-intensity focused ultrasound (HIFU) non-invasive therapy is studied by Liu et al., Tan et al. and Gutierrez et al., concerning respectively the prediction of HIFU propagation in a dispersive medium, the influence of dynamic tissue properties on HIFU hyperthermia and acoustic field of focused ultrasound transducers.

In addition to HIFU, ultrasound waves are presented in several and various applications [7–10]. In particular, Choo et al. proposes a method to estimate the soil depth based on elastic wave velocity. Some underwater applications are presented by Wang et al. and Wang et al., considering respectively underwater acoustic communication and channel modelling and estimation, and underwater acoustic sources estimation. Jin et al. [11] presents an application of hydro-elastic analysis of a submerged floating tunnel under extreme wave and seismic excitations.

The topic of the acoustic ultrasound emissions study for monitoring of fatigue crack growth in mooring chains [12] and for rail defect detection [13] is faced respectively by Angulo et al. and Shi et al., with very interesting results. Dobrzychi et al. applies the acoustic emissions technique to epoxy resin electrical treeing study [14]. Also, Teng et al. [15] proposes the evaluation of cracks in metallic material.

Vibration and vibroacoustic studies are presented in [16–19], by Chatterhee et al., Wu et al. and Qian et al., with applications to parabolic tapered annular circular plate (Chatterhee et al. 2018 and 2019) and sensorized prodder for landmine detection (Wu et al.). Also, Flückiger et al. studies the vibrations, but referred to piano keys and their influence on piano players' perception and performance [20].

Yin et al. [21] and Jiang et al. [22] present their results on transducers, respectively on a 3D model of electromagnetic acoustic transducers and balanced armature receiver optimal design.

Target speech with nonlinear soft masking is presented by Zou et al. [23] while Tronchin et al. [24] reports a study about spatial information on voice generation.

Propagation in a fluid-filled polyethylene pipeline is studied from an experimental point of view by Li et al. in [25].

Bo et al. [26] presents a study on the acoustics of the Syracuse open-air theatre, with experiments and simulations. A room acoustic experiment is proposed by Wang et al. in [27] to investigate fingerprinting acoustic localization indoors.

Noncontact audio recording and a multi-frame Principal Component Analysis (PCA) based stereo audio coding method are proposed respectively by Sato et al. [28] and Wang et al. [29].

Tarrazó-Serrano et al. [30] proposes a material for acoustic lenses compatible with magnetic resonance imaging.

An acoustic detection method for localization is proposed by Yin et al. in [31].

Kirkup [32] proposes a survey of the boundary element method in acoustics.

Sparse impulse response estimation is treated in [33] by Lim et al.

#### **3. Conclusions**

All the researches presented, published in this Special Issue, suggest that the topic of modelling and simulation in acoustic problems is extremely important and popular in the scientific community. The advances proposed by all the authors push further the knowledge in this area and open the way to new and interesting possible evolutions. We believe that this issue could become a reference in the near future of modelling and simulation in acoustics. The new horizons in this research area will be traced starting by state of the art innovations largely represented by the papers of this issue.

**Acknowledgments:** The success of this Special Issue is strongly related to the huge work and the great contributions of all the authors. In addition, we acknowledge the hard work and the professional support of the reviewers and of the editorial team of Applied Sciences. We are extremely grateful to all the reviewers involved in the issue, for their time and their knowledge. We congratulate the assistant editors from MDPI that collaborated with us and thank them for their tireless support. We hope that the editorial process, starting from the submission and focusing on the review, was appreciated by all the authors, despite the final decisions. The real value of the time and the work spent in this process must be traced in the help provided to the authors to improve their papers.

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

#### **References**


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

### *Article* **Acoustic Emission Monitoring of Fatigue Crack Growth in Mooring Chains**

**Ángela Angulo 1,2,\* , Jialin Tang 1, Ali Khadimallah 1, Slim Soua 1, Cristinel Mares <sup>2</sup> and Tat-Hean Gan 1,2**


Received: 11 March 2019; Accepted: 20 May 2019; Published: 28 May 2019

#### **Featured Application: Structural health monitoring experimental methodology for crack initiation and crack growth analysis for damage detection in mooring chains using acoustic emission.**

**Abstract:** Offshore installations are subject to perpetual fatigue loading and are usually very hard to inspect. Close visual inspection from the turret is usually too hazardous for divers and is not possible with remotely operated vehicles (ROVs) because of the limited access. Conventional nondestructive techniques (NDTs) have been used in the past to carry out inspections of mooring chains, floating production storage and offloading systems (FPSOs), and other platforms. Although these have been successful at detecting and assessing fatigue cracks, the hazardous nature of the operations calls for remote techniques that could be applied continuously to identify damage initiation and progress. The aim of the present work is to study the capabilities of acoustic emission (AE) as a monitoring tool to detect fatigue crack initiation and propagation in mooring chains. A 72-day large-scale experiment was designed for this purpose. A detailed analysis of the different AE signal time domain features was not conclusive, possibly due to the high level of noise. However, the frequency content of the AE signals offers a promising indication of fatigue crack growth.

**Keywords:** structural health monitoring; acoustic emission; mooring chain; fatigue crack growth; structural integrity

#### **1. Introduction**

Offshore operators are constantly concerned by the safety and integrity of their assets, due to the high stakes involved. Ageing structures notably pose a significant threat to human lives and can incur exorbitant costs when unplanned shutdowns or catastrophic failures occur. Furthermore, in addition to the structural challenges that onshore structures experience, offshore assets withstand harsh marine environments as a result of severe storms, highly corroding sea water, seaquakes, and cyclic wave loading [1,2]. In this context, structural health monitoring (SHM) tools have been developed over the last few decades to mitigate such risks and offer continuous monitoring solutions that can be used in different industries.

One of the major problems in the design of offshore equipment is fatigue damage accumulation. Although this topic has been extensively studied in the literature, theoretically, numerically, and experimentally [3–9], the available inspection and monitoring technologies developed to date have not been able to fully overcome the severe environmental challenges associated with offshore service activities. Remotely operated vehicles (ROVs) have been widely used since the 1970s but face serious

difficulties, despite technological advances, due to the highly unpredictable operating environment characterised by poor visibility and unstable conditions [10,11]. For fatigue damage detection in structural applications in general, several sensing techniques have been developed [12], including guided ultrasonic waves [13,14], fibre Bragg gratings [15–17], strain gauges [18], and piezoelectric sensors [19]. Acoustic emission (AE) has also been proposed as a potential solution to detect and monitor cracking in structures, such as vessels and pipelines [20] and bridges [21]. Roberts et al. [22] found a correlation between crack evolution rates and the AE count rates for a narrow range of loading in steel specimens subject to tension. Yu et al. [23] predicted the crack growth behaviour in compact tension steel specimens based on AE data. The AE data were filtered using specific techniques and obtained from AE sensors placed around the crack tip. The authors showed that an accurate life prediction model could be established when the absolute energy of the AE signals was analysed. Despite its limitations, it is now recognised that the AE technique is capable of monitoring fatigue crack initiation [24] and propagation [25] in steels and other metals.

Amongst the offshore assets that are vulnerable to corrosion-enhanced fatigue damage, mooring chains are one of the most crucial mooring components used in permanently anchored structures [26]. Despite their importance, limited experimental work on mooring chains exists in the literature. Studies have been carried out on the chains' material microstructural properties [27,28], whereas others have been numerical modelling oriented [29]. Few large-scale testing attempts have been made. Rivera et al. [30] conducted a 4-month feasibility study to establish the AE technique's capabilities in monitoring damage in mooring chain links subjected to stress corrosion cracking in artificial sea water. Although a metallurgical examination was not performed, the authors suggested that the AE technique has promising potential in monitoring fatigue cracking in mooring chains.

The present work is a continuation of a research programme [30–32] aimed at identifying the key AE signal features for the prediction of fatigue crack growth in mooring chain links. The primary goal of this study was to investigate the applicability of using ultrasonic guided waves (UGW) and AE approaches for detecting and monitoring crack initiation, location, and propagation on a mooring chain. The study shows preliminary modelling simulations using different finite element analysis (FEA) methods.

In this paper, in Section 2, a description of the full-scale test rig and the monitoring setup is presented. The results and discussion then follow before a brief conclusion. Due to the complexity of the experiment, a substantial amount of data was collected and analysed. For the purpose of simplicity and brevity, only the final outcomes relevant to the scope of this paper are presented.

#### **2. Experimental Procedure**

#### *2.1. Test Rig Description*

The experiment carried out was part of a large Joint Industry Project (JIP) research programme, which aimed to perform large-scale tests for the analysis of the fatigue performance of mooring chains in seawater. A full-scale fatigue test rig was arranged to perform the AE measurements. The setup was designed to test a short section of a chain made of 7 links in artificial seawater, with a link diameter of 127 mm and a maximum load capacity of 700 tons. The total length and width of the rig was 7.65 m and 2.2 m, respectively.

The chain was initially subjected to a tension of 1000 kN in order to calibrate the load cell. A cycling tensile loading ranging between 3113 kN and 3497 kN at 0.5 Hz frequency was applied during the 72-day experiment (Figure 1).

#### *2.2. Hardware Selection and Sensor Deployment*

The AE data acquisition was performed using a Vallen AMSY-6 (Vallen Systeme GmbH, Icking, Germany) digital multi-channel AE-measurement system (ASIP-2 dual channel acoustic signal processor). Different AE sensors with different resonance frequencies and frequency bandwidths were

tested: Vallen-VS150-WIC-V01 (resonant frequency 150 kHz, bandwidth 100–450 kHz), VS375-WIC-V01 (resonant frequency 375 kHz, bandwidth 250–700 kHz), and VS900-WIC-V01 (resonant frequency 350 kHz, bandwidth 100–900 kHz). It has been found [33,34] that chain failure most likely occurs at the point of the intrados (KT point) and crown positions, due to higher localised stresses in these areas (Figure 2). Therefore, on each link, the sensor was placed 10 cm away from the weld on an accessible side (Figure 3b). Each sensor was equipped with an integrated 34 dB pre-amplification. Four sensors were used on Links 3, 4, 5, and 6 (Figure 3a, and Figure 3c): one VS150, one VS375, and two VS900. A water-based couplant was used to facilitate the transmission of the sound signal between the transducer and the link's surface. Links 1 and 7 were not monitored in the present setup.

**Figure 1.** Test rig illustration: the chain is fixed at one end (right) and the loading (strain) is applied at the other end (left).

**Figure 2.** Chain link potential failure locations.

The sensors used in this study had already been calibrated according to the American Society for Testing and Materials standard, ASTM E1106 [35], prior to the experiment to ensure the reliability of the collected signals. Since the experiments were to be performed with the chain submerged in water, the reproducibility of the AE sensors response was verified in air and underwater, according to ASTM E976 [36], by carefully breaking a 0.5 mm pencil lead against the link's surface (the Pencil Lead Breakage (PLB) test). PLB, also known as the Hsu and Nielsen [37] pencil lead break, is a well-established technique and has long been used as a method to artificially generate reproducible AE signals [36]. Additionally, to verify the sensor coupling, the pulsing function was applied. Each sensor was used as a signal generator that sends signals to be intercepted by the rest of the sensors. These details are briefly summarised in Table 1.



**Figure 3.** Sensor deployment on the chain in (**a**) the full-size test rig. Sensors fixed using magnetic holders 10 cm away from the link's weld (**b**). Illustrative drawing of the chain (**c**).

Figure 4 shows the results of the PLB calibration in air. Sensor 1 (S1), located on Link 3 (L3) and referred to as S1L3, and sensor S3L5 both correspond to the VS900 broadband AE sensor and seemed to receive higher dB levels of ~88 dB. For sensor S2L4 and sensor S4L6, corresponding respectively to VS150 and VS375, the dB levels received were approximately 10 dB lower at around 78 dB. From the figure, it is clear that the AE wave propagated through the links—the dB levels dropped to ~64 dB in the first adjacent link, ~54 dB in the second, and ~45 dB in the third (Figure 4a). Signal reproducibility was also demonstrated using the pulsing function. Each sensor generated four consecutive signals of peak amplitude equal to 81 dB (top rectangle in Figure 4b), which were intercepted by the rest of the sensors. As can be seen in Figure 4b, the wave was attenuated to a minimum of ~54 dB in the furthest sensor.

**Figure 4.** Calibration results for S1L3 (yellow), S2L4 (red), S3L5 (blue), and S4L6 (green) using (**a**) PLB and (**b**) the pulsing technique.

Figure 5 shows the behaviour of a wave generated in Link 3 from a PLB calibration event, where the peak amplitude of the signal ranges from ~90 dB in Link 3 to less than 45 dB in Link 6. The evolution of the other AE signal time domain features as the wave propagates through the chain links is shown in Figure 5. As the wave propagates through the links, the peak amplitude, number of counts, and average frequency decrease, whereas the rise time, duration, and time of arrival increase.

**Figure 5.** A pencil lead calibration (PLC) event generated in Link 3.

It should be noted that the signal reproducibility was also verified when the tank was being filled with water and the links were partially (almost completely) submerged (Figure 6).

In-air calibration was performed in order to understand in detail the propagation and resulting AE patterns across the links. The results obtained from underwater calibration during and after the tank filling was completed are similar to those obtained from the in-air calibration. However, due to the complexity of the setup, the values and conclusions acquired when full access was granted to the sensors' positions (in-air) were considered as the final calibration results.

**Figure 6.** AE reproducibility was verified as the tank was being filled with water (**a**) before the onset of the experiment (**b**).

#### **3. Results and Discussion**

The total duration of the AE monitoring period was 72 days. Taking into account the findings of the reproducibility tests, the AE acquisition threshold was originally set at 45 dB. However, due to the large unsustainable level of events captured, mostly noise, the threshold was initially increased to 55 dB then finally to 60 dB within the first 24 h of the experiment. The data recorded in this experiment included the number of AE hits (number of emissions detected by the sensors) and AE time domain features, including the number of counts, energy, rise time, time of arrival, and duration, in addition to the load and displacement measured by the load cells. Sensors S1L3 and S4L6 failed during the first days of the experiment for unknown reasons and were not replaced. Only the data from S2L4 and S3L5 were considered in the present analysis.

The first time the pre-set limit displacement was exceeded occurred after 4,333,424 cycles. After the water tank was fully drained, a complete inspection of the chain using visual inspection and magnetic particle inspection (MPI) revealed a ~5 mm crack located at the weld area in Link 3 (Figure 7a). No other indication of cracking was found, and the experiment was resumed. At 4,923,552 cycles, close to the fatigue life (5 <sup>×</sup> 10<sup>6</sup> cycles when similar loading conditions are applied), the experiment was interrupted. A MPI inspection revealed that the same crack in Link 3 had grown to ~120 mm (Figure 7b). A complete inspection of the chain did not reveal any other cracks.

**Figure 7.** Crack indications after magnetic particle inspection (MPI) at (**a**) 4,333,424 cycles (5 mm) and (**b**) 4,923,552 cycles (120 mm).

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

The crack of interest was revealed 31 days before the end of the experiment. It should be noted that a few interruptions occurred at the beginning of the experiment. Fortunately, these were during the early stages of the test and the information collected after day 41 can be considered accurate. For these reasons, only the AE data that were acquired during the last 31 days were processed. The key AE signal features were examined in an effort to establish correlations that would link the AE activity to the fatigue crack evolution. Some of the relationships considered were:


The AE signal features detected during the experiment were studied in detail. An example of the peak amplitude vs. load plot is shown in Figure 8. As can be seen in the figure, the events detected by the sensor attached to Link 5 are exclusively located within certain load ranges: 3150–3250 kN, 3700–3780 kN, and 3800–3850 kN. The corresponding events detected in Link 4 cannot be distinguished due to their uniform distribution. The location pattern of the red events without a clear connection to any cracking suggests they are most likely due to a mechanical noise, possibly caused by the test setup or surface friction where the links meet. It is usually accepted that in fatigue cracking, the AE signals generated close to the high end of the loading range are related to crack growth. Roberts et al. [22] showed that any correlation attempt was unsuccessful when the complete loading range was considered. Attributing the top 10 and 5% of the load range to crack growth improved the correlation to AE significantly. Bhuiyan et al. [39] associated different AE signal groups to the crack behaviour based on a spectral analysis of the AE data. The authors discovered that the AE events generated in the top 75–85% of the load range were statistically different from those below 60% and were related to the fatigue crack growth. In the present study, due to the large amount of noise generated by the full-scale experiment, only data collected above 3820 kN are shown to offer good correlations.

**Figure 8.** Peak amplitude (dB) vs. load (kN). Details of the minimum–maximum loading, and the 3820–3873 kN range is highlighted.

A closer look at the frequency content of the signals detected during the experiment, presented in Figure 9, shows a continuous shift of the average frequency towards higher values throughout the monitoring period.

**Figure 9.** Hits vs. average frequency (kHz). Shift of the average frequency.

The shift in the AE signal frequency content may be indicative of a change in the damage mode. For example, in reinforcing steel bars used in reinforced concrete, a decrease of the average AE frequency is attributed to the onset of phase three of the corrosion loss phenomenological model [40]. Additionally, tensile cracking is associated with a higher AE signal frequency range when compared to shear cracking in the same class of materials [41]. In composites, the AE signal frequency can be correlated to the fracture mechanism [42].

Finally, the observation of the cumulative energy during the last 31 days of the experiment did not show any particular change except during the final cycles (yellow circle in Figure 10). The fast rise indicates a rapid release of energy corresponding to rapid crack propagation. No statistical change was noted in the slope of the cumulative energy after the crack was initially detected. A quick increase in the cumulative energy precedes the exposure of the 120 mm crack. The displacement measured by the gauge did not exceed the pre-set limit until 120 mm was reached. This indicates a sharp crack growth, as seen in Figure 10.

**Figure 10.** Cumulative energy (load > 3,820 kN) evolution as measured by the sensor fixed on Link 4. Photograph shows the fracture surface corresponding to the final crack.

#### **4. Conclusions**

The aim of this work was to assess the ability of AE to detect crack initiation and growth in mooring chains under realistic loading and environmental conditions. After careful evaluation of the different AE signal features and all possible correlations, it appears that the frequency content of the AE signals is the most promising parameter. An increase of the average frequency is observed with the growth of the crack in the chain link. However, due to the challenging environment and the high level of noise recorded, a more comprehensive frequency analysis will be needed.

**Author Contributions:** Conceptualization, Á.A.; Data curation, Á.A., J.T.; Formal analysis, Á.A., S.S.; Investigation, Á.A., A.K.; Methodology, Á.A., J.T.; Project administration, Á.A.; Resources, Á.A.; Software, Á.A., J.T.; Writing—original draft preparation, Á.A.; Supervision, C.M., S.S., T.-H.G.; Validation, C.M., S.S.; Writing—review & editing, Á.A., A.K., T.-H.G.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors gratefully acknowledge the opportunity provided by TWI's JIP project (22116) sponsors to make use of the mooring chain full-scale test rig at TWI's facilities.

**Conflicts of Interest:** The JIP sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


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

### *Article* **Using ANN and SVM for the Detection of Acoustic Emission Signals Accompanying Epoxy Resin Electrical Treeing**

### **Arkadiusz Dobrzycki 1,\* , Stanisław Mikulski <sup>1</sup> and Władysław Opydo <sup>2</sup>**


Received: 26 February 2019; Accepted: 8 April 2019; Published: 12 April 2019

**Featured Application: From a practical point of view, the phenomenon of electrical treeing can be an exploitation problem, especially in the elements or places with locally increasing of electric field intensity, because it is an irreversible process. Such places are e.g., connecting points of the power network elements or parts of electrical devices in which there is constant insulation. The results of the carried out analyzes may complement the knowledge of identifying the type of PD occurring based on selected signal parameters.**

**Abstract:** Electrical treeing is one of the effects of partial discharges in the solid insulation of high-voltage electrical insulating systems. The process involves the formation of conductive channels inside the dielectric. Acoustic emission (AE) is a method of partial discharge detection and measurement, which belongs to the group of non-destructive methods. If electrical treeing is detected, the measurement, recording, and analysis of signals, which accompany the phenomenon, become difficult due to the low signal-to-noise ratio and possible multiple signal reflections from the boundaries of the object. That is why only selected signal parameters are used for the detection and analysis of the phenomenon. A detailed analysis of various acoustic emission signals is a complex and time-consuming process. It has inspired the search for new methods of identifying the symptoms related to partial discharge in the recorded signal. Bearing in mind that a similar signal is searched, denoting a signal with similar characteristics, the use of artificial neural networks seems pertinent. The paper presents an effort to automate the process of insulation material condition identification based on neural classifiers. An attempt was made to develop a neural classifier that enables the detection of the symptoms in the recorded acoustic emission signals, which are evidence of treeing. The performed studies assessed the efficiency with which different artificial neural networks (ANN) are able to detect treeing-related signals and the appropriate selection of such input parameters as statistical indicators or analysis windows. The feedforward network revealed the highest classification efficiency among all analyzed networks. Moreover, the use of primary component analysis helps to reduce the teaching data to one variable at a classification efficiency of up to 1%.

**Keywords:** solid dielectrics; acoustic emission; artificial neural networks; electrical treeing; wavelets; non-destructive testing; high-voltage insulating systems

#### **1. Introduction**

The major factor that contributes to the deterioration of the insulation characteristics of dielectrics is partial discharges (PD), which occur under the influence of high-intensity electric fields. Partial discharges occur on the surface of or inside dielectrics, causing deterioration of their electrical insulating characteristics. PD in solid dielectrics often entails electrical treeing, which involves the formation of conductive or semi-conductive channels in the shape of trees inside the dielectric. Discharges in the channels result in further development of trees and finally lead to the low-resistance short-circuit of electrodes. That is why the forecasts for insulation in which treeing occurs are disastrous. Therefore, the information on whether treeing has been initiated in the insulation becomes of supreme importance for the purposes of high voltage (HV) insulation condition diagnostics.

Early studies on partial discharges (PD) occurring in dielectrics were carried out at the beginning of the Twentieth Century by Rayner, who investigated the impact of partial discharges on breakdown [1]. The development of the studies was stimulated by the development of power engineering and the high unreliability of paper and oil insulation used in HV cables, transformers, and generators. The studies were continued in the 1930s by Robinson [2] and Whitehead [3], who demonstrated that partial discharges in paper and oil insulation caused electrical treeing, which involves the formation of channels in the shape of trees or shrubs, or a sight in a dielectric, and are the main cause for the low durability of high-voltage cables with this type of insulation.

The main direction of studies in the second half of the Twentieth Century included the following: electrical treeing initiation [4], studies on the relationship between tree length, voltage impact time, and time to breakthrough [5,6], studies on the properties of cables with solid polyethylene insulation [7], and insulation breakthrough tests [8]. There have also been studies on partial discharge detection and measurement methods. Initially, mainly oscilloscopic methods were proposed by Tykociner [9]. Some studies on the development of insulation ageing models were also carried out [10].

In addition to detailed papers, several publications summarized and systematized the current knowledge status [11,12].

The turn of the Twenty-first Century saw studies related to forecasting treeing in dielectrics. Engineering progress contributed to the use of modern and sophisticated methods of partial discharge testing.

Shimizu [13,14] pointed out the possibility of examining treeing incubation by means of electroluminescence. Papers by Kudo [15] and Dissado [16] discussed the ways to predict treeing propagation using the theory of fractals and deterministic chaos. Since trees have a fractal-like nature, the authors used fractal size to identify the shape of the tree being formed. The theory of deterministic chaos was, in turn, used to identify the probable directions of tree development. One can observe a tendency for the development of non-invasive partial discharge methods, and the acoustic emission method in particular [17–22]. Advanced signal analysis methods, such as wavelet transformation or artificial neural networks, used for analysis parameters of solid insulation systems and identification of the discharge source and type have been employed [10,23–29].

Among the methods used for the detection of partial discharges (e.g., in the transformer tank), the acoustic emission method is common, focusing on the recording and analysis of an acoustic wave propagating in the material as a result of external impact exerted on it. The sources of the impact may include mechanical pressure (testing the stress inside the material) [30–34], electric field (partial discharge testing), or magnetic fields (used for Barkhausen noise analysis).

The purpose of the study was to try to automate the process of identifying the condition of insulation materials by developing a neural classifier allowing the detection of symptoms related to treeing in the recorded acoustic emission signals.

#### **2. Test Setup and the Course of the Experiment**

The first stage of the studies involved measurements of acoustic emission signals in electrically-stressed epoxy resin samples. The samples had cubicoid shapes of 25 mm × 10 mm × 4 mm. One sample surface of 25 mm × 4 mm was ground and coated with varnish conducive to electric current.

Since the process of tree channels forming in a dielectric may last very long, an electrode made of a T10 surgical needle with liquid resin poured inside the sample during sample formation was used, such that the distance between the sample bottom and the electrode end ranged from 1–3 millimeters. The procedure allowed obtaining the electric field intensity between the needle electrode and the electrode applied to the bottom of the sample, which was high enough to reduce the tree-forming duration to a few hours. The needle was connected to the neutral terminal of the transformer. The opposite base of the sample was adjacent to a plane copper electrode. The electrode was connected through a resistor of 0.5 MΩ resistance to a high-voltage terminal of a test transformer with the ratio of 220 V/30 kV and power of 10 kVA. The voltage was measured with an electrostatic voltmeter.

As part of the tests, acoustic signals were recorded for a dozen or so samples in which, under the influence of a high intensity of the electric field, the process of forming an electric tree began. The measurements were made for variable values of the supply voltage with a 50-Hz frequency and different distances between the electrodes.

The elastic waves of acoustic frequencies were emitted from the analyzed sample through a wave-guide made of a steel rod of 2 mm in diameter. One of its ends was put into a hole bored in the sample, while the other was connected to an electroacoustic converter.

As regards measurement time, the studied sample was placed in a methyl polymethacrylate vessel filled with electrical insulating oil. The voltage between the electrodes during the tests ranged from a few to several kV. The AE signals were measured by means of a Physical Acoustic Corporation (PAC) R3α electroacoustic converter and a filtering and amplifying system composed of 2/4/6 type pre-amplifier from PAC, a 20 ÷ 1000 kHz transmission band filter, and PAC AE5A amplifier. Upon amplification, the signal was recorded in the computer memory with an data acquisition (DAQ) NI–USB6251 card, which enabled signal recording with a sampling frequency up to 1 MS/s and 16-bit resolution. The test setup schema, the view of the measuring equipment, and the view a sample holder are presented respectively in Figure 1a–c.

(**a**)

**Figure 1.** *Cont.*

(**b**)

(**c**)

**Figure 1.** The test setup schema (**a**), view of the measuring equipment (**b**), and view of the sample holder (**c**) for measuring acoustic signals accompanying the treeing of solid dielectrics.

To ensure that the recorded signals concerned the electrical treeing while we were preparing the experiment, we also tested and observed under a microscope different states of the system, i.e., without any PDs, with corona, etc. We noticed that each of the PD types was connected with different waveforms of signals, and this can be found also in different researchers' results [35]. After that test, we rebuilt/reconfigured the setup and chose a voltage range to be sure that other types of PDs would not be present. Physically, electrical treeing is a combination of chemical and physical changes inside the specimen, so no one can be sure whether that particular signal is connected with the breaking polymer chain or PD inside the existing channel.

#### **3. Method of Features' Extraction**

#### *3.1. Features' Definition*

In order to extract the teaching data for neural classifiers, the analyzed signal fragment was divided into *x* blocks with *N* length of samples, where:

$$\mathbf{x} = \{\mathbf{x}\_1, \dots, \mathbf{x}\_N, \dots, \mathbf{x}\_N\}. \tag{1}$$

For each *x* signal block, a set of statistical parameters typically used in the analysis of AE signals was identified [36,37], which was further used as input data for the neural network. The following parameters were used: signal energy (*e*), band power (*pBP*), signal upper envelope (*env*), skewness (*skew*), and kurtosis (*kurt*).

The signal energy is identified based on the following relationship:

$$\mathfrak{e} = \sum\_{n=1}^{N} \mathfrak{x}\_{n-1}^{2} \tag{2}$$

Skewness, as the probability distribution asymmetry measure, is identified in the following way:

$$
\kappa \kappa w = \frac{\overline{x} - med(x)}{\sqrt{\sigma}} \tag{3}
$$

where *x* is the mean value of samples, med (*x*) the block median, and σ the sample block variance.

Kurtosis, which is the measure of the results' concentration around the mean, amounts to:

$$kurt = \frac{\mu\_4}{\sigma^4} - 3\tag{4}$$

where μ<sup>4</sup> is the fourth central moment of samples in the block and σ the sample block variance. A teaching data vector *d* was identified for each *x* signal sample block, where:

$$d = \left\{ \varepsilon, p\_{\text{BP}}, env, skew, kurt \right\} \tag{5}$$

In order to obtain training information *y*, which belongs to the set {0,1}, the recorded signal was blasted using the wavelet decomposition with the fifth-order Daubechies wave. Then, using the method of exceeding the trigger threshold, which was set at 20 mV for the presented problem, the information value y was allocated for each training set: 0 for no acoustic emission and 1 for the acoustic emission event (these are the red areas in Figure 2). Figure 2 presents one signal fragment with the sampling frequency of 1 MS/s selected for further studies. The red areas mark acoustic emission events determined on the basis of threshold crossing.

**Figure 2.** Fragment of the recorded signal where dielectric degradation occurred.

#### *3.2. Passband Power*

The selection of the band used for identifying band power pBP was made based on the averaged spectrum for the recorded AE signals. Figure 3 presents an averaged spectrum for 1000 AE signals recorded and the spectrum of a selected signal fragment where no acoustic emission was observed. For both spectrums, the transfer function of the selected detector (R3α), given in the calibration datasheet, was taken into account.

**Figure 3.** Averaged signal spectrum for 1000 recorded acoustic emission signals and measurement noise spectrum.

The presented spectrum was used to select a signal band with a frequency ranging from 15 kHz–30 kHz for further analysis because the power of the signal recorded within the band increased significantly when acoustic emission occurred.

#### *3.3. Block Length*

Selecting the right length of signal blocks to study is among the most important factors affecting classification quality. When performing an analysis using a taught network, we are not able to predict if the analyzed signal block starts before, during, or at the end of the acoustic emission signal, because it is necessary to know the influence of the signal window position against the acoustic signal on the values of selected network teaching features.

In order to identify the aforementioned relationships, window lengths of 10, 60, 200, and 300 samples were used. The windows were applied to a signal containing a single acoustic emission and moved by one sample (the window position against the AE signal was changed this way). A dataset was calculated for each position of the window. Figure 4 presents the values of each teaching feature as a function of the window position, at different window lengths, against the AE signal: kurtosis (a), skewness (b), and energy (c) respectively. Moreover, Figure 4c presents the identified upper envelope of the sample signal.

**Figure 4.** Influence of the position of different window lengths against the AE signal on the teaching properties' values: kurtosis (**a**); skewness (**b**); energy (**c**); upper envelope (**d**); the recorded signal wavelength is in the background.

#### *3.4. Principal Component Analysis*

Principal component analysis (PCA) is a statistical method for factor-based analysis. A collection *N* of *K* variable observations is analyzed as a collection of *N* points distributed in a *K*-dimensional space. Using the distribution of singular values for the deviation covariance matrix, it is possible to develop a new system that maximizes the variance of subsequent coordinates. This way, it is possible to reduce the space size (decrease the number of the analyzed input data) [30] and the size of the neural network, while maintaining the highest possible amount of information regarding the input process.

#### **4. Artificial Neural Network Classifiers Used to Select EA Signals**

#### *4.1. Feedforward Neural Network*

In a forward non-linear neural network (FNN), the flow of information (signals) is unidirectional. Its structure consists of the input layer, two or more hidden layers, and the output layer (Figure 5). All layer inputs can be linked only with the neurons in the preceding layer. Input signals *x* in a neuron are added to weights *w*. Neuron *y*'s output signal is calculated using a non-linear activation function ϕ, according to the relationship [38]:

$$\mathbf{y} = \phi \left(\sum\_{m=1}^{M} w\_m \mathbf{x}\_m \right) \tag{6}$$

The unipolar sigmoid function is among the most commonly-used activation functions:

$$\varphi(z) = \frac{1}{1 + e^{-z}}\tag{7}$$

alongside a bipolar function (hyperbolic tangent):

$$\varphi(z) = \frac{1 - e^{-z}}{1 + e^{-z}} \tag{8}$$

For learning feedforward neural network (FFN), the Levenberg–Marquardt algorithm was used. This method is one of the most effective learning algorithms. It has high convergence when network weights are near the optimal solution (as the Gauss–Newton method) and when the network is far from the optimal solution (as the descent gradient method). Detailed information about the Levenberg–Marquardt algorithm can be found in [39].

**Figure 5.** Sample structure of a feedforward network with three hidden layers.

#### *4.2. Radial Basis Functions*

A radial basis functions (RBF) network, similar to FNN, is a network with an oriented flow of signals. The network neurons have a radial activation function ϕ. The neuron output y is described by the relationship [40]:

$$w\_j(\mathbf{x}) = \sum\_{i}^{m} w\_{ij} \varphi(\|\mathbf{x} - \mathbf{c}\_i\|) \tag{9}$$

where *i* is the hidden layer neuron number, *j* the radial network output number, *wij* the coefficient of the weight, *x* the input data vector, and *c<sup>i</sup>* the center of the radial function for the *i*th neuron of the hidden layer.

In such a neuron, the value of the output signal is not proportional to the scalar product of *x* inputs and *w* neuron weights, but inversely proportional to the distance between *x* and the central point of radial function *c* located at the hyperspace of the network input parameters. The Gaussian function is among the most commonly-used radial functions. Figure 6 presents the structure of an RBF network.

**Figure 6.** RBF network diagram.

#### *4.3. Wavelet Neural Network*

A wavelet neural network (WNN) structure consists of an input, output, and hidden layer of wavelet neurons (Figure 7).

**Figure 7.** Wavelet neural network structure.

The function of wavelet neurons' activation ϕ is described by the relationship [41]:

$$\varphi(\mathbf{x}) = \prod\_{i}^{m} \psi(\frac{\mathbf{x}\_{i} - \zeta\_{i}}{\varepsilon\_{i}}).\tag{10}$$

where *x* is the input signal vector, ψ the wavelet function, ζ the translation parameter, and ε the scale parameter.

The process of such network teaching involves the selection of the scale parameters and shifting each wavelet ψ used inside a wavelet neuron nucleus. For the study, the authors used Mexican hat wavelet function.

Contrary to other network types, drawing of the initial values of the wavelet parameters ζ and ξ may result in the teaching algorithm being stuck in the local minimum. That is why network teaching is preceded by the initialization of initial parameters. A feature selection method developed by Oussar and Dreyfus was used for the study [42]. It consists of the following steps:


Subsequently, fully-initialized WNN was learned with the backpropagation algorithm. A detailed description of the WNN's structure and learning algorithm used for the study can be found in [41].

#### *4.4. Support Vector Machine*

A support vector machine (SVM) was another network used in the studies. The purpose of the SVM network operation is to use a hyperplane spread in the hyperspace of the input parameters, which separates a collection of points (sets of input parameters) belonging to two different classes, with a certain error margin.

It is a type of binary classifier. Its operating principle can be presented on the example of a network with two input parameters *x* = (*x*1, *x2*). Figure 8 presents the operating principle of such a classifier. The points lying in the plane (*x*1, *x2*) are divided into two classes marked in green and viloet. The SVM classifier searches for a straight line (black line in the drawing), which will maximize the error margin limited on both sides of the curve with support vectors, the brown lines in the drawing. The points exceeding the limits of the error margin will be misclassified. The purpose of SVM teaching is to select the coefficients of the straight line separating the points in such a way that the error margin is maximized (distance between the support vectors). Algorithms of non-linear optimization with constraints, and Lagrange's method in particular, were used for the network teaching [43].

**Figure 8.** Principle of support vector machine operation.

#### **5. Results**

Each of the described networks was taught in order to classify the signals into two groups: that which contains and that which does not contain acoustic emission signals accompanying dielectric treeing. The previously-described properties were used to develop teaching data for the blocks with the following lengths: 10, 60, 100, 200, and 300 samples.

In the first part of the results, classification efficiency was specified for a test fragment of the signal other than the network teaching signal. In the second section, the signal classification efficiency was presented for the teaching data reduced by means of PCA.

Additionally, the optimum size of the network was achieved in the cross-validation, which consisted of dividing the set of reference data into equinumerous *k* subsets. The network with a given number of neurons was learned on the basis of *k*-1 subsets, and the unused subset was the validation data. This process was repeated *k* times. Thus, the obtained *k* learning outcomes networks were averaged.

#### *5.1. Classifier E*ffi*ciency Analysis*

Each network taught was subjected to a test developed based on a random signal fragment. Table 1 presents the classification efficiency for each network type at variable block lengths based on which, calculations were made.


**Table 1.** Classification efficiency test results for variable data block lengths. WNN, wavelet neural network.

#### *5.2. Classification Using PCA*

Since the most stable efficiency results for each network were obtained for the data identified based on the block length of 300 samples, the data were subjected to PCA, and PCs were developed from them. Table 2 shows variances for different components and their percent share in the total variance.


**Table 2.** Classification efficiency test results for variable data block lengths.

Based on the results presented in the table above, it can be concluded that the share of the first primary component in the whole signal amounted to over 98%, which denotes the possibility to reduce a system of five decisive variables to only one variable with no significant information loss. Table 3 presents the efficiency of the analyzed networks for teaching when only the first primary component was used.

**Table 3.** Classification efficiency test results for variable data block lengths.


#### **6. Discussion**

The paper proposed the use of ANN for the analysis of acoustic signals accompanying electrical treeing in epoxy resin.

The recorded signals were divided into blocks used to identify the values of the teaching data. Then, the features of the reference signal were analyzed using a moving window of a set length.

The features were divided into two groups whose properties depended on window length and position versus the beginning of the AE pulse. Hence, kurtosis and skewness, due to the high dynamics of changes, can be good indicators for the identification of the beginnings of AE pulses related to treeing. The values of the parameters, however, did not reflect the duration of AE signals, and their use required the application of an analysis window with a minimum length. In the analyzed case, the results were regarded as satisfactory and repeatable for window lengths over 200 samples.

The other group of analyzed parameters, which describes the signal by its power characteristics, i.e., band power and signal energy, demonstrated elevated values during the signal. The steepness of the signal energy increase did not depend on the analysis of the window length, but a relevant selection of the length helped to identify the place in the signal where the signal amplitude dropped significantly. This is seen in the diagram as an inflexion point.

Since AEs in electrical treeing are related to the breaking of polymer chains or the presence of partial discharge in the existing channels, the value and dynamics of signal energy changes during a single pulse can be linked to treeing intensity.

By analyzing the results for each neural network, one can see that FFN demonstrated the highest efficiency. Slightly poorer results were obtained for RBF, WNN, and SVM. Moreover, for FFN and RBF, the results for each window length applied to the signal were comparable. WNN and SVN rendered better results for the greater analysis of window lengths. Moreover, WNN in its realization used only two wavelet neurons, while the structure of FNN required at least 15 neurons.

In the case of networks taught by means of the first primary component, one can state that the efficiency of each network dropped by a value that did not exceed 0.5%.

Summing up, the high efficiency in detecting AE signals accompanying electrical treeing needs to be highlighted. In order, however, to obtain as much information as possible on the AE pulse course, it seems justified to use an algorithm based on the following three stages:


A signal envelope analysis could provide hints for identifying the optimum window length at the discrimination level assumed on an arbitrary basis.

The performed experiments revealed that ANN can be successfully used for signal detection and classification, with the kurtosis, skewness, and energy of a single signal as the classifying parameters.

The effectiveness of the applied method was confirmed by the results obtained for various electrical parameters, i.e., different inter-electrode distances and different intensity of the electric field in this space. Regardless of the values of the above parameters, the results of AE signal identification were characterized by equally high efficiency, which proves the correctness of the assumptions made. The next stage covered a detailed analysis of the signals classified in the previous step with regard to the criteria assumed by the experimenter.

The completed analyses justify the usefulness of further studies on using ANN for partial discharge analysis and electrical treeing for solid dielectrics.

This knowledge could be of particular value to persons making diagnoses on the condition that the proposed method could be used online. That is why future efforts should focus on developing algorithms that would satisfy this demand.

**Author Contributions:** A.D. and S.M. developed a concept for the use of neural networks for the analysis of AE signals, analyzed the results of measurements, edited text and prepare conclusions. A.D. author of the methodology, designed and built a measurement stand, planned a research program, carried out preliminary tests. S.M developed computational algorithms, program code and performed numerical calculations. W.O. analyzed the results of measurements and gave some valuable suggestions, improved the text of article.

**Funding:** This research was funded by Polish Government, Ministry of Science and Higher Education.

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

#### **References**


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

### *Article* **An Ultrasonic Guided Wave Mode Selection and Excitation Method in Rail Defect Detection**

#### **Hongmei Shi 1,2, Lu Zhuang 1, Xining Xu 1,2,\*, Zujun Yu 1,2 and Liqiang Zhu 1,2**


**\*** Correspondence: xnxu@bjtu.edu.cn

Received: 27 February 2019; Accepted: 14 March 2019; Published: 19 March 2019

**Abstract:** Different guided wave mode has different sensitivity to the defects of rail head, rail web and rail base in the detection of rail defects using ultrasonic guided wave. A novel guided wave mode selection and excitation method is proposed, which is effective for detection and positioning of the three parts of rail defects. Firstly, the mode shape data in a CHN60 rail is obtained at the frequency of 35 kHz based on SAFE method. The guided wave modes are selected, combining the strain energy distribution diagrams with the phase velocity dispersion curves of modes, which are sensitive to the defects of the rail head, rail web and rail base. Then, the optimal excitation direction and excitation node of the modes are calculated with the mode shape matrix. Phase control and time delay technology are employed to achieve the expected modes enhancement and interferential modes suppression. Finally, ANSYS is used to excite the specific modes and detect defects in different rail parts to validate the proposed methods. The results show that the expected modes are well acquired. The selected specific modes are sensitive to the defects of different positions and the positioning error is small enough for the maintenance staff to accept.

**Keywords:** ultrasonic guided waves; SAFE; rail defect detection; mode excitation

#### **1. Introduction**

It is of great significance to discover the internal defects of rails in time for maintaining and ensuring the safety of trains. At present, the hand-push detection vehicle and large-scale detection vehicle are the main tools to detect the defect detection of Continuous Welded Rails [1–3]. These two kinds of detection vehicles are running in the maintenance time, generally from 12 a.m. to 4 a.m., and cannot be used for real-time monitoring of rail internal defects. Guided waves can propagate a long distance in the rail [4–6]. Based on the pulse-echo method [7–9], the on-line monitoring of rail internal defects can be implemented by detecting the echo signals generated at the rail defects during the propagation of the guided waves.

In the research of rail defect detection based on guided waves, some scholars selected the suitable modes for rail defect detection by a mode classification method. Hayashi [10] divided the modes of the rail base into three categories. By analyzing the dispersion characteristics of guided waves, the concept of dominant modes is proposed. In the frequency range of 60~200 kHz, the transverse and vertical vibration modes with smooth dispersion curves of rail base are more suitable for defect detection than the longitudinal vibration modes [11]. A rail base defect of 20 mm in length was successfully detected by exciting the vertical modes. Subsequently, Sheng Huaji [12] obtained the optimal excitation frequency and pulse period of the propagating modes at the rail base through simulation, and tested oblique cracks at different angles by numerical simulation. The results show that the detection effect of the vertical vibration modes is better than the transverse vibration modes. Lu Chao et al. [13] excited the vertical vibration modes on the rail head with a kind of mode force hammer. The experimental results show that the vertical vibration modes can effectively detect the transverse crack of the rail head. In addition, some scholars chose the modes to detect the damage of different parts of the rail based on the mode energy distribution. Gharaibeh et al. [14] selected the mode with energy concentrated on the rail head to detect the rail head defect. And the signal was excited according to the experimental method. The transverse defect of rail head at a distance of 9 m from the excitation location was successfully identified. In order to distinguish large transverse cracks with a diameter of 30 mm from 6 mm welded seam, Long and Loveday [15] compared the symmetric mode 1 and the antisymmetric mode 2 focusing energy at the rail head, and the symmetric mode 1 successfully distinguished the crack and welded seam at the rail head. Mode 3 and mode 4 focused with energy at the rail web, mode 4 was suitable for detecting the welded seam as well as damage in the rail web. It can be seen that in the research of the excitation mode applied to the detection of rail defects, most of the scholars use the experimental methods to excite the expected modes without specific theoretical guidance and mathematical model. The excitation signal usually contains multiple modes, which also makes it difficult for echo signal analysis and defect location.

In view of the above problems, this paper will focus on the guided wave mode selection and excitation methods based on a mathematical model in rail defect detection, by which the modes can be selected and excited with different sensitivities to the defects of the rail head, rail web and rail base, and accordingly the defect position can be rapidly found in different parts of the rail. The paper is organized as follows: A mode selection method is described in Section 1, which is used to select the modes having the concentrated strain energy and small attenuation corresponded to the rail head, rail web and rail base. The mode excitation method is presented in Section 2, to acquire the optimum excitation direction and excitation node of guided wave modes. Section 3 introduces the results of the target modes excitation simulation and verification. The simulation and detection results of the rail defects are demonstrated and analyzed in Section 4.

#### **2. Mode Selection Based on Semi-Analytical Finite Element (SAFE) Method**

In order to study the mode selection method of guided waves in rail defect detection, it is necessary to obtain the mode shape data. Many studies have proved that SAFE method is suitable to solve the propagation characteristics of guided waves in waveguide, such as plates [16], cylinders [17], rails [18], etc. SAFE method is also used in this paper to solve the mode shape data of CHN60 rail, which is a 60 kg/m China rail. The CHN60 rail coordinate system is shown in Figure 1.

**Figure 1.** Coordinate system of a rail.

The displacement *u*, stress *σ* and strain *ε* of each particle in the rail can be expressed as Equation (1):

$$\begin{aligned} \boldsymbol{u} &= \begin{bmatrix} \boldsymbol{u}\_{\boldsymbol{x}} & \boldsymbol{u}\_{\boldsymbol{y}} & \boldsymbol{u}\_{\boldsymbol{z}} \end{bmatrix}^{T} \\ \boldsymbol{\sigma} &= \begin{bmatrix} \boldsymbol{\sigma}\_{\boldsymbol{x}} & \boldsymbol{\sigma}\_{\boldsymbol{y}} & \boldsymbol{\sigma}\_{\boldsymbol{z}} & \boldsymbol{\sigma}\_{\boldsymbol{y}\boldsymbol{z}} & \boldsymbol{\sigma}\_{\boldsymbol{x}\boldsymbol{z}} & \boldsymbol{\sigma}\_{\boldsymbol{xy}} \end{bmatrix}^{T} \\ \boldsymbol{\varepsilon} &= \begin{bmatrix} \boldsymbol{\varepsilon}\_{\boldsymbol{x}} & \boldsymbol{\varepsilon}\_{\boldsymbol{y}} & \boldsymbol{\varepsilon}\_{\boldsymbol{z}} & \boldsymbol{\gamma}\_{\boldsymbol{y}\boldsymbol{z}} & \boldsymbol{\gamma}\_{\boldsymbol{x}\boldsymbol{z}} & \boldsymbol{\gamma}\_{\boldsymbol{xy}} \end{bmatrix}^{T} \end{aligned} \tag{1}$$

ε

Based on SAFE theory, it is assumed that the displacement field of the guided wave propagation direction is harmonic [19]. The cross section of the rail is discretized into 255 triangular elements and 177 nodes, as shown in Figure 2.

**Figure 2.** Discretization of triangle elements.

After discretization, the strain energy *φ* of each triangular element can be written as Equation (2):

$$\Phi = \int\_{V\_{\varepsilon}} \delta(\varepsilon^{(\varepsilon)T}) \mathbf{C} \varepsilon^{(\varepsilon)} dV\_{\varepsilon} = \delta q^{(\varepsilon)^{T}} (\mathbf{K}\_{1} - \mathrm{i}\_{5}^{\pi} \mathbf{K}\_{2} + \mathrm{i}\_{5}^{\pi} \mathbf{K}\_{3}) q^{(\varepsilon)} \tag{2}$$

where *q*(*e*) is the nodal unknown displacement of each triangular element; *ξ* is the wavenumber; *K*1, *K*2, and *K*<sup>3</sup> are stiffness matrices; *C* is the elastic constant matrix of the rail.

By substituting strain energy and potential energy at any node in the rail cross section into Hamilton's formula, the general homogeneous wave equation of guided waves can be derived as Equation (3) [19]:

$$\left[\mathbf{K}\_1 + i\xi\mathbf{K}\_2 + \xi^2\mathbf{K}\_3 - \omega^2\mathbf{M}\right]\_M \mathbf{U} = 0\tag{3}$$

where *M* is the mass matrix; *U* is the nodal displacements; *ω* is frequency.

By solving the eigenproblem of Equation (3), the corresponding wavenumber *ξ* can be obtained for a given frequency *ω*. The triangular element nodal displacements *U* are calculated from the eigenvalue problem. And the strain fields are reconstructed from Equation (2).

At the frequency of 35 kHz, the nodal displacement *q*(*e*), stiffness matrix *K*1, *K*2, and *K*3, and wavenumber *ξ* are substituted into the Equation (2), and the strain energy of each mode in each triangular element can be acquired. The strain energy distribution diagram of each mode is normalized as shown in Figure 3.

At the frequency of 35 kHz, there are 20 kinds of guided wave modes propagating in the rail. Our expected mode is that the strain energy is concentrated in a single part of the rail. From Figure 3, we can see that mode 1, mode 2, mode 3, mode 7, mode 10, and mode 20 have such characteristics while strain energy of the other modes is distributed at three rail parts.

In Figure 4a,b, the strain energy of mode 7 and mode 20 is concentrated on the rail head. In Figure 4c, it can be seen that the strain energy of mode 3 is focused on the rail web. In Figure 4d, Figure 4e,f, the strain energy of mode 1, mode 2, and mode 10 is concentrated on the rail base. Therefore, mode 7 and mode 20 can be selected to detect rail head defects, similarly, mode 3 for rail web and mode 1, mode 2 and mode 10 for rail base.

**Figure 3.** Strain energy distribution diagram of modes.

**Figure 4.** *Cont*.

**Figure 4.** Normalized strain energy amplitude of modes. (**a**) mode 7; (**b**) mode 20; (**c**) mode 3; (**d**) mode 1; (**e**) mode 2; (**f**) mode 10.

The dispersion curves of the above modes are shown in Figure 5. When the frequency dispersion occurs, the wave packet composed of different frequency components will be dispersed, causing the waveform of the received signal to be distorted relative to the transmitted signal. Therefore, it is necessary to select modes with good non-dispersion characteristics. At the frequency of 35 kHz, the dispersion curves of mode 7 are gradual than that of mode 20. This means that the non-dispersion characteristic of mode 7 is better than that of mode 20. Hence mode 7 is selected to detect the rail head defect. Similarly, mode 3 is selected to detect the rail web defect. Mode 1, mode 2, and mode 10 are suitable to detect rail base defects.

**Figure 5.** Group velocity dispersion curves.

#### **3. Mode Excitation Based on Mode Matrix**

In order to excite the expected modes, it is necessary to determine the optimal excitation direction and excitation node. The detailed excitation method is described in this section.

#### *3.1. Determine Excitation Direction*

Here Euclidean distance is used to describe the true distance between the vibration vectors of modes in the rail. The distances of displacement vector between mode *m* and other modes in *x*, *y* and *z* directions are calculated, and the direction with the largest Euclidean distance is chosen as the optimum excitation direction of mode *m*.

There are *p* nodes in the profile of rail cross section. Each node has 3 degrees of freedom. The displacements are represented as *xi*, *yi*, and *zi*, respectively. The Euclidean distances *Xmn*, *Ymn*, and *Zmn* of the vibration vectors of mode *m* and mode *n* in *x*, *y*, and *z* directions can be defined as follows:

$$X\_{mn} = \sqrt{\sum\_{i=1}^{p} (x\_{im} - x\_{in})^2} \tag{4}$$

$$Y\_{mn} = \sqrt{\sum\_{i=1}^{p} (y\_{im} - y\_{in})^2} \tag{5}$$

$$Z\_{mn} = \sqrt{\sum\_{i=1}^{p} (z\_{im} - z\_{in})^2} \tag{6}$$

At the frequency of 35 kHz, there are 20 kinds of guided wave modes in rails. Assuming the expected mode is mode *m*. Firstly, we can calculate the Euclidean distance between mode *m* and all other modes in *x*, *y*, and *z* directions respectively at the same frequency. Then 3 matrices of 20 × 20 can be obtained and the column m or the row m of the matrix represents the Euclidean distance between the mode *m* and other modes. Next, we calculate the mean of the distance between a given mode to all others in *x*, *y*, and *z* directions and call the value *xm*, *ym*, and *zm*. This is done separately for all three directions. By comparing the values of *xm*, *ym*, and *zm*, the direction with the largest value is selected as the optimum excitation direction.

According to the above method, the mode shape data of 20 modes at the frequency of 35 kHz are substituted into Equations (4)~(6). Then 3 matrices of 20 × 20 can be obtained and each column or each row of a matrix represents the Euclidean distance between the corresponding mode and other modes. This paper selects mode 1 to detect the rail base defect. Then we can calculate *x*1, *y*1, and *z*<sup>1</sup> of mode 1. By comparing the 3 values of *x*1, *y*1, and *z*1, the *z* direction with the largest value is selected as the optimum excitation direction. The calculation process of mode 3 and mode 7 is the same and *x*, *y*, and *z* of mode 1, mode 3, and mode 7 in *x*, *y*, and *z* directions are shown in Table 1.


**Table 1.** The values of *x*, *y*, and *z*.

According to the data in Table 1, the optimal excitation direction of mode 1, mode 3 and mode 7 is *z*, *y*, and *z* direction, respectively.

#### *3.2. Determine Excitation Nodes*

Covariance is employed to analyze the vibration displacement deviation of each node on the rail profile from the average displacement of all nodes. Then the optimal excitation point is selected according to the overall vibration displacement of the nodes and the node with the largest covariance is selected as the optimum positive excitation node in mode *m*. If there is a negative covariance value, the node with the minimum covariance value is the best inverse excitation node. Equation (7) describes the covariance of node *t* and node *e* of mode *m* in *x*, *y*, and *z* directions at the rail base.

$$RB\\_COV\_m(t,\varepsilon) = \frac{\sum\_{s=1}^{3} \left(t\_s - \overline{t}\right) \left(\varepsilon\_s - \overline{\varepsilon}\right)}{3 - 1} \tag{7}$$

where *ts* is the vibration displacement of node *t* in mode *m*. *t* is the average value of the vibration displacement of node *t*. *es* is the vibration displacement of node *e* of mode *m*. *e* is the average value of the vibration displacement of node *e*.

Similarly, Equations (8) and (9) respectively describes the covariance of node *t* and node *e* of mode *m* in *x*, *y*, and *z* directions at the rail web and rail head.

$$RW\\_COV\_m(t,\varepsilon) = \frac{\sum\_{s=1}^{3} \left(t\_s - \overline{t}\right) \left(\varepsilon\_s - \overline{\varepsilon}\right)}{3 - 1} \tag{8}$$

$$RH\\_COV\_m(t,\varepsilon) = \frac{\sum\_{s=1}^{3} \left(t\_s - \overline{t}\right) \left(\varepsilon\_s - \overline{\varepsilon}\right)}{3 - 1} \tag{9}$$

The mode shape data of mode 1, mode 3, and mode 7 is substituted into Equations (7)~(9), respectively. The optimal excitation nodes of each mode are calculated and the excitation method is selected according to mode shape characteristics. Figure 6 shows the optimal excitation nodes.

**Figure 6.** Excitation nodes of modes. (**a**) Excitation node of mode 1; (**b**) excitation nodes of mode 3; (**c**) excitation nodes of mode 7.

In summary, mode 1 selects the two-side symmetric excitation, which is excited on node 49 and node 33 in *z* positive direction. Node 33 receives the signal in *z* direction. Mode 3 selects the one-side symmetric excitation on node 73 and node 74 in *y* positive direction. Node 74 is selected to receive the signal in *y* direction. Mode 7 selects the symmetric excitation on nodes 43, node 39, node 10 and node 23 in *z* direction. Then node 43 receives the signal in *z* direction. The detailed selection process of excitation node can be referred to our previous work [20].

#### **4. Mode Excitation Simulation and Verification**

Mode excitation is simulated with ANSYS to verify the above proposed method. The excitation signal is a sine wave modulated by the Hanning window with a frequency of 35 kHz as shown in Figure 7.

**Figure 7.** Excitation signal.

In the simulation process, the excitation signal is applied on the rail at 4 m from the right end of the rail as shown in Figure 8 in order to avoid the influence of the rail end echo on the simulation results. Between 0.8 and 2.7 m from the excitation node, a set of data acquisition array is set at intervals of 5 mm. There are 380 data acquisition nodes.

The two-dimensional Fast Fourier Transform(2D-FFT) is used on the acquired data to identify the phase velocity of the mode. The simulation results are compared with the frequency wavenumber dispersion curves calculated by SAFE method as shown in Figure 9.

**Figure 8.** *Cont*.

**Figure 8.** Node excitation simulation diagram. (**a**) Node excitation simulation diagram of mode 7; (**b**) node excitation simulation diagram of mode 3; (**c**) node excitation simulation diagram of mode 1.

**Figure 9.** *Cont*.

**Figure 9.** Frequency wavenumber curves. (**a**) 2D-FFT result of mode 7; (**b**) 2D-FFT result of mode 3; (**c**) 2D-FFT result of mode 1.

In Figure 9, the result of 2D-FFT coincides well with the frequency wavenumber curve of the corresponding mode. The ordinate value of the wavenumber *ξ* can be acquired from Figure 9 when the abscissa *f* is 35 kHz. The relationship between wavenumber *ξ* and frequency *f* is shown in Equation (10):

$$C\_p = \frac{2\pi f}{\xi} \tag{10}$$

Then the frequency *f* and wavenumber *ξ* are substituted into the Equation (10) to calculate the simulation phase velocity of each mode. The calculation results are shown in Table 2.


**Table 2.** Phase velocity of modes.

It is found from Table 2 that the phase velocity errors of mode 7, mode 3 and mode 1 are 2.49%, 0.16%, 0.07% respectively. It is proved that the method of node excitation is effective, but there are still some interference modes. In order to excite a relatively single mode, it is necessary to study the method of mode enhancement and interference mode suppression.

#### **5. Mode Enhancement and Defect Detection Simulation**

#### *5.1. Mode Enhancement Theory Based on Phase Control and Time Delay Technology*

According to the theory of Rose [21], the installation interval and the excitation sequence of transducer array are set based on the expected guided wave mode, and the excited guided waves have the same phase angle. The excitation time interval is set as periodic *T*. According to the principle of wave superposition, the synthetic amplitude of two coherent waves with the same phase is the largest when the interval of transducers is offered based on the expected mode wavelength. As shown in Figure 10, 5 transducers are installed on the rail.

**Figure 10.** Phased array excitation diagram.

Transducer 0# is an excitation signal without time delay. The displacement function of each particle is *u*0(*x*, *y*, *z*, *t*).

$$u\_0(\mathbf{x}, \mathbf{y}, z, t) = \begin{bmatrix} u\_x(\mathbf{x}, \mathbf{y}, z, t) \\ u\_y(\mathbf{x}, \mathbf{y}, z, t) \\ u\_z(\mathbf{x}, \mathbf{y}, z, t) \end{bmatrix} = \begin{bmatrix} \mathcal{U}\_x(\mathbf{y}, z) \\ \mathcal{U}\_y(\mathbf{y}, z) \\ \mathcal{U}\_z(\mathbf{y}, z) \end{bmatrix} e^{i(\xi \mathbf{x} - \omega t)} \tag{11}$$

where *ξ* is the wavenumber.

The distance between transducer 1# and transducer 0# is *λ*, and the time delay is *T*. The displacement function of each particle is *u*1(*x*, *y*, *z*, *t*).

$$u\_{\mathbb{L}}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z},t) = \begin{bmatrix} u\_{\mathbb{L}}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z},t) \\ u\_{\mathbb{Y}}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z},t) \\ u\_{\mathbb{Z}}(\boldsymbol{x},\boldsymbol{y},\boldsymbol{z},t) \end{bmatrix} = \begin{bmatrix} \mathcal{U}\_{\mathbb{L}}(\boldsymbol{y},\boldsymbol{z}) \\ \mathcal{U}\_{\mathbb{Y}}(\boldsymbol{y},\boldsymbol{z}) \\ \mathcal{U}\_{\mathbb{Z}}(\boldsymbol{y},\boldsymbol{z}) \end{bmatrix} e^{\int\_{\mathbb{R}}^{\mathbb{L}}(\boldsymbol{x}-\boldsymbol{\Lambda}\boldsymbol{\!}) - \omega(t-\boldsymbol{\Lambda}t)]} = \begin{bmatrix} \mathcal{U}\_{\mathbb{L}}(\boldsymbol{y},\boldsymbol{z}) \\ \mathcal{U}\_{\mathbb{Y}}(\boldsymbol{y},\boldsymbol{z}) \\ \mathcal{U}\_{\mathbb{Z}}(\boldsymbol{y},\boldsymbol{z}) \end{bmatrix} e^{\int\_{\mathbb{R}}^{\mathbb{L}}(\boldsymbol{x}-\boldsymbol{\lambda}) - \omega(t-\boldsymbol{\Lambda})} \tag{12}$$

As we know:

$$
\xi \lambda = 2\pi \tag{13}
$$

$$
\omega T = 2\pi \tag{14}
$$

Substitute the Equations (13) and (14) into the Equation (12) to obtain the Equation (15).

$$u\_1(\mathbf{x}, y, z, t) = \begin{bmatrix} u\_x(x, y, z, t) \\ u\_\mathcal{Y}(x, y, z, t) \\ u\_\mathcal{Z}(x, y, z, t) \end{bmatrix} = \begin{bmatrix} \mathcal{U}\_x(y, z) \\ \mathcal{U}\_\mathcal{Y}(y, z) \\ \mathcal{U}\_\mathcal{Z}(y, z) \end{bmatrix} e^{i[\frac{\mathcal{Y}}{\mathcal{Y}}(x - \lambda) - \omega(t - T)]} = \begin{bmatrix} \mathcal{U}\_x(y, z) \\ \mathcal{U}\_\mathcal{Y}(y, z) \\ \mathcal{U}\_\mathcal{Z}(y, z) \end{bmatrix} e^{i[\frac{\mathcal{Y}}{\mathcal{Y}} - \omega t]} = \mathfrak{u}(x, y, z, t) \tag{15}$$

Similarly, the displacement and phase of transducer array 0#, 1#, 2#, 3#, and 4# can be derived as Equation (16):

$$\mathfrak{u}\_0(\mathbf{x}, y, z, t) = \mathfrak{u}\_1(\mathbf{x}, y, z, t) = \mathfrak{u}\_2(\mathbf{x}, y, z, t) = \mathfrak{u}\_3(\mathbf{x}, y, z, t) = \mathfrak{u}\_4(\mathbf{x}, y, z, t) \tag{16}$$

When the guided wave propagates along x direction, the expected mode is enhanced and the interference mode is suppressed. The transducer intervals of each mode are shown in Table 3.


**Table 3.** Transducer intervals.

#### *5.2. Mode Enhancement Simulation Results*

The simulation process of mode enhancement is the same as Section 3, only the excitation method is changed to transducer array. The simulation results are compared with the frequency wavenumber dispersion curves as shown in Figure 11.

**Figure 11.** Frequency wavenumber curves. (**a**) 2D-FFT result of mode 7; (**b**) 2D-FFT result of mode 3; (**c**) 2D-FFT result of mode 1.

From Figure 11 we can see, the result of 2D-FFT coincides well with the frequency wavenumber curve of the corresponding mode and there is no interference mode. The wavenumber *ξ* and frequency *f* are substituted into Equation (10). The calculation results of simulated phase velocities are shown in Table 4, and compared with theoretical values. It can be seen that the maximum error is 0.68%.


**Table 4.** Phase velocities of modes.

The above discussion is about phase velocity, then the process of the group velocity calculation is described as following. The simulation data at 1 m and 2.2 m from the excitation nodes are extracted and the simulation results of mode 7, mode 3, and mode 1 are shown in Figure 12.

**Figure 12.** *Cont*.

**Figure 12.** Mode simulation results. (**a**) Mode 7 simulation result; (**b**) Mode 3 simulation result; (**c**) Mode 1 simulation result.

The simulation group velocity can be calculated by calculating the peak time difference between 1 m and 2.2 m and the results are shown in Table 5.


**Table 5.** Group velocities of modes.

As can be seen from Table 5, the group velocity error of each mode is relatively small. The results show that the method of array excitation can enhance the expected mode and suppress the interference mode, which indicates a relatively single mode is successfully excited.

#### *5.3. Defect Detection Simulation*

At present, rail head transverse defects are the most dangerous damage for rail operation. Figure 13a presents the simulation of detecting an internal transverse defect. The shape of the defect on the rail cross section is a circle with a diameter of 23 mm and a length of 5 mm in x direction. In Figure 13b, the rail web defect is assumed as a vertical crack. The depth of the crack in z direction is 30 mm and 5 mm in x direction. In Figure 13c, the rail base defect is a transverse defect, the length of which is 25 mm in the y direction and 10 mm in x direction. In defect detection simulation, the modes excitation and data acquisition method are the same as Section 5.2 and the rail defects are 3 m from the left end of the rail.

**Figure 13.** Defect detection simulation diagram. (**a**) diagram of array excitation mode 7 defect detection; (**b**) diagram of array excitation mode 3 defect detection; (**c**) diagram of array excitation mode 1 defect detection.

It is necessary to verify whether the received defect echo at the receiving nodes has occurred in mode conversion before calculating the defect location. The defect echoes extracted separately are transformed by 2D-FFT as shown in Figure 14.

**Figure 14.** Frequency wavenumber curves. (**a**) mode 1 defect reflected echo 2D-FFT results; (**b**) mode 3 defect reflected echo 2D-FFT results; (**c**) mode 7 defect reflected echo 2D-FFT results.

Similarly, the phase velocity of each mode can be calculated from the frequency wavenumber curve. Table 6 shows the phase velocities of reflected echo modes.


**Table 6.** Phase velocities of reflected echo modes.

As can be seen from Table 6 that the error between the simulated phase velocity and the theoretical phase velocity of each defect echo can be accepted for application. The reflected echoes received at the corresponding nodes have not occurred mode conversion.

During the process of simulation, the receiving node is 1 m away from the defect. The position of the rail defect can be obtained by calculating the time difference between the received excitation signal and the reflected echo at the receiving node. As shown in Figure 15, the defect is 3 m from the left end of the rail. The receiving node is 1 m from the defect and 2 m from the excitation signal.

**Figure 15.** Diagram of the rail head defect detection.

When the defect is in the rail head, the detection results of mode 7, mode 3 and mode 1 are as shown in Figure 16. Only mode 7 has reflected echo. The time difference between the received excitation signal and the reflected echo is 0.6137 ms.

(**a**)

**Figure 16.** *Cont*.

**Figure 16.** Rail head defect detection. (**a**) Mode 7 rail head defect detection; (**b**) mode 3 rail head defect detection; (**c**) mode 1 rail head defect detection.

When the defect is in the rail web, only mode 3 has reflected echo as shown in Figure 17 and the time difference is 0.6478 ms.

**Figure 17.** *Cont*.

**Figure 17.** Rail web defect detection. (**a**) Mode 7 rail web defect detection; (**b**) mode 3 rail web defect detection; (**c**) mode 1 rail web defect detection.

(**c**)

Similarly, when the defect is in the rail base, only mode 1 has reflected echo as shown in Figure 18 and the time difference between is 0.6913 ms.

**Figure 18.** *Cont*.

(**c**)

**Figure 18.** Rail base defect detection. (**a**) Mode 7 rail base defect detection; (**b**) mode 3 rail base defect detection; (**c**) mode 1 rail base defect detection.

The defect position can be calculated based on the time difference. Table 7 shows the detect results of the defect position. The results can be seen from the Table 7 that the defect positioning error of mode 7, mode 3 and mode 1 is 1.34%, 2.19%, and 1.40% respectively.



#### **6. Conclusions**

In this paper, a novel method for mode selection and excitation in rail defect detection is proposed based on ultrasonic guided waves. The mode shape data in the CHN60 rail is obtained at the frequency of 35 kHz by using SAFE method. And the modes with good non-dispersion characteristics, which

are respectively sensitive to the defects of the rail head, rail web and rail base, are selected combining the strain energy distribution diagrams with group velocity dispersion curves. The optimal excitation direction of the expected modes is solved by the Euclidean distance, and the optimal excitation node is calculated by covariance matrix. Then, phase control and time delay technology is applied to achieve the enhancement of expected modes and suppression of interferential modes. After that, ANSYS is used to simulate the expected modes excitation and defects detection in the different rail parts with phase velocity and group velocity to validate the proposed methods. The simulation results present an acceptable error and demonstrate the effectiveness of the guided wave mode selection and excitation method, which proposes a feasible scheme for the field test and application to rail defect detection.

**Author Contributions:** This paper is a result of the full collaboration of all the authors; H.S. and L.Z. (Lu Zhuang) conceived and proposed the main idea of the paper; X.X. combined with phase control and time delay technology to improve the mode excitation method; L.Z. (Liang Zhu) simulated the defect detection and analyzed the data; L.Z. (Liang Zhu) and Z.Y. put forward their views on data processing; H.S. and L.Z. (Lu Zhuang) wrote the paper.

**Acknowledgments:** The work is supported by the National Key Research and Development Program of China (2016YFB1200401), Supported by Foundation of Key Laboratory of Vehicle Advanced Manufacturing, Measuring, and Control Technology (Beijing Jiaotong University), Ministry of Education, China.

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

#### **References**


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

### *Article* **Development of a Vibroacoustic Stochastic Finite Element Prediction Tool for a CLT Floor**

#### **Cheng Qian 1,\*, Sylvain Ménard 1, Delphine Bard <sup>2</sup> and Juan Negreira <sup>3</sup>**


Received: 31 January 2019; Accepted: 8 March 2019; Published: 15 March 2019

**Abstract:** Low frequency impact sound insulation is a challenging task in wooden buildings. Low frequency prediction tools are needed to access the dynamic behavior of a wooden floor in an early design phase to ultimately reduce the low frequency impact noise. However, due to the complexity of wood and different structural details, accurate vibration predictions of wood structures are difficult to attain. Meanwhile, a deterministic model cannot properly represent the real case due to the uncertainties coming from the material properties and geometrical changes. The stochastic approach introduced in this paper aims at quantifying the uncertainties induced by material properties and proposing an alternative calibration method to obtain a relative accurate result instead of the conventional manual calibration. In addition, 100 simulations were calculated in different excitation positions to assess the uncertainties induced by material properties of cross-laminated-timber A comparison between the simulated and measured results was made in order to extract the best combination of Young's moduli and shear moduli in different directions of the CLT panel.

**Keywords:** wooden constructions; acoustics; low frequency noise; modelling

#### **1. Introduction**

Multi-story wooden constructions have increased in the market during the last two decades. However, this kind of constructions still face the challenge coming from low-frequency sound insulation. Although these types of buildings comply with the present regulations, subjective ratings of inhabitants have shown complaints [1–6] due to low-frequency impact noise. One of the important reasons resulting in this discrepancy between the subjective and objective evaluations is that the current standards are initially designed for heavy constructions, i.e., concrete constructions, but without an appropriate modification, they are directly applied to the wooden constructions. On the other hand, the low frequency impact noise is neglected in the standardized impact sound insulation evaluations. From recent research [3], it was shown that the subjective ratings are correlated better with the objective ratings with the help of the adaptation term *CI*, 50−<sup>3000</sup> following the standard ISO 717-2 [7]. But the correlation between the inhabitants' satisfaction (reported by means of questionnaires) and objective ratings can be further improved by taking the measurement from 20 Hz in comparison with the measurements performed from 50 Hz instead. Furthermore, the dynamic behavior of wood (as a natural material) is hard to predict due to its inhomogeneity. As a consequence, product development in the wooden industry is still based on empirical models and experimental tests, which could lead to an over-designed and expensive acoustical improvement solution. To end that, more knowledge

about wooden construction is needed and accurate and handy prediction tools are called for in order to enable acoustic comfort in wooden constructions, especially in a low-frequency range.

The finite element (FE) method is a widely employed approach to develop numerical prediction models in wooden industry. By performing numerical simulations, experimental acoustical tests can be reduced, and parametric studies can also be carried out to investigate the influence of specific geometrical changes in construction as well as the influence of the variations/uncertainties in material properties, which always have a markable influence on the results. In Reference [8], experimental tests were conducted on a full-scale cross-laminated timber (CLT) floor. The material properties of the CLT were collected from the literature and then put into the established FE model to compare the simulation results with the measured ones. A better correlation between the testing results and the modelling results was attained after tuning the collected material properties of CLT. The latter points out the importance of knowing the material properties if a proper calculation of dynamic properties needs to be achieved. A similar conclusion was drawn in References [9–11].

However, one should know that the variety of the wood species, the variation in theoretical identical wood structure, as well as the lack of material properties data base, make the reliable material properties of wood difficult to obtain. Consequently, the calibration of a model is always tedious and time-consuming. Moreover, even with the calibrated model to predict the theoretical identical wooden structure, the prediction results may be different from the realistic case due to the uncertainties induced by wooden material properties, workmanship and geometrical details in structures. As a result, the deterministic model may not be representative enough in a realistic situation. Therefore, it is a necessary step to quantify the variations in the model's outputs in order to ultimately establish an accurate prediction tool. To achieve that, uncertainties of the wood material properties can be addressed in a model by introducing probabilistic parameters [6]. In Reference [12], a generalized probabilistic model was constructed to take into account the statistical fluctuation associated with the elastic properties in the model. The uncertainties of the mechanical constants of a wooden structure were also investigated in Reference [13]. A small data base of the Young's modulus in the longitudinal direction was established by means of vibration measurements. Then, by sampling the Young's modulus distribution established by the obtained data base, Monte Carlo simulations were performed in the FE model development to quantify the uncertainties induced by the elastic constants in modelling the vibro-acoustic behavior of wooden buildings in a low-frequency range. However, not only does the Young's modulus in the longitudinal direction have a great impact on the dynamic response of a wooden structure, but the Young's moduli in the other two directions as well as the shear moduli play an important role in its vibroacoustic performance. Regarding the uncertainties in structural dynamics, Shannon's maximum entropy principle [14] is an optimal choice to model the random data and the uncertainties [15]. In Reference [16,17], a probabilistic model was proposed to construct the probability distribution in high-dimension of a vector-valued random variable using the maximum entropy principle. This proposed probabilistic model was then extended to a tensor-valued coefficient (stiffness tensor) with different symmetric levels [15,18,19].

The research reported in this case aims at developing a stochastic FE prediction tool of a CLT slab in order to quantify the uncertainties induced by material properties and, subsequently, acquire the reliable material properties of the under investigated structure to accurately predict the vibroacoustic behavior of CLT in a low-frequency range. To achieve that, different mechanical constant distributions without available material property data base are constructed by means of the maximum entropy principle. The random elastic constants generated following the established distributions are considered as the inputs for the FE model to calibrate CLT in a low frequency range by comparing the simulated results with its dynamic responses obtained from experimental modal analysis (EMA) method. The best combination of the material properties of CLT panel is selected by minimizing different error metrics. The ultimate objective is to acquire the knowledge about the modeling method of CLT and to propose a generalized method to quantify the uncertainties coming from the orthotropic

level material properties and to accurately calibrate the corresponding model by avoiding the repetitive manual calibration and, subsequently, to increase the prediction accuracy.

#### **2. Measurements**

#### *2.1. Experimental Structure*

The test structure was a 4 × 1.5 m<sup>2</sup> 5-ply 175 mm thick CLT panel made of Canadian black spruce of machine stress rated grade 1950f-1.7E in parallel layers and visual grade No. 3 in perpendicular layers with density of 520 kg/m3. At the initial stage, two CLT slabs were nailed together by a thin wooden panel (c.f. Figure 1a) and they were placed on a standardized sound insulation measurement console. However, this weak connection between two CLT panels resulted in a low stiffness coupling between two slabs and the imperfect simply supported boundary conditions of the measurement console makes the CLT floor easily "jump up" when there is an external excitation. In a real building construction, these problems can be resolved by adding a top layer on the CLT bare floor, which can add mass on CLT and, subsequently, enforce the connection between two CLT leaves. Consequently, due to the weak connection between two assembled CLT panels and the imperfect boundary conditions in the sound insulation laboratory, only the first eigen-mode can be extracted from the measurements, which cannot provide enough information to calibrate the established FE model. As a result, in order to easily extract more meaningful information, only one leaf of CLT floor was under investigation.

**Figure 1.** (**a**) Two CLT panels connected by a thin wooden lath. (**b**) One CLT panel simply supported on two I-steel beams.

#### *2.2. Measurement Procedure*

#### 2.2.1. Measurement Setup

EMA was performed on the CLT slab to characterize the dynamic properties of the CLT slab. The CLT panel was set free at two opposite long sides and simply supported (SFSF boundary condition) with two steel I-beams at two shorter edges (shown in Figure 1b) in order to simplify the boundary conditions. A predefined mesh was drawn on the CLT surface to determine the excitation positions and the measurement positions. In that way, the CLT was divided into five parts in the long direction and three parts in the short direction, which gives a total number of 16 excitation points (the nodes on the short edges were not excited). The size of mesh was decided based upon the highest frequency of interest, i.e., 200 Hz, to avoid the hammer exciting on the nodes of the eigen-modes.

A single input multiple output system (SIMO) can provide redundant data to better identify the eigen-frequencies and the local modes of complex structures based on different frequency response function (FRF) matrices [20,21]. Meanwhile, these different FRFs can be used as references to validate the FE model by comparing different simulated FRFs with different measured FRFs. Concerning the complexity and the inhomogeneity of CLT, the SIMO system was employed to characterize the dynamic properties of the CLT panel via different FRF matrices at different measurement positions. Five uni-axial accelerometers (Brüel & Kjær Accelerometer Type 4507 001, Brüel & Kjær Sound & Vibration Measurement A/S, Nærum, Denmark) were attached on the different nodes (red dots in Figure 2) with the help of Faber-Castell Tack-It Multipurpose Adhesive. An impact roving hammer (Brüel & Kjær Impact hammer Type 8208 serial No. 51994, Brüel & Kjær Sound & Vibration Measurement A/S, Nærum, Denmark) was used as excitation.

**Figure 2.** Mesh on the CLT panel. Five accelerometers (red dots) were placed at points 10, 11, 13, 17, and 24), whereas the hammer was moved around all nodes, except on the short edges.

Depending on the impulse shape and the force spectrum of the impact hammer shown in Figure 3, the medium hard hammer tip was selected for this measurement to ensure most of the input energy localizing within the frequency range of interest (up to 200 Hz). The roving hammer approach was employed. The accelerometers were kept fixing on the selected positions and the hammer were roved over the predefined mesh grid except for the points on boundaries. The accelerometers and the impact hammer were connected to Brüel & Kjær multi-channel front-end LAN-XI Type 3050. Then, the acceleration signals and the impact force signals were recorded by using the software Brüel & Kjær PULSE Labshop. There were three averages for each excitation position. The resolution of the FRF was 0.625 Hz.

**Figure 3.** (**a**) Impulse shapes of the hammers showing the shape as a function of the used impact tip. (**b**) Force spectra of the hammers showing the frequency response as a function of used impact tip [22] Copyright © Brüel & Kjær.

#### 2.2.2. Modal Parameters Extraction

The eigen-frequencies, the mode shapes, and the damping ratios were extracted by applying rational fraction polynomial—Z algorithm in Brüel & Kjær PULSE Connect. The frequency band of FRFs was divided into several parts, which means that identification of poles was restricted in a narrow frequency band each time in order to easily select the stable poles (eigen-frequencies) from the FRFs in each divided frequency range. An example of the stabilization diagram is shown in Figure 4.

**Figure 4.** Stabilization diagram with rational fraction polynomial—Z algorithm.

Iteration times of rational fraction polynomial—Z algorithm was 40. The synthesized FRF is shown in Figure 5. The corresponding eigen-frequencies and the corresponding damping ratios are shown in Table 1.

**Figure 5.** Synthesized FRF, measured FRF, and relative errors between the synthesized FRF and the measured FRF.


**Table 1.** Measured and simulated Eigen frequencies of the bending modes and the measured corresponding damping ratios.

Since the measurement system is the SIMO system. Different FRFs (measurement point 11/excitation point 11, measurement point 13/excitation point 13, measurement point 17/excitation point 17, measurement point 24/excitation point 24) can also be saved as references to calibrate the FE model.

#### **3. Finite-Element Model**

#### *3.1. Model Description*

A numerical model of the CLT slab was created in the commercial FE software Abaqus [23]. Five different layers were modelled by assigning different oriented material properties in different layers of the model. Therefore, there is no relative displacement between different layers in this model, which means that the CLT model is a complete model and each layer is fully tied together. The same material was assigned to each layer, except that the in-plan material orientation assignments were 90◦ oriented from the adjacent layer to mimic the cross-laminate layers of CLT panel. The material properties of CLT were gathered from the literature [24,25], reported in Table 2. The meshes of 20-node quadratic brick, reduced integration (C3D20R) quadratic type were assigned to the entire model. Since the shape of each layer is a simple rectangular so that there is no discontinuity between different layers. Details of mesh are shown in Figure 6. The mesh size was 0.1 to ensure the accuracy of the highest frequency interest. The eigen-frequencies and the eigen-modes were calculated by the linear perturbation frequency step. The FRF of the CLT was obtained by the Steady-state dynamics, Modal step. The damping extracted from the measurement was included in the model by means of the direct modal damping (c.f. Table 1). In this framework, the FRFs of CLT were first calculated to quantify the uncertainties of material properties. Afterward, the best FRF justified by different criterions was selected to extract the material properties of this under investigated CLT.

**Table 2.** Material Properties of CLT collected from the literature. Stiffness parameters have the unit of MPa, Poisson ratios are dimensionless, and the density is given in kg/m3.


In general, it is a challenge task to create appropriate constraints to describe the boundary conditions. Since dynamic responses of the structure are sensitive to the boundary conditions, slight changes in the FE model can lead to a big variation of eigen-frequencies and mode order. In order to mimic the simply supported boundary condition, all the displacements in three directions at one shorter edge were constrained and, at the other shorter edge, all the displacements were also constrained except the displacement in a vertical direction (U1 in Abaqus). This boundary set-up of the FE model was kept throughout the entire modelling process. More explanations about the choice of the boundary condition can be found in Section 3.3.

**Figure 6.** Meshes of FE CLT model.

#### *3.2. Model Validation Criterion*

The established model needs to be validated by two different metrics. One is the normalized relative frequency difference (NRFD), which characterizes the discrepancies between the simulated and the measured eigen-frequencies, defined by the equation below.

$$NRFD\_i = \frac{\left|f\_{ref\_i} - f\_i\right|}{f\_{ref\_i}},\tag{1}$$

where *frefi* is the measured resonance frequency and *fi* is the simulated resonance frequency.

The other one is the modal assurance criterion (MAC), which quantifies the similarity of simulated and measured mode shapes, defined by the formula below.

$$\text{MAC} = \frac{|\left(\Phi\_i^{sim}\right)^T \left(\Phi\_j^{exp}\right)|^2}{\left(\Phi\_i^{sim}\right)^T \left(\Phi\_i^{sim}\right) \left(\Phi\_j^{exp}\right)^T \left(\Phi\_j^{exp}\right)}\tag{2}$$

where Φ*sim <sup>i</sup>* is the *<sup>i</sup>*-th simulated mode shape and <sup>Φ</sup>*exp <sup>j</sup>* is the *j*-th measured mode shape. The range of the MAC number is from 0 to 1. When the MAC number equals 1, the simulated eigen-mode perfectly correlates to the measured one. However, MAC equals 0, which implies irrelevant simulated and measured eigen-modes. Different from the conventional calibration procedure, stochastic process calibration begins with comparing eigen-frequencies (NRFDs) between the simulated FRFs and measured FRFs. Therefore, several FRFs of different positions are needed to make sure the calibrated model validating at different positions, which can increase the credibility of the model. However, the FRFs calibration can only ensure consistency of simulated and measured eigen-frequency but not the mode order. The simulated eigen-frequency may be the same as the measured eigen-frequency, but they could have different mode shapes. Therefore, both indictors are needed to ensure the simulated eigen-frequencies correlated with the measured ones to keep the simulated mode order corresponding to the measured one.

#### *3.3. Preliminary Sensitivity Analysis*

Wood as a kind of orthotropic material has nine different variables (three Young's moduli, three shear moduli, and three Poisson's ratios) to be calibrated. In order to decrease the complexity of calibration, a sensitivity analysis should be performed to investigate the effect of different elastic constants on simulated eigen-frequencies before the stochastic process is introduced to the FE model. In this section, Young's moduli in different directions were increased or decreased 25% when compared to the Young's moduli given by Table 2. Shear moduli in different directions were increased or

decreased 15% when compared to the shear moduli given by Table 2. The Poisson ratios, *ν*12, *ν*13, *ν*<sup>23</sup> were set to be 0.25, 0.25, and 0.35, whereas the initial values were 0.3, 0.3, and 0.4. The reference eigen-frequencies employed in NRFDs were calculated, according to the elastic constants reported in Table 2. The measured eigen-frequencies were not selected as a reference since the objective of sensitivity analysis aims to investigate how different elastic constants affect the eigen-frequencies of an FE model.

From the NRFDs shown in Figures 7 and 8, it can be seen that the influence of Young's moduli and shear moduli on eigen-frequencies cannot be ignored. Among all the elastic constants, Young's modulus in the longitudinal direction has the most important influence on eigen-frequency changes. However, Young's modulus in a vertical direction barely changes the eigen-frequencies. Therefore, the variation of Young's modulus in the vertical direction was not reported in Figure 7. When looking at Figure 9, we could find that all the NRFDs are lower than 0.5%, which indicates that the influences of Poisson's ratios on eigen-frequencies are negligible. From this sensitivity analysis, it can be concluded that the Young's moduli and shear moduli have a more significant influence on eigen-frequency calculations than Poisson's ratios. As a result, the calibration of material properties can be reduced to six variables.

**Figure 7.** NRFDs of Young's moduli.

**Figure 8.** NRFDs of shear moduli.

**Figure 9.** NRFDs of Poisson's ratios.

However, from this sensitivity analysis, it was also noticed that the first four simulated eigen-frequencies are always higher than the first four measured eigen-frequencies even decreasing different elastic constants. This may be caused by the over-stiffened boundary condition. Since the tested CLT is only placed on the top two steel I beams. It is difficult to have a perfect simply supported boundary conditions in reality. Therefore, restricting all the displacements at the boundaries of the FE model can create over stiffened boundary conditions, which results in over-estimated eigen-frequencies. Consequently, the displacement in a vertical direction at one boundary of the FE model is released to try to mimic the real boundary conditions.

#### **4. Stochastic Process**

Uncertainties of material properties are always assumed to follow the Gaussian distribution because of its simplicity and the lack of relevant experimental data, even though most material property distributions are non-Gaussian in nature [26,27]. The theory introduced in this case is about the probabilistic modeling of a random elasticity tensor in an orthotropic symmetric level within the framework of the maximum entropy principle under the constraint of the available information [18,19,28]. The established random elasticity tensor is considered as the inputs in the FE model to quantify the uncertainties induced by the CLT material properties and to seek the best combination of CLT material properties to calibrate the CLT model.

In this section, the elastic tensor is first decomposed in terms of random coefficients and tensor basis so that the fluctuation of different elastic constants can be characterized by the probability distribution functions (PDF). Next, construction of PDFs of different elastic constants in high-dimension [16] is shortly introduced. Lagrange multipliers associated with the explicit PDFs of random variables in high dimensions is estimated with the help of the Itô differential equation. The established PDFs is sampled by Metropolis-Hastings algorithm to obtain the random data to construct a random elasticity matrix [18,19] in order to derive the corresponding random combinations of elasticity constants to quantify the uncertainties of material properties and to calibrate the CLT model. A flow chart of the application of the stochastic procedure is shown in Figure 10.

**Figure 10.** Flow chart of the application of the stochastic process.

#### *4.1. Decomposition of the Random Elastic Tensor*

Construction of PDF of a random value in a high-dimension approach can be applied to any arbitrary material symmetry class [15,18,19], such as isotropic symmetry, cubic symmetry, and transversely isotropic symmetry. In this work, this stochastic approach aims to seek a reasonable probability distribution of the random elastic tensor of the target material (CLT). Therefore, only the orthotropic symmetry case is considered in this section. The dimension of the random variable *N* is limited to 9.

Let *C* be a fourth-order random elastic tensor, which could be decomposed by using the equation below.

$$\mathbf{C} = \sum\_{i=1}^{N} c\_i \mathbf{E}\_{i\prime} \tag{3}$$

where *ci* is a set of random coefficients that can be described by its PDFs and *Ei* is the tensor basis of the random elastic tensor *C*, based on Walpole's derivation [29]. The CLT slab is modelled as orthotropic material in Abaqus so that the tensor basis *Ei* of the orthotropic symmetric elastic tensor are defined in the following form [18].

$$\begin{cases} \begin{aligned} \mathsf{E}^{11} &= \mathsf{a}\otimes\mathsf{a}\otimes\mathsf{a}\otimes\mathsf{a}, \mathsf{E}^{12} = \mathsf{a}\otimes\mathsf{a}\otimes\mathsf{b}\otimes\mathsf{b}, \mathsf{E}^{13} = \mathsf{a}\otimes\mathsf{a}\otimes\mathsf{c}\otimes\mathsf{c}, \\ \mathsf{E}^{21} &= \mathsf{b}\otimes\mathsf{b}\otimes\mathsf{a}\otimes\mathsf{a}, \mathsf{E}^{22} = \mathsf{b}\otimes\mathsf{b}\otimes\mathsf{b}\otimes\mathsf{b}, \mathsf{E}^{23} = \mathsf{b}\otimes\mathsf{b}\otimes\mathsf{c}\otimes\mathsf{c}, \\ \mathsf{E}^{31} &= \mathsf{c}\otimes\mathsf{c}\otimes\mathsf{a}\otimes\mathsf{a}, \mathsf{E}^{32} = \mathsf{c}\otimes\mathsf{c}\otimes\mathsf{b}\otimes\mathsf{b}, \mathsf{E}^{33} = \mathsf{c}\otimes\mathsf{c}\otimes\mathsf{c}\otimes\mathsf{c}, \\ &\quad \mathsf{E}^{4}\_{ijkl} = \left(a\_{i}b\_{j} + b\_{i}a\_{j}\right)\left(a\_{k}b\_{l} + b\_{k}a\_{l}\right)/2, \\ &\quad \mathsf{E}^{5}\_{ijkl} = \left(b\_{i}c\_{j} + c\_{i}b\_{j}\right)\left(b\_{k}c\_{l} + c\_{k}b\_{l}\right)/2, \\ &\quad \mathsf{E}^{6}\_{ijkl} = \left(c\_{i}a\_{j} + a\_{i}c\_{j}\right)\left(c\_{k}a\_{l} + a\_{k}c\_{l}\right)/2, \end{aligned} \end{cases} \tag{4}$$

where *a*, *b*, and *c* are the unit orthogonal vectors, ⊗ is the Kronecker product.

The fourth-order elastic tensor *C* is decomposed as:

$$\mathbf{C} = c\_1 \mathbf{E}^{\mathbf{1}\mathbf{1}} + c\_2 \mathbf{E}^{\mathbf{2}\mathbf{2}} + c\_3 \mathbf{E}^{\mathbf{3}\mathbf{3}} + c\_4 \left(\mathbf{E}^{\mathbf{12}} + \mathbf{E}^{\mathbf{21}}\right) + c\_5 \left(\mathbf{E}^{\mathbf{23}} + \mathbf{E}^{\mathbf{32}}\right) + c\_6 \left(\mathbf{E}^{\mathbf{31}} + \mathbf{E}^{\mathbf{13}}\right) \tag{5}$$
 
$$+ c\_7 \mathbf{E}^{\mathbf{4}} + c\_8 \mathbf{E}^{\mathbf{5}} + c\_9 \mathbf{E}^{\mathbf{6}}.\tag{6}$$

#### *4.2. Construction of Probability Distribution Function in High-Dimension Using the Maximum Entropy Principle*

The objective of this section is to establish the PDFs of the random coefficients *ci*, which control the statistical fluctuation of the fourth-order random elastic tensor. Let *c* = (*c*1,..., *cN*) be a vector in R*N*-valued second order random variable, which obeys certain unknown probability distribution *Pc*(*dc*) with the Lebesgue measure *dc* = *dc*<sup>1</sup> ... *dcN*. The element in vector *c* = (*c*1,..., *cN*) represent the random coefficient of the random elasticity matrix in the previous section (Equation (3)). The unknown probability distribution *Pc*(*dc*) of the vector *c* can be presented by a probability density function *pc*(*c*), which satisfies the following normalization condition.

$$\int p\_{\mathfrak{c}}(\mathfrak{c})d\mathfrak{c} = 1.\tag{6}$$

The Maximum Entropy Principle applied here aims to construct the unknown probability distribution *Pc*(*dc*) with the help of the available information. In this case, the probability density function (PDF) *pc* could be written using the equation below.

$$p\_{\mathfrak{C}} = \arg\max \mathcal{S}(p),\tag{7}$$

where the entropy *S*(*p*) is defined by the equation below.

$$S(p) = -\int p(\mathbf{c}) \log(p(\mathbf{c})) d\mathbf{c}.\tag{8}$$

In order to find an explicit probability density function *pc*(*c*), several constraints should be set as available information: (1) the mean values of the variables, (2) the log condition, and (3) the normalization condition.

$$E\{\mathcal{C}\} = \overline{\mathfrak{c}}, \text{ with } \overline{\mathfrak{c}} = (\overline{\mathfrak{c}}\_1, \dots, \overline{\mathfrak{c}}\_\mathcal{O}) \tag{9}$$

$$E\left\{\log\left(\det\left(\sum\_{i=1}^{N}c\_i E\right)\right)\right\} = \nu\_{\mathbb{C}'} \text{ with } |\nu\_{\mathbb{C}}| < +\infty \tag{10}$$

and

$$\int p\_{\mathfrak{c}}(\mathfrak{c})d\mathfrak{c} = 1.\tag{11}$$

Equation (9) indicates that the mean values of variables are supposed to be known and Equation (10) ensures that both the *C* and *C*−<sup>1</sup> are second-order random variables. This equation also creates the statistical dependence between the different random variables.

To optimize the problem defined by Equation (7), Lagrange multipliers associated with Equations (9)–(11) are introduced. Let *<sup>λ</sup>*<sup>0</sup> <sup>∈</sup> <sup>R</sup>+, *<sup>λ</sup>***<sup>1</sup>** <sup>∈</sup> <sup>A</sup>*λ*<sup>1</sup> , and *<sup>λ</sup>*<sup>2</sup> <sup>∈</sup> <sup>A</sup>*λ*<sup>2</sup> be Lagrange multipliers associated with *Appl. Sci.* **2019**, *9*, 1106

the constraints defined by Equations (9)–(11). It could be proven that the optimized Equation (7) could be written as [28]:

$$p\_{\mathcal{C}}(\mathbf{c}) = k\_0 \exp\{- < \lambda\_{\text{sol}}, \mathbf{g}(\mathbf{c}) >\_{\mathbb{R}^{N+1}} \}, \forall \mathbf{c} \in \mathbb{R}^N,\tag{12}$$

where *k*<sup>0</sup> = exp −*λ*<sup>0</sup> is the normalizing constant, the operator , is the Euclidean inner product, the *<sup>c</sup>* <sup>→</sup> *<sup>g</sup>*(*c*) is the mapping defined on *<sup>S</sup>* <sup>×</sup> <sup>R</sup>, with the values in <sup>R</sup>*N*<sup>+</sup>1. *<sup>g</sup>*(*c*) is defined by Equation (13) below.

$$\log(\mathfrak{c}) = (\mathfrak{c}, \mathfrak{q}(\mathfrak{c})), \text{ with } \mathfrak{q}(\mathfrak{c}) = \log\left(\det\left(\sum\_{i=1}^{N} c\_i E\right)\right),\tag{13}$$

where det(∑*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *ciE*) > 0 is the support of Equation (13).

An R*N*-valued random variable *<sup>B</sup><sup>λ</sup>* parameterized by *<sup>λ</sup>* should be introduced to identify the Lagrange multipliers. Supposing the probability density function *b* → *pB<sup>λ</sup>* (*b*, *λ*) of the random variable *B<sup>λ</sup>* is written as:

$$\mathfrak{p}\_{\mathbb{B}\_{\lambda}}(\mathfrak{b},\lambda) = k\_{\lambda} \exp \{- < \lambda, \mathfrak{g}(\mathfrak{b}) >\_{\mathbb{R}^{N+1}} \}, \forall \mathfrak{b} \in \mathbb{R}^{N}, \tag{14}$$

where *k<sup>λ</sup>* is the normalization constant parameterized by *λ*. Taking *k*<sup>0</sup> = *kλ*, from Equations (12) and (14), it can be deduced that:

$$p\_{\mathcal{C}}(\mathfrak{c}) = p\_{B\_{\lambda}}(\mathfrak{b}, \lambda). \tag{15}$$

According to Equations (9), (10), (13), and (15), the calculation of Lagrange multipliers can be derived by evaluating the expectation of *g*(*bλ*):

$$E\{\mathbf{g}(\mathbf{b}\_{\lambda})\} = (\overline{\mathbf{c}}, \boldsymbol{\nu}\_{\mathbb{C}}).\tag{16}$$

As a result, the problem of the calculation of the Lagrange multipliers converts into generating the independent realizations of the random variable *<sup>B</sup><sup>λ</sup>* defined over R*<sup>N</sup>* and then evaluating the left-hand side of Equation (16).

#### 4.2.1. Calculation of Lagrange Multipliers

To derive Lagrange multipliers introduced in the previous section, there are several different methods to generate the independent the random variable *B<sup>λ</sup>* with respect to the corresponding probability density function (Equation (15)), such as the Metropolis-Hastings method [30,31], Gibbs method [32]. However, it should be noticed that the Metropolis-Hastings method demands an appropriate proposal distribution, which is sometimes difficult to choose, and the Gibbs method requires us to know the conditional distributions. As a consequence, it could be intricate to give a robust calculation without an adequate initial guess, especially for a high-dimension case. Therefore, another alternative algorithm is introduced in Reference [16] to generate the independent random variable *Bλ*.

#### Random Number Generator

Let *u* → Φ(*u*, *λ*) be the potential function defined as:

$$\Phi(\mathfrak{u}, \lambda) = <\lambda, \mathfrak{g}(\mathfrak{u})>\_{\mathbb{R}^{N+1}}.\tag{17}$$

Let {(*U*(*r*),*V*(*r*))},*<sup>r</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> be the Markov stochastic process defined on the probability space (Θ, *<sup>T</sup>*, *<sup>P</sup>*) indexed by <sup>R</sup><sup>+</sup> with values in <sup>R</sup><sup>+</sup> <sup>×</sup> <sup>R</sup>+, for *<sup>r</sup>* <sup>&</sup>gt; 0, satisfying the following It<sup>ô</sup> stochastic differential equation (ISDE).

$$\begin{cases} \begin{array}{c} d\mathcal{U}(r) = \mathcal{V}(r)dr, \\ d\mathcal{V}(r) = -\nabla\_{\mathcal{U}}\Phi(\mathcal{U}(r),\lambda)dr - \frac{1}{2}f\_{0}\mathcal{V}(r)dr + \sqrt{f\_{0}}d\mathcal{W}(r), \end{array} \tag{18}$$

where *W*(*r*) is the normalized Wiener process defined on (Θ, *T*,) indexed by R<sup>+</sup> and with values in R*N*. The probability distribution of the initial condition *U*(0) and *V*(0) are supposed to be known. The parameter *f*<sup>0</sup> is a real positive number, which could dissipate the transition part of the response generated by the initial condition and ensures a reasonable fast convergence of the stationary solution corresponding to the invariant measure.

When *r* tends to infinity, the solution *U*(*r*) of ISDE converges to the probability distribution of the random variable *Bλ*.

$$\lim\_{r \to \infty} (\mathcal{U}(r)) = \mathcal{B}\_{\lambda}. \tag{19}$$

*ns* independent realization of *B<sup>λ</sup>* is denoted as *Bλ*(*θ*1), ... , *Bλ*(*θns* ). Let *ro* be the iteration step that the solution of ISDE converge. When *rk* ≥ *r*0, the ISDE could be written as the equation below.

$$\begin{cases} \begin{array}{c} d\mathbf{U}(r\_k) = \mathbf{V}(r\_k) dr\_k, \\ d\mathbf{V}(r\_k) = -\nabla\_\mathbf{u} \Phi(\mathbf{U}(r\_k), \lambda) dr\_k - \frac{1}{2} f\_0 \mathbf{V}(r\_k) dr\_k + \sqrt{f\_0} d\mathbf{W}(r\_k). \end{array} \end{cases} \tag{20}$$

Therefore, the independent realizations *Bλ*(*θl*) can be presented by the stationary solution of ISDE.

$$\mathcal{B}\_{\lambda}(\theta\_l) = \mathcal{U}(r\_k, \theta\_l). \tag{21}$$

It is worthwhile to mention that the Itô stochastic differential equation defined by Equation (18) can be discretized by the Explicit Euler scheme to obtain an approximate solution.

$$k = 1, \dots, M - 1, \left\{ \begin{array}{c} \mathbf{U}^{k+1} = \mathbf{U}^{k} + \Delta r \mathbf{V}^{k}, \\ \mathbf{V}^{k+1} = \left( 1 - \frac{f\_0}{2} \Delta r \right) \mathbf{V}^{k} + \Delta r \mathbf{L}^{k} + \sqrt{f\_0} \Delta \mathbf{W}^{k+1}, \end{array} \tag{22}$$

with the initial conditions:

$$\mathbf{U}^1 = \mathfrak{u}\_{0\prime} \ \mathbf{V}^1 = \mathfrak{v}\_{0\prime} \tag{23}$$

where <sup>Δ</sup>*<sup>r</sup>* is the iteration step and <sup>Δ</sup>*Wk*+<sup>1</sup> <sup>=</sup> *<sup>W</sup>k*+<sup>1</sup> <sup>−</sup>*W<sup>k</sup>* is a second-order Gaussian centered <sup>R</sup>*N*-valued random variable with a covariance matrix *E* Δ*Wk*+<sup>1</sup> Δ*Wk*+<sup>1</sup> *T* <sup>=</sup> <sup>Δ</sup>*r*{*IN*}, where *<sup>W</sup>*<sup>1</sup> <sup>=</sup> <sup>0</sup>*N*. In Equation (22), *L<sup>k</sup>* is an R*N*-valued random variable, which is the partial derivative of Φ(*u*, *λ*) defined by the equation below.

$$L^k\_j \cong -\frac{\Phi\left(\Delta \mathcal{U}^{k,j}, \lambda\right) - \Phi\left(\mathcal{U}^k, \lambda\right)}{\mathcal{U}^{k+1}\_j - \mathcal{U}^k\_j},\tag{24}$$

with

$$
\Delta \mathbf{U}^{k,j} = \left( \mathcal{U}^{k}\_{1}, \dots, \mathcal{U}^{k}\_{j-1}, \mathcal{U}^{k}\_{j} + \Delta \mathcal{U}^{k+1}\_{j}, \mathcal{U}^{k}\_{j+1}, \dots, \mathcal{U}^{k}\_{N} \right), \\
\Delta \mathcal{U}^{k+1}\_{j} = \mathcal{U}^{k+1}\_{j} - \mathcal{U}^{k}\_{j}. \tag{25}
$$

#### Mathematical Expectation Estimation

After the random number generator has been established and the ISDE has been discretized to obtain the random numbers, the expectation of these independent random numbers should be calculated to derive Lagrange multipliers. The mathematical expectation of the random variable *B<sup>λ</sup>* can be estimated by using the Monte Carlo method. The evaluation of the mathematical expectation of the random variable *B<sup>λ</sup>* is given by the equation below.

$$E\{\mathbf{g}(\mathcal{B}\_{\lambda})\} \cong \frac{1}{n\_s} \sum\_{s=1}^{n\_s} \mathbf{g}(\mathcal{B}\_{\lambda}(\theta\_l)).\tag{26}$$

After Lagrange multipliers are derived by evaluating the mathematical expectation of the random variable *Bλ*, the explicit PDFs of different elastic elements in the random elasticity matrix can be established. Depending on Equation (12), the PDF of the elastic tensor for the orthotropic symmetric class material could be defined by using the equation below.

$$p\_{\mathcal{C}}(\mathbf{c}) = p\_{\mathcal{C}\_1, \dots, \mathcal{C}\_6}(\mathbf{c}\_1, \dots, \mathbf{c}\_6) p\_{\mathcal{C}\_7}(\mathbf{c}\_7) p\_{\mathcal{C}\_8}(\mathbf{c}\_8) p\_{\mathcal{C}\_9}(\mathbf{c}\_9), \tag{27}$$

with

$$p\_{\mathbb{C}\_1,\dots,\mathbb{C}\_6}(\mathbf{c}\_1,\dots,\mathbf{c}\_6) = k \det(\text{Mat}(\mathbf{c}\_1,\dots,\mathbf{c}\_6)) \exp\left(-\sum\_{i=1}^6 \lambda\_i^{(1)} \mathbf{c}\_i\right),\tag{28}$$

with

$$\text{Mat}(\mathfrak{c}\_1, \dots, \mathfrak{c}\_6) = \left( \begin{array}{cccc} \mathfrak{c}\_1 & \mathfrak{c}\_4 & \mathfrak{c}\_6 \\ \mathfrak{c}\_4 & \mathfrak{c}\_2 & \mathfrak{c}\_5 \\ \mathfrak{c}\_6 & \mathfrak{c}\_5 & \mathfrak{c}\_3 \end{array} \right) \tag{29}$$

and

$$p\_{\mathbb{C}\_j}(c\_j) = k\_j \exp\left(-\lambda\_i^{(1)} c\_i\right) c\_i^{-\lambda^{(2)}}, \; j = 7, 8, 9. \tag{30}$$

**Remarks 1.** *The random variables* (*C*1,..., *C*6)*, C*7*, C*8*, and C*<sup>9</sup> *are mutually independent of each other. C*7*, C*8*, and C*<sup>9</sup> *are Gamma-distributed and the k and kj are the normalization constants.*

#### *4.3. Numerical Application of the Orthotropic Symmetric Material (CLT)*

The material properties of CLT gathered from the literature (c.f. Table 2) are regarded as an initial starting point (mean value) to procced with the stochastic approach presented in the previous sections.

In Abaqus, the orthotropic materials compliance matrix can be defined by the engineering constants:

$$\begin{Bmatrix} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \varepsilon\_{33} \\ \gamma\_{12} \\ \gamma\_{13} \\ \gamma\_{23} \end{Bmatrix} = \begin{bmatrix} \frac{1}{E\_1} & -\frac{\nu\_{21}}{E\_2} & -\frac{\nu\_{31}}{E\_3} & 0 & 0 & 0 \\ -\frac{\nu\_{13}}{E\_1} & \frac{1}{E\_2} & -\frac{\nu\_{32}}{E\_3} & 0 & 0 & 0 \\ -\frac{\nu\_{13}}{E\_1} & -\frac{\nu\_{32}}{E\_2} & \frac{1}{E\_3} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G\_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G\_{13}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G\_{23}} \end{bmatrix} \begin{Bmatrix} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{33} \\ \sigma\_{12} \\ \sigma\_{13} \\ \sigma\_{23} \end{Bmatrix} . \tag{31}$$

The compliance matrix of the CLT panel can be derived depending on the elastic constants reported in Table 2 and Equation (31). The stiffness matrix is the inverse of the compliance matrix.

$$\mathbf{C} = \begin{bmatrix} 10.58 & 2.3 & 2.3 & 0 & 0 & 0\\ 2.3 & 5.2619 & 2.4028 & 0 & 0 & 0\\ 2.3 & 2.4028 & 5.2619 & 0 & 0 & 0\\ 0 & 0 & 0 & 0.9 & 0 & 0\\ 0 & 0 & 0 & 0 & 0.09 & 0\\ 0 & 0 & 0 & 0 & 0 & 0.063 \end{bmatrix} \times 10^9. \tag{32}$$

From Equation (32) and the corresponding orthotropic symmetric matrix basis (Equation (4)), the mean value of *ci* defined in Equation (5) can be deduced by using the formulas below.

$$(\overline{\varepsilon}\_1, \dots, \overline{\varepsilon}\_9) = (10.58, 5.2619, 5.2619, 2.3, 2.4048, 2.3, 0.9, 0.09, 0.063) \times 10^9. \tag{33}$$

Let *f target* = (*c*, *νC*) be the target vector to compute the Lagrange multipliers.

$$f^{\text{target}} = (10.58, 5.2619, 5.2619, 2.3, 2.4048, 2.3, 0.9, 0.09, 0.063, 5.3059). \tag{34}$$

And let *f est*(*λ*) = (*c*(*λ*), *νC*(*λ*)) be the estimated vector to compare with the target vector so that the optimal solution of *λ* can be obtained by solving the following optimization function.

$$J(\lambda^{opt}) = a\gamma \text{s/min} \left( (1 - a)(\mathbb{E} - \mathbb{E}(\lambda))^2 + a(\nu\_{\mathcal{C}} - \nu\_{\mathcal{C}}(\lambda))^2 \right), \tag{35}$$

where *α* ∈ [0, 1] is a free parameter. In this scenario, *α* is set to be 0.5 to give a robust estimation.

Lagrange multipliers associated with *λ*(**1**) can be further expressed by means of *λ*(2) [15].

$$\begin{cases} \begin{array}{c} \lambda\_1^{(1)} = -\lambda^2 \frac{c\_2 c\_3 - c\_5^2}{\triangle}, \lambda\_2^{(1)} = -\lambda^{(2)} \frac{c\_1 c\_3 - c\_5^2}{\triangle}, \lambda\_3^{(1)} = -\lambda^{(2)} \frac{c\_1 c\_2 - c\_4^2}{\triangle}, \\\ \lambda\_4^{(1)} = -\lambda^{(2)} \frac{(\frac{c\_5 c\_6 - c\_3 c\_4}{\triangle})}{\triangle}, \lambda\_5^{(1)} = -\lambda^{(2)} \frac{(\frac{c\_4 c\_6 - c\_1 c\_5}{\triangle})}{\triangle}, \lambda\_6^{(1)} = -\lambda^{(2)} \frac{(\frac{c\_4 c\_5 - c\_2 c\_6}{\triangle})}{\triangle}, \\\ \lambda\_7^{(1)} = -\lambda^{(2)} \frac{1}{c\_7}, \lambda\_8^{(1)} = -\lambda^{(2)} \frac{1}{c\_8}, \lambda\_9^{(1)} = -\lambda^{(2)} \frac{1}{c\_9}, \end{array} \tag{36}$$

with

$$
\triangle = c\_1 c\_2 c\_3 + 2c\_4 c\_5 c\_6 - c\_1 c\_5^2 - c\_2 c\_6^2 - c\_3 c\_4^2. \tag{37}
$$

From the sensitivity analysis, the initial guess of *<sup>λ</sup>*(2) is set to be <sup>−</sup>2. The target optimization function (Equation (35)) is evaluated by using the interior-point method (*fmincon* function) in Matlab [33].

From Figure 11, it could be observed that the optimization algorithm converges fast to a small value and the corresponding optimal values of the Lagrange multipliers are found to be:

$$
\lambda^{(1)} = (0.1263, 0.2904, 0.2905, -0.0758, -0.2324, -0.0758, 12.9089, 18.4413, 1.2909), \\
\lambda^{(2)} = -1.1618.
$$

**Figure 11.** Convergence of the optimization algorithm.

The estimated vector *f est*(*λ*) is evaluated by the Monte Carlo method.

$$\text{conMC}(n\_s) = \left| n\_s^{-1} \sum\_{i=1}^{n\_s} \mathbf{g}(\mathbf{C}(\theta\_i)) \right|. \tag{39}$$

Therefore, the estimated vector *f est*(*λ*) yields the following.

*f est* = (10.5882, 5.2677, 5.2675, 2.3035, 2.4121, 2.2860, 0.8997, 0.0904, 0.0636, 5.3088). (40)

The cost function of the target vector *J λopt* is 1.048211 × <sup>10</sup>−3, which implies good agreement of the estimated values with the reference values.

#### *4.4. Sampling the Defined Probability Distribution Function by Metropolis-Hastings Alforithm*

Following the process introduced in the previous section, the PDFs of the random elasticity tensor of the CLT were constructed. The objective of this section is to generate the random numbers that obey the defined PDFs. The corresponding random stiffness matrix of CLT can be constructed following the generated data. Subsequently, the compliance matrix of CLT can be derived by inversing the stiffness matrix. The random elastic constants (engineer constants in Abaqus) of CLT can be determined according to Equation (31). Lastly, the generated elastic constants should be imported into Abaqus to analyze the dynamic response of the CLT panel. The Markov Chain Monte Carlo method is wildly used to sample the high-dimension PDFs. A specific algorithm, called the Metropolis-Hastings algorithm (MHA), is used in this case to sample the target the function. The proposed distribution is a conventional multivariate Gaussian distribution. The mathematical support of *pC*1,...,*C*<sup>6</sup> (*c*1,..., *c*6) is *det*( *N* ∑ *i*=1 *ciE*) > 0 and the mathematical support of *pCj cj* is *cj* > 0, *j* = 7, 8, 9. When sampling the target function, the generated data should stay in the supports of the sampled function. A total of 50,000 combinations of (*C*1, ... , *C*9) are realized by performing the MHA and by obeying the mathematical constraints (supports) of the target functions. The marginal distributions of the different mechanical constants reconstructed by *ksdensity* function in Matlab are shown in Figure 12.

**Figure 12.** (**a**) Joint probability density function of random variables *C*<sup>11</sup> and *C*22. (**b**) Joint probability density function of random variables *C*<sup>22</sup> and *C*33.

#### **5. Implementation of Stochastic Data in Abaqus**

A total of 50,000 generated random numbers (*C*1, ... , *C*9) have been generated. The corresponding random stiffness tensor could be determined. The random compliance tensor can also be derived by inversing the random stiffness tensor. Since Poisson's ratios have a very slight influence on the eigen-frequencies and the mode shape order of the CLT numerical model, the variation of the Poisson ratios will not be taken into account in the FE model. 50,000 generated random elastic constants are only satisfied with the mathematical constraints. The constraints associated with the physical meaning are not considered when generating the random numbers. In order to ensure the date implemented in the FE model has physical meanings, several physical constraints are set: (1) Variation of Young's moduli and shear moduli should be in a reasonable range of wood. The ranges of Young's modulus and the Shear modulus are assumed in ±50% of the respective mean values.

$$E\_1 \in [4600 MPa, \ 13800 MPa], \ E\_2 \in [2000 MPa, \ 6000 MPa], \\ E\_3 \in [2000 MPa, \ 6000 MPa], \tag{41}$$

$$G\_{12} \in \left[450 MPa, 1350 MPa\right], \ G\_{13} \in \left[45 MPa, 135 MPa\right], \ G\_{23} \in \left[45 MPa, 135 MPa\right]. \tag{42}$$

(2) Young's modulus (Shear modulus) in a principle direction should be larger than the Young's moduli (Shear moduli) in the other two directions.

$$E\_1 > E\_2; \ E\_1 > E\_3;\tag{43}$$

$$G\_{12} \geqslant G\_{13} \colon G\_{12} \geqslant G\_{23}.\tag{44}$$

(3) Young's modulus (Shear modulus) in direction 2 should be larger or equal to Young's modulus (Shear modulus) in direction 3.

$$E\_2 \ge E\_3; \ G\_{13} \ge G\_{23}.\tag{45}$$

Only the generated random elastic constants fulfilled with the above requirements (from Equation (41) to Equation (45)) could be imported in Abaqus. In the work reported in this case, due to limited time and limited computer calculation capacity, only 100 different combinations of elastic constants were selected to calibrate the model and to investigate the influence of material properties on the dynamic response of CLT. This large sum of Abaqus input files with different input mechanical constants were realized with the help of Python scripts.

#### **6. Results and Discussion**

#### *6.1. Quantification of Uncertainties*

One of the objectives of this article is to investigate the effects of uncertainties induced by material properties on the model output and to ultimately obtain a reliable model to predict the vibro-acoustic behavior of different designs. To achieve that, 100 steady-state simulations with different combinations of material properties were carried out in Abaqus.

In Figure 13, each subfigure has 100 FRF simulations and the measurement results are shown in blue. Figure 13 shows that there is an obvious envelope overlapped around the first four peaks lower than 100 Hz. The large variation envelope range is due to varying five mechanical constants in one time since five elastic constants variations have a larger impact on the dynamic response of CLT than only changing one elastic constant in one time.

On the contrary of the frequency range lower than 100 Hz, the simulated FRFs begin to scatter above 100 Hz. No clear envelope peaks can be found around the measured resonances higher than 100 Hz. One possible reason of these scattering curves in a relative higher frequency range is the complexity of the mode shapes. It is known that the lower the eigen-frequency is, the simpler the mode shape is and the longer the wave length is. Long wavelength are not sensitive to the small details in the CLT panel, such as the non-uniform air gaps throughout the laminas and the edge bonding (c.f. Figure 14). It implies that the dynamic behavior of CLT in a low frequency range can be mimicked by a simplified homogeneous orthotropic laminated FE model. However, when the mode shapes become more complex in a higher frequency range, the wavelength becomes smaller. As a consequence, small details of the CLT panel begin to affect the vibration of the CLT panel. In this case, the homogeneous orthotropic laminated FE model could not properly describe the dynamic response of the CLT panel in the higher frequency range. In order to increase the accuracy of the FE in a high frequency range, non-homogeneous laminate layer, such as an account for the irregular air gap in laminas, should be modeled. Nevertheless, it should be aware that the calculation time will become longer, when more details are taken into account in the model. The stochastic method needs a large number of calculations to quantify the uncertainties induced by the material properties so that a compromise should be carefully made between the accuracy of the model and the calculation time.

**Figure 13.** Measured (blue) and simulated (red) FRFs at point 11, 13, 17, and 24.

**Figure 14.** (**a**) Air gaps in the laminate layers of CLT. (**b**) No edge-bonding of CLT.

#### *6.2. Calibration of the CLT Panel*

In this section, the best combination of elastic constants of CLT should be identified by selecting the best NRFDs and MACs among 100 simulations. To select the best combination of mechanical constants, the NRFDs of the first four resonances at point 11, point 13, point 17, and point 24 were calculated. The simulations with the smallest NRFDs of the first four resonances at each point were selected from 100 simulations at each point. The NRFDs of the simulated and measured eigen-frequencies at four excitation positions are shown in Figure 15. It can be seen from the NRFDs of each excitation position that the NRFDs of the first four resonances are lower than 5%. However, the NRFDs of the 5th and 6th resonances are relatively high when compared to the first four resonances. This result emphasizes that more structural details should be involved in the FE model to calibrate the dynamic behavior of CLT in the frequency range higher than 100 Hz. However, only NRFD values are not enough to justify the best combination of elastic constants. Since the NRFDs can only represent the simulated eigen-frequency shifts when compared with the experimental results, the mode order can be different even with a low NRFD. Therefore, MAC numbers are needed to validate the model by ensuring the modes in the same order with reference even if there are low NRFDs. The simulated eigen-frequencies and mode shape are reported in Table 1 and Figure 17. The corresponding material properties of CLT are reported in Table 3. Figure 18 shows that the first six simulated modes are in the same order with the measured ones. However, there exist two extra modes in the simulation results.

**Figure 15.** NRFDs of the first six resonances at points 11, 13, 17, and 24.

(**a1**) Measured 1st mode. (**a2**) Simulated 1st mode.

(**b1**) Measured 2nd mode. (**b2**) Simulated 2nd mode.

(**c1**) Measured 3rd mode. (**c2**) Simulated 3rd mode.

(**d1**) Measured 4th mode. (**d2**) Simulated 4th mode.

(**e1**) Measured 5th mode. (**e2**) Simulated 5th mode.

**Figure 17.** *Cont.*

(**f1**) Measured 6th mode. (**f2**) Simulated 6th mode.

(**g1**) Measured 7th mode. (**g2**) Simulated 7th mode.

**Figure 17.** Measured and simulated modes.

**Figure 18.** Cross-MAC.

**Table 3.** Material Properties of CLT used in the calibrated FE model. Stiffness parameters have the unit of MPa and the density is given in kg/m3.


The same results can also be seen in the mobility of different excitation points of the lowest NRFD values (c.f. Figure 19). The simulated FRFs at these four different excitation points correlate better with the measured ones, while there are extra peaks and eigen-frequency shifts in the simulated FRFs in the frequency range from 110 Hz to 170 Hz. We suspect that these discrepancies are higher than the 110 Hz result from the over-simplified homogenous laminated FE model, which ignores the geometrical details contained in the real CLT panel. Yet, the boundary condition set-up in the model could not describe the real measurement boundary conditions.

**Figure 19.** Magnitude of the complex mobility in the vertical direction of point 11, 13, 17, and 24. Simulated FRFs in red, measured FRFs in blue.

The dynamic properties of wooden structures highly depend on the material properties of the structure, the geometry details, and the workmanship. Consequently, the deterministic model may not be able to represent the dynamic response of the wooden structures in a realistic way. A calibrated model may not be able to accurately predict the dynamic behavior of the theoretical identical wooden structure due to the uncertainties. The stochastic method is applied in this case to quantify the uncertainties induced by material properties so that this model can estimate the dynamic response of a class of wooden structures, instead of only one structure. Moreover, the influence of material properties on the vibration of CLT is the coupling effect of Young's moduli and shear moduli in all directions, so that calibration is always time consuming and tedious work to find the appropriate combination of elastic constants in different directions. The calibration employing the stochastic approach could start from the material properties collected from the literature and set the mathematical and physical constraints to generate the input data to find the best combination of the mechanical constants of the under investigated structure. This method could automate the calibration step to avoid the repetitive manual calibrations. However, we should pay attention to the mathematical and physical constraints before generating the input elastic constant data. Because the generated elastic constants should be in a reasonable range of the under investigated material. Otherwise, the input elastic constants may not

have a physical meaning, even though the calibration results fit well with the reference. The stochastic method uses a large amount of data to describe an unknown problem (the data base of CLT in our case). More calculations are made and more accurate calibration can be achieved. However, a trade-off between the calculation time and the accuracy of the result should be made in order to keep the calculation time in a reasonable range. This method can not only be applied to CLT but also can be employed to calibrate the other wooden structures whose stiffness constants are difficult to obtain.

Furthermore, one of the objectives of this stochastic approach is to calibrate the material properties of the target structure. It would be better to decrease the influence of other influence factors, such as boundary conditions. Therefore, it is suggested to hang up the under investigated structure (free-free boundary condition) or fix the structure boundary to the ground (perfect simply supported condition) to eliminate the influence of boundary conditions as far as possible. In the work reported in this case, due to a lack of support materials, the CLT panel just laid on top of the I-steel beam and it was not screwed into the ground. Consequently, when the CLT is excited, the deformation of the I-steel beam can affect the vibration of the CLT slab. Furthermore, the laboratory boundary conditions are always different from the in-situ boundary conditions [34]. Thus, it would be necessary to investigate the dynamic response of CLT in a real building. To achieve that, the FRFs could be first measured from a CLT bare floor in real mounting conditions. Then, the same CLT bare floor could be set in the simplified laboratory conditions to compare the relative differences between different FRFs under different boundary conditions.

From the FE CLT modelling perspective, the model validation criterions (NRFD and MAC) and the simulated FRFs suggest that dynamic behavior of the CLT panel can be modelled by the homogenous orthotropic laminated FE model in the frequency range lower than 100 Hz. In a higher frequency range, as the inhomogeneity of the laminated layers of CLT slab begins to pronounce in the vibration of CLT panel, more geometrical details in the CLT panel should be taken into account in the FE CLT model to obtain more accurate results.

#### **7. Conclusions**

Low frequency sound insulation is always a challenge for the wooden constructions, especially for the multi-family dwellings. Even though the wooden constructions are satisfied with the standards in force, acoustic comfort is not always met. Since the evaluation frequency range even with the adaptation term of the current standards is from 50 Hz to 3150 Hz, however, the first few resonance frequencies of the wooden floor, which are believed to cause most annoyances, are left out of the evaluation scope. Low frequency prediction tools are needed to access the vibratory performance of wooden buildings at the early design stage due to complaints coming from the inhabitants in wooden buildings. Accessing an accurate low frequency prediction tool requires involving the structure details in the model. Moreover, material properties are another important factor, which can induce a remarkable change in the modelling output.

In this paper, we introduced the stochastic process into the FE model to quantify the uncertainties generated by the material properties. By performing Monte Carlo simulations, variation of Young's moduli and shear moduli in different directions were taken into account in FE model to investigate the coupling effect of different elastic moduli on the dynamic response of structure. Furthermore, 100 simulations were calculated at four different driving points. Clear envelopes can be observed from the simulations lower than 100 Hz. However, the simulations begin to scatter in the frequency range higher than 100 Hz. The best combination of material properties is selected from 100 different combinations of elastic constants to calibrate the FE CLT model. It was noticed that the simulated dynamic response that was lower than 100 Hz was correlated better with the measured dynamic response of CLT. From the promising results, it was concluded that the stochastic method can be applied to a deterministic model (FE model) to quantify the uncertainties of the structures. Furthermore, this method can be employed to calibrate the FE model to acquire the material properties of the under-investigated structure.

**Author Contributions:** Data curation, C.Q. Investigation, C.Q. Methodology, C.Q. and J.N. Supervision, S.M. and D.B. Validation, C.Q. Visualization, C.Q. Writing—original draft, C.Q. Writing—review & editing, C.Q. and J.N.

**Funding:** The authors are grateful to FPInnovation and Nordic Structure and Industrial Chair on Eco-responsible Wood Construction (CIRCERB).

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

#### **References**


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

#### *Article*

### **Nonuniform Bessel-Based Radiation Distributions on A Spherically Curved Boundary for Modeling the Acoustic Field of Focused Ultrasound Transducers**

#### **Mario Ibrahin Gutierrez <sup>1</sup> , Antonio Ramos 2, Josefina Gutierrez 3,\*, Arturo Vera <sup>4</sup> and Lorenzo Leija <sup>4</sup>**


Received: 14 December 2018; Accepted: 19 February 2019; Published: 4 March 2019

**Abstract:** Therapeutic focused ultrasound is a technique that can be used with different intensities depending on the application. For instance, low intensities are required in nonthermal therapies, such as drug delivering, gene therapy, etc.; high intensity ultrasound is used for either thermal therapy or instantaneous tissue destruction, for example, in oncologic therapy with hyperthermia and tumor ablation. When an adequate therapy planning is desired, the acoustic field models of curve radiators should be improved in terms of simplicity and congruence at the prefocal zone. Traditional ideal models using uniform vibration distributions usually do not produce adequate results for clamped unbacked curved radiators. In this paper, it is proposed the use of a Bessel-based nonuniform radiation distribution at the surface of a curved radiator to model the field produced by real focused transducers. This proposal is based on the observed complex vibration of curved transducers modified by Lamb waves, which have a non-negligible effect in the acoustic field. The use of Bessel-based functions to approximate the measured vibration instead of using plain measurements simplifies the rationale and expands the applicability of this modeling approach, for example, when the determination of the effects of ultrasound in tissues is required.

**Keywords:** focused transducer; acoustic field; nonuniform radiation distribution; Bessel radiation distribution; spherically curved uniform radiator; rim radiation; Lamb waves; finite element modeling

#### **1. Introduction**

In recent years, the use of focused ultrasounds (FUS) has been increased in biological applications for both high intensity and low intensity modalities [1–5]. A high-intensity focused ultrasound is used for the rapid destruction of tissues by thermal ablation [2,3,6], for example in oncology, while low-intensity applications are based on producing midterm hyperthermia and nonthermal ultrasonic therapy [3,7], with multiple possible applications [4,8,9]. Among the non-ablating FUS applications reported in literature can be mentioned the low-intensity pulsed ultrasound (LIPUS) using focused

transducers [8], drug delivery in deep tissues (as blood–brain barrier disruption) [10], gene transfer therapy [11], and sonothrombolysis [12,13], among others. In all these applications, the control of the dose in the target volume should be precise during long periods of time to avoid cell death in non-treated zones [12,14]. Noninvasive (and non-expensive) technologies to monitor the temperature in the treated zone [12], more precise and simpler calibration techniques for FUS transducers to determine effective radiating parameters [15,16], in conjunction with accurate and simple computational models capable of effectively representing the acoustic fields of real FUS transducers are required. This could provide more information to adequately study the produced effects along the ultrasound pathway and to quantify undesired consequences in surrounding tissues, previous the application of the therapy. However, in medical applications, the use of simple, but inefficient, ideal models for the acoustic field of FUS transducers is a common and not very questioned practice [17–19].

One of the first approaches to calculate the acoustic field of curved radiators is the O'Neil solution proposed in 1949 [20]. This solution is based on the Rayleigh integral, and it assumes a spherically curved uniform radiator (SCUR) oscillating with a uniform velocity distribution. Although this approximation could be appropriate for many modern transducers at the focus, in the regions where the pressure amplitudes are highly affected by diffraction, e.g., before the focus, the discrepancies are very evident [17]. These variations between "ideal" theory and experiments occur because the assumption of a normal uniform velocity distribution is very conservative and, usually, unreal [17,21,22]. The rather complex vibration of piezoelectric plates not only is composed by a thickness-extensional (TE) vibration mode but also includes contributions of radial modes [23,24], edge waves [17], and Lamb waves [21,25], whose effects are more noticeable under a continuous regime [17]; this vibration can be more complicated if we consider that the piezoelectric plate is not vibrating freely, but it is somehow clamped by its edge producing a higher vibration amplitude at the center of the plate with an attenuated vibration at the edge [26,27]. This complex vibration occurs in both planar and concave plates, and it should be accounted for when producing models of acoustic fields for more accurate results.

In this paper, we are proposing an approach to model the acoustic field produced by FUS transducers in a low-intensity regime (considering linear propagation) using polynomial-Bessel based functions as nonuniform radiating distributions on a curved surface. The reason of proposing these functions is based on the reported vibration patterns produced in piezoelectric curved disks composed of a main vibrating thickness-extensional (TE) mode generating a wave in the thickness direction and a second component of Lamb waves propagating radially [17,28,29]. These two components are modified when the disk is fixed to the transducer case, which produced a combined vibration pattern that can be approximated by a polynomial-Bessel based function, in accordance with the measurements for any specific transducer. With this approach, we got better results than the widely used SCUR, which can be comparable with the results using the intuitive approach for modeling the acoustic field using the velocity distribution measured on the radiating surface [29,30]. Using analytical functions instead of measurements will permit to propose mathematical algorithms to produce realistic acoustic fields of actual FUS transducers; however, these demonstrations are beyond the scope of this paper. The future applications of these new models are open, since the assumptions taken for this work are not particular.

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

In this paper, the acoustic field of a FUS transducer has been modeled. A 2 MHz spherically focused transducer with a 20 mm nominal focal length and a 20 mm of nominal aperture (Onda Corporation, Sunnyvale, CA, USA) was used for the experiments; these values usually differ from the measured ones as it will be discussed later [15]. This kind of transducer, designed for high-intensity applications but used here for low-intensity measurements, has its negative terminal exposed, i.e., it does not have a nonconductive front-layer. This aspect will be useful for the explained below curvature measurements.

#### *2.1. Measurement Setup*

Three different sets of measurements were carried out for this work using the setup shown in Figure 1. For the first set, the curvature of the FUS transducer was determined by taking advantage of its exposed negative terminal on the radiating surface; for this, a multimeter (Fluke 289, Fluke Corporation, Everett, WA, USA) and a 3D positioning system (Onda Corporation, Sunnyvale, CA, USA) driven in manual mode were used. The transducer was placed at the final position for acoustic field measurements (to be made after) in which the transducer will emit in positive *z*-direction; this first measurement was made in the air. Then, the multimeter was set to measure electric continuity, and one of its probes was attached to the positioning system while the other was connected to the negative terminal of the transducer. The measurement probe was moved step-by-step (step resolution of 6.36 μm) towards the transducer surface, down in the *z*-direction (relative to the direction of the ultrasound emission), for a fixed *x*-coordinate until the probe touched the transducer conductive surface, indicated by a "beep" from the multimeter. This procedure was repeated diametrically for every *x*-coordinate from −14.573 mm to 14.573 mm with a step resolution of 0.3312 mm, covering the transducer case and the curved piezoelectric plate. The used step-resolution was smaller than a half of a wavelength of 2 MHz ultrasound in water at 20 ◦C (approx. 0.38 mm). The determined *z*-coordinates for all the "*x*" were saved for further use. These final coordinates were used for the second below explained measurement.

**Figure 1.** The setup for the curvature/acoustic pressure measurements.

For the second set, the transducer was immersed in a tank filled with degassed distilled water carefully keeping intact the transducer setup used for the first set (see Figure 1). This second set of measurements was carried out using the previous obtained coordinates for the curvature to determine the radiation pattern "very close to the transducer surface". The radiation measured at this region would closely represent the vibration distribution of the curved piezoelectric disk, assuming the effect of other regions of the radiating surface is negligible. For this, a wideband needle hydrophone PZTZ44-0400 (Onda Corporation, Sunnyvale, CA, USA) with a 40 μm aperture and sensitivity of −260 dB referred to 1 V/μPa was mounted on the positioning system, replacing the multimeter probe. The measurements were made at the positions determined previously but 1 mm in front of the transducer surface (in *z*-direction). The transducer was driven with a wave generator (Array 3400, Array Electronic Co., Taiwan, China) using s 20 Vpp sine tone-burst; the received data were recorded in a PC. The ultrasound signals were recorded for 1, 5, 10, and 20 sine cycles to determine the radiation pattern under different excitation conditions but, more specifically, to know in which number of cycles the vibration presents a quasi-stationary pattern. The data were post-processed in MATLAB (R2017a, MathWorks, Natick, MA, USA) to determine the peak-to-peak values at each spatial point and the full radiation distribution for each excitation condition.

The third set of measurements was performed to determine the full acoustic fields of the focused transducer. This was measured using the previously described setup. The acoustic fields were obtained by saving only the peak-to-peak acoustic pressure at each point in the radiated volume. It was recorded on an XZ plane covering the transducer dimensions starting at *z* = 2 mm and finishing at *z* = 50 mm from the transducer case, with a step resolution of 0.3312 mm in the *x*-direction and 0.5000 mm in the *z*-direction. The XY planes were measured at two depths, 0.2 cm and 1.7 cm (at the focus), using a step resolution of 0.3312 mm in both directions; these planes were obtained to verify the symmetry of the acoustic radiation. As formerly, the acoustic fields were captured for 1, 5, 10, and 20 sine cycles to determine the conditions for a nearly continuous emission pattern. The data were postprocessed in MATLAB to reconstruct the fields.

#### *2.2. Acoustic Field Modeling Using FEM*

The acoustic field was modeled using the finite element method (FEM) based on the geometry shown in Figure 2. The software used for the FEM processing was COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden) working in a PC of 8-core 3 GHz microprocessor and 64 RAM (Dell, Round Rock, TX, USA). Based on the cylindrical symmetry of the transducer, the problem was assumed axisymmetric. The validity of this assumption was verified with acoustic field measurements in which the profile of the radiation along the azimuthal coordinate is similar for any angle (data not shown). The mesh in the rectangular part of Figure 2 consisted of 10 square elements per wavelength, i.e., more than 530,000 elements; triangular elements were used only for the zone created with the boundaries 1 and 5 to simplify meshing. The requirements to mesh using squares are very specific, and this kind of element cannot be used in certain curved geometries; triangular elements are sometimes the only option for these complicated zones. The mesh convergence was verified by increasing the mesh resolution, which indicated an error of 0.01% at the focus amplitude between meshing using 9 and 10 elements per wavelength. Using square elements instead of triangular permitted to increase the spatial resolution with the same interpolation functions because this kind of mesh has a larger number of nodes (and therefore more degrees of freedom) than the triangular one with the same number of elements (more nodes per element area); this can be noticed with the reduction of the solution time compared to the time required for solving the problem using the same number of triangular elements. The largest solution time registered for our main model, which was barely the most demanding of computational resources, was 45 s.

**Figure 2.** The finite element method (FEM) axisymmetric geometry for modeling the focused acoustic field.

It was assumed the linear ultrasound propagation (i.e., rather low acoustic intensity) was in an attenuation-free homogeneous media (degassed distilled water). The ultrasound propagation is determined with FEM based on the homogeneous Helmholtz wave equation for harmonic radiation, assuming a purely harmonic source and a no-frequency dispersion. This equation can be written as

$$
\nabla^2 p + k^2 p = 0 \tag{1}
$$

where the wavenumber is *k* = *ω*/*c*, *ω* is the angular frequency (rad), and *c* is the speed of sound in water (m/s) assumed constant. Boundary 1 was set with the harmonic normal acceleration. The amplitude of this normal acceleration was adjusted for the different radiation conditions presented in this paper. For instance, in the uniform radiation model of the next section, the amplitude of the harmonic acceleration was set constant along the radius. Usually, the radiation of a boundary is expressed in terms of particle velocity [20,26,31]. In harmonic conditions [32], the particle acceleration *a*<sup>0</sup> on the radiator surface can be obtained by time-deriving the constant-amplitude particle velocity *v*0*ejω<sup>t</sup>* as

$$a\_0 = \frac{d}{dt} \left( \upsilon \omega^{j\omega t} \right) = \upsilon \upsilon \omega e^{j(\omega t + \frac{\omega}{2})} = \frac{\omega}{\rho c} p \upsilon e^{j(\omega t + \frac{\omega}{2})},\tag{2}$$

where *ρ* is the medium density (kg/cm3). The last term was determined by considering *p*<sup>0</sup> = *ρcv*0, which is true in the field very close to the radiator surface [27]. The term *<sup>π</sup>* <sup>2</sup> in Equation (2) is a time phase shift between the velocity and the acceleration, and this is not related to the spatial profile along the transducer surface [32]. Then, under harmonic simulations, using either acceleration or velocity for the radiation distributions of the plate is not relevant if adequate amplitude considerations are taken.

Boundaries 2–4 were configured to match the acoustic impedance of water and to reduce the ultrasound reflections at the walls [27]. Then, *Z* = *ρc* = 1.5 MRayls at 25 ◦C, where *Z* is the acoustic impedance (MRayls) given by the product of the media density and the speed of sound; the walls were considered to be perfectly flat. However, the dimensions of the FEM geometry (10 cm × 3 cm) were big enough to not have residual reflected waves affecting the region of interest for this application, i.e., a region of 5.0 cm depth and 1.5 cm width after the transducer. Boundary 5 was set with a continuity condition, and it was used only to simplify the meshing. Boundary 6 represents the transducer rim, which was set, accordingly to the measurements, to uniformly radiate a relative pressure of 8% of the average pressure radiated by the curved surface in the effective radiating area; this area was determined as the area producing 95% of the transducer's radiation, as defined in other applications for planar radiators [27]. The use of this effective area instead of the nominal area improved the model results as it will be explained later; this was already proposed by other means in Reference [15]. The rim radiation was included because it represents an important contribution in the field for this transducer, more evident in the post-focus field but with some little effects in the pressure amplitude on the propagation axis (along *z*-axis) in the prefocus zone.

#### 2.2.1. Radiator with Uniform Vibration Distribution

Conventionally, the acoustic field produced by planar ultrasound transducers has been represented as the product of the radiation coming from a uniform vibrating surface (the piston in a Baffle and the Rayleigh equation) [33]. This assumption produces adequate results for wideband transducers, in which the vibration is usually damped by the backing material and, thus, produces quasi-uniform displacements along the transducer radiating surface [34–36]. However, narrowband transducers often do not have a backing material (air-backed) which makes the vibration of their piezoelectric components less uniform [27]. These nonuniformities have an important effect in the field.

The acoustic field of focused transducers has been historically modeled following the same supposition of planar plates, in which it is assumed the plate is a slightly curved uniform radiator (SCUR) [20]. For this, the curvature of the piezoelectric plate of the transducer is usually assumed spherical; conversely, this curvature is not often reported in datasheets. The most well-known theoretical proposal to determine the acoustic field produced by curved surfaces was made by O'Neil in 1949 [20]. This model works adequately for low-power wideband transducers with a backing material, and it has produced poor results when used to model the acoustic field of air-backed narrowband radiators [17,36]. In order to compare our proposal with the most used ideal approach of curved transducers, the acoustic field produced by a SCUR was determined. For this, a constant amplitude harmonic normal acceleration was used in the curved boundary 1, in Figure 2.

#### 2.2.2. Radiator with "Classical" Nonuniform Vibration

Nonuniform vibration distributions have been proposed to generate consistent acoustic field models of some planar radiators but with a limited range of applicability. These are based on the assumption that, under certain conditions, a radiator can behave not only as a piston but also as a membrane and a clamped plate [26,31,37]. The general equation of "classical" nonuniform radiation distributions include the two more common theoretical conditions for fixing the edge of a rigid piezoelectric plate [38]: 1) a plate with simply supported edges that restrict edge movement in any direction but allow rotation by the edge (simply supported radiator) and 2) a plate with clamped edges that restricts the movement and rotation in any direction (clamped radiator). The general equation for the nonuniform radial acceleration *a*0(*r*) of the surface of a plate radiator can be expressed by [26,31,37]

$$a\_0(r) = \sum\_{mn} A\_{mn} \left[ 1 - \left(\frac{r}{R}\right)^{2m} \right]^{n+1},\tag{3}$$

where *R* is the fixed radius of the plate (m), *Amn* is the adjustable amplitude of the acceleration, and constants *m* and *n* will determine the shape of the distribution. However, only two cases have been reported to represent an actual meaning in Equation (3), specifically, when *n* = 0, the radiation distribution corresponds to the simply supported radiator (SSR) for *m* = 1, 2, 3, ...; when *m* = 1, the radiation distribution is that of a clamped radiator (CR) for *n* > 0. This radial acceleration is set in boundary 1 shown in Figure 2. The summation in Equation (3) can be reduced after assuming that the first vibration mode dominates [31].

#### 2.2.3. Radiator with Bessel Acceleration

The vibration distribution of a piezoelectric plate is closer to a combination of Bessels than a uniform distribution [38–40]. This behavior depends on the material's piezoelectric properties, the constraint conditions, the coupling layers, and the shape. When the plate is concave, the vibration of the concave-shape piezoelectric plates can also be composed of Bessel-like vibration distributions due to the Lamb waves [17]. The components affecting the vibration of an ultrasound transducer are vast and difficult to quantify. Then, for the proposal of this paper, the vibration was determined by measuring the acoustic pressure very close to the radiating surface [27]. If we suppose the ultrasound behaves as a plane wave in the region very close to the radiating surface (when *r* → 0), the acoustic pressure is related with the acceleration as shown in Equation (2).

Thus, the proposal of this paper is to use, as the nonuniform radiation distribution, a Bessel-based function composed by two sub-functions [27]. For this, the acceleration on the curved radiating surface will be expressed by

$$a\_0(r) = A\_0 f(r) \cdot g(r),\tag{4}$$

where *r* is the variable radius related to each point in the curved boundary and *A*<sup>0</sup> is the amplitude of the emitted acoustic pressure, provided that the amplitude of the product of *f*(*r*) and *g*(*r*) is 1 at *r* = 0. The function *f*(*r*) is given by the acceleration of Equation (3), and it represents the effect of attaching the piezoelectric plate to the transducer case (edge clamping), then *<sup>f</sup>*(*r*) = <sup>∑</sup>*mn Amn*- 1 − (*r*/*R*) <sup>2</sup>*mn*+<sup>1</sup> . This function can be determined by adjusting the parameters to get a correct representation of the measured profile, which can be evaluated by correlation.

The function *g*(*r*) represents the radial measured peak distribution (MPD) shown in Figure 3. This was approximated using the Bessel-based function

$$g(r) = \mathbb{C}\_1 l\_0 \left(\beta\_{2N} \frac{r}{R}\right) + \mathbb{C}\_2. \tag{5}$$

Equation (5) was proposed after considering the reports and simulations about Bessel patterns in the vibration distribution on some transducer surfaces [17,28,29]. For the profile measured for the transducer used for this paper, we used a pure negative Bessel *J*<sup>0</sup> mounted over either an SSR or CR function to approximate the experimental MPD. Then, the Bessel-SSR and Bessel-CR distributions were obtained by combining Equations (5) and (3) into Equation (4), with adequate constants. Here, *C*<sup>1</sup> and *C*<sup>2</sup> are the coefficients that control the amplitude and the offset of Equation (5), respectively, and *β*2*<sup>N</sup>* is the 2*N* zero of *J*0, where *N* is the number of peaks in the MPD. When *C*<sup>1</sup> is close to zero, the acoustic field is close to the ideal uniform radiation (SCUR) for *C*<sup>2</sup> = 0; for practical concerns, the average amplitude into the effective radiating radius of the function *a*0(*r*), after selecting adequate functions *f*(*r*) and *g*(*r*), should be 1, which can be adjusting with an adequate value of *Amn*. The effective radius was determined by calculating the effective area containing 95% of the total radiation. Then, the amplitude of the emitted pressure can be controlled with *A*0. In Figure 3, *C*<sup>1</sup> = −0.3, *C*<sup>2</sup> = 0.9, *R* = 1 cm, *m* = 1, *n* = 0.3, *N* = 5, *Amn* = 11.7, and *β*2*<sup>N</sup>* = 30.67, which is the tenth zero of Bessel *J*0. The values of *m* and *n* of Equation (3) were determined by the values producing the maximum correlation with the measured data. The Bessel-SSR of Figure 3 requires *n* = 0; Bessel-MOD (MOD, modified Bessel-SSR) have the same constants as the Bessel-SSR, but *n* = 0.3. These values are highly dependent on the transducer construction, not only on the piezoelectric constants and disk geometry; they are dependent on the way the disk was attached to the case, the thickness and properties of extra layers (electrodes, glue, etc.), and the mechanical properties of the case. Because of these reasons, it would not be possible to specify any rule to determine the parameters of Equations (3) and (5), neither by simple transducer inspection nor even after knowing the transducers' electrical characteristics. Actual measurements of the emitted field very close to the transducer should be made.

**Figure 3.** The measured acoustic pressure at 1 mm from the radiating concave surface along its radius compared with the uniform and proposed Bessel distributions. The graphs are normalized with the average radiated pressure at the radiator surface.

#### **3. Results**

In this paper, we propose a Bessel-based nonuniform vibration distribution for the surface of a spherically curved radiator to get a modeled acoustic field closer to a real measured field. Three classical approaches based on uniform (SCUR) and nonuniform distributions are included to contrast our main results. Figure 3 shows the measured pressure profile at 1 mm from the piezoelectric curved plate that is composed of local peaks distributed along the radius. This distribution can be represented with Equation (4), combining an adequate Bessel function (Equation (5)) and a function representing the radiator clamping using the traditional nonuniform approaches (SSR or CR with Equation (3)). Using a radiating function instead of the measured data for determining the acoustic field could permit the generalization of the radiation coming from this kind of transducer and eventually simplify the calculations by proposing analytic solutions for certain conditions. Although SSR and CR nonuniform distributions have been suggested since the 1960s decade for planar radiators, these have not been formally proposed for curved radiators [26,37]; their inclusion in this paper is not just for comparison purposes but as two valid alternatives for producing certain types of acoustic fields. Then, the use of SSR/CR nonuniform distribution for curved radiators is not discarded by this work, so this can be an effective approach if the operation conditions of the transducer produce these kinds of vibration profiles.

In Figure 4, it is shown the relative pressure profile along the propagation axis (*z*-direction) of classical SCUR and two "traditional" nonuniform vibration distributions compared with the measured field. The improvement in the fields using nonuniform distributions on the radiator can be noticed. When the vibration of a curved radiator is uniform, the field has very sharp local peaks along the propagation axis that are smoothed when the emitted radiation is gradually reduced at the plate rim, as in the proposed nonuniform approaches for *r* → *R*. The concordances of these fields and the measurements are much better, except in the prefocus zone in which the diffraction patterns differ. In both nonuniform approaches, the post-focus zone is smoother, following the same tendency as the measured profile. There is an after-lobe at about 3 cm, which is well-represented by both nonuniform approaches but only when those included the rim radiation. Without this component, this lobe disappears [41], keeping the most other components of the field. This rim radiation was also included in the main results of this paper.

**Figure 4.** The amplitude of the relative acoustic pressure along the propagation axis of the spherically curved radiator with uniform distribution (SCUR), simply supported radiator (SSR), and clamped radiator (CR) conditions vs. measured field. The pressure is relative to the average pressure on the radiator surface of each condition.

The results of using the Bessel-based radiation distributions are shown in Figure 5. For these proposals, Equations (3) to (5) were used to set the acceleration of the curved radiating boundary; also, the radiation from a 3.5 mm rim was included (boundary 6). For the Bessel-SSR approach, it can be noticed that the acoustic field at the prefocal zone has a peak distribution very close to the measured field, and the post-focal zone has the after-lobe produced by the rim radiation with an amplitude equivalent to that of the measurements. The Bessel-CR approach did not provide a good result in the prefocal zone, preserving only one peak of the three located at the measured field. The after-lobe in the far field was not present for this approximation, even when it included the rim radiation, which could indicate this approach is useless for this particular transducer. However, controlling the value of *n* in Equation (3) can reduce the amplitude of the focus of the pure SSR condition to make the model fit with the measurements if correlation does not provide an acceptable result. For this paper, correlation was useful, and the value of *m* and *n* were easily found by this method, which permitted to find the constants for the modified SSR profile (Bessel-MOD) in Figure 5, with a more accurate result comparable with the measurements.

**Figure 5.** The amplitude of the relative pressure on the propagation axis using the Bessel-based nonuniform approximations with SSR and CR conditions on the curved radiating surface. Bessel-MOD is the Bessel-SSR with *n* = 0.3, which emulates correctly the measured field distribution. The pressure is relative to the average emitted pressure. The three models include an 8% of rim radiation, equivalent to the radiation measured at that zone.

The graphs of Figure 6 show the radial profiles at different depths. Important regions were chosen to plot these graphs: Figure 6a is the graph obtained at the first depth measured for the XZ plane; Figure 6b,c is at two equally space depths before the focus; Figure 6d is at the focus; Figure 6e is at the first minimum after the focus; and Figure 6f is at the after-lobe (after the focus). In these figures, the radial profiles at different *z*-distances for the Bessel-MOD are closer to the measurements than the uniform approach, not only in amplitude but also in shape. Other approaches presented in this paper were omitted from this figure for clarity, even if they also would have shown an improvement compared with the uniform approach. For instance, the pure SSR produced in the prefocal zone a profile like the Bessel-MOD without the peaks (a smoother profile), while after the focus, the SSR profiles were the same as those presented for the Bessel-MOD in Figure 6e,f.

**Figure 6.** The radial relative pressure at different depths of the measured acoustic field compared with the SCUR distribution and the Bessel-MOD distribution, i.e., Bessel-SSR with *n* = 0.3 in Equation (3). Other analyzed distributions did not fully match the measured profile, and they were not included in these graphs for clarity, namely, SSR, CR, and Bessel-CR. The depths are (**a**) 0.2 cm, (**b**) 0.5 cm, (**c**) 1.0 cm, (**d**) 1.7 cm (focus); (**e**) 2.7 cm (down peak); and (**f**) 3.1 cm (lobe after focus). All the graphs have the same type of lines as in (**a**) and the same axis labels. The pressure is relative to the average pressure of Figure 3 of each distribution.

In Figure 7, the full acoustic fields of the more representative models of this work are shown. The Bessel-based approaches with modified SSR (SSR with *n* = 0.3) and CR components are presented in Figure 7c,d, respectively, in order to be compared with the measured field of Figure 7a and the most used classical SCUR distribution of Figure 7b. This ideal SCUR distribution is composed of a very characteristic diffraction pattern before the focus, which is not present in the measured profile. The Bessel-MOD proposal provides a better result in the prefocal zone, with a similar diffraction pattern to the measured field. The Bessel-CR still has a similar diffraction pattern but with a narrower acoustic field and larger focus size, probably because of the reduced effective radiating area [41]. The focus locations in the Bessel-based models were more in agreement with the measured field than the SCUR; this was probably because of the reduction of the effective radiating areas, which were determined from each data set (models and measurements) to normalize the average emitted radiation that was the reference of the relative pressure used in the field comparisons.

**Figure 7.** The normalized acoustic field of the focused transducer at a central longitudinal plane: the transducer would be at the bottom of each figure. The linear color scale: white = 1 and black = 0. (**a**) The measured acoustic field; (**b**) the modeled field using SCUR; (**c**) the modeled field using a modified Bessel-SSR (*m* = 1, *n* = 0.3); and (**d**) the modeled field using Bessel-CR (*m* = 1, *n* = 1).

#### **4. Discussion**

A proposal to effectively model the acoustic field of focused transducers has been presented. This was based on another previously proposed nonuniform distribution for planar radiators named Bessel-based nonuniform radiation distributions (shown in Figure 3). In this work, two other more general proposals are also used as comparison: the simply supported radiator (SSR) and clamped radiator (CR). The three proposals have shown that can be used to generate acoustic fields in good agreement with the field produced by focused radiators. The distributions were set in a curved boundary, and the results were obtained with a FEM commercial software. Only low-power linear behavior was considered in this work, but the main rationale presented here could be applied for high-power nonlinear measurements. The emission was considered to be purely harmonic operating at the transducer's nominal frequency of 2 MHz; the measured transducer emission bandwidth was 100 kHz, which represents about 5% of the central operation frequency. Because of that narrow bandwidth, the effect of this parameter in the acoustic field was considered negligible. The use of an analytical expression in the radiating boundary instead of the raw measurements, for acoustic field processing, would open the possibility not only to reduce computation but also to find a simpler way to calculate analytically the acoustic field in a closed form under certain conditions, for instance using the Hankel transform [42] or the direct integration of the Rayleigh integral [22,24]. However, this process is beyond the scope of this paper.

Figure 4 shows the axial acoustic fields for the SCUR approach and the two classical nonuniform distributions SSR and CR already proposed for planar radiators. It can be noticed that the latter produce good results in the overall field amplitude when they are set in this curved radiating boundary but with poor agreement in the prefocal zone. At this region, the diffraction effects are quite different between the measured profiles in both the SSR and CR models. Actually, the diffraction profile of pure SSR and CR fields are similar to the SCUR field but with a reduced amplitude and an "offset". This indicates that reducing the emitted field amplitude at the edges of the radiator provokes also a smooth overall acoustic field. The SSR and CR fields are closer in amplitude to the measured field, and this could occur because the piezoelectric plate of the transducer is fixed at the edges, which produces an emitted radiation closer to a SSR/CR distribution than a SCUR, as seen in the measured field in Figure 3. This is also noticeable in the focus locations, which were closer to the measurements in the

SSR and CR distributions than the SCUR; this was improved probably because of the reduction of the effective radiating area in the nonuniform distributions [15]

Although these classical nonuniform distributions are closer in amplitude to the measured field, their results are not fully satisfactory. Then, it is proposed the use of Bessel-based distributions is closer to the real vibration distribution produced in concave transducers. The results of combining a Bessel with SSR and CR edge conditions are shown in Figure 5. From those results, it can be noticed that the prefocal zone of the acoustic field is better modeled after including the Bessel behavior in the radiation distribution. The peak-distribution in the modeled prefocal zone using Bessel-SSR has the same number of peaks as the measured field and almost at the same position. The focus amplitude was improved when adapting the equation with *n* = 0.3, for the Bessel-MOD, which produced a field more congruent with the measurements. The focus locations of these representations were also more in agreement with the measurements than the classical models presented before, probably due to the reduction of the effective radiating areas. In the post-focus zone, the model can also produce the after-lobe observed in the measurements but with a slightly lower amplitude. This after-lobe was only present in the Bessel-SSR and Bessel-MOD simulations with 8% of rim radiation, which was based on the measurements (see rim radiation in Figure 3). The Bessel-CR condition produced a similar profile at the prefocal zone but with less pronounced peaks. In this case, the after-lobe in the post-focal zone was practically unnoticeable, even with the inclusion of the rim radiation, probably because of the total dominance of the acoustic intensities coming from the main radiating surface.

Figure 6 shows the radial distributions at different depths of the three more significant cases for our comparison purposes, namely, the SCUR, the modified Bessel-SSR (Bessel-MOD), and the measurements. Before the focus, in Figure 6a–c, the radial profile using the Bessel-MOD follows adequately the measured field with some variations in amplitude in the central peak. In spite of the complexity of adequately simulating this region, the Bessel-MOD approach correctly produces acceptable results. A similar agreement was obtained in other less complicated regions as the focus and after the focus in both amplitude and field width. In the former, shown in Figure 6d, the Bessel-MOD produced almost the same amplitude as the measurements, with a 0.4% relative error; the focus with the SCUR is 40% larger than the measured one. After the focus (Figure 6e,f), our model still has a similar tendency as the measured field, with a good match of the graphs in Figure 6e and little variation in the central peak in Figure 6f. From these graphs, the SCUR does not represent an adequate model for this kind of radiators, i.e., those with nonuniform Bessel-based vibration patterns.

The acoustic fields presented in Figure 7 show the differences among each approach in a detailed manner using a linear color scale. The prefocal zones of both Bessel approaches are similar but with a larger focus for the Bessel-CR approach. Focus locations were also improved using our proposed models. The SCUR approach differs at practically any region from the measured field. For this transducer, the Bessel-MOD presents better results than the other nonuniform distributions studied here, with a clear good agreement in most of the graphs shown in this paper. However, this does not mean that the applicability of this model is universal for curved radiators. Some transducers could still behave as ideal pistons, simply supported curved disks, or clamped curved disks that could require the use of any of the other approaches (SCUR, SSR, CR, and Bessel-CR). The use of any of these models should depend on the real radiation profile measured very close to the radiating surface. This paper has presented a new modeling approach to effectively simulate the acoustic field of focused transducers that can be adapted to most of the devices used in different medical applications. Eventually, using the proposed functions as radiation distributions of this kind of transducers could potentially permit to find more analytical solutions of this kind of radiators that could better match the measurements.

#### **5. Conclusions**

In this paper, we have presented four proposals of nonuniform vibration distributions on the radiating surface of a curved transducer to obtain more realistic simulated acoustic fields very congruent with field measurements. From the results here shown, it was possible to conclude that these models provide better approximations of the vibration distribution capable of producing accurate representations of acoustic field for focused applications. When using Bessel-based functions for the vibration of the curved radiator, the prefocal zone of the transducer was correctly simulated. For the post-focal zone, our Bessel proposal had to include the radiation coming from the transducer rim, which allowed the incorporation of the central "after-lobe" observed in the measured field. In the other proposed approaches, the amplitude of the focus significantly varied with respect to the measurements. This happens possibly because of the differences in the proportion of the emitted average radiation used for determining the relative pressure and because of the effective radiating area in each condition, which is a parameter rarely considered for this kind of transducers [15].

Having a model to correctly simulate the acoustic field in the prefocal zone for focused radiators is a very important improvement in the field. This will permit to increase the accuracy in therapy planning, to improve the prediction of thermal increments in tissues outside the focus, to produce better thermal models in hyperthermia, and in general, to have a better dose control. The use of more realistic but still simple models of acoustic field for focused radiators will help to control the undesired effects out of the treatment zone, i.e., before and after the focus, and to easily incorporate this model proposal to therapy planning. In this work, it was proven that our models represent better alternatives for focused radiators than the widely used ideal uniform approaches.

**Author Contributions:** Conceptualization, M.I.G. and A.R.; data curation, M.I.G.; formal analysis, M.I.G.; funding acquisition, M.I.G., A.R., A.V., and L.L.; investigation, M.I.G.; methodology, M.I.G. and A.R.; project administration, M.I.G., A.V., and L.L.; resources, M.I.G., J.G., A.V., and L.L.; software, M.I.G.; supervision, M.I.G., A.R., J.G., A.V., and L.L.; validation, M.I.G.; visualization, M.I.G.; writing—original draft, M.I.G.; writing—review and editing, M.I.G. and A.R.

**Funding:** This research was funded by CONACyT, grant number 257966; CSIC, grant number COOPB20166; ERAnet-EMHE CSIC, grant number 200022; Spanish P.N RETOS, grant number DPI2017-90147-R; and CYTED-Ditecrod Network Ref. 218RT0545.

**Acknowledgments:** The authors would like to thank Rubén Pérez Valladares for his technical support during the acoustic field measurements.

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

#### **References**


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

*Article*

### **Comparison of Multi-Physical Coupling Analysis of a Balanced Armature Receiver between the Lumped Parameter Method and the Finite Element/Boundary Element Method**

### **Yuan-Wu Jiang 1, Dan-Ping Xu 2, Zhi-Xiong Jiang 1, Jun-Hyung Kim <sup>1</sup> and Sang-Moon Hwang 1,\***


Received: 24 January 2019; Accepted: 25 February 2019; Published: 27 February 2019

**Abstract:** The balanced armature receiver (BAR) is a product based on multiphysics that enables coupling between the electromagnetic, mechanical, and acoustic domains. The three domains were modeled using the lumped parameter method (LPM) that takes advantage of an equivalent circuit. In addition, the combined finite element method (FEM) and boundary element method (BEM) was also applied to analyze the BAR. Both simulation results were verified against experimental results. The proposed LPM can predict the sound pressure level (SPL) by making use of the BAR parts dimension and material property. In addition, the previous analysis method, FEM/BEM, took 36 h, while the proposed LPM takes 1 h. So the proposed LPM can be used to check the BAR parts' dimension and material property influence on the SPL and develop the BAR product efficiently.

**Keywords:** balanced armature receiver; lumped parameter method; finite element method and Boundary element method

#### **1. Introduction**

In terms of acoustic transducers, there are the MEMS receiver [1,2], the dynamic receiver [3], the speaker box with passive radiator [4], and the balanced armature receiver [5–20]. The Balance armature structure was devized by Olsen [5]. Three types of structure are described: the unpolarized armature, the polarized reed, and the polarized balanced armature. A new magnetic circuit balanced armature structure transducer has been developed for use in hearing aids [6,7]. In the balanced armature structure transducer, the coil is moved outside the magnetic structure. Stationary gaps are located at the two side legs. Based on the proposed polarized balanced armature microactuator, an implantable hearing device was developed [8]. The balanced electromagnetic separation transducer used in a bone-anchored hearing aid is presented in ref. [9]. To minimize harmonic distortions, quadratic distortion forces and static forces are counterbalanced. A closed loop armature is proposed in the structure design of BAR [10,11].

Nowadays, the balanced armature receiver (BAR) is widely used in hearing aids and earphones because of its small size and high sensitivity. The BAR is a product that follows the principles of multiphysics. When current is passed through a coil, the flux density in the upper and lower air gap is different, which contributes to the generated force on the armature's end. With the input force, the armature and pin vibrate. The vibration of the pin contributes to the diaphragm's vibration, which leads to sound radiation through the spout. The BAR consists of electromagnetic, mechanical, and acoustic physical domains, which are coupled with each other.

To construct the electromagnetic mathematic modeling of the BAR, the lumped parameter method (LPM) was proposed in previous research [12,13] which does not investigate the acoustic domain and cannot predict the SPL result, the key performance criteria of BAR. The multimode of BAR is investigated by using the lumped parameter in refs. [14,15]. The electromagnetic domain and acoustic domain is not involved. To obtain the SPL result, the acoustic domain is considered [16–18] while the input electromagnetic parameter (resistance, inductance, force factor) and mechanical parameters (mass, stiffness, force factor) come from experiment, which means the input parameter cannot be obtained in the analysis if there is no sample or experiment.

This paper is organized in the following structure. In the first part, to predict the SPL and develop a new structure BAR product, LPM is proposed to analyze the performance of the BAR by modeling the electromagnetic, mechanical, and acoustic domains along with their coupling effect according to the modeling dimension and material property, which is defined as the proposed method. Second, BAR is analyzed by FEM (electromagnetic, mechanical domain) and BEM (acoustic domain) which is defined as the previous method [19,20]. Consequently, the samples were manufactured according to the analyzed modeling. Based on the samples, the SPL experiment result was obtained. Finally, the result shows that LPM was verified by experiment, as the previous method (FEM and BEM). The difference is that the proposed method (LPM) takes 1 hour, while the previous method (FEM and BEM) took 36 hours. In conclusion, the contribution of the paper is that the proposed LPM method makes use of the BAR modeling dimension and material property to predict the SPL result without parameter identification by experiment and is more efficient than the previous method. The proposed LPM method can be used to design a new structure BAR and shorten the development cycle of the BAR. The method detail is listed in Table 1.

**Table 1.** Comparison of simulation methods in the paper.


#### **2. Analysis by LPM**

#### *2.1. Electromagnetic Analysis*

The electromagnetic modeling can be described by an equivalent circuit, which is demonstrated in Figure 1. By solving Kirchhoff's voltage and current laws, the flux density in the armature can be expressed as a dimension of the electromagnetic circuit and material property, which includes the B-H curve and coercive force. The B-H curve indicates that the permeability changes with flux density in the armature and magnet house. When the position of the armature end changes, the reluctance of the air gap between armature and permanent magnet is changed, which means the different position contributes to a different flux density and force on the armature. Consequently, the electromagnetic characteristic such as inductance and force factor are nonlinear. Hence, to handle the nonlinear characteristic property of the B-H curve, the under-relaxed Newton–Raphson method is adopted to solve the Kirchhoff's current and voltage laws [21–23].

By post processing flux density, the cogging force, force factor, and inductance can be obtained which are depicted in Figure 2.

**Figure 1.** The equivalent circuit of the electromagnetic part.

**Figure 2.** Electromagnetic parameters vs displacement.

The cogging force, force factor, and inductance expressions are listed in the following equations:

$$\begin{aligned} F\_{\text{clog},\text{sing}}(\mathbf{x}) &= 2\mu\_0 A\_{\mathcal{S}} F\_M^2 \frac{\left(\mathbf{D}\_{eff2} + \mathbf{D}\_{eff1} + \mathbf{4D}\_{eff}\right) \left(\mathbf{D}\_{eff2} - \mathbf{D}\_{eff1} + \mathbf{2x}\right)}{\left[\left(\mathbf{D}\_{eff1} - \mathbf{x}\right)\left(\mathbf{D}\_{eff2} + \mathbf{x}\right) + \mathbf{D}\_{eff}\left(\mathbf{D}\_{eff1} + \mathbf{D}\_{eff2}\right)\right]^2} \\ L\_E(\mathbf{x}) &= \mu\_0 A\_{\mathcal{S}} N^2 \frac{\left(\mathbf{D}\_{eff2} + \mathbf{D}\_{eff1}\right)}{\left(\mathbf{D}\_{eff1} - \mathbf{x}\right)\left(\mathbf{D}\_{eff2} + \mathbf{x}\right) + \mathbf{D}\_{eff}\left(\mathbf{D}\_{eff2} + \mathbf{D}\_{eff1}\right)} \\ B(\mathbf{x}) &= \frac{\mu\_0 A\_{\mathcal{S}}}{2} \frac{F\_M N \left[\left(\mathbf{D}\_{eff2} + \mathbf{D}\_{eff1} + \mathbf{4D}\_{eff}\right)\left(\mathbf{D}\_{eff2} + \mathbf{D}\_{eff1}\right) + \left(\mathbf{D}\_{eff2} - \mathbf{D}\_{eff1} + \mathbf{2x}\right)^2\right]}{\left[\left(\mathbf{D}\_{eff1} - \mathbf{x}\right)\left(\mathbf{D}\_{eff2} + \mathbf{x}\right) + \mathbf{D}\_{eff1}\left(\mathbf{D}\_{eff2} + \mathbf{D}\_{eff1}\right)\right]^2} \end{aligned} \tag{1}$$

where

$$\begin{aligned} D\_{eff1} &= \mu\_0 A\_\% (R\_M + R\_{Hi1}) + D \\ D\_{eff2} &= \mu\_0 A\_\% \left( R\_M + R\_{Hi2} + \frac{1}{2} R\_H \right) + D \\ D\_{eff} &= \mu\_0 A\_\% (R\_A + R\_i + R\_{ii}) \end{aligned}$$

where *Fcogging*, *Bl* and *LE* are cogging force, force factor, and voice coil inductance, respectively, which are expressed by the armature position *x*. *RA*, *RH* are the reluctance of the armature and magnet house. *FM*, *Ag*, *N* and *D* are the magnetomotive force of the magnet, the area of the magnet, the coil turns, and the air gap width respectively.

The mathematic equation for the electrical part is given in the following equation:

$$Z(s) = R\_E + sL\_E \tag{2}$$

where *RE* is the electrical voice coil resistance at DC and *LE* is the voice coil inductance.

#### *2.2. Mechanical Analysis*

The modeling of the mechanical simulation is demonstrated in Figure 3. The structure contains the armature, pin, and diaphragm. Magnetic force is generated on the armature end. The diaphragm vibrates through the pin connection. The edge of the diaphragm is fixed.

**Figure 3.** Mechanical simulation modeling.

To describe the mechanical system, one degree-of-freedom of the vibration system's governing equation is adopted as follows:

$$F = M\_{\rm ms}\ddot{\mathbf{x}} + R\_{\rm ms}\dot{\mathbf{x}} + \frac{1}{\mathcal{C}\_{\rm ms}}\mathbf{x} \tag{3}$$

where *Mms* is the mechanical mass of the driver diaphragm, *Rms* denotes the mechanical resistance of the total driver losses, and *Cms* denotes the mechanical compliance of the driver suspension.

In the modeling, if the input force is 0.01 N, there will be 6.38 × <sup>10</sup>−<sup>3</sup> mm displacement on the output point. If the input force is 0.02 N, the displacement becomes 12.76 × <sup>10</sup>−<sup>3</sup> mm. Hence, the stiffness of the mechanical system is calculated by dividing the difference of the force with the displacement, which is 1570 N/m and defined as original stiffness (*Koriginal*). The cogging force stiffness is 720 N/m, which can be treated as negative stiffness. Therefore, the modified stiffness is

850 N/m, which means *Cms* is 1.176 m/N. By conducting a modal analysis of the mechanical system, the resonance frequency is obtained, which is 3169.9 Hz. According to the following equation:

$$f\_0 = \frac{1}{2\pi} \sqrt{\frac{K\_{original}}{M\_{ms}}}\tag{4}$$

The mass is calculated as 3.95 × <sup>10</sup>−<sup>6</sup> kg.

#### *2.3. Acoustic Analysis*

The acoustic domain is modeled as described in the following sections:

There are four tubes in the acoustic modeling which are depicted in Figures 4 and 5. These tubes can be treated as transmission line models which are depicted in Figure 6.

**Figure 4.** Acoustic modeling in the BAR.

**Figure 5.** Acoustic modeling in the test jig.

**Figure 6.** Tube model.

The governing equation is listed below:

$$
\begin{pmatrix} p\_1 \\ q\_1 \end{pmatrix} = \begin{bmatrix} \cos(kl) & jZ\_w \sin(kl) \\ (j/Z\_W)\sin(kl) & \cos(kl) \end{bmatrix} \begin{pmatrix} p\_2 \\ q\_2 \end{pmatrix} \tag{5}
$$

and

$$\begin{aligned} p\_1 &= \cos(kl) \times p\_2 + jZ\_w \sin(kl) \times q\_2\\ q\_1 &= (j/Z\_w)\sin(kl) \times p\_2 + \cos(kl) \times q\_2 \end{aligned} \tag{6}$$

The parameters are *k* = *ω*/*c*, where *ω* = 2*πf*; *l* = tube length; and *Zw* = *ρc*/*S*, where *ρ* is the air density and *c* is the speed of sound propagation.

The 2-cc coupler is modeled as the acoustic capacity, which is shown in Figure 7.

**Figure 7.** Equivalent circuit of cavity.

The related equation defining the acoustic capacity is given as follows:

$$\mathbf{C}\_{a} = \frac{V}{\rho c^{2}}\tag{7}$$

The SPL (Sound Pressure Level) can be calculated as below:

$$SPL = 20\log\left(p/20 \times 10^{-6}\right) \tag{8}$$

#### *2.4. Electromagnetic-Mechanical–Acoustic Coupling Factors*

The force factor is the electromagnetic–mechanical coupling factor and treated as a gyrator which generates a back EMF in the electromagnetic domain and a driving force in the mechanical domain. The effective area is the mechanical–acoustic coupling factor and is treated as the transformer which changes the vibration velocity of the diaphragm into the volume velocity of the air motion. The effective area of the diaphragm is calculated as half of the area of the diaphragm, because it vibrates just like a cantilever beam. The total simulation tool is depicted in Figure 8 which contains the electromagnetic, mechanical, and acoustic domains along with the coupling factor.

**Figure 8.** Equivalent circuit of BAR.

#### **3. Analysis by FEM and BEM**

To check the efficacy of LPM, the BAR is also analyzed by FEM and BEM. In the electromagnetic simulation part, the BAR is modeled with a different deformation, i.e., different displacement. Varying current is then input to the coil. The magnetic vector potential is defined as zero on the surface of the surrounded air box which is shown in Figure 9. Finally, the flux density is solved through FEM simulation. The electromagnetic FEM simulation governing equation is demonstrated by the following equation:

$$\nabla \times \left(\frac{1}{\mu\_r} \nabla \times \mathbf{A}\right) = \mathbf{J}\mu\_0 + \mu\_0 \nabla \times \mathbf{H}\_c \tag{9}$$

where **A** is the magnetic vector potential, **J** is the current density, **H***c* is the permanent magnetic intensity, *μ*<sup>0</sup> is the permeability in the vacuum. After solving the governing equation, the flux density of every node in the modeling can be obtained. By post process, the force factor, inductance, and cogging force are obtained, which become the input of the vibro-acoustic simulation.

**Figure 9.** Modeling of electromagnetic FEM simulation.

The governing equation of the mechanical FEM simulation is listed in equation

$$[\mathcal{M}]\ddot{u} + [\mathcal{C}]\dot{u} + [\mathcal{K}]u = [F\_i] \tag{10}$$

where [*M*], [*C*], [*K*], {*Fi*}, and u denote the mass matrix, damping matrix, stiffness matrix, vector of current force, and displacement. The surround part of the diaphragm is fixed which means that the displacement is zero. With a given input force, the displacement of every node can be solved.

In the acoustic domain, the velocity of the tube surface is defined as zero, which means a rigid wall and is depicted in Figure 10. The sound pressure on the test point is solved by solving the Helmholtz governing equation.

$$
\nabla^2 p - k^2 p = -j\rho\_0 w q \tag{11}
$$

where *p*, *k*, *ρ*0, *w*, and *q* are sound pressure, wave number air density, angular frequency, and volume velocity. The volume velocity is from the displacement result obtained in the mechanical FEM simulation. The details are outlined in Tables 2 and 3.

**Figure 10.** Modeling of acoustic BEM simulation.




Table 2 shows the simulation detail. ANSYS is taken advantage of in the electromagnetic and mechanical system and Virtual Lab is adopted to simulate the acoustic system. In order to obtain an accurate result, the numbers of nodes in the three systems are 146817, 525006, and 1329. Based on the same computer and the listed node number in Table 2, the computation time by FEM and BEM is 36 h while, the time by the LPM method is 1 hour. So the computation time is improved.

#### **4. Experiment**

Furthermore, according to 3D modeling, the balanced armature sample is manufactured and assembled. The sample is shown in Figure 11. The experimental condition is demonstrated in Figure 12. The sweep input source ranges from 100 Hz to 20,000 Hz. The SPL of the BAR is tested using a microphone through the 2-cc coupler jig. In the experimental result, the SPL can maintain 100 dB because the BAR is tested in an enclosed tube and chamber. The first SPL peak is due to resonance of the mechanical structure. The second SPL peak is due to the tube and front chamber. Figure 13 shows the comparison of SPL between the experiment and simulation. The LPM simulation results were verified through experimental results, using FEM and BEM.

**Figure 11.** Parts of BAR.

**Figure 12.** Experimental condition.

**Figure 13.** Comparison of experiment and simulation.

#### **5. Conclusions**

The BAR is a product based on multiphysics that considers electromagnetic, mechanical, and acoustic domains. These physical domains are coupled with each other. This study derived LPM simulation, which modeled the multiphysics characteristic and considered the coupling effect. The SPL performance of the BAR was predicted through LPM, which provided the same validation as FEM and BEM. Additionally, it was found that LPM was more efficient and could be used to develop the BAR. In the future, the proposed LPM method will be discussed in relation to the BAR dimension and material property influence on SPL. For example, in the electromagnetic domain, the sensitivity analysis of the magnetomotive force of the magnet will be done to check its influence on SPL. Finally, according to the design target, every dimension can be determined.

**Author Contributions:** Conceptualization, Y.-W.J. and D.-P.X.; methodology, Y.-W.J.; software, Y.-W.J.; validation, Y.-W.J., Z.-X.J. and J.-H.K.; formal analysis, Y.-W.J.; D.-P.X. investigation, Y.-W.J.; D.-P.X. resources, Y.-W.J.; D.-P.X. data curation, Y.-W.J.; D.-P.X. writing—original draft preparation, Y.-W.J.; writing—review and editing, Y.-W.J., Z.-X.J. and J.-H.K.; visualization, Y.-W.J.; supervision, S.-M.H.; project administration, S.-M.H.

**Funding:** This research received no external funding.

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

*Appl. Sci.* **2019**, *9*, 839

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


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

#### *Article*
