Simulation of detecting contact nonlinearity in carbon fibre polymer using ultrasonic nonlinear delayed time reversal

A finite element method simulation of a carbon fibre reinforced polymer block is used to analyse the nonlinearities arising from a contacting delamination gap inside the material. The ultrasonic signal is amplified and nonlinearities are analysed by delayed Time Reversal -- Nonlinear Elastic Wave Spectroscopy signal processing method. This signal processing method allows to focus the wave energy onto the receiving transducer and to modify the focused wave shape, allowing to use several different methods, including pulse inversion, for detecting the nonlinear signature of the damage. It is found that the small crack with contacting acoustic nonlinearity produces a noticeable nonlinear signature when using pulse inversion signal processing, and even higher signature with delayed time reversal, without requiring any baseline information from an undamaged medium.


Introduction
In the past, the use of carbon fibre reinforced polymer (CFRP) has been limited to non-structural parts of high-tech aeronautical products. In recent times, due to the effort of weight reduction and product lifetime enhancement, the application areas of CFRP have widened to the load-bearing parts of the aeronautical, automotive and civil engineering products. Due to the increased demands on the strength of the CFRP products and possible complex failure mechanisms, the Non-Destructive Testing (NDT) methods of CFRP have been an important applied and academic problem.
The complex failure mechanisms of CFRP include microcracking and delamination. Microcracking can occur at lower loads or due to aging and can be difficult to examine using ultrasonic NDT. With increased loading, the damage can evolve to delaminations, a very fine cracking between the layers of the CFRP. These damages are difficult to detect us-ing ultrasonic methods due to their small thicknesses. The damage can exhibit itself as a contact acoustical nonlinearity [1]. A statistical distribution of microcracks or delamination damage in the material could also be described by hysteresis in a continuum material model [2,3,4]. This can also be applicable for other materials than CFRP, for example biological tissues [5,6]. Nonlinear ultrasonic methods have been in development for detecting and localizing fatigue and micro-crack damage by their nonlinear effects [7,8]. The detection of harmonic overtones is one of the simplest measures of nonlinearities [9]. Many nonlinear analysis methods not requiring filtering have been developed, for example scaling subtraction method [10,11] or pulse inversion (PI) with its generalizations [12,13], and applications of time reversal using scattering as new sources. This paper proposes a delayed TR-NEWS signal processing method [14] for detecting the nonlinear signature of a single small crack in CFRP as contact acoustical nonlinearity. In the Finite Element Method (FEM) simulation, the CFRP is modelled as anisotropic, layered medium. The ultrasonic signal is focused by TR-NEWS to the region of the material with the defect. The nonlinear signature of the crack is analysed by PI and compared with the delayed TR-NEWS method, which allows to create arbitrary wave envelope at the focusing region of TR-NEWS. It is used here to create an interaction of waves near the damage. The signature of the damage appears as the nonlinear effect of the wave interaction on the contacting crack. This signal processing requires only one transmitting and one receiving transducer. The effectiveness of the delayed TR-NEWS method has been shown in the previous work by physical experiments and simulations in undamaged and linear materials [14]. In this paper, the FEM simulation model is advanced further by including absorbing boundary conditions and the contacting crack defect in the material.

Mathematical and simulation model
This section describes the simulation which is based on a physical experiment, and describes the differences and similarities between the simulation and the experiment. It shows some important points about the mathematical model, the delayed TR-NEWS signal processing and the FEM simulation. Detailed information about the derivation of mathematical and FEM model is available at [15].

Mathematical model
The test object is a CFRP block consisting of 144 layers (Fig. 1). It is composed of fabric woven from yarns of fibre and impregnated with epoxy. The crosssection of the yarns have elliptical shape (Fig. 2) and the material has inclusions of pure epoxy, so a wave propagating through the material will encounter yarns (fibres with epoxy) and areas of pure epoxy. The simulation is in time domain, since the TR-NEWS procedure relies on transient echoes and complex wave motion for the wave energy focusing process. Due to the heavy computational cost of time domain simulation, a simple laminate model is used where: i) the material consists of homogeneous layers, ii) each layer has its own elasticity properties, and iii) dispersion arises due to the periodical discontinuity of the material properties. It consists of CFRP layers with 90 • /0 • weave, 45 • /45 • weave and epoxy layer. The thicknesses of the layers are given by random variable functions which reflect the actual structure of the material. The random variable distribution, describing the CFRP structure, is measured from a close-up image of the CFRP test object [15]. This links the distribution of the microstructure inside the actual material with the thicknesses of the layers in the laminate model. It should enable a more realistic simulated material having dispersion effects due to discontinuities.
The three different kind of layers have the following mechanical properties: i) isotropic pure epoxy: E =  The boundary conditions of the model (Fig. 3) include Lysmer-Kuhlemeyer absorbing boundary conditions [16] so the wave energy would pass through the simulation region. Four degrees of freedom (DOFs) are fixed, the rest are free. The simulation model includes a contacting delamination defect in the material near the receiving transducer (Fig. 4). The transmitting shear wave transducer can send maximum 50 kPa pulse at 70 • degree angle.

TR-NEWS signal processing
This subsection repeats the signal processing method that is applied for this problem and has been published in the previous work [14]. It is included for a self-contained discourse in this paper.
In the physical experiments, on which the simulation is based on, the CFRP block ( Fig. 1) was studied using TR-NEWS NDT equipment and signal processing methods [14]. The 2D FEM simulations reflect  The roles of the transducers are not changed during the experiment: the focusing of the ultrasonic wave relies on the TR-NEWS signal processing. This is a two-pass method where the receiving and transmitting transducers do not change their roles. In this sense the "Time Reversal" describes the signal processing method which accounts for internal reflections of the material as virtual transducers, used for focusing the wave in the second pass of the wave transmission. The placement of the transducers is not important from the signal processing standpoint: in NDT investigation they could be placed arbitrarily and they do not have to be in line with each other, but the configuration must remain fixed during the complete TR-NEWS procedure. Figure 5 outlines the TR-NEWS signal processing steps. The simulation uses the same signal processing steps as are usually applied to physical experiments. Firstly the chirp-coded excitation c(t) is transmitted through the medium.
where ψ(t) is linearly changing instantaneous phase. In this work, a linear sweep from 0 to 2 MHz was used. Then the chirp-coded coda response y(t) with a time duration T is recorded at the receiver where h(t − t , T ) is the impulse response of the medium. The y(t, T ) is the direct response from the receiving transducer when the chirp excitation c(t) is transmitted through medium. Next the correlation Γ(t) between the received response y(t, T ) and chirpcoded excitation c(t) is computed during some time period ∆t . Therefore the actual correlation Γ(t) ∼ h(t) contains information about the wave propagation paths in complex media. (2) Additional echoes coming from interfaces and scatterers in its response R x could be associated to a virtual source T (2) x .
(3) Applying reciprocity and TR process to R x . (4) The time reversed new excitation T x = R x (−t) produces a new response R x (the TR-NEWS coda y T R (t)) with a spatio-temporal focusing at z = 0; y = 0; t = t f and symmetric side lobes with respect to the focusing.
Time reversing the correlation Γ(t) from the previous step results in Γ(−t) used as a new input signal. Re-propagating Γ(−t) in the same configuration and direction as the initial chirp yields where y T R ∼ δ(t − T ) is now the focused signal under receiving transducer where the focusing takes place at time T . This is because Γ(t) contains information about the internal reflections of the complex media, and transmitting its time reversed version Γ(T − t) will eliminate these reflection delays by the time signal reaches the receiver, resulting in the focused signal y T R (Eq. (4)). The test configuration must remain constant during all of these steps, otherwise the focusing is lost. The steps of this focusing process in a physical experiment are shown in Fig. 6.  Figure 6: Bi-layered aluminium experimental chirpcoded TR-NEWS signal processing steps: (1) chirp excitation, (2) output recorded at Rx, (3) crosscorrelation between input and output, (4) focusing resulting from taking time-reversed cross-correlation as new input [14].
PI is an established method for detecting nonlinearities [12]. The procedure used here involves conducting TR-NEWS measurements with positive and negative sign for A in Eq. (1) and comparing the focused signals. Differences could indicate the presence of nonlinearities.
Delayed TR-NEWS signal processing considers a single y T R focusing wave as a new basis which can be used to build arbitrary wave shapes at the focusing. This is done by time-delaying and superimposing n time-reversed correlation Γ(T − t) signals ( Fig. 7 left column) where a i is the i-th amplitude coefficient and τ i the i-th time delay. In case of uniform time delay the ∆τ is the time delay between samples. Upon propagating this Γ s (t − T ) through the media according to the last step of TR-NEWS, a delayed and scaled shape of signal at the focusing point can be created. The delayed TR-NEWS signal processing optimization can be used for amplitude modulation, signal improvement and sidelobe reduction [14].
It is possible to predict what the delayed TR-NEWS focusing output would be in a linear material (Fig. 7 right column): The purpose of the prediction is twofold. Firstly it can be used to figure out optimal delay and ampli- tude parameters a i and τ i beforehand for the delayed TR-NEWS experiment, using the original focusing peak y T R . Secondly it could be possible to analyse the differences between the measured delayed TR-NEWS result and its prediction, which acts as a baseline for comparison. The difference could indicate the magnitude of nonlinearity, because the prediction relies on the applicability of linear superposition and is found to be accurate in experiments with linear material [14].

The FEM simulation model
The simulation program considers 2D wave propagation in a solid material with linear elasticity. The nonlinearity comes from an internal defect, a crack in the computational region (Fig. 4) which can come into contact with itself. This contacting nonlinearity has asymmetric stiffness and is therefore nonclassically nonlinear. Since the CFRP is a complex material, then in this work it is modelled as a laminate with anisotropic layers arranged in a periodic manner, described in Section 2.1. Because the physical experiment was conducted on the corner of a large CFRP block, the simulation is also in a semi-infinite quarterspace. The region has two free surfaces for reflection and two absorbing boundaries for the wave energy to escape.
The constitutive equation of the material itself is linear (although anisotropic). The linear plane strain elastodynamics problem is solved where ρ is material density, u i is displacement com-ponent, σ ij is stress component and, b i is body force component [17]. Einstein summation convention is used and comma in index denotes spatial derivative. The constitutive equation in the variational formulation is where ε ij is strain and t i is traction component on boundary. In our case the region Ω is a 2D space and boundary Γ surrounding it is a 1D line. The body forces are zero in this simulation. Strain is assumed to be small. The matrix formulation of the finite element model with damping is where M is mass matrix, K is stiffness matrix, F is external forcing and ∆ is displacement vector [15].
The damping matrix C is used to apply the Lysmer-Kuhlemeyer absorbing boundary conditions [16] as a diagonal matrix, allowing to take advantage of the explicit solution scheme. The element matrices are where C e is here the constitutive matrix for the plane strain elasticity. Linear triangular three-node elements (T3), also known as constant strain triangles [17], were chosen for this problem for the following reasons. Firstly because the epoxy layers in the laminate model can be very small, therefore small elements are needed anyway, with T3 being computationally cheapest. Secondly, linear elements are well suited for nonlinear problems: since the strain is constant throughout the element, the computation of nonlinear constitutive relations would also be simple. In this simulation, the material itself is linear but future work might include nonlinearity or hysteresis.
The T3 element lumped mass matrix [18] is where I 6 is 6 × 6 identity matrix and A e is the area of the element. The element stiffness matrix is where matrix B is and with x i and y i being the node coordinates [17], then The external distributed force is simply divided into relevant nodes. The Lysmer-Kuhlemeyer absorbing boundary conditions are applied as viscous stresses on the boundaries, which means that they can be applied on DOF basis, making the damping matrix C diagonal. The viscous stresses on the boundary DOFs are where Γ is the boundary portion of the element [16].
In this work the scaling parameters are a = 1 and b = 1. The wave velocities used for these boundary conditions are V p = 2972 m/s and V s = 1956 m/s [19]. Equation (9) is solved for each timestep ∆t = 5 · 10 −10 s by explicit central difference scheme This scheme is solved by dividing the equation by the term in the first parentheses, which is simple if M and C are diagonal matrices. Each simulation considers a 60 µs time window.

Contact gap treatment
There is a single source of nonlinearity in this simulation: the contacting crack fully inside the material (Fig. 4). If the material is at rest, then the crack is small and straight. In this work, there is neither a preload nor an initial gap in the contacting crack. This simple material defect results in a localised nonclassical nonlinearity, which can be analysed by various signal processing methods. It is known that frictional contact problems can be sensitive to timestep length and loading path [20]. In this work, it is assumed that the small timestep length and relatively small forces involved keep the error small. Therefore an explicit solution method scheme is utilized, similarly to [21]. A more precise solution could be expected from an implicit scheme, but that is left for the future. Further refinements could include thermoelastic contribution to the constitutive equation at the frictional contact gap [3].
The node-to-node contact model is used [22] with Coulomb friction. The defect is horizontal, simplifying the calculation of normal gap between the nodes. If the position of a node on a slave (lower) surface is (n s x , n s y ) and on master (higher) (n m x , n m y ), then the normal contact gap is g N = n s y − n m y and the tangential gap (offset) is g T = n s x − n m x . In case of normal penetration of one surface into another, then g N > 0. If there is no penetration, then g N ≤ 0. The coefficient of friction is µ = 0.6, and the solution aims to satisfy the Kuhn-Tucker conditions on the crack surface: where λ N is the normal force on crack, σ is stress and n is the normal vector of the surface. The penalty plus Lagrange multiplier method is used for normal contact and the penalty method for friction [23]. The contact logic for the node pairs can be summarized by following steps.
• Vector gap functions are found: g N = n s y − n m y and g T = n s x − n m x .
• Normal forcing is updated λ N = λ N + g N b where b is some big penalty value and λ N ≥ 0.
• Logic diverges to 3 paths: No force is applied in case of no contact.
Only normal force is applied if preceding step had no contact or had contact with tangential slip.
Normal and tangential forces are applied if previous iteration had non-slip contact.
• The normal contact condition is verified by setting the penetration value g P = g N where g N ≥ 0. Then the L 2 -norm of penetration is evaluated g P |g P < ε where ε is the limiting value for the error due to contact penetration. If the condition is not fulfilled, the iteration is repeated, otherwise new timestep is taken.
A more thorough explanation of this contact gap logic is available at [15].

Results
The signal analysis of the time domain simulation results of the damaged and undamaged medium are compared, describing some analysis measures which could allow to detect the presence of damage as nonlinearity. The simulation follows ultrasonic TR-NEWS NDT procedures where the transducer data is available as time-series, measured at some specific location. The signals are low-pass filtered to keep only the ultrasonic component. Here five measurement points are analysed at various distances from the crack damage and transmitting transducer (Fig. 4). A video of the displacement fields for TR-NEWS focusing to point 3 in cracked medium is available at [24]. Figure 8 shows the undamaged CFRP TR-NEWS focusing for the receiver positions 1 to 5 (Fig. 4). It is an ordinary TR-NEWS focusing where at the middle of the signal (30 µs) is the focusing, surrounded by the sidelobes. There are two aspects to note about this is figure Figure 9 shows the TR-NEWS results of the cracked CFRP test object simulation for the receiving transducer positions 1 to 5 (Fig. 4). Here the PI signal processing is also applied and it shows the nonlinearity as difference between results from initial chirp signals with positive and negative sign. This nonlinear, damaged case exhibits nonlinearity particularly strongly in receiving position 3 (near the middle of the crack). Also, the sidelobes are unsymmetrical in respect to the main lobe. Figure 10 shows the envelopes of the PI measure of nonlinearity across the measuring points. The non-  2. Comparing the amplitudes of the positions 2 and 4, at far and near side of the crack end respective to transmitter: the farther position has larger focusing amplitude than the nearer position. Since the simulation region has two absorbing boundaries, the wave propagation is mostly in one direction, therefore the defect between pos. 2 and 4 must be capturing the wave energy and the TR-NEWS signal processing is using that energy as a new "virtual transducer" for the pos. 2 focusing. This could be further analysed in future works from the correlation signals which generate these focused signals.

TR-NEWS with pulse inversion
3. Amplitudes from the measurement positions 1 and 5 are "right way" around: the nearer measurement point has larger focusing amplitude than the farther. Figure 12 shows a snapshot of the simulation u 2 displacement at a time moment t = 32.6 µs, just after the focusing. The defect in material is acting as a source of new excitation after TR-NEWS focusing. Wave energy is captured between the damage and outside wall of the material and emitted as a wave.

Delayed TR-NEWS analysis
Section 2.2 describes the delayed TR-NEWS signal processing method which allows to create arbitrary envelope wave at the focusing (Eq. (5)), instead of the simple peak of the TR-NEWS. Equation (6) shows that in linear material, the outcome of the delayed TR-NEWS process can be predicted. Since this method with prediction works very well in physical NDT measurements of linear materials [14], it is now tested in simulation with the nonlinearity, supposing Figure 12: (Colour online) Displacement u 2 field at time t = 32.6 µs with a wave emission coming from the damaged region. Video available at [24] that the difference between the simulation result and the linear prediction (Eq. (6)) is due to nonlinear interaction of waves in the presence of nonlinearities or damage. Figure 13 shows the comparison between the linear superposition prediction and the simulation result of a simple delayed TR-NEWS process where two focusing peaks are at superposition with 1 µs time delay. The difference between the prediction and the simulation is large and obvious, indicating the presence of nonlinearity. This measure of nonlinearity seems to be stronger than the measure calculated from PI ( Fig. 9), making it a good candidate for further investigation. The delayed TR-NEWS signal processing could also be used for activating the contacting gap as the energy pocket. This could be done by creating a new focusing wave envelope which would have the resonant frequency of the defect, permitting higher amplitude waves near the damaged region, which would enhance the extraction of the nonlinear signature. This study is left for the future.

Conclusion
This paper investigated nonlinear NDT by using a simple FEM simulation model for a crack nonlinearity in CFRP. In the laminate model, the damage is a simple horizontal contacting crack near the receiving transducer. The signal processing uses TR-NEWS method for focusing the available wave energy near the receiving transducer. The magnitude of nonlinearity due to the damage is measured firstly with PI, secondly with the proposed delayed TR-NEWS signal processing procedure. While PI indicates the presence of the nonlinearity, the simple delayed TR-NEWS procedure shows it even more strongly and is promising for future investigations and further development due to its signal processing flexibility.
Since the delayed TR-NEWS procedure allows to generate a wave at the focusing with arbitrary envelope, it could be used in the future to excite the crack damage by its resonance frequencies, using the damage as an energy pocket. Other perspectives include a more detailed simulation model for the CFRP in order to take more of its microstructure geometry into account to have stronger focusing. Additionally, the damage could be modelled either by a collection of various cracks at various angles or by hysteresis. Moreover, heating from the frictional forces at the damage could be considered for a more precise simulation model.