ultrasonics for the detection of closed delaminations in composites

Accepted Manuscript Finite element modelling of nonlinear acoustics/ultrasonics for the detection of closed delaminations in composites Ashish Kumar Singh, Bo-Yang Chen, Vincent B.C. Tan, Tong-Earn Tay, HeowPueh Lee PII: DOI: Reference:

S0041-624X(16)30199-8 http://dx.doi.org/10.1016/j.ultras.2016.09.019 ULTRAS 5383

To appear in:

Ultrasonics

Received Date: Revised Date: Accepted Date:

16 June 2016 3 August 2016 25 September 2016

Please cite this article as: A.K. Singh, B-Y. Chen, V.B.C. Tan, T-E. Tay, H-P. Lee, Finite element modelling of nonlinear acoustics/ultrasonics for the detection of closed delaminations in composites, Ultrasonics (2016), doi: http://dx.doi.org/10.1016/j.ultras.2016.09.019

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Finite element modelling of nonlinear acoustics/ultrasonics for the detection of closed delaminations in composites Ashish Kumar Singha , Bo-Yang Chena , Vincent B. C. Tana,1,∗, Tong-Earn Taya , Heow-Pueh Leea a

Department of Mechanical Engineering, National University of Singapore. 9 Engineering Drive 1, Singapore 117575

Abstract Linear ultrasonics methods based on the principle of reflection, transmission, dissipation of sound waves have been traditionally used to detect delaminations in composite structures. However, when the delamination is in very early stages such that it is almost closed, or closed due to a compressive load, the linear methods may fail to detect such cases of delaminations. Nonlinear acoustics/ultrasonics have shown potential to identify damages in composite structures which are difficult to detect using conventional linear ultrasonic methods. The nonlinear method involves exciting the structure with a sinusoidal signal of certain (or multiple) frequency and observing the vibrations of the structure. The vibrations of the damage region differ significantly from intact regions and can be used to identify the damage. However due to the complex and varying nature of the nonlinear phenomena created by the interaction between the exciting signal and the damage, there are ∗

Corresponding author Email address: [email protected] (Vincent B. C. Tan)

Preprint submitted to Ultrasonics

September 26, 2016

many variables at play which can lead to success or failure of the method. While experiments lead to the establishment of the method to be used as a damage detection technique, numerical simulations can help to explain the various phenomena associated with nonlinearity. This work presents a quick approach to model the nonlinear behavior caused by closed delaminations. The model is validated with a previously available approach for nonlinear vibrations modelling and a comparison is made between the two. The local nature of the nonlinearity enables to map out the area of damage in the structure. Additionally, a few parametric studies are performed to study the effect of various parameters related to the nonlinear phenomenon. Keywords: nonlinear acoustics/ultrasonics, modelling and simulation, contact acoustic nonlinearity, closed delaminations, kissing bonds 1. Introduction Composite materials are being increasingly used in aerospace industry due to their high strength to weight ratios leading to highly efficient structures. However, due to layer by layer composition of composite materials, internal defects like delaminations can easily develop, which can reduce the strength of the structure severely. The delaminations in composite structures can be identified with the help of linear ultrasonic methods which employ the principles of reflection, transmission, and absorption of sound waves to locate the delamination. However, the linear ultrasonic techniques may fail to detect defects in certain cases like when the delamination is in very early stages or it is closed due to a compressive load. Such type of defects form what is known as kissing bond defects and are encountered in other type of structures as 2

well. For example, kissing bonds are known to occur in adhesive joints and have been studied extensively [1–6]. For the detection of kissing bond defects, nonlinear acoustics/ultrasonics methods have shown great potential [7–15]. Nonlinear methods involve exciting the structure with a signal of certain frequency or multiple frequencies and observing the output response of the structure. The region which has the defect vibrates nonlinearly and in case of a single frequency excitation input signal, the output frequency spectrum consists of the exciting frequency (linear response) along with harmonics and subharmonics of the exciting frequency (nonlinear response). In case of a delamination with cleanly separated surfaces the nonlinearity is due to the opening and closing of the delamination gap commonly known as contact acoustic nonlinearity or breathing of the crack. An introduction to the nonlinear acoustics for damage detection can be found in [16–18]. A vast amount of experimental studies have shown the potential of the nonlinear acoustic methods to detect different types of defects in structures. However, there have been few studies to explain the complex nonlinear phenomenon exhibited in experiments, and hence a need arises for more numerical simulations. Numerical simulations provide a method to explain the nonlinear behavior and provide detailed insight which can help to optimize the technique further. In case of a delamination with cleanly separated surfaces, the nonlinear behavior can be modellled by employing nonlinear spring dampers between the delamination surfaces [19–21] or using contact between the delamination surfaces [22–29]. The nonlinear spring damper models can be tedious to implement and require a number of parameters related to the spring stiffness and damping coefficients, the estimates of which 3

Input point

Output point Delamination

0.3 mm

20 mm 40 mm

40 mm

Figure 1: Geometry of the model. Left: Isometric view; Right: Cross-sectional view

are not available readily. A recent review of various nonlinear phenomenon and their impact on the damage detection process can be found in [30, 31]. This paper discusses the approach to model the nonlinear behavior by employing the contact boundary conditions between the delamination surfaces. The model is validated with the nonlinear spring damper model presented in [20]. A critical comparison between the two models is also discussed. The local nature of nonlinear behavior has been used to map out the area of damage. Parametric studies which discuss the effect of delamination diameter and delamination depth were also performed. 2. Model description 2.1. Geometry and material properties The model consists of a square plate with 40 mm x 40 mm dimensions and a thickness of 2 mm [20]. The plate consists of a closed circular delamination of diameter 20 mm at a depth of 0.3 mm from the top surface as can 4

Table 1: Material properties.

E1 (GPa)

161.0

E2 = E3 (GPa)

11.38

G12 = G13 (GPa)

5.17

G23 (GPa)

3.98

ν12 = ν13

0.32

ν23

0.436

ρ (kg/m3 )

1800

be seen in cross-sectional view of Figure 1. The material is a unidirectional carbon fibre reinforced composite the properties of which are listed in Table 1 [20]. To implement material damping in the transient model, Rayleigh damping has been used. It uses the mass damping parameter (α) and stiffness damping parameter (β) to characterise the damping behaviour of the material. The Rayleigh damping parameters are related to the damping ratio (ξ) and the circular excitation frequency (ω) by Equation 1. In the current model the damping ratio ξ value has been taken as 0.1. Equal contributions are considered to the damping ratio from the stiffness damping (β) and the mass damping(α).

ξ=

α βω + 2ω 2

(1)

2.2. Boundary conditions A sinusoidal excitation of frequency 26165 Hz (f0 ) with an amplitude of 1 µm is provided to the plate at the input point (Figure 1). The chosen 5

frequency is one of the natural frequencies of an intact plate with same dimensions and material properties but without the delamination defect. This particular frequency was chosen as it had a mode shape which facilitates the opening and closing of the delamination if there were to be a defect at the centre. All the edges of the plate are kept free. Contact boundary conditions [32, 33] have been defined between the two surfaces of the delamination. Contact simulations are highly nonlinear and require robust solvers to avoid convergence problems. Therefore, it was decided to use the commercial finite element package Abaqus/Standard [34] to build the 3D dynamic implicit model. It is a well established software with significantly advanced stabilization controls for contact modelling. In contact modelling it is recommended to select the surface with higher stiffness as the master surface and the surface with lower stiffness as the slave surface. Since the upper part of the delamination is thinner than the lower part of the delamination, it is less stiff and hence, the upper surface of delamination has been defined as the slave surface. Similarly, the bottom part of the delamination is thicker than the upper part of the delamination and therefore the bottom surface of the delamination has been defined as the master surface. Any standard contact formulation in a numerical framework should address the following aspects; constraints formation, pressure-overclosure relationship, constraints enforcement, and constraints evolution upon sliding. Abaqus provides a number of possible options for these aspects. Therefore in order to have a well defined contact model an appropriate selection should be made for each of the aspects. The upcoming discussion would be about 6

the possible options available in the these aspects and the selection criteria for the chosen particular option. The details presented in this discussion are from the documentation of Abaqus [34]. Also, the focus will be on the implication of the selection rather than the mathematical implementation of the selection. For details about the mathematical implementation Abaqus documentation [34] can be referred to. For modelling of contact between two surfaces, two constraint formation schemes are available namely, node to surface discretisation and surface to surface discretisation. Node to surface discretisation involves node against a surface while surface to surface discretisation involves a constraint which is based on an integral over the region surrounding a slave node. In this simulation a surface to surface contact discretisation has been used since it is more accurate than the node to surface discretisation method. Also, the surface to surface discretisation is significantly less dependent on the selection of master and slave surfaces. The pressure-overclosure relationship can be defined as the hard contact or softened contact. The current model uses the hard contact conditions as it minimizes the penetration of the slave surface into the master surface. Also, it is known that for implicit models involving collisions, the soft contact modelling approach is computationally expensive in comparison to hard contact modelling approach. The contact pressure (p) and overclosure (h) are related by Equation 2(also see Figure 2a) in hard contact relationship. It is to be noted that a negative overclosure is a positive clearance.

p = 0 f or h < 0 (open), and h = 0 f or p > 0 (closed) 7

(2)

(a) Ideal presuure-overclosure relation- (b) Linear and nonlinear penalty method ship in hard contact modelling

for constraint enforcement

Figure 2: Hard contact modelling

The hard contact conditions have the following implications: 1) The contact pressure is zero unless the master and slave surfaces are in contact with each other. 2) Penetration is not allowed between the master and slave surfaces. However, depending on the enforcement method this condition will either be strictly satisfied or approximated. 3) When the surfaces are in contact with each other there is no limit on the value of contact pressure. To enforce the constraints in hard contact modelling, Abaqus provides three different enforcement methods, namely the direct method; which strictly enforces a given pressure-overclosure behaviour per constraint, the penalty method; which approximates pressure-overclosure behaviour, and the augmented Lagrange method; which uses same kind of approximation as the penalty method, but also uses augmentation iterations to improve accuracy. In the current model a penalty method has been used for the enforcement of contact constraints. In the penalty method, the contact force is a function of

8

the penetration distance and some degree of penetration is involved (usually insignificant). The relationship between the contact force and the penetration distance can be defined to be linear or nonlinear (see Figure 2b). The current model assumes a linear relationship between the contact pressure and overclosure with a contact stiffness of ten times the material stiffness. The penalty method was chosen because it offers numerical softening which is beneficial for the computational efficiency as it provides improved convergence. The effect of other constraint enforcement methods (direct enforcement and augmented Lagrange) in hard contact modelling and nonlinear penalty stiffness on the results will be discussed in section 3 at a later stage. To model the evolution of constraints upon sliding two options are available i.e. small sliding and finite sliding. As the name suggests small sliding assumes that the relative sliding between the two surfaces in contact is smaller. The current model uses the finite sliding scheme which provides for any arbitrary motion between the contacting surfaces. The mathematical details of the two schemes can be found in [34]. Additionally, the tangential behaviour is modelled as frictionless. The tree below summarises the various options available for contact modelling in Abaqus/Standard.

Contact modelling Constraint formation Node to surface discretisation Surface to surface discretisation Pressure-overclosure relationship Hard contact Direct enforcement Penalty method enforcement 9

Linear Non-linear Augmented Lagrange enforcement Softened contact (direct enforcement) Exponential relationship Linear relationship Tabular piecewise-linear relationship Constraint evolution upon sliding Small sliding Finite sliding

The current model assumes that the delamination surfaces are smooth and there exists no bonding force as well as no friction between the two contacting surfaces. However, in experiments the type of nonlinearity may not be known beforehand and there may even exist a combination of different types of nonlinearities. For example in a defective adhesive joint some part of the damaged surfaces maybe be completely separated while some part may still be weakly bonded. To model this type of defects additional options need to be incorporated into the contact model. To facilitate such modelling, the contact regions can have additional options to include springs, dampers, friction, and cohesive behaviour as well. However, these options are not implemented in the current model as the current defect is considered to be smooth cleanly separated delamination surfaces. 2.3. Meshing details The mesh used for the model is shown in Figure 3. A combination of quadratic tetrahedral and wedge elements was used for the analysis. The mesh resolution has been kept higher near the defect region so that accurate 10

Figure 3: Mesh used for the model. Left: Isometric view; Middle: Top view; Right: Bottom view

modelling of the contact is possible. Also, as a general rule of thumb for contact modelling the size of mesh elements on the slave surface (delamination upper surface) should be twice (or more) the size of mesh elements on the master surface (delamination lower surface). The remaining mesh is still refined enough such that in one wavelength of the excitation frequency there exist enough number of mesh elements for accurate wave propagation modelling. 3. Results and discussions The results obtained from the contact model are presented in this section. It is to be noted that similar results have already been reported using the nonlinear spring damper model in [20]. Later in subsection 3.3 the results between the two models will be compared and it will also be validated that the results produced by the new contact model are similar to the results obtained using the nonlinear spring damper model. 11

T/6

2T/6

3T/6

4T/6

5T/6

6T/6

7T/6

8T/6

9T/6

10T/6

11T/6

12T/6

Figure 4: Opening and closing of the delamination at the cross-sections during two time cycles. The deformations are scaled by a factor of 10000 for better visualization.

Figure 4 shows the opening and closing of the delamination at two perpendicular cross sectional vertical planes (x-y plane and y-z plane, Figure 1) of the plate for two cycles of the excitation frequency. The deformations have been scaled by a factor of 10000 for better visualization. It can be observed that the part of the plate above the delamination acts like a membrane and the vibration pattern is quite different from that of the global plate. The upper and lower portions of the delamination come in contact with each other a couple of times during the two cycles. The delamination region near the output point (Figure 1) experiences an impact around time 4T /6 and around time 11T /6. However, the impact around time 4T /6 is more severe and leads to a large opening of the delamination region. Similar to a study described in [20], the y-displacement time history at five points is extracted and is shown in Figure 5. The solid blue line is

12

Normalised FFT amplitude

Normalised y-displacement

Measurement locations

x z

t/T

f/f0

Figure 5: Left: Measurement location of points; Middle: y-displacement time history at the measured points on the upper surface (solid blue line) and lower surface (dashed red line) of the delamination; Right: Frequency domain results of the measurement points on the upper surface of the delamination

13

for y-displacement at a point just above the delamination and the dashed red line is for a point just below the delamination. The displacement has been normalised by the peak to peak amplitude of the point just above the delamination and in the centre of it (this point experiences the maximum vibration amplitude). From the figure it can be observed that the points above the delamination vibrate radically different from the points below the delamination. The points below the delamination have more or less sinusoidal profile at each of the points. This is due to the fact that the lower part having more inertia due to it’s mass does not change it’s direction of movement upon the impact with the upper part of the delamination. To find out the frequency content of the vibrations of the upper surface of the delamination, the displacement time history is converted to frequency domain using the Fast Fourier transform (FFT). In the frequency spectrum, the presence of harmonics of the excitation frequency (2f0 , 3f0 , 4f0 ...) can easily be identified. Also present is the subharmonic of the excitation frequency (f0 /2) and it’s higher harmonics (3f0 /2, 5f0 /2, ..). The presence of harmonics and subharmonics is consistent with the observations reported in a number of experimental and numerical works [7–15, 19–21, 35–40]. The amplitude of the subharmonic (f0 /2) is highly dependent on the position of the point. At some points the amplitude of the subharmonic (f0 /2) even exceeds the amplitude of the linear frequency. This nonlinear frequency content (harmonics and subharmonics) is a result of the nonlinear interaction of the contact surfaces. The hard contact model with constraint enforcement through penalty method (linear behaviour) was compared with penalty method (nonlinear 14

Figure 6: y-displacement time history at the output point (Figure 1) for different constraint enforcement methods of hard contact modelling

behaviour), augmented Lagrange method, and direct stiffness method shown in Figure 6. The results are shown at the output point (Figure 1). The pressure overclosure relationship for hard contact modelling with constraint enforcement using linear or nonlinear penalty method has been shown in Figure 2b. It can be seen that the results from all the constraint methods were found equivalent to each other. Therefore, it can be said that when hard contact modelling approach has been used, any of the constraint enforcement method can be used without any considerable effect on the output results. 3.1. Cause of the harmonics and subharmonics It is known that if a sinusoidal signal has a variable amplitude, it’s frequency spectrum consists of the harmonics of the frequency of the sinusoidal signal. The same effect is happening in the current model as well. Due to the sinusoidal loading, the delamination will repetitively experience the tension 15

and compression cycles. When the delamination surfaces are in tension, it causes the delamination to open, changing the amplitude of the signal. For the compression part, the amplitude does not change because the delamination surfaces get in touch with each other behaving like an intact surface. This difference of the amplitude between the compression and tension cycles lead to the presence of higher harmonics in the frequency spectrum of points near the delamination region. The two collisions between the upper and lower surfaces of the delamination around time instants 4T /6 and 11T /6 can be observed in Figure 4 and Figure 5. The collision which occurs at 4T /6 time causes the upper part of the delamination to effectively change it’s direction of motion and there is a large opening of the delamination gap after the impact event. The lower part of the delamination being thick does not deviate much from its original path due to this collision. However, the upper part of the delamination being thin changes the direction of its original path. The larger is the impact of this collision more are the chances of a change in the direction of velocity. This opening of the delamination gap has an effect on the time period of the vibrations (and thus the frequency) of the upper part of the delamination. It can be seen from Figure 5 that at output points 1, 4, and 5 the time period of the y-displacement time history of points on the upper part is roughly double the time period of the y-displacement time history of points on the lower part. The doubling of the time period effectively leads to a significant presence of subharmonic (f0 /2) in the frequency spectrum. However, it is to be noted that if the impact energy is not enough to cause this change in the direction of the upper half of the delamination, the sub16

harmonics are absent from the frequency spectrum. This happens when the delamination is deep inside the surface of the plate or the delamination size is small in comparison to the size of the plate. The effect of the delamination depth and the size of the delamination will be discussed in subsection 3.4. 3.2. Harmonic imaging of the damage In order to obtain a global map of the damage area the output y-displacement was extracted at all the nodes on the top surface of the plate. The nonlinear behavior observed is local in nature and therefore the area of the damage can be mapped out by measuring the amplitude of the nonlinear frequencies (harmonics and subharmonics) in the frequency spectrum of all the nodes. The ratio of the amplitude of the second harmonic (2f0 ) to the amplitude of the linear frequency (f0 ) is plotted in Figure 7a. It can be seen that the ratio is more pronounced near the defect region giving a qualitative idea about the location of the damage. Similarly, the ratio of the amplitude of subharmonic (f0 /2), to the amplitude of the excitation frequency (f0 ) is also shown in Figure 7b. The ratio is larger near the defect region. The nonlinear harmonic imaging technique for the detection of defects has been applied experimentally in a number of works [35–40]. However, in both the plots there are some points which are not above the delamination but still have a larger magnitude of the ratio, giving a false indication of the damage in that region. This happens because these points lie on the modal pattern where the vibration amplitudes are very small. This is one of the limitations of the nonlinear imaging methods for damage detection. Since, different frequencies will create different modal patterns the position of these points will change if some other frequency is 17

(a)

(b)

Figure 7: (a) Ratio of the second harmonic (2f0 ) amplitude to excitation frequency (f0 ) amplitude (b) Ratio of the subharmonic (f0 /2) amplitude to excitation frequency (f0 ) amplitude

used. Therefore, a possible alternative is to conduct the investigations at more than one frequency. For example in [20] sweep excitation signals have been studied to overcome the dependency on excitation frequency. 3.3. Model validation The nonlinear interaction of the delamination surfaces can also be modelled by assuming nonlinear spring dampers between the contacting surfaces. The current model using contact conditions have been compared to the earlier nonlinear spring damper model presented in [20]. The details of the nonlinearity modelling using springs and dampers from [20] are also desribed below. The interaction between the delamination surfaces have been defined by forces on the upper and lower surfaces of the delamination. The spring force 18

Figure 8: Illustration of the virtual spring damper elements introduced for the implementation of the delamination boundary conditions [20]

. per unit area on the top surface (Fst ) has been defined as:

F st =

         

k1 (Z0 − ∆z)

if

∆z < Z0

k2 (Z0 − ∆z)

if

Z0 < ∆z < aZ0

        

k3 (bZ0 − ∆z)

if

aZ0 < ∆z < bZ0

if

∆z > bZ0

0

(3)

The separation distance between the interfaces (∆z) is the difference between the displacement of the top (uzt ) and bottom (uzb ) surfaces.

∆z = uzt − uzb

(4)

The value of the parameters involved in Equation 3 is given in Table 2. The parameter k3 is defined as:

k3 = k2

1−a b−a

19

(5)

Table 2: Spring force parameters.

k1 (N/m3 )

5 · 1012

k2 (N/m3 )

3 · 1012

Z0 (nm)

0.1

a

4

b

5

The damping force per unit area on the top surface (Fdt ) is a function of the velocity of the top surface (vzt ) and the bottom surface (vzb ) and is defined in Equation 6. The damping coefficient(γ) has a value of 106 N s/m3 .

Fdt =

  

−γ(vzt − vzb )

if

∆z < bZ0

 

0

if

∆z > bZ0

(6)

The spring and damping force per unit area on the lower surface (Fsb and Fdb respectively) have the same magnitude as that of the spring and damping force per unit area on the top surface (Fst and Fdt respectively) but act in the opposite directions. Implementing the above procedure which involves a force dependent on the opening and closing of the delamination would require a new user element to be defined in Abaqus. Due to the complexity involved with the implementation of a UEL subroutine it was instead decided to use the Comsol Multiphysics software [41] itself to build the model as was used in [20]. To implement forces dependent on the deformations, Comsol does not require any external subroutine and the implementation is quite straightforward. The above nonlinear spring damper model has been implemented in a num20

ber of works [42–45]. An equivalent of the contact model was built with same geometry, mesh, and boundary conditions. In place of defining contact between the delamination surfaces, the spring damper elements are implemented in Comsol by defining four forces, Fst ; Fdt ; Fsb and Fdb respectively acting on both the top (subscript t) and bottom (subscript b) interface of the delamination as shown in Figure 8. The contact model and nonlinear spring damper model were compared by plotting the time displacement history at the measurement points (Figure 5) on the top surface of the delamination as shown in Figure 9. The solid blue line is for the results obtained from the contact model while the red dashed line is for the results obtained from the nonlinear spring damper model. The results from both the models are similar to each other along with a minor differences which are to be expected since the two methods are not exactly equivalent to each other. Although, the vibration patterns at the points more or less resembles between the two models, it can be seen that the spring damper model overestimates the opening of the delamination upon the impact when compared to the contact model. A major limitation of the nonlinear spring damper model is that the parameters for spring stiffness (Table 2) and damping coefficients are not available readily. This set of parameters are suited to the current simulation. However if there is a major change in the size of the specimen or the material properties of the specimen, the parameters need to be altered. This is not a problem with contact modelling since the required parameters are established internally by the commercial software based on then stiffness of 21

Normalised FFT amplitude

Normalised y-displacement

Measurement locations

x z

t/T

f/f0

Figure 9: Left: Measurement location of points; Middle: y-displacement time history at the measured points on the upper surface of the delamination for contact model (solid blue line) and spring damper model (dashed red line); Right: Frequency domain results of the measurement points on the upper surface of the delamination for contact model (solid blue line) and spring damper model (dashed red line)

22

the underlying material. Usually, for hard contact modelling with a linear penalty approach a contact stiffness which is ten times the stiffness of the material can be used. A very small contact stiffness increases the interpenetration between the surfaces and a very large contact stiffness may lead to convergence issues. Both of these conditions are undesirable. Another advantage of the contact model is that it is much more computationally efficient than the spring damper model. But this claim is difficult to prove since the two models were solved in two different commercial softwares and the difference between the solvers of both models need to be considered. For a realistic comparison of time efficiency the two models need to be implemented in a single package or platform. The contact model was solved using Abaqus while the spring damper model was solved using Comsol Multiphysics software. The opposite combination is not suitable since time dependent contact formulations are tedious to implement in Comsol, and deformation dependent force formulations are tedious to implement in Abaqus. 3.4. Parametric studies This section presents the effect of the delamination depth and the diameter of the delamination on the response of the plate. It is to be noted that the parametric studies are now carried out using the contact model only. The analysis was performed at two more depths of the delamination and the results at the output point (see Figure 1) are shown in Figure 10. The time domain results are shown for four cycles of the excitation frequency. While the harmonics of the excitation frequency are present for all three values of the depth, the subharmonics of the excitation frequency are absent 23

(a) Time domain

(b) Frequency domain

Figure 10: Effect of change in the depth of the delamination.

(a) Time domain

(b) Frequency domain

Figure 11: Effect of change in the diameter of the delamination.

for delamination depth of 0.4 mm. It can be seen that when the delamination distance from the surface increases the subharmonics (f0 /2, 3f0 /2, ...) disappear from the frequency spectrum. This is because for larger delamination depth, the impact of upper surface with the lower surface does not cause a sudden change in the direction of the movement of the upper part, resulting in the loss of the subharmonic component. The change in amplitude of the sinusoidal signal is still present leading to harmonics in the frequency spectrum. The phenomena of disappearance of subharmonics with the delamination depth was also reported in [20]. The simulation was carried out for two more values of the delamination

24

diameter and the results are shown in Figure 11. The effect of an increase in the diameter of the delamination has a similar effect to a decrease of the delamination depth. When the diameter of the delamination decreases the impact between the lower and upper surface does not cause the movement of the upper surface to change its direction resulting in the loss of the subharmonic component. 4. Conclusions and summary A 3D finite element model was developed to study the nonlinear vibrations created by the opening and closing of a closed delamination. The current nonlinear acoustic method involves exciting the structure at a certain location with a signal of certain frequency (f0 ) and extract the time displacement history at other locations of the structure. The time displacement history when converted to frequency domain shows the presence of nonlinear frequencies in the form of harmonics (2f0 , 3f0 , ..) and subharmonics (f0 /2, 3f0 /2, 5f0 /2, ..). It was found out that the presence of the harmonics is due to the opening of the delamination when the surfaces experience a tension force. The cause of subharmonics was established as the impact between the upper and lower surfaces of the delamination which causes a sudden change in the direction of motion of the upper part of the delamination leading to a large opening of the delamination gap and effectively changing the time period of vibrations (and frequency) of the upper part of the delamination. The current model uses contact between the surfaces of the delamination to model the nonlinear behavior. This nonlinearity can also be modelled by inserting nonlinear spring and damper elements between the contact surfaces 25

[20]. A model was built with this approach and the results of the two models were compared. It was found out that the spring damper model shows more opening of the delamination after the impact in comparison to the contact model. Overall, the results were similar to each other for both the models. The presence of harmonics and subharmonics in frequency spectrum signifies damage in the structure. The nonlinearity due to the defect is local in nature and is larger in the vicinity of the crack. This property can be utilized to map out the area of damage in a structure by measuring the output displacement at a number of points. The ratio of the amplitude of the subharmonic (f0 /2) and first harmonic (2f0 ) to the amplitude of excitation frequency (f0 ) was plotted on the upper surface of the 3d plate. The ratio near the defect region was found to be significantly higher. Parametric studies were performed to determine the effect of delamination diameter and delamination depth on the response of the structure. It was noted that as the time displacement of the delamination decreases, or as the depth of delamination increases, the impact between the top and bottom surfaces of the delamination does not change the direction of movement of the upper part of the delamination and the subharmonic component of the frequency spectrum disappears. 5. Acknowledgement The support in the form of research scholarship for the first author and resources from grant no. R-265-000-463-112 of NUS are gratefully acknowledged.

26

6. References [1] P. B. Nagy, “Ultrasonic-detection of kissing bonds at adhesive interfaces,” Journal of Adhesion Science and Technology, vol. 5, no. 8, pp. 619–630, 1991. [2] P. B. Nagy, “Ultrasonic classification of imperfect interfaces,” Journal of Nondestructive evaluation, vol. 11, no. 3-4, pp. 127–139, 1992. [3] T. Kundu, A. Maji, T. Ghosh, and K. Maslov, “Detection of kissing bonds by Lamb waves,” Ultrasonics, vol. 35, no. 8, pp. 573–580, 1998. [4] C. J. Brotherhood, B. W. Drinkwater, and F. J. Guild, “The effect of compressive loading on the ultrasonic detectability of kissing bonds in adhesive joints,” Journal of Nondestructive Evaluation, vol. 21, no. 3, pp. 95–104, 2002. [5] C. J. Brotherhood, B. W. Drinkwater, and S. Dixon, “The detectability of kissing bonds in adhesive joints using ultrasonic techniques,” Ultrasonics, vol. 41, no. 7, pp. 521–529, 2003. [6] S. L. Poveromo and J. C. Earthman, “Analysis of “Kiss” Bonds Between Composite Laminates,” JOM, vol. 66, no. 6, pp. 970–978, 2014. [7] P. B. Nagy, P. McGowan, and L. Adler, “Acoustic nonlinearities in adhesive joints,” in Review of progress in quantitative nondestructive evaluation, pp. 1685–1692, Springer, 1990. [8] M. Rothenfusser, M. Mayr, and J. Baumann, “Acoustic nonlinearities in adhesive joints,” Ultrasonics, vol. 38, no. 1, pp. 322–326, 2000. 27

[9] D. Donskoy, A. Sutin, and A. Ekimov, “Nonlinear acoustic interaction on contact interfaces and its use for nondestructive testing,” NDT & E International, vol. 34, no. 4, pp. 231–238, 2001. [10] M. Meo and G. Zumpano, “Nonlinear elastic wave spectroscopy identification of impact damage on a sandwich plate,” Composite Structures, vol. 71, no. 3, pp. 469–474, 2005. [11] Y. Ohara, T. Mihara, and K. Yamanaka, “Effect of adhesion force between crack planes on subharmonic and DC responses in nonlinear ultrasound,” Ultrasonics, vol. 44, no. 2, pp. 194–199, 2006. [12] I. Solodov and G. Busse, “Nonlinear air-coupled emission: The signature to reveal and image microdamage in solid materials,” Applied Physics Letters, vol. 91, no. 25, p. 251910, 2007. [13] D. Yan, B. W. Drinkwater, and S. A. Neild, “Measurement of the ultrasonic nonlinearity of kissing bonds in adhesive joints ,” NDT & E International, vol. 42, no. 5, pp. 459 – 466, 2009. [14] D. Yan, S. A. Neild, and B. W. Drinkwater, “Modelling and measurement of the nonlinear behaviour of kissing bonds in adhesive joints,” NDT & E International, vol. 47, pp. 18–25, 2012. [15] B. Y. Chen, S. K. Soh, H. P. Lee, T. E. Tay, and V. B. C. Tan, “A vibro-acoustic modulation method for the detection of delamination and kissing bond in composites,” Journal of Composite Materials, p. 0021998315615652, 2015.

28

[16] I. Y. Solodov, N. Krohn, and G. Busse, “CAN: an example of nonclassical acoustic nonlinearity in solids,” Ultrasonics, vol. 40, no. 1, pp. 621– 625, 2002. [17] I. Solodov, N. Krohn, and G. Busse, “Nonlinear Ultrasonic NDT for Early Defect Recognition and Imaging,” in European Conf. on NDT (ECNDT), Moscow, Citeseer, 2010. [18] A. Klepka, T. Stepinski, T. Uhl, and W. Staszewski, “Nonlinear acoustics,” Advanced Structural Damage Detection: From Theory to Engineering Applications, pp. 73–107, 2013. [19] B. Sarens, B. Verstraeten, C. Glorieux, and G. Kalogiannakis, “Investigation of contact acoustic nonlinearity in delaminations by shearographic imaging, laser doppler vibrometeric scanning and finite difference modeling,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency control, vol. 57, no. 6, pp. 1383–1395, 2010. [20] S. Delrue and K. Van Den Abeele, “Three-dimensional finite element simulation of closed delaminations in composite materials,” Ultrasonics, vol. 52, no. 2, pp. 315–324, 2012. [21] S. Biwa and Y. Ishii, “Second-harmonic generation in an infinite layered structure with nonlinear spring-type interfaces,” Wave Motion, vol. 63, pp. 55–67, 2016. [22] F. Ju, H. P. Lee, and K. H. Lee, “Dynamic response of delaminated composite beams with intermittent contact in delaminated segments,” Composites Engineering, vol. 4, no. 12, pp. 1211–1224, 1994. 29

˙ [23] A. Zak, M. Krawczuk, and W. Ostachowicz, “Vibration of a laminated composite plate with closing delamination,” Journal of Intelligent Material Systems and Structures, vol. 12, no. 8, pp. 545–551, 2001. [24] S. Biwa, S. Nakajima, and N. Ohno, “On the acoustic nonlinearity of solid-solid contact with pressure-dependent interface stiffness,” Journal of applied mechanics, vol. 71, no. 4, pp. 508–515, 2004. [25] M. K¨ogl, S. Hurlebaus, and L. Gaul, “Finite element simulation of nondestructive damage detection with higher harmonics,” NDT & E International, vol. 37, no. 3, pp. 195–205, 2004. [26] U. Andreaus, P. Casini, and F. Vestroni, “Non-linear dynamics of a cracked cantilever beam under harmonic excitation,” International journal of non-linear mechanics, vol. 42, no. 3, pp. 566–575, 2007. [27] V. N. Burlayenko and T. Sadowski, “Finite element nonlinear dynamic analysis of sandwich plates with partially detached facesheet and core,” Finite Elements in Analysis and Design, vol. 62, pp. 49–64, 2012. [28] D. Garcia, R. Palazzetti, I. Trendafilova, C. Fiorini, and A. Zucchelli, “Vibration-based delamination diagnosis and modelling for composite laminate plates,” Composite Structures, vol. 130, pp. 155–162, 2015. [29] M. Yuan, T. Lee, T. Kang, J. Zhang, S. J. Song, and H. J. Kim, “Absolute measurement of ultrasonic non-linearity parameter at contact interface,” Nondestructive Testing and Evaluation, vol. 30, no. 4, pp. 356– 372, 2015.

30

[30] D. Broda, W. J. Staszewski, A. Martowicz, T. Uhl, and V. V. Silberschmidt, “Modelling of nonlinear crack-wave interactions for damage detection based on ultrasound–A review,” Journal of Sound and Vibration , vol. 333, no. 4, pp. 1097 – 1118, 2014. [31] A. Bovsunovsky and C. Surace, “Non-linearities in the vibrations of elastic structures with a closing crack: A state of the art review,” Mechanical Systems and Signal Processing, vol. 62, pp. 129–148, 2015. [32] K. L. Johnson, Contact mechanics. Cambridge university press, 1987. [33] P. Wriggers and T. A. Laursen, Computational contact mechanics. Springer, 2006. [34] Abaqus Documentation, “Version 6.14-2,” Dassault Syst´emes Simulia Corp., Providence, RI. [35] N. Krohn, R. Stoessel, and G. Busse, “Acoustic non-linearity for defect selective imaging,” Ultrasonics, vol. 40, no. 1, pp. 633–637, 2002. [36] I. Solodov, J. Wackerl, K. Pfleiderer, and G. Busse, “Nonlinear selfmodulation and subharmonic acoustic spectroscopyfor damage detection and location,” Applied physics letters, vol. 84, no. 26, pp. 5386–5388, 2004. [37] K. Kawashima, M. Murase, R. Yamada, M. Matsushima, M. Uematsu, and F. Fujita, “Nonlinear ultrasonic imaging of imperfectly bonded interfaces,” Ultrasonics, vol. 44, pp. e1329–e1333, 2006.

31

[38] U. Polimeno, M. Meo, D. P. Almond, and S. L. Angioni, “Detecting low velocity impact damage in composite plate using nonlinear acoustic/ultrasound methods,” Applied Composite Materials, vol. 17, no. 5, pp. 481–488, 2010. [39] I. Solodov, D. D¨oring, and G. Busse, “New opportunities for NDT using non-linear interaction of elastic waves with defects,” Strojniˇski vestnikJournal of Mechanical Engineering, vol. 57, no. 3, pp. 169–182, 2011. [40] L. Pieczonka, A. Klepka, W. J. Staszewski, and T. Uhl, “Nonlinear acoustic imaging of structural damages in laminated composites,” in EWSHM-7th European Workshop on Structural Health Monitoring 2014, 2014. [41] COMSOL Multiphysics Documentation, “Version 5.1,” COMSOL Inc. [42] S. Delrue and K. Van Den Abeele, “Simulations of nonlinear air-coupled emission (NACE),” in International Congress on Ultrasonics, pp. 510– 515, 2013. [43] S. Delrue and K. Van Den Abeele, “A simulation study of Local Defect Resonances (LDR),” in Emerging Technologies in Non-Destructive Testing VI: Proceedings of the 6th International Conference on Emerging Technologies in Non-Destructive Testing (Brussels, Belgium, 27-29 May 2015), p. 167, CRC Press, 2015. [44] S. Delrue and K. Van Den Abeele, “Detection of defect parameters using nonlinear air-coupled emission by ultrasonic guided waves at contact acoustic nonlinearities,” Ultrasonics, vol. 63, pp. 147–154, 2015. 32

[45] S. Delrue, M. Tabatabaeipour, J. Hettler, and K. Van Den Abeele, “Applying a nonlinear, pitch-catch, ultrasonic technique for the detection of kissing bonds in friction stir welds,” Ultrasonics, vol. 68, pp. 71–79, 2016.

7. Appendix The following python script can be run to produce the contact model in Abaqus/Standard version 6.14. from abaqus import ∗ from abaqusConstants import ∗ from caeModules import ∗ from d r i v e r U t i l s import executeOnCaeStartup import math executeOnCaeStartup ( ) Mdb( )

#=============================Parameters============================ j0

= ’3 d defect ’

# F i l e name

c0

= 12

# Number o f cpu c o r e s

ht

= 2 . 0 e−3

# Height of the p l a t e

lt

= 4 0 . 0 e−3

# Length o f t h e p l a t e

wt

= 4 0 . 0 e−3

# Width 33

of the plate

dd

= 0 . 3 0 e−3

# Depth

dr

= 1 0 . 0 e−3

# Radius o f d e l a m i n a t i o n

m0

= 0 . 6 e−3

# Mesh s i z e

f0

= 26165.0

# Excitation frequency

A0

= 1 . 0 e−6

# Excitation amplitude

T0

= 1.0/ f0

# Time p e r i o d

w0

= 2 . 0 ∗ math . p i ∗ f 0

# Circular excitation frequency

z0

= 0.1

# Damping r a t i o

a0

= 2 . 0 ∗ ( 0 . 5 ∗ z0 ) ∗w0

# Mass p r o p o r t i o n a l damping c o e f f i c i e n t

b0

= 2 . 0 ∗ ( 0 . 5 ∗ z0 ) /w0

# S t i f f n e s s p r o p o r t i o n a l damping c o e f f i c i e n t

N0

= 45.0

# Number o f time p e r i o d s

n0

= 32.0

# R e s o l u t i o n o f one time p e r i o d

of delamination

Tt0 = N0∗T0 ;

# T o t a l s i m u l a t i o n time

t s 0 = T0/n0

# S t e p time

E1

= 1 6 1 . 0 e9

# Young ’ s modulus a l o n g z d i r e c t i o n

E2

= 1 1 . 3 8 e9

# Young ’ s modulus a l o n g x d i r e c t i o n

E3

= E2

# Young ’ s modulus a l o n g y d i r e c t i o n

G12 = 5 . 1 7 e9

# Shear modulus

G23 = 3 . 9 8 e9

# Shear modulus

G13 = G12

# Shear modulus

v12 = 0 . 3 2

# Poisson ’ s r a t i o

v23 = 0 . 4 3 6

# Poisson ’ s r a t i o

v13 = v12

# Poisson ’ s r a t i o

rh0 = 1 . 8 e3

# Density

#=================================================================== 34

#=========================Geometry================================== model = mdb . models [ ’ Model−1 ’ ] # Create Part −1 s = model . C o n s t r a i n e d S k e t c h ( name= ’

profile

’ , s h e e t S i z e =200.0)

g , v , d , c = s . geometry , s . v e r t i c e s , s . dimensions , s . c o n s t r a i n t s s . s e t P r i m a r y O b j e c t ( o p t i o n=STANDALONE) s . C o n s t r u c t i o n L i n e ( p o i n t 1 = ( 0 . 0 , −100.0) , p o i n t 2 = ( 0 . 0 , 1 0 0 . 0 ) ) s . r e c t a n g l e ( p o i n t 1 = ( 0 . 0 , ht−dd ) , p o i n t 2 =(dr , ht ) ) p1 = model . Part ( name= ’ Part −1 ’ , d i m e n s i o n a l i t y=THREE D, type=DEFORMABLE BODY) p1 . B a s e S o l i d R e v o l v e ( s k e t c h=s , a n g l e =360.0 , f l i p R e v o l v e D i r e c t i o n=OFF) s . uns etP rim ary Obj ect ( ) del model . s k e t c h e s [ ’

profile

’]

# Create Part −2 s = model . C o n s t r a i n e d S k e t c h ( name= ’

profile

’ , s h e e t S i z e =200.0)

g , v , d , c = s . geometry , s . v e r t i c e s , s . dimensions , s . c o n s t r a i n t s s . s e t P r i m a r y O b j e c t ( o p t i o n=STANDALONE) s . C o n s t r u c t i o n L i n e ( p o i n t 1 = ( 0 . 0 , −100.0) , p o i n t 2 = ( 0 . 0 , 1 0 0 . 0 ) ) s . F i x e d C o n s t r a i n t ( e n t i t y=g . f i n d A t ( ( 0 . 0 , 0 . 0 ) ) ) s . r e c t a n g l e ( p o i n t 1 = ( 0 . 0 , 0 . 0 ) , p o i n t 2 =(dr , ht−dd ) ) p2 = model . Part ( name= ’ Part −2 ’ , d i m e n s i o n a l i t y=THREE D, type=DEFORMABLE BODY) p2 . B a s e S o l i d R e v o l v e ( s k e t c h=s , a n g l e =360.0 , f l i p R e v o l v e D i r e c t i o n=OFF) s . uns etP rim ary Obj ect ( ) 35

del model . s k e t c h e s [ ’

profile

’]

# Create Part −3 s = model . C o n s t r a i n e d S k e t c h ( name= ’

profile

’ , s h e e t S i z e =200.0)

g , v , d , c = s . geometry , s . v e r t i c e s , s . dimensions , s . c o n s t r a i n t s s . s e t P r i m a r y O b j e c t ( o p t i o n=STANDALONE) s . r e c t a n g l e ( p o i n t 1=(−wt / 2 . 0 , − l t / 2 . 0 ) , p o i n t 2 =(wt / 2 . 0 , l t / 2 . 0 ) ) s . C i r c l e B y C e n t e r P e r i m e t e r ( c e n t e r = ( 0 . 0 , 0 . 0 ) , p o i n t 1 =(dr , 0 . 0 ) ) p3 = model . Part ( name= ’ Part −3 ’ , d i m e n s i o n a l i t y=THREE D, type=DEFORMABLE BODY) p3 . B a s e S o l i d E x t r u d e ( s k e t c h=s , depth=ht ) s . uns etP rim ary Obj ect ( ) del model . s k e t c h e s [ ’

profile

’]

# Make a s s e m b l y a = model . rootAssembly a . DatumCsysByDefault (CARTESIAN) a1 = a . I n s t a n c e ( name= ’ Part −1−1 ’ , p a r t=p1 , dependent=ON) a2 = a . I n s t a n c e ( name= ’ Part −2−1 ’ , p a r t=p2 , dependent=ON) a3 = a . I n s t a n c e ( name= ’ Part −3−1 ’ , p a r t=p3 , dependent=ON) a . r o t a t e ( i n s t a n c e L i s t =( ’ Part −3−1 ’ , ) , a x i s P o i n t = ( 0 . 0 , 0 . 0 , 0 . 0 ) , a x i s D i r e c t i o n = ( 1 . 0 , 0 . 0 , 0 . 0 ) , a n g l e =90.0) a . t r a n s l a t e ( i n s t a n c e L i s t =( ’ Part −3−1 ’ , ) , v e c t o r = ( 0 . 0 , ht , 0 . 0 ) ) #===================================================================

#=====================M a t e r i a l and s e c t i o n s========================= m a t e r i a l 1 = model . M a t e r i a l ( name= ’ M a t e r i a l −1 ’ ) 36

m a t e r i a l 1 . E l a s t i c ( type=ENGINEERING CONSTANTS, t a b l e =((E1 , E2 , E3 , v12 , v13 , v23 , G12 , G13 , G23 ) , ) ) m a t e r i a l 1 . D e n s i t y ( t a b l e =(( rh0 , ) , ) ) m a t e r i a l 1 . Damping ( be ta=b0 , alpha=a0 ) # Create s e c t i o n model . H o m o g e n e o u s S o l i d S e c t i o n ( name= ’ S e c t i o n −1 ’ , m a t e r i a l= ’ M a t e r i a l −1 ’ , t h i c k n e s s=None ) # Assign s e c t i o n Part −1 datum1 = p1 . DatumCsysByThreePoints ( name= ’Datum c s y s −1 ’ , coordSysType=CARTESIAN, o r i g i n = ( 0 . 0 , 0 . 0 , 0 . 0 ) , point1 =(0.0 , 0.0 , 1 . 0 ) , point2 =(1.0 , 0.0 , 1 . 0 ) ) c = p1 . c e l l s c e l l s = c . f i n d A t ( ( ( 0 . 0 , ht , 0 . 0 ) , ) ) r e g i o n = r e g i o n T o o l s e t . Region ( c e l l s = c e l l s ) p1 . S e c t i o n A s s i g n m e n t ( r e g i o n=r e g i o n , sectionName= ’ S e c t i o n −1 ’ , o f f s e t =0.0 , o f f s e t T y p e=MIDDLE SURFACE, o f f s e t F i e l d= ’ ’ , t h i c k n e s s A s s i g n m e n t=FROM SECTION) o r i e n t a t i o n = p1 . datums [ datum1 . id ] p1 . M a t e r i a l O r i e n t a t i o n ( r e g i o n=r e g i o n , o r i e n t a t i o n T y p e=SYSTEM, a x i s=AXIS 3 , l o c a l C s y s=o r i e n t a t i o n , fieldName= ’ ’ , a d d i t i o n a l R o t a t i o n T y p e=ROTATION NONE, a n g l e =0.0 , a d d i t i o n a l R o t a t i o n F i e l d= ’ ’ , s t a c k D i r e c t i o n=STACK 3) # Assign s e c t i o n Part −2 datum2 = p2 . DatumCsysByThreePoints ( name= ’Datum c s y s −1 ’ , coordSysType=CARTESIAN, o r i g i n = ( 0 . 0 , 0 . 0 , 0 . 0 ) , 37

point1 =(0.0 , 0.0 , 1 . 0 ) , point2 =(1.0 , 0.0 , 1 . 0 ) ) c = p2 . c e l l s c e l l s = c . findAt ( ( ( 0 . 0 , 0.0 , 0 . 0 ) , )) r e g i o n = r e g i o n T o o l s e t . Region ( c e l l s = c e l l s ) p2 . S e c t i o n A s s i g n m e n t ( r e g i o n=r e g i o n , sectionName= ’ S e c t i o n −1 ’ , o f f s e t =0.0 , o f f s e t T y p e=MIDDLE SURFACE, o f f s e t F i e l d= ’ ’ , t h i c k n e s s A s s i g n m e n t=FROM SECTION) o r i e n t a t i o n = p2 . datums [ datum2 . id ] p2 . M a t e r i a l O r i e n t a t i o n ( r e g i o n=r e g i o n , o r i e n t a t i o n T y p e=SYSTEM, a x i s=AXIS 3 , l o c a l C s y s=o r i e n t a t i o n , fieldName= ’ ’ , a d d i t i o n a l R o t a t i o n T y p e=ROTATION NONE, a n g l e =0.0 , a d d i t i o n a l R o t a t i o n F i e l d= ’ ’ , s t a c k D i r e c t i o n=STACK 3) # Assign s e c t i o n Part −3 datum3 = p3 . DatumCsysByThreePoints ( name= ’Datum c s y s −1 ’ , coordSysType=CARTESIAN, o r i g i n = ( 0 . 0 , 0 . 0 , 0 . 0 ) , point1 =(0.0 , 1.0 , 0 . 0 ) , point2 =(1.0 , 1.0 , 0 . 0 ) ) c = p3 . c e l l s c e l l s = c . f i n d A t ( ( ( wt / 2 . 0 , l t / 2 . 0 , 0 . 0 ) , ) ) r e g i o n = r e g i o n T o o l s e t . Region ( c e l l s = c e l l s ) p3 . S e c t i o n A s s i g n m e n t ( r e g i o n=r e g i o n , sectionName= ’ S e c t i o n −1 ’ , o f f s e t =0.0 , o f f s e t T y p e=MIDDLE SURFACE, o f f s e t F i e l d= ’ ’ , t h i c k n e s s A s s i g n m e n t=FROM SECTION) o r i e n t a t i o n = p3 . datums [ 2 ] p3 . M a t e r i a l O r i e n t a t i o n ( r e g i o n=r e g i o n , o r i e n t a t i o n T y p e=SYSTEM, a x i s=AXIS 3 , l o c a l C s y s=o r i e n t a t i o n , fieldName= ’ ’ , 38

a d d i t i o n a l R o t a t i o n T y p e=ROTATION NONE, a n g l e =0.0 , a d d i t i o n a l R o t a t i o n F i e l d= ’ ’ , s t a c k D i r e c t i o n=STACK 3)

#===========================Meshing================================= elemType1 = mesh . ElemType ( elemCode=C3D20R , e l e m L i b r a r y=STANDARD) elemType2 = mesh . ElemType ( elemCode=C3D15 , e l e m L i b r a r y=STANDARD) elemType3 = mesh . ElemType ( elemCode=C3D10 , e l e m L i b r a r y=STANDARD) # Mesh Part −1 p1 . s e e d P a r t ( s i z e=m0, d e v i a t i o n F a c t o r =0.1 , m i n S i z e F a c t o r =0.1) c = p1 . c e l l s c e l l s = c . f i n d A t ( ( ( 0 . 0 , ht , 0 . 0 ) , ) ) p1 . s e t M e s h C o n t r o l s ( r e g i o n s=c e l l s , elemShape=WEDGE) p i c k e d R e g i o n s =( c e l l s , ) p1 . setElementType ( r e g i o n s=p i c k e d R e g i o n s , elemTypes=(elemType1 , elemType2 , elemType3 ) ) p1 . generateMesh ( ) # Mesh Part −2 p2 . s e e d P a r t ( s i z e =2.0∗m0, d e v i a t i o n F a c t o r =0.1 , m i n S i z e F a c t o r =0.1) c = p2 . c e l l s c e l l s = c . findAt ( ( ( 0 . 0 , 0.0 , 0 . 0 ) , )) p2 . s e t M e s h C o n t r o l s ( r e g i o n s=c e l l s , elemShape=WEDGE) p i c k e d R e g i o n s =( c e l l s , ) p2 . setElementType ( r e g i o n s=p i c k e d R e g i o n s , elemTypes=(elemType1 , elemType2 , elemType3 ) ) p2 . generateMesh ( ) 39

# Mesh Part −3 p3 . s e e d P a r t ( s i z e =2.5∗m0, d e v i a t i o n F a c t o r =0.1 , m i n S i z e F a c t o r =0.1) c = p3 . c e l l s p i c k e d R e g i o n s = c . f i n d A t ( ( ( wt / 2 . 0 , l t / 2 . 0 , 0 . 0 ) , ) ) p3 . s e t M e s h C o n t r o l s ( r e g i o n s=p i c k e d R e g i o n s , elemShape=TET, t e c h n i q u e=FREE) c e l l s = c . f i n d A t ( ( ( wt / 2 . 0 , l t / 2 . 0 , 0 . 0 ) , ) ) p i c k e d R e g i o n s =( c e l l s , ) p3 . setElementType ( r e g i o n s=p i c k e d R e g i o n s , elemTypes=(elemType1 , elemType2 , elemType3 ) ) e = p3 . e d g e s p ickedEdges = e . f i n d A t ( ( ( dr , 0 . 0 , ht ) , ) ) p3 . seedEdgeBySize ( e d g e s=pickedEdges , s i z e =2.0∗m0, d e v i a t i o n F a c t o r =0.1 , m i n S i z e F a c t o r =0.1 , c o n s t r a i n t=FIXED) p ickedEdges = e . f i n d A t ( ( ( dr , 0 . 0 , 0 . 0 ) , ) ) p3 . seedEdgeBySize ( e d g e s=pickedEdges , s i z e =1.0∗m0, d e v i a t i o n F a c t o r =0.1 , m i n S i z e F a c t o r =0.1 , c o n s t r a i n t=FIXED) p3 . generateMesh ( ) #===================================================================

#=========================S e l e c t i o n s================================ s 1 = a3 . f a c e s s i d e 1 F a c e s 1 = s1 . f i n d A t ( ( ( dr , ht / 2 . 0 , 0 . 0 ) , ) ) a . S u r f a c e ( s i d e 1 F a c e s=s i d e 1 F a c e s 1 , name= ’ Part−3−s u r f a c e ’ ) v1 = a3 . v e r t i c e s 40

v e r t s 1 = v1 . f i n d A t ((( − wt / 2 . 0 , 0 . 0 , l t / 2 . 0 ) , ) ) a . S et ( v e r t i c e s=v e r t s 1 , name= ’ Part−3−s e n s o r −t i e ’ ) s 1 = a2 . f a c e s s i d e 1 F a c e s 1 = s1 . f i n d A t ( ( ( dr , ( ht−dd ) / 2 . 0 , 0 . 0 ) , ) ) a . S u r f a c e ( s i d e 1 F a c e s=s i d e 1 F a c e s 1 , name= ’ Part−2−s u r f a c e ’ ) s i d e 1 F a c e s 1 = s1 . f i n d A t ( ( ( 0 . 0 , ht−dd , 0 . 0 ) , ) ) a . S u r f a c e ( s i d e 1 F a c e s=s i d e 1 F a c e s 1 , name= ’ Part−2−c o n t a c t −s u r f a c e ’ ) s 1 = a1 . f a c e s s i d e 1 F a c e s 1 = s1 . f i n d A t ( ( ( dr , ht−dd / 2 . 0 , 0 . 0 ) , ) ) a . S u r f a c e ( s i d e 1 F a c e s=s i d e 1 F a c e s 1 , name= ’ Part−1−s u r f a c e ’ ) s i d e 1 F a c e s 1 = s1 . f i n d A t ( ( ( 0 . 0 , ht−dd , 0 . 0 ) , ) ) a . S u r f a c e ( s i d e 1 F a c e s=s i d e 1 F a c e s 1 , name= ’ Part−1−c o n t a c t −s u r f a c e ’ ) #===================================================================

#===================Study and BCs=================================== # Study s t e p model . I m p l i c i t D y n a m i c s S t e p ( name= ’ Step −1 ’ , p r e v i o u s= ’ I n i t i a l ’ , t i m e P e r i o d=Tt0 , maxNumInc=1000000 , i n i t i a l I n c=t s 0 , minInc=1e −09 , maxInc=t s 0 ) model . f i e l d O u t p u t R e q u e s t s [ ’F−Output−1 ’ ] . s e t V a l u e s ( f r e q u e n c y =1) # Define t i e c o n s t r a i n t s r e g i o n 1=a . s u r f a c e s [ ’ Part−3−s u r f a c e ’ ] r e g i o n 2=a . s u r f a c e s [ ’ Part−2−s u r f a c e ’ ] model . Tie ( name= ’ C o n s t r a i n t −1 ’ , master=r e g i o n 2 , s l a v e=r e g i o n 1 , p o s i t i o n T o l e r a n c e M e t h o d=COMPUTED, a d j u s t=ON, t i e R o t a t i o n s=ON, 41

t h i c k n e s s=ON) r e g i o n 1=a . s u r f a c e s [ ’ Part−3−s u r f a c e ’ ] r e g i o n 2=a . s u r f a c e s [ ’ Part−1−s u r f a c e ’ ] model . Tie ( name= ’ C o n s t r a i n t −2 ’ , master=r e g i o n 1 , s l a v e=r e g i o n 2 , p o s i t i o n T o l e r a n c e M e t h o d=COMPUTED, a d j u s t=ON, t i e R o t a t i o n s=ON, t h i c k n e s s=ON) # D e f i n e c o n t a c t BC r e g i o n 1=a . s u r f a c e s [ ’ Part−2−c o n t a c t −s u r f a c e ’ ] r e g i o n 2=a . s u r f a c e s [ ’ Part−1−c o n t a c t −s u r f a c e ’ ] i n t p r o p = model . ContactProperty ( ’ IntProp −1 ’ ) i n t p r o p . T a n g e n t i a l B e h a v i o r ( f o r m u l a t i o n=FRICTIONLESS) i n t p r o p . NormalBehavior ( p r e s s u r e O v e r c l o s u r e=HARD, a l l o w S e p a r a t i o n=ON, c o n s t r a i n t E n f o r c e m e n t M e t h o d=DEFAULT) model . S u r f a c e T o S u r f a c e C o n t a c t S t d ( name= ’ Int −1 ’ , createStepName= ’ Step −1 ’ , master=r e g i o n 1 , s l a v e=r e g i o n 2 , s l i d i n g=FINITE , t h i c k n e s s=ON, i n t e r a c t i o n P r o p e r t y= ’ IntProp −1 ’ , adjustMethod=NONE, i n i t i a l C l e a r a n c e=OMIT, datumAxis=None , c l e a r a n c e R e g i o n=None ) # D e f i n e e x c i t a t i o n BC model . P e r i o d i c A m p l i t u d e ( name= ’Amp−1 ’ , timeSpan=STEP, f r e q u e n c y=w0 , s t a r t =0.0 , a 0 =0.0 , data = ( ( 0 . 0 , 1 . 0 ) , ) ) r e g i o n=a . s e t s [ ’ Part−3−s e n s o r −t i e ’ ] model . DisplacementBC ( name= ’BC−2 ’ , createStepName= ’ Step −1 ’ , r e g i o n=r e g i o n , u1=UNSET, u2=A0 , u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, amplitude= ’Amp−1 ’ , f i x e d=OFF, 42

d i s t r i b u t i o n T y p e=UNIFORM, fieldName= ’ ’ , l o c a l C s y s=None ) # Create and s u b m i t j o b mdb . Job ( name=j0 , model= ’ Model−1 ’ , d e s c r i p t i o n= ’ ’ , type=ANALYSIS, atTime=None , waitMinutes =0, waitHours =0, queue=None , memory=90 , memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True , e x p l i c i t P r e c i s i o n=SINGLE , n o d a l O u t p u t P r e c i s i o n=FULL, e c h o P r i n t=OFF, modelPrint=OFF, c o n t a c t P r i n t=OFF, h i s t o r y P r i n t=OFF, u s e r S u b r o u t i n e= ’ ’ , s c r a t c h= ’ ’ , r e s u l t s F o r m a t=ODB, m u l t i p r o c e s s i n g M o d e=THREADS, numCpus=c0 , numGPUs=0, numDomains=c0 ) mdb . saveAs ( pathName=j 0 ) mdb . j o b s [ j 0 ] . submit ( c o n s i s t e n c y C h e c k i n g=OFF) #===================================================================

43

Highlights 1) A simple easily replicable contact model for nonlinear acoustics. 2) The modelling can be done without the need of any additional nonlinearity parameters. 3) The model provides explanation of harmonics and subharmonics in the frequency spectrum. 4) The results compare well with nonlinear spring damper model available in literature.