A Preliminary Numerical Study on the Interactions Between Nonlinear Ultrasonic Guided Waves and a Single Crack in Bone Materials With Motivation to the Evaluation of Micro Cracks in Long Bones

Fatigue damage in a bone occurs in the form of micro-scale cracks with the lengths as small as to a few microns. The evaluation of cracks in bones has recently been a hot topic. However, the current most frequently used method is based on traditional linear ultrasound, which is just sensitive to gross damages rather than micro cracks. Nonlinear ultrasonic technique, which is capable of detecting micro-scale damages, has been widely used in metallic structures. However, few study has been directed to employing second harmonic generation of nonlinear ultrasound to evaluate fatigue damages in bones. In this study, a preliminary study is conducted on the interactions between nonlinear guided waves and a single crack in bone materials motivating to the evaluation of micro cracks in long bones. Considering the symmetry, asymmetry, location, orientation, length and width of the crack in bone materials, their influences on second harmonic responses are discussed in detail. Our results are presented as follows. Firstly, not only the primary S0 but also A0 mode Lamb wave generates the second harmonic component of S0 mode after interacting with a symmetric crack. Secondly, along the thickness direction, the primary S0 mode Lamb wave possesses almost the same detection sensitivity, no matter where the crack is located. Thirdly, the amplitudes of S0 mode second harmonic component are increased slightly when the oblique angle of the crack is set from 0°to 45°, while the amplitudes increase dramatically when the oblique angle changes to 67.5° and 90°. Lastly, the second harmonic amplitude and the relative acoustical nonlinear parameter are increased with the length of the crack, while they possess a monotonically decreasing relationship with the width of the crack. Our preliminary studies will provide priori knowledge when using nonlinear ultrasonic guided waves for detecting micro-scale cracks in long bones in future clinical inspections.


I. INTRODUCTION
Fatigue damage in a bone occurs in the forms of micro-scale cracks with the lengths of 5-500 µm [1], [2].This micro crack damage contributes to the formation of stress fractures and acts as a stimulus for bone remodeling [3].Bones, therefore, The associate editor coordinating the review of this manuscript and approving it for publication was Jun Shi .
have an outstanding advantage over most engineering structures in that they possess an inherent ability to repair damage.However, if this crack damage accumulates at a rate that the capacity for repair is exceeded, stress fractures result.These fractures occurred commonly among the people who undertake intensity and repetitive activities, such as athletes, soldiers and gymnasts [4].If this crack accumulates at normal rates but the bone's repair mechanism is deficient, fragility VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/factures result, which occur commonly in ageing bone [5].Consequently, the evaluation and detection of cracks in bones prior to fracture is of great urgency and importance.Conventional radiography, e.g., dual X-ray absorptiometry (DXA), which has been established as a reliable means of measuring bone density, is one of the most common method to evaluate bone fracture [6], [7].Ultrasonic testing technique has some advantages of non-ionizing radiation, portability and low cost, which holds a promise to become a viable alternative for radiography [8], [9].In the past decade, the ultrasonic axial transmission technique has been received considerable interest for the assessment of bones [10]- [12].The so-called first arriving signal is one contribution to the ultrasonic axial transmission technique.It has conventionally been used to measure the parameters of bones, such as mineral density and bone geometry [13]- [15].Besides of the first arriving signal, bone structures also support the propagation of ultrasonic guided waves.It has been demonstrated that ultrasonic guided waves propagating along the bones provides more comprehensive information on bone material and structural characteristics [16], [17].The characteristics of ultrasonic guided waves propagating in bones have been clearly described using theoretical analysis, numerical and experimental studies [18]- [21].Furthermore, ultrasonic guided wave technique has been applied to assess the quality of bones with significant results [22]- [25].
Almost all of the above studies have been restricted to use linear ultrasonic technique to evaluate the status of bones.It is well known that linear ultrasound is mainly dependent on measuring particular parameters, such as sound velocity, attenuation, and/or the transmission and reflection coefficient of propagating waves.These parameters are just sensitive to gross and macro-scale damages within the bone structures.Consequently, ultrasonic methods based on linear theory are not suitable to detect micro-scale cracks in bones.
Nonlinear ultrasonic method is much more sensitive to micro-scale damages than linear technique [26]- [30].Recently, nonlinear ultrasonic techniques have been applied to evaluate bone damage status including nonlinear resonant ultrasound spectroscopy [1], [31], nonlinear dynamic response [32] and nonlinear modulation measurements [33], [34].These studies are mainly based on using bulk waves.Furthermore, in addition to the above mentioned nonlinear effects, higher harmonics generation (normally second harmonic generation) as one of the most frequently nonlinear behaviors in NDE (Nondestructive Evaluation) community [35]- [37] has seldom been used for the evaluation of damages in bones.Zhang et al. [38] experimentally studied nonlinear Lamb waves propagating in long bones and clearly observed accumulated second harmonic signals from the primary Lamb waves.The nonlinear ultrasonic guided waves technique combines the high sensitivity to micro damages of nonlinear method and the expeditious inspection of guided waves.It has thus great potential for the assessment of bone structures.Up to now, according to our literature review, few study has been directed to employing second harmonic generation of nonlinear ultrasonic method to evaluate fatigue damages in bone structures.
The objective of this study is to investigate the feasibility and the effectiveness of using nonlinear ultrasonic guided waves for the evaluation of micro-scale cracks in bone materials.The characteristics of second harmonic generation from the primary S0 and A0 mode Lamb waves interacting with a single crack in bone materials are investigated.The influences of the symmetry, asymmetry, location, orientation, length and width of the crack in bones on the second harmonic generation are studied, respectively.Our study provides motivation to the evaluation of micro cracks in long bones.The rest of this article is organized as follows.Theoretical fundamentals are briefly introduced in section II.Numerical setups are presented in section III.Results and discussions are illustrated in section IV. Conclusions and future work are illustrated at last in section V.

II. THEORETICAL FUNDAMENTALS A. DISPERSION CURVES IN A PLATE WITH BONE MATERIALS
Previous studies [16], [17], [24], [38] have demonstrated that bone structures support the propagation of ultrasonic guided waves.Long bones are treated as circular hollow pipe structures and ultrasonic guided waves propagating in long bones are studied by Ta et al. [16], [17].It is found [16] that when the ratio of inner radius to thickness of a circular hollow pipe is increased, the propagation characteristics of a pipe are becoming similar to those of a plate-like structure.Recently, more and more researchers [14], [25], [38]- [40] have directed to use plate-like structure model for numerical study of bones instead of circular tube structure.In this study, two dimensional plate structure model with bone materials is employed for numerical simulations.
Ultrasonic guided waves propagating in a solid plate are known as Lamb waves.They are commonly grouped as symmetric and asymmetric modes according to the different vibrations features, governed by the Rayleigh-Lamb frequency equations [41]- [43].From the Rayleigh-Lamb equations, phase and group velocities dispersion curves can be derived.
The density, the velocities of longitudinal and shear waves of a cortical bone are 1500 kg/m 3 , 4060 m/s and 1840 m/s, respectively [39].The obtained dispersion curves of a bone plate with the thickness of 4 mm are presented in Figure 1.It is found that Lamb waves are composed of two modes (S0 and A0 modes) within the low frequency range (0-200 kHz).It is know that [44] the zero-order modes (S0 and A0 modes) carry more energy than high-order modes, with a smaller energy attenuation during the propagation distance compared to high-order modes.Therefore, S0 and A0 modes are selected as the primary waves to investigate second harmonic generation effect from the interactions between the Lamb waves and a crack.

B. ACOUSTICAL NONLINEAR PARAMETER
A large number of previous studies have shown that nonlinear ultrasonic technique based on second harmonic generation is effective to detect the microstructural changes in metallic materials.The second harmonic waves are generated when the fundamental ultrasonic waves pass through a material with micro-scale damages.The micro-scale damage can be quantitatively evaluated by acoustic nonlinear parameter which is involved by the amplitude of the second harmonic waves.
Physically, the phenomenon of second harmonic generation is related to nonlinearity in the elastic behavior of the material, which indicates that the relationship between stress σ and strain ε is nonlinear.For one dimensional case, the nonlinear stress-strain relationship can be expressed as where E is the linear elastic modulus of materials and β denotes the second-order classical nonlinearity parameter.By combining the wave motion equation with Eq. ( 1), one dimensional nonlinear wave motion equation of longitudinal waves can be obtained.The second-order nonlinear parameter in terms of the amplitude of second harmonic waves can be derived and expressed as [35] where A 1 refers to the amplitude of fundamental waves, A 2 is the amplitude of the second harmonic waves and k denotes the wave number of the fundamental waves.
Considering the propagation of Lamb waves, a relative acoustical nonlinearity parameter β is used to characterize the damage state of a wave guide.In the previous studies, two different definitions are employed.The first definition considers the slope of the amplitude's ratio with the propagation distance and the second considers the amplitude's ratio at a fixed propagation distance at a specimen

C. CONTACT ACOUSTICAL NONLINEARITY CAOUSED BY CRACKS
When an ultrasonic wave excited by a large amplitude is incident to an imperfect interface, higher harmonic waves are generated.This phenomenon is known as CAN (Contact Acoustic Nonlinearity), and has attracted increasing amounts of attention for its potential to characterize closed cracks or imperfect bond interfaces.
The basic physical mechanism of CAN is that a crack driven by longitudinal acoustic traction causes clapping of the crack interface.This clapping nonlinearity originates from asymmetrical dynamics of the contact stiffness which is higher in the compression phase than in the tensile phase.As a result, the compressional part of the waves can penetrate it, but their tensile part cannot.Therefore, after penetrating the interface, the waves exhibit half-wave rectification, which means that they have obvious nonlinearity.This nonlinearity can then be detected by second harmonics [35].
The CAN is generated at the crack.Unlike the material nonlinearity which is distributed across the whole waveguide, it is a kind of localized nonlinearity.Therefore, the Eq. ( 4) which does not include the propagation distance is employed to quantify the CAN.

III. NUMERICAL SETUPS A. THE ALID TECHNIQUE
In this study, long bone is simplified to two dimensional plate model as discussed in section II.Lamb waves propagating in long bones and interacting with micro cracks to generate nonlinear effect of second harmonic are studied.In order to eliminate unwanted boundary reflections, ALID (Absorbing Layer using Increasing Damping) technique is applied to two dimensional bone plate finite element models.In addition to ALID, the infinite elements method is another approach to VOLUME 8, 2020 remove the unwanted reflections.However, previous studies [45], [46] have concluded that infinite elements are not suitable for high accuracy removal of unwanted reflection of bulk waves and ultrasonic guided waves.The ALID technique is thus employed in this study.
The ALID is an absorbing layer that is made of a material with the same properties with those of the area of study except for having a gradually increasing damping.
The equation of dynamic equilibrium in the time domain is expressed as with [M], [C] and [K] denoting the mass, damping and stiffness matrices.
Stiffness or mass proportional damping can be introduced in time domain finite element models and it is generally termed as Rayleigh damping.Consequently, the damping matrix [C] can be expressed as where α and β denote the mass and stiffness proportional damping coefficients.
In an ALID with a boundary perpendicular to the x axis, the value of α and β are gradually increased in the x direction.The following formulations are set.
where α max and β max are positive real numbers and X (x) varies from 0 at the interface between the ALID and the area of study to 1 at the end of the ALID following a power law whose order is defined by m.
It is noted that the introduction of damping decreases with the value of the stable time increment when solving the finite element model with central difference explicit scheme [47].The damping value at the end of an ALID is usually very large compared to the values commonly used in the structures.A high value of α causes a relatively small decrease in the stable increment whereas a value of β usually has a very strong effect leading to a great loss in computational efficiency.Therefore, it is preferable to avoid using β to define ALID with an explicit scheme.In this article, we only have α for numerical studies.The Eq. ( 7) changes to the following formulation It is obvious to see that proper definition of the layer parameters, i.e., the length of the layer L a , variation of the attenuation parameter α and the power law m is essential to achieve an efficient and accurate model.In Finite element models, as the space is discretized, the gradual increase of α occurs by steps.An ALID is defined as a series of sub layers having the same material properties but different values of α.It is preferable to minimize the change of α between two adjacent sub layers.It is recommended to have one element thick sub layers.
Suppose an ALID with the length L a has n sub layers, and the length of each sub layer is l a , the attenuation parameter of i th sub layer α (i) is defined as.
where i varies from 1 to n.The value of i equals to 1 corresponding to the sub layer next to the interface and n corresponding to the sub layer at the end of the ALID.According to the reference [45], in order to achieve an efficient and accurate model, these parameters can be selected as follows.
Normally, α max is selected to larger than 10 f 0 , where f 0 denotes the excitation frequency.L a is set to larger than 2λ, in which λ refers to the wavelength.To eliminate the reflections of second harmonics, α max is set to larger than 20 f 0 and L a to 4λ, The length of the sub layer l a equals to element size.The value of m is set to 2 or 3.

B. FINITE ELEMENT MODEL
Numerical simulations for studying nonlinear ultrasonic Lamb waves interacting with a single crack in a plate bone structure are performed by using ABAQUS software.The schematic of two dimensional finite element model is illustrated in Figure 2 (a).The length and thickness of the bone plate is set to 500 mm and 4 mm, respectively.Two ALID regions are applied at both ends with a length of L a .In this study, the value of L a is set to more than 4λ, where λ referring the length wave of the excitation signal.The crack region is located at distance of 100 mm from the excitation.Vertical, horizontal and oblique crack models are established and integrated in the finite element model.They are presented at Figures 2 (b), (c) and (d), respectively.The shape of the crack is modeled as an ellipse [37], [48].The surfaces of the crack are simulated by hard contact with a frictionless model.The major and minor axis of the ellipse are defined as the length l c and width w c of the crack, respectively.It is noted that in order to acquire obvious second harmonic generation, the crack with millimeter level of the length and micron grade of width is introduced to the finite element model.The traction excitation is located at 50 mm from the left end of the bone plate.Symmetric and anti-symmetric excitations are applied to generate the primary S0 and A0 mode Lamb waves, respectively.They are clearly illustrated in Figures 2(e) and 2(f), respectively.

C. ELEMENT SIZE AND TIME STEP
In general, a higher-order element type, a denser mesh and smaller time step will result in a more accurate result, but will also cost more in terms of calculation time and computer resources.In order to obtain adequate accuracy and high efficiency, a second-order rectangular element type is used to discrete the bone plate and the maximum element size and time step is adopted according to the references [49], [50].
where I is the element size and t the time step; λ min and f max are shortest wavelength and highest frequency of interest, respectively.
To ensure the adequate accurate accuracy of second harmonic responses, f max is referred to the upper limit of the frequency bandwidth of the second harmonic.In this study, the primary centre frequency of excitation signal f 0 is set to 100 kHz.According to Eq. ( 11), the element size setting to 0.2 mm for the normal region of the model is sufficient to ensure accuracy.The value of time step setting to 5e-8 s is also adequate.The crack zone is more densely meshed, with much smaller elements to accommodate the complicated mechanical response.The element size at the crack is set to 0.02 mm.The meshing results at the normal and cracked regions are presented in Figure 3.

D. SCHMETIC DIAGRAM FOR DERIVING THE RELATIVE ACOUSTICAL NONLINEARITY PARMATER β
In this study, the primary Lamb wave modes are S0 and A0 at low frequency range.Due to their dispersive characteristics, there is always a difference of the velocities (phase and group velocities) between the fundamental and the second harmonic components.There will be a time flight gap between the fundamental and second harmonic components in the received time domain waveforms.Applying FFT (Fast Fourier Transform) directly to the received temporal waveforms to derive the amplitudes of the fundamental (A 1 ) and second harmonic (A 2 ) components will result in a poor accuracy of the relative acoustical nonlinearity parameter β .Therefore, a schematic diagram for deriving accurate relative acoustical nonlinearity parameter β is proposed.The proposed idea is that we firstly extract time domain waveforms of the fundamental and second harmonic components from the received signal, and then apply FFT to the temporal waveforms of fundamental and second harmonic components to get the corresponding amplitude, respectively.The schematic flowchart is shown in Figure 4 and explained in detail as following steps.
Step 1: Numerical simulations are conducted under the excitations with both positive and negative polarities of the tone burst as illustrated in Figures 4 (a) and (b).The corresponding received waveforms denoted by x (t, F 0 ) and x (t, −F 0 ) are shown in Figures 4 (c) and (d), respectively.For the cracked plate bone structure, in addition to the fundamental components, the received signals also include the waveforms of nonlinear responses which are the static displacement component [29], [51] and high harmonic components (mainly second harmonic component).
Step 2: The equation ) is used to remove the fundamental component of the received signal, and derive the signal which just includes the static and second harmonic components [52], [53].An example illustrating the superimposition of the static and second harmonic components is presented in Figure 4 (e).The equation ) is employed to eliminate the static and second harmonic components, and the fundamental component in the received signal is obtained.An example of the obtained fundamental component is shown in Figure 4 (f).x 02 (t) refers to the derived waveforms of the superimposition of static displacement and second harmonic component.x 1 (t) denotes the extracted fundamental component.
Step 3: The EMD (Empirical Mode Decomposition) method [54], [55] is used to decompose the signal of x 02 (t) into a series of IMFs (Intrinsic Mode Function) and a residue.The residue represents the temporal waveform of the static component.The red curve in Figure 4 (g) demonstrating an example of extracted static displacement component.The sum of the IMFs is the results of the extracted time domain waveform of the second harmonic component with an example shown in Figure 4 (i).
Step 4: By applying FFT to the extracted time domain waveforms of the fundamental and second harmonic components, respectively, the frequency spectrums of the fundamental and second harmonic components are obtained.
The corresponding examples are illustrated in Figures 4 (h) and (j).From these two spectrums, the amplitudes of the fundamental (A1) and second harmonic (A2) components can be derived.Step 5: By using the formulation β = A 2 A 2 1 , the relative acoustical nonlinearity parameter β can be derived.

IV. RESULTS AND DISCUSSIONS
In this section, nonlinear effect of second harmonic generation from the primary Lamb waves interacting with a crack in a plate with bone material is investigated.The dependence of second harmonic generation on the symmetry, asymmetry, location, orientation, length and width of a crack in the bone plate is discussed, respectively.The proposed schematic diagram in Figure 4 is used to extract the time domain waveof the fundamental and second harmonic components in the received waveforms and to derive the relative acoustical nonlinearity parameter.

A. SECOND HARMONIC GENERATION FROM THE PRIMARY LAMB WAVES INTERACTING WITH A SYMMETRIC CRACK
In this subsection, the primary S0 and A0 mode Lamb waves at the center frequency of 100 kHz interacting with a vertical 169174 VOLUME 8, 2020 symmetric crack in a bone plate are studied.The length and width of the crack are set to 2 mm and 0.1 µm, respectively.As illustrated in Figure 2 (a), along the thickness direction, the center of the crack is at the middle of the bone plate.It means that the center of the crack is located at 2 mm from the top surface.The crack is thus symmetric.Along the length direction, the crack locates between the excitation and the monitoring points.The monitoring point at a distance of 150 mm from the crack is used to acquire the received waveforms.
Figures 5 (a) and (b) illustrate the extracted temporal waveforms of the fundamental and second harmonic components from the received signal after the primary S0 and A0 mode Lamb waves interacting with a symmetric vertical crack, respectively.The blue and red curves denote the corresponding fundamental and second harmonic components, respectively.In Figure 5 (a), it is clearly observed that the time domain wave packet of second harmonic lags a little behind the fundamental component.The fundamental component wave packet appears at 7.99e-5 s, and it propagates 250 mm.Therefore, the calculated group velocity is 250 mm / 7.99e-5 s = 3129 m/s.It is quite close to the theoretical group velocity of S0 mode at 100 kHz with the value of 3190 m/s as illustrated in Figure 1 (b).The primary S0 mode Lamb wave at the center frequency of 100 kHz propagates 100 mm in the long bone and encounters a vertical crack.It then interacts with the crack.After interaction, the second harmonic component is generated.
The second harmonic appears at 8.55e-5 s in the time domain waveform.Suppose the group velocity of the second harmonic component is V 2g1 , 100 mm / 3190 m/s + 150 mm / V 2g1 = 8.55e-5 s.By solving this formulation, V 2g1 equals to 2770 m/s, which is in accordance with the theoretical group velocity of S0 mode at 200 kHz with the value of 2808 m/s as also shown in Figure 1 (b).The error is relative quite small and acceptable.It is reasonable and thus verified that the primary S0 Lamb wave generates the mode second harmonic component.
In Figure 5 (b), it is obviously found that the wave packet of second harmonic component appears ahead of the wave packet of the fundamental component.The fundamental component is appeared at 1.34 e-4 s.Its group velocity is thus calculated as 250 mm / 1.34e-4 s = 1866 m/s, which is close to the theoretical group velocity of A0 mode at the frequency of 100 kHz with the value of 1808 m/s.Similar to the primary S0 mode, the mode of the second harmonic generated from the primary A0 mode interacting with the crack is identified as follows.The wave packet of second harmonic component appears at 1.013e-4 s, and its group velocity is supposed to be V 2g2 .The calculated formulation is V 2g2 = 150 mm / (1.08e-4 s -100 mm / 1808 m/s) = 2847 m/s.This value is also accord with the theoretical group velocity of S0 mode at the frequency of 200 kHz with the value of 2808 m/s.The error is also quite small and acceptable.Therefore, it is verified that the primary A0 mode Lamb waves interacting a symmetric crack in a bone plate generates S0 mode second harmonic component.This finding accords to the results reported in previous published reference [56] that second harmonic of S0 mode is generated when the primary A0 mode Lamb waves propagating in a carbon steel plate with considering evenly distributed material nonlinearity.

B. SECOND HARMONIC GENERATION FROM THE PRIMARY LAMB WAVES INTERACTING WITH AN ASYMMETRIC CRACK
In this subsection, the primary S0 and A0 mode Lamb waves at the center frequency of 100 kHz interacting with an asymmetric crack in a bone plate are investigated.The length and width of the crack are 1.6 mm and 0.1 µm, respectively.As illustrated in Figure 2 (c), along the thickness direction, the center of the crack is located at 1 mm from the top surface.Therefore, the crack is thus asymmetric.Along the length direction, the monitoring point locating at a distance of 300 mm from the crack is employed to acquire the received waveforms.
Figures 6 (a) and (b) show the extracted time domain waveforms of the fundamental and second harmonic components from the received signal after the primary S0 and A0 mode Lamb waves interacting with an asymmetric vertical crack, respectively.In Figure 6 (a), there are two wave packets for time domain waveforms of both the fundamental and second harmonic components.For the fundamental component, the first wave packet appears at 1.253e-4 s.The propagation distance is 400 mm.Therefore, the calculated group velocity is 400 mm / 1.253e-4 s = 3192 m/s, which is in accordance with the theoretical group velocity value of S0 mode at 100 kHz.It is thus proved that the first wave packet of the fundamental component is the direct S0 mode.Similarly, the second wave packet with a small amplitude of the fundamental component is identified as the A0 mode, which is converted from the primary S0 mode due to the asymmetry of the crack.Regarding to the time domain waveform of the second harmonic component, there are also two wave packets.The first wave packet lags a little behind the first corresponding wave packet of the fundamental component.Its calculated group velocity is close to the theoretical group velocity of S0 mode at 200 kHz.Therefore, the first wave packet of second harmonic is S0 mode.It is generated due to the primary S0 mode interacting with the asymmetric crack.As the group velocity of S0 mode Lamb waves at 100 kHz is larger than at 200 kHz, the second harmonic wave packet lags a little behind the fundamental wave packet in the time domain waveforms.Furthermore, the second wave packet of the second harmonic component almost coincides with the corresponding second wave packet of the fundamental component.It is identified that the second wave packet of the second harmonic component is A0 mode.It is converted from the second harmonic of S0 mode due to the asymmetry of the crack.
In Figure 6 (b), for the fundamental component, the first and second wave packets are identified as the S0 and A0 mode, respectively.The presence of S0 mode wave packet of the fundamental component is due to the conversion from the primary A0 mode.As the group velocity of S0 mode at 100 kHz is much larger than A0 mode at 100 kHz, the S0 mode wave packet appears head of A0 mode wave packet.The first and second wave packets of the second harmonic components are also identified as S0 and A0 mode, respectively.The primary A0 mode Lamb waves interacting with the crack generates the second harmonic S0 mode Lamb waves.Meanwhile, due to the asymmetry of the crack, the phenomenon of mode conversion takes place.Part of the second harmonic of S0 mode converts to A0 mode.Consequently, the wave packet of A0 mode Lamb waves are in presence in the extracted time domain waveform of the second harmonic component.
From the above analysis, it is inferred that the asymmetry of the -scale crack in a bone plate will result in the mode conversion of both the fundamental and second harmonic components.The presence of mode conversion phenomenon in the time domain waveforms for both the fundamental and second harmonic components denotes the asymmetry of the crack.

C. SECOND HARMONIC GENERATION FROM THE PRIMARY LAMB WAVES INTERACTING WITH A HORIZONTAL CRACK AT DIFFERENT LOCATIONS ALONG THE THICKNESS DIRECTION
In this subsection, the primary S0 and A0 mode Lamb waves at the center frequency of 100 kHz interacting with a horizontal crack at different locations along the thickness direction is studied.The length and width of the horizontal crack are 2 mm and 0.1 µm, respectively.Along the thickness direction, the distance from the center of the crack to middle plane of the bone structure h 0 is set to 0, 0.3 mm, 0.6 mm, 0.9 mm, 1.2 mm, 1.5 mm and 1.8 mm, respectively.The degree of the asymmetry of the crack is increased.Along the length direction, the monitoring point locating at a distance of 300 mm from the crack is used to record the received waveforms.
Figures 7 (a) and (b) show the extracted time domain waveforms of the second harmonic component from the received signal after the primary S0 and A0 mode Lamb waves interacting with a horizontal crack at different locations along the thickness direction, respectively.The blue, red, green and black waveforms correspond to horizontal crack located at h 0 equal to 0, 0.6 mm, 1.2 mm and 1.8 mm, respectively.h 0 equal to 0 represents that the crack is symmetric, while h 0 equal to other values denotes the asymmetric crack.In Figure 7 (a), the low right part shows the local enlargement of the second wave packets with small amplitudes.The first wave packets are the second harmonic of S0 mode generated from the primary S0 mode.It is observed that the generated S0 mode second harmonic displacements are close when h 0 is equal to 0, 0.6 mm, 1.2 mm and 1.8 mm.From the local enlargement, it is obviously found that the second wave packets appear when h 0 equals to 0.6 mm, 1.2 mm and 1.8 mm.They are the second harmonic components of A0 mode, which are converted from the second harmonic of S0 mode due to the asymmetry of the crack.
In Figure 7 (b), the first wave packets denote the extracted second harmonic of S0 mode generated the primary A0 mode Lamb waves when h 0 equals to 0, 0.6 mm, 1.2 mm and 1.8 mm.Similar to the primary S0 mode, the second wave packets are the converted A0 mode from the second harmonic of S0 mode.It is illustrated that when h 0 is equal to 0, there is no converted second harmonic of A0 mode.The converted second harmonic of A0 mode are negligible when h 0 equals to 0.6 mm.However, when h 0 equals to 1.2 mm and 1.8 mm, the displacements of the converted A0 mode second harmonic components are greatly increased.Maximum displacement of second harmonic component of S0 mode generated the primary S0 and A0 mode Lamb waves with the horizontal crack at different locations are derived.Figure 8 illustrates the trend of the obtained maximum displacement with the h 0 with varied values from 0 to 1.8 mm with a step of 0.3 mm.It is shown that the maximum displacement of S0 mode second harmonic generated from the primary S0 mode is much larger than that generated from the primary A0 mode.This result is in accordance with the time domain waveforms presented in Figure 7. Furthermore, it is also found that the maximum displacement of S0 mode second harmonic component from the primary S0 mode is slightly declined with increased of h 0 , while the maximum displacement from the primary A0 mode is increased as h 0 increases.This finding is reasonable.For the primary S0 mode Lamb waves, the horizontal displacement field across the thickness direction is almost uniform as shown in Figure 9 (a).The displacement amplitude of the generated second harmonic component is thus dependent on the wave structure of S0 mode at 200 kHz.The slightly declined trend of the maximum displacement with h 0 coincides with the lightly deceased horizontal wave structure of S0 mode Lamb waves at 200 kHz as illustrated in Figure 9 (c).The quite slightly change in the displacement amplitude of the second harmonic component can always be negligible.Therefore, the horizontal amplitudes of the second harmonic S0 mode generated from the primary S0 mode Lamb waves interacting with cracks at different locations along thickness direction can be regarded the same.For the primary A0 mode Lamb waves, as illustrated in Figure 9 (b), the increasing trend of the wave structure at 100 kHz from the middle plane to the upper surface contributes to the trend that the maximum displacement from the primary A0 mode increases with h 0 .
From the above analysis, it is inferred that the primary S0 mode Lamb wave possesses almost the same detection sensitivity, no matter where the crack is located along the thickness direction.Regarding that the maximum displacement of second harmonic generated from the primary S0 mode is much larger than that generated from the primary A0 mode, it is recommended that using the nonlinear S0 mode Lamb waves for the detection of cracks in bone structure is preferable.

D. SECOND HARMONIC GENERATION FROM THE PRIMARY LAMB WAVES INTERACTING WITH A CRACK WITH DIFFERENT ORIENTATIONS
In the previous subsections, horizontal and vertical cracks are considered.In this subsection, the primary S0 mode Lamb waves at the center frequency of 100 kHz interacting with a crack with different orientations are investigated.The oblique angle of the crack θ is set to 0 • , 22.5 • , 45 • , 67.5 • and 90 • .θ setting to 0 • and 90 • refers to the horizontal and vertical crack, respectively.The length and width of the horizontal crack are 2 mm and 0.1 µm, respectively.Along the length direction, the monitoring point locating at a distance of 300 mm from the crack is employed to acquire the received waveforms.
Figure 10 illustrates the extracted temporal waveforms of the second harmonic component from the received signal after the primary S0 mode Lamb waves interacting with a crack with different orientations.The blue, red, green, yellow and black curves the crack with the oblique angle of 0 • , 22.5 • , 45 • , 67.5 • and 90 • , respectively.In Figure 10, when θ is set to 0 • and 90 • , the crack is symmetric.There is one wave packet of S0 mode second harmonic in the waveform.On the other hand, there are two wave packets in the waveforms when θ is equal to 22.5 • , 45 • and 67.5 • .The first and second wave packets denote the second harmonic components of S0 and A0 mode, respectively.The second harmonic components of A0 mode are converted from the second harmonic of S0 mode due to the asymmetry of the crack.It is found that the amplitudes of S0 mode second components are increased slightly when θ is set from 0 • to 45 • , while the amplitudes increase dramatically when θ changes to 67.5 • and 90 • .Furthermore, the amplitudes of the A0 mode second harmonic components are also increased with θ .This results from the dominant horizontal displacement filed of S0 mode Lamb waves.

E. THE INFLUENCE OF THE LENGTH OF THE CRACK ON THE SECOND HARMONIC GENERATION
In this subsection, the dependence of the second harmonic generation on the length of the crack is investigated.The primary S0 mode Lamb waves at 100 kHz interacting with a symmetric vertical crack with varied length and constant width are studied.The width of the crack is set to 0.1 µm.The length of the crack is set to 0.5 mm, 1 mm, 1.5 mm, 2 mm, 2.5 mm and 3 mm.Along the length direction, the monitoring point locating at distance of 300 mm from the crack is used to acquire the received waveforms.Following the schematic diagram presented in Figure 4, the amplitude of the second harmonic component and the relative acoustical nonlinear parameter β' are derived.
The amplitudes of the second harmonic components generated from the primary S0 mode Lamb waves interacting with the -scale crack with different length and constant width are illustrated in Figure 11 (a).It is found that the second harmonic amplitude has a monotonically increasing relationship with the crack's length.The longer the crack, the larger amplitude of the second harmonic component becomes.Our findings are consistent with the results reported in the previous reference [57], in which nonlinear longitudinal waves were employed to detect cracks in steel structures.Furthermore, our results are also in accordance with the observation presented in the reference [37], in which nonlinear Lamb waves were used for the detection buried cracks in metallic structures.A straightforward explanation is that second harmonic component is generated due to the contact acoustical nonlinearity which is increased with the length of the crack.As the crack is increased, the contact stiffness of the interface of the crack is declined.According to the theory reported in the previous references [58]- [60], as the contact stiffness of the interface is decreased, that is, as the crack becomes longer, the acoustical nonlinearity is increased.Consequently, the increased second harmonic amplitude with the length of the crack conforms to this theory.The derived relative acoustical nonlinear parameter β' with the length of the crack is shown in Figure 11 (b).It is clearly observed that the relative acoustical nonlinear parameter β' possesses a similar monotonically increasing trend with the crack's length.This results from the quite small changes of fundamental amplitude at different crack's length.

F. THE INFLUENCE OF THE WIDTH OF THE CRACK ON THE SECOND HARMONIC GENERATION
In this subsection, the influence of the crack's width on the second harmonic generation is discussed.The primary S0 mode Lamb waves at 100 kHz interacting with a symmetric vertical crack with varied width and constant length are investigated.The length of the crack is selected as 2 mm.The width of the crack is set to 0.1 µm, 1 µm, 3 µm, 5 µm, 8 µm and 10 µm, respectively.The monitoring point locating at a distance of 300 mm from the crack is used to record the received waveforms.Similarly, following the schematic diagram illustrated in Figure 4, the amplitude of the second harmonic component and the relative acoustical nonlinear parameter β' with the width variation are obtained..The amplitude of the second harmonic component and the relative acoustical nonlinear parameter β' with the width of a crack are illustrated in Figures 12 (a) and (b), respectively.Unlike the positive relationship between the second harmonic amplitude and the length of the crack, the second harmonic amplitude and the relative acoustical nonlinear parameter β' both have a monotonically decreasing relationship with the crack's width.The wider the crack is, the smaller second harmonic amplitude becomes.As the crack becomes wider, the gap between the two interfaces of the crack becomes larger.Consequently, some area of the interfaces may not be in contact during the compressional phase of incident wave, and the contact area is reduced [37].The amplitude of the second harmonic component is thus reduced with the width of the crack.It is noted that the trends of the second harmonic amplitude and β' with the crack's width are similar.The reason is that the fundamental amplitude is almost unchanged with the increment of the crack's width.

V. CONCLUSIONS AND FUTURE WORK
In this study, the second harmonic generation from the interactions between the ultrasonic guided waves and the crack in bone material is studied.The influences of the symmetry, asymmetry, location, orientation, length and width of the crack in bones on the second harmonic generation are investigated, respectively.Finite element model of a bone plate with ALID regions is built.Schematic flow chart to extract the time domain waveforms of fundamental and second harmonic components and to derive the relative acoustical nonlinear parameter based on EMD method is proposed.Several conclusions are drawn and stated as follows.
Firstly, second harmonic responses from the primary S0 and A0 mode Lamb waves interacting with a symmetric and asymmetric cracks are discussed.It is found that not only the primary S0 but also the primary A0 mode Lamb wave generates the second harmonic component of S0 mode after interacting with the symmetric crack.This finding is in accordance with the previous published results that second harmonic of S0 mode is generated when the primary S0 and A0 mode Lamb waves propagating in a carbon steel VOLUME 8, 2020 plate with considering evenly distributed material nonlinearity.Furthermore, when the primary S0 and A0 mode Lamb waves interact with an asymmetric crack, part of the generated S0 mode second harmonic component will converted to A0 mode due the asymmetry of the crack.The wave packet of A0 mode is thus appeared in the extracted time domain waveform of the second harmonic component.
Secondly, second harmonic generation from the primary S0 and A0 mode Lamb waves interacting with a crack at different locations along the thickness direction is investigated.It is observed that along the thickness direction, the primary S0 mode Lamb wave possesses almost the same detection sensitivity, no matter where the crack is located.Further regarding that the maximum displacement of second harmonic generated from the primary S0 mode is much larger than that generated from the primary A0 mode, it is preferred to use the nonlinear S0 mode Lamb waves for the detection of cracks in bone structures.
Thirdly, the primary S0 mode Lamb wave interacting with a -scale crack at different orientations is studied.It is shown that when the oblique angle is set to 22.5 • , 45 • and 67.5 • , A0 mode wave packet appears in the temporal waveforms of the second harmonic component.It is converted form the S0 mode second harmonic component.Furthermore, it is illustrated that the amplitudes of S0 mode second components are increased slightly when the oblique angle is set from 0 • to 45 • , while the amplitudes increase dramatically when the oblique angle changes to 67.5 • and 90 • .
Lastly, the dependences of second harmonic generation on the crack's length and width are explored.It is illustrated that the second harmonic amplitude and the relative acoustical nonlinear parameter are increased with the length of the crack, while they possess a monotonically decreasing relationship with the width of the crack.
Nonlinear effect of second harmonic component generation from the primary ultrasonic Lamb waves interacting with a crack under different circumstances, i.e., the symmetry and asymmetry, the location, the orientations, the length and the width are investigated in detail.Our study results have shown that the use of nonlinear ultrasonic Lamb waves for the detection of cracks in bone materials is feasible and effective.Our results will provide some priori knowledge when using nonlinear ultrasonic guided waves for the evaluation of micro cracks in bones in future clinical inspection.In future work, more accurate model of long bone will be built.Platforms for detecting fatigue damages in bones by using nonlinear ultrasonic guided waves will also be set up.The related experimental studies and practical inspections will be conducted.

FIGURE 1 .
FIGURE 1. Dispersion curves of a bone plate with the thickness of 4 mm: (a) Phase velocity dispersion curves; and (b) Group velocity dispersion curves.

FIGURE 3 .
FIGURE 3. Meshing results in normal and cracked regions.

FIGURE 5 .
FIGURE 5. Extracted time domain waveforms of the fundamental and second harmonic components from the received signal after the primary Lamb waves interacting with a symmetric vertical crack: (a) the primary S0 mode Lamb waves; and (b) primary A0 mode Lamb waves.

FIGURE 6 .
FIGURE 6. Extracted temporal waveforms of the fundamental and second harmonic components from the received signal after the primary Lamb waves interacting with an asymmetric vertical crack: (a) the primary S0 mode Lamb waves; and (b) the primary A0 mode Lamb waves.

FIGURE 7 .
FIGURE 7. Extracted temporal waveforms of second harmonic components from the received signals after the primary Lamb waves interacting with a horizontal crack at different locations along the thickness: (a) the primary S0 mode Lamb waves; and (b) the primary A0 mode Lamb waves.

FIGURE 8 .
FIGURE 8.The trend of maximum displacement of second harmonic component of S0 mode generated the primary S0 and A0 mode Lamb waves with h 0 .

FIGURE 9 .
FIGURE 9. Wave structures of Lamb waves propagating in a bone plate: (a) S0 mode at 100 kHz; (b) A0 mode at 100 kHz; and (c) S0 mode at 200 kHz.

FIGURE 10 .
FIGURE 10.Wave structures of Lamb waves propagating in a bone plate: (a) S0 mode at 100 kHz; (b) A0 mode at 100 kHz; and (c) S0 mode at 200 kHz.

FIGURE 11 .
FIGURE 11.The influence of the length of the crack on the second harmonic generation from the primary S0 mode Lamb waves: (a) the amplitude of the second harmonic components with the length of the crack; and (b) the relative acoustical nonlinear parameter β' with the length of the crack.

FIGURE 12 .
FIGURE 12.The influence of the width of the crack on the second harmonic generation from the primary S0 mode Lamb waves: (a) the amplitude of the second harmonic components with the width of the crack; and (b) the relative acoustical nonlinear parameter β' with the width of the crack.