Nonlinear Modelling and Simulation of Impact Events and Validation with Experimental Test

STOCKHOLM, SWEDEN 2018

Nonlinear Modelling and

Simulation of Impact Events and

Validation with Experimental Test

ANDRÁS PÉTER SZÉKELY

NICOLE AL HANNA

1

Nonlinear Modelling and

Simulation of Impact Events

and Validation with

Experimental Test

NICOLE AL HANNA

ANDRÁS PÉTER SZÉKELY

TRITA-SCI-GRU 2019:042 Master Thesis, 2019

KTH Royal Institute of Technology Engineering Sciences

3

Acknowledgements

A big thank you to our supervisor at Volvo Car Corporation, Mohsen Bayani Khaknejad, our co-supervisor, Henrik Viktorsson, and everyone at the Solidity Department who helped accomplish this thesis work.

We would also like to thank our examiner, Per Wennhage, our university, KTH Royal Institute of Technology, and all the professors in the Aeronautical and Vehicle Engineering Department for providing a professional and friendly environment suitable for knowledge and self-growth.

4

Abstract

5

Sammanfattning

6

Contents

Current CAE rattle predicting methods 10

Impact modelling 10

Linearity and non-linearity in FEM 11

Explicit and implicit solvers 11

Method 12

Physical tests 13

Experimental setup 13

Modal measurements 14

Time domain measurements 15

Signal processing 15

Mean removal 18

Butterworth filter 18

Moving average filter 19

Omega arithmetic 20

Setbacks and Solutions 22

FEA model 23

Mesh 24

Boundary conditions 24

Excitation 25

Material and property cards 27

Contact modelling 28

Solver Settings 29

Model validation 30

Modal correlation 30

Time domain correlation before impact event 31

Gap Size Determination 33

Time domain correlation at impact 35

Sensitivity Analysis 35

Nominated parameters and the ranges 36

Evaluating the Fundamental Rattle Criteria 38 Results of the sensitivity analysis and guidelines 41 The Effect of Contact Surface Definition 43

The Effect of Friction 44

The Effect of Contact Mechanics 47

The Effect of Material Damping 49

Guidelines 52

Conclusion 56

Findings and Lessons Learned 56

Future work 56

Distribution of labour 58

Appendix 59

7

Nomenclature

α Mass-proportional Rayleigh damping coefficient [Hz]

β Stiffness-proportional Rayleigh damping coefficient [s]

[C] Damping matrix [N/(m/s)] [K] Stiffness matrix [N/m] [M] Mass matrix [kg] δ Logarithmic decrement [-] 𝜁 Damping ratio [-] n Number of periods T Period [s] x(t) Displacement at time t [m]

µ0 Critical damping fraction [-]

∆t Time increment [s]

ω Frequency [Rad/s]

f Frequency [Hz]

t0 Time at which the start curve is connected to the rest of the excitation signal [s]

8

Abbreviations

VCC Volvo Car Corporation

NVH Noise, Vibration, Harshness

CAE Computer Aided Engineering

FEA Finite Element Analysis

FEM Finite Element Modelling

FFT Fast Fourier Transformation

S&R Squeak and Rattle

BC(s) Boundary Condition(s)

DoF Degrees of Freedom

MPC Multi-Point Constraint St-St Steel-Steel P-St Polypropylene-Steel P-P Polypropylene-Polypropylene CP Contact pair Gen General Kin Kinematic Pen Penalty

Mdamp Mass proportional material damping

Sdamp Stiffness proportional material damping

El-el Surfaces of the material pair modelled as element-based

Node-el Surfaces of the material pair are modelled one as element based and the other as node-based

Introduction

At Volvo Car Corporation (VCC), the Solidity department is concerned with Squeak & Rattle (S&R), two types of unwanted sounds that could occur while driving, and which the user may consider an indication of low product quality. The ideal scenario for the engineers at VCC would be to be able to predict the location, intensity and annoyance level of those sounds, and take the necessary measures to eliminate them, before the manufacturing starts. That way the development process is made faster and more efficient.

In this thesis work, rattle is the phenomenon under investigation. It is an impact induced noise defined as “a rapid succession of short, sharp, sounds” (1) and occurs in vehicles when two parts interact in an oscillating contact-no-contact manner. The engineers’ primary tools are CAE softwares, such as ABAQUS, NASTRAN, and LS-DYNA, which manifest powerful prediction abilities in the automotive world. Rattle prediction, however, is still not fully possible due to its complicated nature, adding to the fact that it has only recently become of bigger importance with the boom of car electrification.

9 Background

VCC has invested a lot of effort in research to find powerful tools that are able to predict S&R before the first prototypes are built. The E-LINE method (2), which is currently used at VCC for rattle prediction, evaluates the clearances between parts and, in the predicted contact areas, statistically evaluates the relative displacements in time domain. However, without modelling the impact events, little information is given on the variables most significant in rattle sound generation, which makes the assessment of these predicted contact areas less reliable.

Current efforts are put into investigating the connection between sound annoyance levels and the fundamental variables of different rattle events, and creating a simulation-based methodology with a reliable assessment method of the predicted rattle. This work is carried out as a doctoral thesis, split up into different parts, one of them being this current master thesis. This part of the study is to create an accurate method to model the rattle events in FEA, with special attention towards the dependency of the fundamental rattle variables on the parameters and settings of the FEA model.

The base of this work is a previous thesis project (3), where a physical test rig was designed, built and used for rattle measurements. The main aim of these measurements was study the impact of ambient conditions and variations in boundary conditions on the generated rattle sounds, based on objective and subjective rattle sound evaluation. However, these measurements also captured acceleration and force data so the annoyance levels, rattle variables and simulation parameters can be associated directly. Therefore, the findings of the whole study have a higher accuracy, since the sound evaluation, and the FEA modelling are based on the same measurement.

Objectives

The deliverables expected by July 7, 2018 were: 1. Build an FEA model of the test rig

2. Correlate the CAE rattle simulations with the experimental measurements

3. Conduct a sensitivity analysis to achieve a better understanding of the factors affecting rattle modelling

4. Write a set of guidelines for modelling rattle with CAE

Limitations

10

way to see it in the same condition as it was used to make the measurements, so it was not possible to check any undocumented details that could have an effect on the behaviour of the setup. This limitation became an even bigger problem when signal processing details contradicted the given specifications that were supposed to serve as a foundation for the correlations. All these missing details and resulting setbacks are documented in section 6 where approaches used to solve them are explained in detail.

More limitations include the fact that the experiment was already conducted. Therefore, only the captured data could be used for the correlation part, with errors (described in the section 5) and without any means of using feedback from the simulations. Additionally, when an attempt was made to recreate the rattle experiment (explained in section 6), restrictions in the signal generating software prevented the use of time domain input, and the creation of a more suitable measurement for the correlation.

CAE software were limited by what is used at VCC. ANSA 17.1.1 was used for pre-processing, META 17.1.1 for post-pre-processing, and ABAQUS 6.14 and 2017 as solver. There was a possibility to use Nastran and LS-Dyna solvers, but, due to the limitations in time, this did not happen.

Finally, the biggest limitation was time. However, we did manage to solve most of the problems, get a correlation, and perform the sensitivity analysis that would be of use for future rattle simulations and validation on car parts.

Literature study

Current CAE rattle predicting methods

Two methods are used in the industry today to predict rattle. Both are linear, use standard NVH models, and measure relative displacement between the critical parts. One of them, the SAR-LINETM _{or E-line (2), operates in time domain and utilizes} user defined 3D lines to make the relative displacement calculations. The other method operates in frequency domain and evaluates the displacement by using connector elements (in ABAQUS terminology) referred to as “virtual sensors” (4). In both methods, rattle is detected if the gap size and tolerance are exceeded.

The limitation of the methods is that they can only predict the possibility of rattle, rather than its certainty. That is because contact does not necessarily mean a sound will be produced. Therefore, a method that could detect only the impacts that will for sure create rattle sounds is sought after in the current study.

Impact modelling

11

experiment, and showed the significance of non-linear time domain modelling of the vibration problem, as linear motion transmissibility could not capture the vibration phenomenon and predict acceleration levels in the parts. The authors established a soft constrained penalty-based contact with static and dynamic friction coefficients and viscous damping. They also give a good way to reduce the computational costs of long duration simulations, while maintaining accurate contact modelling, by employing a dynamic sub-structuring strategy. This strategy is based on dividing the structure to parts that behave linearly and nonlinearly, and then reduce the linear substructures’ degree of freedom to speed up the computation.

In (6), a study of several parameters available in RADIOSS for contact modelling is performed. The model consists of a table and two billiard balls, all made of linear elastic materials. The table is fixed in all degrees of freedom at its bottom, and the balls are made to impact each other, defining a non-linear behavior. The parameters varied are the contact interfaces between the balls (types 7, 16 and 17 of RADIOSS), the computation method upon node-element penetration between the table and the balls (tied or sliding definition), and the algorithm for constraint calculation upon contact (penalty or Lagrange multiplier method). After running the simulations and comparing the kinetic energies obtained with analytical results, it is found that the sliding option and type 16 interface give the smallest error.

Linearity and non-linearity in FEM

The linearity in FEA is defined by the mathematical stiffness and load matrices of the modelled structure. In linear FEA, these matrices remain constant during the simulation (7). There is no need for the software to update stiffness and load components when deformation occurs. In nonlinear analyses, which may result from geometry, material properties, contact, the stiffness or load matrices vary constantly. A time increment is then needed, since the stiffness and force components have to be updated at every step of the structure. This is due to the fact that the nonlinear phenomena present in the event make the previously defined matrices unable to describe the structure accurately.

The current study deals with nonlinearities resulting from the contact between the parts coming together to create rattle. In finite element tools, contact modelling can be defined by preventing (kinematic method) or eliminating (penalty method) the penetration of the parts involved using different methods. These methods are further described in 9.1, and the difference between them is of interest in the sensitivity analysis.

Explicit and implicit solvers

12

The analytical difference between the two approaches lies in the number of unknowns to be evaluated at every time increment (∆t). The implicit method, which is an iterative technique, calculates for information about the current and the later states of the system, so coupled equations have to be solved at every time increment. On the other hand, the explicit method utilizes de-coupled equations since it only calculates for the future state based on the current one.

The implicit method is unconditionally stable, meaning it converges irrespective of ∆t. The explicit method, due to the nature of its equations, is stable only if ∆t is below a certain critical value. It is best used for short time simulations because of the ∆t requirement and one of its main advantages is that it is less computationally expensive than an implicit simulation with the same number of time increments. Increasing ∆t means less computations for the implicit method, but it also hinders the accuracy, which could be of importance depending on the behaviour under analysis.

In ABAQUS, two solvers are available: ABAQUS/Standard which calculates using the implicit method, and ABAQUS/Explicit.

Method

As mentioned in section 1.1, prior to this thesis, a physical setup was designed and built by previous thesis students to test for rattle (3). The current thesis aimed to build an FEA model of this physical setup, achieve a good correlation in time domain simulations and use this model for a sensitivity study on the FEA parameters. The correlation was reached in three different steps:

1. Modal measurements were performed in anticipation of the current work. That is, for the validation of the FEA model. This experiment was performed again during this thesis to serve as a better base for correlation work (further explained in section 6). In the repeated experiment, the beam was excited by imposing an initial displacement, and then letting it vibrate freely. The acquired eigen frequency was used calibrate the mass and stiffness properties of the FEA model, which are the most fundamental factors of the structure’s dynamic behaviour. This intermediate step was helpful because it was also easy to interpret any error in the results, and feedback it for a better correlation.

2. Time domain explicit analysis results were correlated to measurement data without impact events. This step was used to validate solver settings and calibrate model parameters like material damping and stiffness of shaker connector rod connections.

3. Finally, time domain explicit analysis results were correlated to measurement data with impact events. This step was used to check the validity of the model during impact events, as a necessary step before the sensitivity study.

13

Finally, the findings were discussed and concluded as preliminary guidelines, along with possible future work.

Physical tests

The current thesis is based on measurements conducted on a rattle producing test rig (3). Although this test rig was not designed or built during the current project, it is essential to describe it in details since the FEA model is aiming to accurately recreate it, and all the measurements used for correlations were captured with this setup.

Experimental setup

The rattle test setup comprised of a cantilever beam with a rectangular cross section. At one end, it was clamped with bolts to a stand. On the other end, it was drilled such that a bolt, with a material sample glued to its head, could be inserted and fastened. The material sample was placed in front of another material sample, similarly glued to the head of another bolt which was fixed to an impact hammer. The material samples could be changed by changing the bolts, carrying different samples. The gap between the two material samples could be adjusted from 10 mm gap to a pretention state that corresponded to -5 mm initial displacement. The beam was connected to an electrodynamic silent shaker by means of a threaded cylindrical rod which was fastened to both, the beam and the shaker.

Figure 1 shows a photograph of the setup, and Table 1 displays the characteristics of the relevant parts.

Part Material Dimensions

Cantilever beam Aluminum 6061 320×20×6 mm

Material sample, beam (st) Steel ASTM-A36 15.3 mm diameter, 1.6 mm thickness Material sample, hammer (st) Steel ASTM-A36 15.4 mm diameter, 1.5 mm thickness Material sample, beam (p1) Daplen EE188HP 15.9 mm diameter, 2.6 mm thickness Material sample, hammer (p1) Daplen EE188HP 15.9 mm diameter, 2.3 mm thickness Shaker connector rod Steel grade 8.8 M6 diameter, 129 mm length

Table 1: Name, material and dimension of main parts of the experimental setup

The test rig was used in two experimental setups:

- One simplified setup for the modal analysis, which was used to measure the eigen frequency and the damping of the structure. In this setup, the electrodynamic shaker was disconnected, and the shaker connector and counter sample were removed. This way the beam could vibrate freely. In this setup, only accelerations were captured.

14

The experiments were performed in a semi-anechoic chamber at room temperature. Acceleration data was acquired from two accelerometers (Bruel & Kjaer TYPE 4524B triaxial CCLD piezoelectric accelerometer) glued to the beam. The shaker point accelerometer was fixed on the connection surface of the electrodynamic shaker. The impact point accelerometer was fixed to the cantilever beam next to the sample carrying bolt. The positions of both accelerometers are marked on Figure 1.

Figure 1: Annotated photograph of the experimental setup. Simplified model for the modal analysis (left). Main setup for the rattle testing (right).

The impact hammer, which served as a fixture for the counter sample, captured impact forces. Sound was recorded with a binaural headset (HEAD Acoustics BHS II) placed on an artificial head (HEAD Acoustics HMS IV). The data acquisition was done at a sampling frequency of 48 kHz using SQuadriga II frontend from Head Acoustics.

The electrodynamic shaker was driven by an MB Dynamics SL500VCF amplifier for which the signal was generated by MB Win2K5 vibration control system.

Modal measurements

15 Figure 2: Experimental eigenvalue analysis result

Time domain measurements

As explained in section 3, the time domain correlation was split into two phases, so it required two control signals to correlate to. However, it was possible to use the same measurement to extract both of these signals. That is because the signal with which the beam was excited is a 10 Hz sine wave that amplifies gradually from stationary state until the acceleration attains the pre-set target amplitude of the measuring sequence (the used sequence had a target acceleration of 0.05 g, see the appendix about file naming conventions). Therefore, impact between the counterparts was not attained until the beam tip displacement was equal to the gap size (the used sequence had a gap size of ~2 mm, see section 8.3). This made it possible to obtain the required excitation curves from the same measuring sequence. The acceleration data was measured near the impact point, and near the shaker connector rod, as stated previously. Impact forces were measured with an impact hammer, similarly explained in the previous section.

As mentioned in section 1.1 , the present work is a part of a doctoral thesis aiming to also investigate the connection between rattle variables and sound annoyance levels. Therefore, sound recordings were performed during these tests. These recordings were not used in the correlation part, neither in the sensitivity analysis. However, they helped to understand the behaviour of rattle events, and even gave some guidelines for the frequency range that should be analysed in the results. This use of the sound data is expounded in section 9.1.

Signal processing

16

velocity and displacement curves as a result of accumulated errors of the numerical integration and noises present in the measured data.

Figure 3: Unfiltered displacement-time graph of impact point. Showing the dominance of drift compared to the magnified part with real vibrations

The signals were analysed with the aim to choose filters and frequencies to filter. Assuming that the presence of noise is constant throughout the measurement, the pre-impact part of the signals proved to be the best to analyse, since they only contained a clear 10 Hz excitation (for the chosen measurement, 5 and 7 Hz were also measured) and the noise to be filtered. FFT was applied to the pre-impact part to analyse characteristic frequencies of the noise. The possible sources of these errors and noises were identified as:

- Ground loop of data acquisition device and signal amplifier (~50 Hz)

- Electromagnetic field of signal amplifier, electrodynamic shaker, data acquisition device, power sources (~50 Hz)

- DC offset due to not accurate zeroing of the data acquisition device (Not visible)

17 Figure 4: Acceleration-time graph for a typical measurement (top), pre-impact part analysed with

FFT (bottom)

These errors require different types of filtering depending on their nature or frequency. The objectives of the signal processing are the following to result a good baseline signal thus enabling further work:

- Eliminate drift

- Filter the identified noises which have an effect on the result

- Keep or improve the shape of the curves near impact events

- Eliminate or avoid phase errors during integration and filtering

The following subsections describe the applied or tested filters addressing the identified problems. All the signal processing was done in MATLAB R2015b.

18

Mean removal

The DC offset of the raw signal was the most straightforward to filter. The signal was simply translated along the y-axis. The extent of the translation was calculated as the mean value of the raw signal, which was based on the nature of the experiment. The cantilever beam is fixed to the base of the rig, therefore the sum displacement and the final velocity of the structure is zero during the experiment. This implies that the mean value of the acceleration should also be zero, to keep the initial zero value of both displacement and velocity.

Despite the fact that this error was small compared to the measured acceleration values, it is still essential to filter, as the signals are relatively long (between 10-30 seconds) and even these small errors can accumulate due to the double integration.

Butterworth filter

Since filtering out the identified noises does not completely eliminate the drift in the displacement-time curves due to other factors – such as the accumulating numerical error – it was necessary to filter the integrated curves. Therefore a high-pass filter was needed. Butterworth was chosen as a popular filter for similar applications. Both 2nd _{and 3}rd _{order high-pass filters were tried, and cut-off frequency was also varied.} The results with cut-off frequency of 1 Hz showed that the frequency content of the drift is not low enough to get good results. Increasing the cut-off frequency to 8 Hz eliminated the drift as hoped, but its proximity to the excitation frequency (10 Hz in this case and all presented figures) was undesirable, especially since other excitation frequencies are below this value (5 and 7 Hz).

Figure 5: Different Butterworth high-pass filter settings applied on curve introduced in Figure 3. Zoomed on the whole length of the curves to assess drift elimination.

19 Figure 6: Different Butterworth high-pass filter settings applied on curve introduced in Figure 3.

Zoomed on individual impact events to assess distortion of the different filters.

The figure above shows two big problem with this kind of filtering. First depending on the proximity of the cut-off and the excitation frequency and the order, there is a phase error in the filtered curves which has different effect on the excitation and the impact part of the curves. This results that the impact seems to slide along the sine curve of the excitation and impact seemingly happens after the minimum value of displacement which is impossible in real life.

The second problem is that as the order of the filtering is increased and the cut-off frequency is closer to the characteristic frequencies there is a distortion in the signal, mostly visible on the upper peak of the sine wave.

Both these problems make classical RF filters unsuitable for this task, therefore they were not used for the final signal processing codes. In the next sections other approaches are described.

Moving average filter

The previous subsection showed that the ideal filter for this application should have zero phase shifting and no distortion on the signal. At the same time its main purpose is to filter the low frequency drift, while knowing relatively lot about the characteristics of the real-life displacement-time curves, and therefore the desired characteristics of the filtered displacement-time curves, such as:

- Due to the periodicity and the constant zero average velocity of the signal the displacement between any time moments one period apart is zero

20

These characteristics impose that the moving average of the real-life curve should be zero, therefore the moving average of the measured signal could be used for filtering. This was done by first calculating a moving average curve with a window size of one period of excitation (that is 0.1 s for the 10 Hz measurement) and then subtract this calculated curve from the unfiltered one. The theoretical result of these operations is a curve with constant and zero moving average.

Figure 7: Removed moving average filter applied on curve introduced in Figure 3. Zoomed on individual impact events to assess distortion of the filter.

On Figure 7 it can be seen that the filter did not removed the drift perfectly, but gave similarly good results to the 3rd _{order Butterworth filter with cut-off frequency of 8} Hz. However, the distortion of it is much smaller compared to all previous methods.

Omega arithmetic

As a result of further literature study (8), and with the help of the examiner, one more non-traditional filtering method was tried. It is based on a FFT algorithm which makes possible to integrate and filter the signal in the frequency domain, and after an inverse FFT to obtain the filtered signal in time domain. This filtering method is called ‘omega arithmetic’ by Prosig the publisher of the referenced paper, where they explain this definition as:

21

components does not have effect on the neighbouring components or causes phase shift in the inversed signal.

These properties made a more sophisticated filtering possible than what was presented in previous sections. Since the proximity of the filtered and unfiltered frequency components does not cause any distortion in the final signal the drift could be filtered by eliminating components between 0 and 8 Hz. Due to the same reasons it was possible to filter frequencies freely on the middle of the frequency range and therefore eliminate the noises identified in the beginning of this section. After a few trial runs filtering out 28 to 60 Hz proved to be the most efficient while keeping the important details of the impact events.

Figure 8: Omega arithmetic filter applied on curve introduced in Figure 3. Zoomed on individual impact events to assess distortion of the filter.

The results show that omega arithmetic is able to filter all the drift while maintaining the shape of the curves. Moreover, it also makes it possible to filter the noises described in the introduction of this section. Since the part of the signal presented in Figure 8 was mainly chosen to check the drift filtering and the distortion of different filters, it is not the best choice to examine rattle events as it consists of one single and small impact. Therefore it worth to check and discuss other parts of the signal, which later will be used for correlation work, to see if the filtering distorted them.

22

Figure 9: Omega arithmetic filtered signal showing bounces during impact event.

Figure 9 shows that the filtering fulfilled all the objectives determined in the beginning of this section.

- Completely eliminated drift

- Filtered the identified noises which have an effect on the result

- Improved the shape of the curves near impact events, bringing local minimums of impact moments to a straight line parallel to time axis

- Performed the whole filtering and integration process without introducing any phase error

The final MATLAB script is attached in the appendix, including the optional use of filters introduced in earlier subsections, and comments explaining the script and its use. The processed signal serves as a good baseline for the following correlation work.

Table 2: Comparison of different filtering methods

Setbacks and Solutions

23

they were overcome, and helps understand why the initial work-plan was distorted, and another one, based on prioritization and limitations, had to be established. The main issues faced are the following five.

1. Rattle was tested in the aim of getting qualitative sound results. So, although the acceleration was recorded, they were not checked for errors because they were not used. Upon checking those and integrating to get the required displacement curve, a shape distortion in the shape of the curve, along with a big drift, were found. This resulted in a 3 month period of trying out filtering techniques to fix the shape of the signals, attempting to recreate the test to check for possible sources of error, considering re-tests, and contacting people who have possibly encountered this problem. Finally, a few online sources presented the Omega Arithmetic as a solution, which was used to overcome the signal issues. All the signal processing details are in section 5.

2. The datasheet of the accelerometer with which the measurements were recorded specified the range of frequencies with which it gives accurate results. The measured frequencies, however, exceeded the accelerometer specifications.

3. The signals had another problem. There was a factor of 10 difference in the units which was assumed to be a wrong setting in the recording software during measurement. The unit in the measurement files was m/s2 _{but, knowing the} distances, especially for the gap, chosen in the setup, the numbers appeared to must have been recorded in g. This was an uncertainty which had to be assumed and overcome. Upon trying to re-test and check Artemis, the recording software, it appeared that the default acceleration unit setting was g, which further enforced the assumption.

4. The signals used to excite the shaker were named with the gap size used during each experiment. Those sizes, however, were not as precise as was needed for the simulations. For example, a gap size of 2 mm was in the file name, but the recordings at the impact point accelerometer showed a bigger gap size. The determination of the actual gap size is detailed in section 8.3.

5. The first eigen frequency that was supposedly recorded during the modal measurements was not matched when the CAE modal analysis was performed. Trying to overcome this problem led to the recreation of a setup which did not include the shaker and which is described in section 4.1.

FEA model

24

Mesh

The finite element model was built in ANSA software with a mixture of solid, shell and beam elements, all of them being first order. The mesh size varied between 1 and 2 mm. The type and size of elements were chosen based on similar models from the field of NVH. An aim was set to only change these parameters, if they were required for the sake of correlation with physical measurements. The model totally consisted of 2742 solid elements and either 1040 or 662 shell elements, depending on the samples’ meshes (explained in section 9.1). The model was prepared for running the solutions with Abaqus solver. Below is a detailed explanation of the elements used for modelling different parts of the model.

- The beam was modelled with 2D shell elements (S4R). Using 3D elements gave no higher accuracy, while increasing the CPU time.

- The bolt and nut pair, holding the samples and the samples themselves were modelled with 3D brick (C3D8) and Penta (C3D6) elements. Since small details of these parts were neglected, the density was scaled up to match the original mass.

- The threaded rod that connected the beam to the shaker was modelled with 1D beam elements (Abaqus B31).

- Small parts with negligible contributing geometry, like accelerometers and nuts of the shaker connector rod were modelled as point masses. The equivalent mass values were based on measurements and they were validated within the first steps of the model correlation to check if they properly represented the parts.

- The connection between the shaker connector rod and the beam was modelled with rigid multi point constraints (MPCs). The mass point representing the accelerometer was also connected to the beam with rigid MPCs.

Boundary conditions

The boundary conditions in the model were defined accordingly to the experimental setup.

- The end of the beam at the support was constrained in all 6 DoF since it was rigidly clamped in the test setup.

- The nodes of the counterpart sample on the hammer side were constrained in all 3 DoF (3D elements were used for these parts). The effect of elasticity coming from the hammer support was investigated and was neglected due to its minor effect on the results.

25

direction in the FEA model. The Z direction had a prescribed displacement as the input signal that will be discussed later.

Figure 10: Finite element model used for the correlation and sensitivity analysis of rattle simulations Excitation

In order to make the simulations computationally efficient (keeping in mind that over 100 different runs were required) only small sections of the signals were used. Ideally, only the impact events would be simulated, starting the runs just before the samples come into contact. Unfortunately, this was hard to achieve since the initial conditions were unknown and calculating them was not straightforward. The displacement, velocity and acceleration values of the beam’s nodes are dependent on both the excitation and the free vibration of the beam. The most convenient way to calculate these values is to run a preceding simulation, which excites the beam from stationary state and then capture the nodal values to be used as initial conditions. Since this project had more computational capacity than time, the simulations simply included this preceding part in every run. The time that was spent solving the unforeseen signal processing problem had to be deducted from this task, where learning how to use the restart files was required. However an option was left to start the simulations from a restart file just before the impact, in case further batch runs are required and CPU time becomes a bottleneck.

26

Although the signals to be used as excitation were already processed, they could not be applied directly to the FEA model as a boundary condition, since that would result in a jerk. That is because the section of the signal gives either an initial velocity or displacement depending on the starting point along the curve. A solution could be to use initial conditions for either displacement or velocities, but the problem with this solution was already discussed in the first paragraph. To overcome the problem of initial conditions, a “start curve” was introduced and attached to the actual signal. This start curve had to be generated so it starts with zero displacement, velocity and acceleration and connects to the signal curve with C2 continuity to avoid any jerk at the connection point (Figure 11). In order to avoid any undesired vibration of the beam, the curve had to have a certain minimum curvature due to the inertia of the beam. A sudden change in displacement would bring the beam into vibration even if the curve itself is continuous.

Mathematically the simplest function that could easily fulfil all these boundary conditions was a sine function. The fitting was ensured by solving the following system of equations: { 𝑎 ∗ 𝑠𝑖𝑛(𝑐 ∗ 𝑡₀+ 𝑐 ∗ 𝑏) + 𝑎 = 𝑠𝑖𝑔𝑛𝑎𝑙(𝑡₀) 𝑑 𝑑𝑡(𝑎 ∗ 𝑠𝑖𝑛(𝑐 ∗ 𝑥 + 𝑐 ∗ 𝑏) + 𝑎) = 𝑑 𝑑𝑡𝑠𝑖𝑔𝑛𝑎𝑙(𝑡0) 𝑑2 𝑑𝑡2(𝑎 ∗ 𝑠𝑖𝑛(𝑐 ∗ 𝑥 + 𝑐 ∗ 𝑏) + 𝑎) = 𝑑2 𝑑𝑡2𝑠𝑖𝑔𝑛𝑎𝑙(𝑡0)

27 Figure 11: Position-time graph for test data showing the measured impact point displacement, shaker point displacement (used for excitation), initial start curve and scaled start curve. The annotated sections mark the simulation phases, not the measurement

The curvature was adjusted by scaling this start curve along the time axis (Figure 11) so that at the connection point it is unaffected to keep the C2 continuity and increasing the scaling factor with the square of the distance from the connection point.

Material and property cards

The finite element model consisted of three types of property cards. One shell property for the cantilever beam, describing its thickness of 6 [mm], solid properties with the corresponding material cards, and the beam property (B31) for the shaker connector rod. For the beam elements of the rod, the material and cross sectional properties were also included in the property as: 7.85e-9 [t/mm3_{] density, 2.5 [mm]} radius of circular cross section, 210 000 [MPa] elasticity modulus and 80 796 [MPa] shear modulus.

28

For the sensitivity analysis, the sample material was changed to polypropylene from “Borealis Daplen™ EE188HP” which has a modulus of 1850 [MPa], 0.42 [-] Poisson’s ratio, and 1.02e-9 [t/mm3_{] density.}

Contact modelling

To model contact, four settings are important to define: the contact type, the contact surfaces, the mechanical constraint, and the interaction.

Two types of contact exist in ABAQUS. They can be used individually or together. Selecting “contact pair” requires the user to manually select the surfaces that are expected to come in contact. On the other hand, the “general contact” formulation allows the definition of an area of, or the whole, model and locates any contact within that region as the simulation runs. In this thesis’s FEA model, the two surfaces that come in contact are known and clear, so the contact pair type was picked for the baseline simulations which had St-St, Shell, Pen, CP, el-el, no material damping, 0.03 contact damping, 0.2 friction settings (see Abbreviations for the meaning and the following sections for the values).

Next, the two surfaces that come in contact, denoted as “master” and “slave” in CAE, have to be defined. The master can only be chosen to be element-based, but the slave can be either element or node-based. Using either of the two is dependent on the problem under investigation (9). Both surfaces were chosen to be element-based for the baseline simulations as it is the default setting in ABAQUS.

The way ABAQUS/Explicit calculates and measures the contact is chosen in the “mechanical constraints” field. One of the two methods can be selected. The hard kinematic contact does not allow any penetrations by using a kinematic predictor/corrector contact algorithm. The penalty contact method searches for penetrations and applies forces, calculated as a function of penetrations, to the slave nodes to eliminate penetration. ABAQUS/Explicit allows the definition of either of the two constraints with the contact pair type, but with general contact, only the penalty enforcement is possible. For the first simulations, the kinematic method, which is the default with contact pair definition, was used.

Finally, a predefined “interaction” must be specified for the contact, where properties such as friction (using the Coulomb model) and contact damping can be edited. For the first simulations, the default ABAQUS/Explicit values were kept. Therefore, there was no friction between the surfaces, and the contact damping coefficient, which is given as a fraction of the critical damping coefficient, was 0.03. A relevant note is that contact damping is not available for hard kinematic contact in ABAQUS/Explicit.

29 Solver Settings

As written in the introduction of the section, two solvers were used during the thesis. Abaqus/Standard was used for the validation step with modal analysis, using the Lanczos method and extracting the first 10 modes. Similarly to the physical tests the BCs at the shaker and the shaker connector rod itself were disabled for modal analysis. Modal and nodal outputs were requested to obtain eigen frequency values and also the shapes of different modes.

The dynamic simulations were run in Abaqus/Explicit solver. The decision of using explicit solver over implicit was based on preliminary assumptions about the impact velocity, that was higher than the suggested limit for implicit solvers (10). The stable time increments were chosen automatically based on the critical element in each increment. No mass scaling was used since the high velocities would result in poor accuracy. All other settings were left default in the solver.

The explicit runs were split into two time steps in order to help eliminate any undesired vibration in the beam due to the excitation. The problem was already explained in 7.3. That part focused on the excitation signal, while this subsection aims to help the problem with introducing material damping in the first time step: that is the start curve and excitation sections marked on Figure 11. A mass proportional damping was applied to the beam material with a coefficient value of 43.35, which corresponds to 0.5 times the critical damping coefficient at the first eigen frequency of the beam (31.27 Hz). The calculations for damping coefficients are explained in 9.1.2. The contact interfaces were disabled for the first time step to reduce the CPU time of the runs.

The second time step contained the impact section marked on Figure 11 (starting at 0.2035 s in simulation time and ending at 0.26 s). In this step, the material damping of the beam material was adjusted to the value studied in the given run. At the same time, contact interfaces were enabled.

The material damping in the first time step eliminated the vibrations in the beam. The variation of material damping is not available through any pre-processors, only through input files as a keyword only option. It required the definition of a tabular damping coefficient input that was dependent on a field variable, changed at the beginning of the second time step. Since the use of this function is not straightforward, the simplified input file of the runs is attached in the Appendix. During the explicit runs different outputs were requested described in Table 3. The outputs were written in either history or field outputs depending on their nature. Energy and contact outputs were requested but not used in this paper. Energy outputs aimed to investigate the dissipated energy during the impacts but proved to be insignificant as a rattle criteria. The contact outputs aimed to analyse the contact outputs, but the reaction force measured on the boundary conditions proved to be a more similar solution to the experimental setup.

30

Table 3: Output requirements by categories for both time steps

Model validation

The method of the model validation used is summarized here and explained in depth in the following subsections.

1. Perform a modal analysis of the beam to ensure the CAE model is a good representation of the experimental setup.

2. Perform a dynamic analysis of the beam without impact, to study its response to the shaker excitation.

3. Perform a dynamic analysis of the beam with a single impact to minimize the amount of variables while still having a well correlated contact model.

Modal correlation

As a first step of ensuring the validity of the model before the planned impact correlations, the mass and dimensions of each part were matched with the FEA model, and material properties were set according to the material datasheets. Then the Young’s modulus and Poisson's ratio of the beam were modified in reasonable range to improve the correlation with the empirical data. A modal analysis was finally run in ABAQUS/Standard and gave a first eigen frequency of 31.334 Hz, as shown in Figure 12. Compared with the measurement result (Figure 2), the error is 0.2%, which meant a valid model was obtained, and the second part of the correlation could start.

Nodal output Time step #1 Time step #2

Acceleration (A) Yes Yes

Velocity (V) Yes Yes

Displacement (U) Yes Yes

Reaction force (RF) Yes Yes

Element output Time step #1 Time step #2

Strain for non-linear analysis (LE) Yes Yes

All stress components (S) Yes Yes

Energy output Time step #1 Time step #2

"Artificial" strain energy (ALLAE) Yes Yes

Total strain energy (ALLIE) Yes Yes

Kinetic energy (ALLKE) Yes Yes

Total energy balance (ETOTAL) Yes Yes

Contact output Time step #1 Time step #2

Contact normal force (CFORCE) No Yes

Contact pressure (CSTRESS) No Yes

Force components from contact pressure(CFN) No Yes

Force magnitude from contact pressure (CFNM) No Yes

Force components from contact pressure and friction (CFT) No Yes

31 Figure 12: Eigenfrequency vs. time curve showing the first eigenfrequency of the CAE model

Time domain correlation before impact event

32

Figure 13: Curves showing the variation of the distance between the beam and the MPC control node, in comparison with the measurement curve

Then, the end point of the start curve (explained in 7.3) was iterated to smoothen the observed breakpoints. The best iteration found was when the end point was chosen towards the first peak of the input signal (Figure 14).

Figure 14: Curves showing the effect of the variation of the end point of the start curve (used on the input signal)

33

the right gap size, a process that is explained in section 8.3. The final correlation before impact is in Figure 15. As visible there is a time shift between the measurement and simulation results. This shift was not investigated further since the displacement curves showed a good correlation near the impact events presented in the following sections. Moreover the difference is decreasing towards the second local minimum thanks to the mass proportional damping described in section 7.6. As visible in Figure 17 the time shift during the impact event decreased to ~0,001 s.

Figure 15: Displacement vs. time curve showing the correlation between the measurement and CAE simulation results before impact

Gap Size Determination

One requirement for the correlation is to know all the parameters of the physical test setup, since the aim of the correlation work is to simulate the same setup and obtain the results of the experiments. As written in section 6, although the gap size was documented during the measurements, it did not have the desired precision. Since knowing the accurate gap size was crucial for the next section, since it has a significant effect on the results, a method was needed to obtain it from the existing data.

34

measured accelerations, was not measured at the impact point, therefore the displacement of this measured point (accelerometer displacement) is less than the gap size (or impact point displacement) when the sample is touching the counter sample. The problem is illustrated on Figure 16.

Figure 16: Illustrating the displacement difference between accelerometer and impact point displacements. Green model represents initial position of the beam, red represents the bended position of the beam. Shells are displayed with their thickness, and displacements are scaled up 10 times to better illustrate the problem.

The problem was solved using the simulations without impacts already introduced in the previous subsection. Since in the simulation results any node’s displacement can be measured it is perfect to convert the accelerometer displacements to impact point displacements. Therefore, the local maximum of measured accelerometer displacement just before the first impact could be converted to impact point displacement that should be equal to the gap size. Keeping in mind that the amplitude is increased by ~0.01 mm per wavelength and the local maximum can be measured just half wavelength before the first impact, due to the symmetric signal, this method theoretically gives a 0.005 mm accuracy. Although this accuracy is probably reduced by undiscovered errors and errors from signal processing, it is still assumed to be more accurate than the documented gap size precision.

35 Time domain correlation at impact

Figure 17: Displacement vs. time curve showing the correlation between the measurement and the chosen baseline CAE simulation results at impact

The time domain correlation at impact did not involve any editing more than what was explained in sections 8.1 and 8.2. The shapes of the measurement and the CAE simulation result curves displayed a good correlation that completed the third objective of the project, and allowed the start of the parameter variation and sensitivity analysis.

Sensitivity Analysis

After the first part of the project was done by establishing the desired correlations, a sensitivity analysis was performed. The use of the word “parameter” in what follows is restricted to the independent specifications that were varied in ABAQUS. The word “criteria” will refer to the dependent variables whose responses to parameter changes were analysed. The aim was to understand the effect of material and model parameters on the criteria that characterize rattle sounds.

The parameters investigated in the sensitivity analysis can be categorized as either 'structural and material properties' or 'CAE-specific options'. The CAE-specific options are explained below using ABAQUS terminology (9), but equivalent terms can be found in similar solvers. The parameters and their explanations are given in section 9.1.

The criteria studied play an influential role in the generation of sound according to (11). They are presented and explained in section 9.2.

36

steel material samples (St-St), two polypropylene plastic material samples (P-P) and a polypropylene and a steel material sample (P-St) was simulated in ABAQUS. Further details of the material properties are available in Table 1 and section 7.4. The material appears in the comparison tables and discussions although it is not a parameter.

Nominated parameters and the ranges

Element type for material samples

The difference between the use of solid and shell elements in contact interactions was of interest. That is why either the whole of both material samples was meshed with 3D brick and penta elements (C3D8, C3D6), referred to as "solid" in the following sections, or the surface was meshed with 2D shell elements (S4) while the rest was in solids. Figure 18 shows the two different meshes. In the result tables (Table 6 onward), this parameter is referred to as “Sample Mesh” for convenience.

Figure 18: Meshing the material samples with (left) solid elements, and (right) shell elements only on the surface of the samples

Material Damping

Rayleigh damping was defined for the material card of the beam. It is a superposition of a mass and a stiffness damping as shown in the equation:

[C] = α[M] + β[K] (1)

where the stiffness proportional coefficient damps out systems with high frequencies, and the mass proportional coefficient damps out systems with low frequencies, as can be seen in the equation below.

ζ = α 2ω+

β

2ω (2)

37

The mass-proportional damping coefficient (α) was calculated from the damping ratio recorded during the eigenvalue test described in section 4.2. Using the logarithmic decrement δ to get the damping ratio 𝜁, α is calculated from equation 2 to a value of 0.8418 with β set to zero.

𝛿 =1 𝑛ln ( 𝑥(𝑡) 𝑥(𝑡 + 𝑛𝑇)) = 1 8ln ( 0.0007954 0.0007142) = 0.01346 𝜁 = 𝛿 2𝜋= 0.0021423

The stiffness proportional damping (β) was calculated as 5% of the target time step, based on a suggestion in the literature (12), giving a value of 6.627e-9.

The remaining parameters, “Contact type definition”, “Contact mechanical constraint”, “Contact damping”, “Contact surfaces” and “Friction”, are explained in more detail in section 7.5, and briefly stated below with their abbreviations.

Contact type definition

Simulation with both contact type definitions, “contact pair” and “general contact”, were run to study their effects. In the result tables (Table 6 onward), this parameter is referred to as “Contact Definition” for convenience.

Contact mechanical constraint

The two contact methods are the kinematic and the penalty. In the result tables (Table 6 onward), this parameter is referred to as “Contact Mechanics” for convenience.

Contact Damping

Contact damping was given as a fraction of the critical damping coefficient and was varied from the default value of 0.03 to the values 0.05, 0.1, and 0.2.

Contact Surfaces

The two variations of contact surfaces investigated in this study are denoted “el-el” and “node-el”. El-el refers to the case where both sample surfaces are defined as element-based, and node-el refers to when the slave surface is defined as node-based.

Friction

38

Evaluating the Fundamental Rattle Criteria

The criteria studied in the sensitivity analysis are: impact velocity, height of the biggest bounce upon impact, duration of the bounce (and frequency), peak impact force, and the duration of impact.

After running the simulations with different parameter combinations, the results were post-processed in META. The criteria of interest were extracted from the relevant plots at the relevant nodes, shown in Figure 19.

39

The impact velocity, height of biggest bounce, and the duration (and frequency) of the biggest bounce were extracted from the plots displaying the displacement at point B vs. time. The reason for choosing this point is that the accelerometer was placed at this point in the experiment, therefore this point of the FEA model has to be investigated to make the findings comparable to the experiment. The impact force and duration of impact were extracted from the plots displaying the force at point A vs. time. This point was also chosen based on the experiment since the acoustic hammer was measuring forces at this point, and the force data was used to determine duration of impact for the experimental data.

Figure 20: Extracted points marked on plotted displacement in Z direction at point B (blue) and hammer force in Z direction at point A (orange). The curves belong to v64_stst_1struct_p248_expl026_nodrift14693_scaled12539_shell_pen_gen_elel_nom atdamp_contdampoc999_nofric_nosslide.odb

The points are extracted and moved to an excel sheet where the values for the criteria are calculated, and the all data is organized and coded such that every parameter and rattle criterion can be filtered and rearranged according to alphabetical order and/or increasing/decreasing order. This kind of arrangement helps analyse the results and deduce the desired conclusions which are presented in section 10.

40

41

Results of the sensitivity analysis and guidelines

42

Table 4: CAE simulations run for the sensitivity analysis

Table 4 displays all the simulation runs in ABAQUS/Explicit for the sensitivity analysis. Those results are analysed and grouped separately in this section to explain the effect of different parameters on the criteria.

Although different material pairs were simulated, conclusions should not be deduced from the differences between these results as the input signals for the shaker point were also different. The inputs signals were based on the measurements and unfortunately the gap size between these measurements varied. Since the signal applied on the electrodynamic shaker was an amplitude sweep, different gap sizes resulted in a different input amplitude at the moment of impact events.

43 The Effect of Contact Surface Definition

Table 5: Results showing the effect of the el-el and node-el contact surface definition on the studied criteria

All models in Table 5, have a CP, Kin, no material damping, default (0.03) contact damping, no friction configuration. The aim is to study the difference between using the node-el and the el-el contact surface parameter. It is visible from the results that there is a significant difference when comparing the values for the models where the sample surfaces are meshed as shell, but barely any difference for samples meshed as solids. The height of the biggest bounce, duration of the biggest bounce, ratio of the height of bounce to impact velocity and peak impact force are higher for the P-P and P-P-St pairs where the node-el surfaces are used. The duration of impact for those pairs is lower with node-el. On the other hand, for the St-St material pair, the opposite is true: the height of the biggest bounce, duration of the biggest bounce, ratio of height of bounce to impact velocity and peak impact force are lower where the node-el surfaces are used, and the duration of impact is higher.

44

The Effect of Friction

In this subsection, the value for contact damping is the default ABAQUS value which is 0.03.

St-St material pair

Table 6: Results showing the effect of friction on the criteria for the St-St material pair

The use of shell or solid elements to model the St-St samples makes a significant difference on how friction affects the height of the biggest bounce, the duration of the bounce, and the ratio of the height of bounce to the impact velocity. As the results for the material pair St-St in Table 6 show, the use of friction:

- Decreases the height of the biggest bounce and the ratio of the height of bounce to the impact velocity when the samples are modelled as shell elements and increases them when the samples are modelled as solid elements. Looking at the numbers and comparing Shell-CP-Pen with Solid-CP-Pen models, friction has a bigger effect when the elements are solid.

- Decreases the duration of the bounce when the samples are modelled as shell elements and increases it when the samples are modelled as solid elements. However, the effect on CP-Pen models is barely visible.

- Increases the peak impact force irrespective of the element type

- Decreases the duration of impact, irrespective of the element type. However, when the Kinematic (Kin) method is in action, friction has a negligible effect on the duration compared with no-friction models.

One more find that can be extracted from this table is that general contact definition has a bigger effect than the contact pair definition. The numbers show that the percent difference between models with and without friction is the biggest for the Solid-Gen models compared with the remaining comparable models.

Sample mesh

Contact definition

Contact

mechanics Damping Friction Impact velocity [mm/s] Percent Difference [%] Height of biggest bounce [mm] Percent Difference [%] Duration of biggest bounce [s] Percent Difference [%] Ratio of height of bounce to impact velocity Percent Difference [%] Peak impact force [N] Percent Difference [%] Duration of impact [s] Percent Difference [%] Shell CP Kin Mdamp Friction -91.07 0.0593 0.00477 0.000651 104.42 0.00111 Shell CP Kin Mdamp No friction -91.07 0.0734 0.00499 0.000806 92.30 0.00113 Shell CP Kin No Friction -90.83 0.0593 0.00476 0.000653 107.67 0.00113 Shell CP Kin No No friction -90.83 0.0717 0.00495 0.000789 94.51 0.00112 Shell CP Pen Mdamp Friction -91.07 0.0529 0.00469 0.000581 63.37 0.00120 Shell CP Pen Mdamp No friction -91.07 0.0564 0.00478 0.000619 59.21 0.00142 Shell CP Pen No Friction -90.83 0.0533 0.00472 0.000587 65.95 0.00118 Shell CP Pen No No friction -90.83 0.0575 0.00479 0.000633 61.50 0.00142 Solid CP Pen No Friction -90.89 0.0463 0.00493 0.000509 87.15 0.00125 Solid CP Pen No No friction -90.89 0.0376 0.00490 0.000414 78.51 0.00151 Solid CP Pen Sdamp Friction -90.89 0.0462 0.00494 0.000508 87.11 0.00125 Solid CP Pen Sdamp No friction -90.89 0.0381 0.00492 0.000419 78.44 0.00151 Solid CP Pen Mdamp Friction -91.13 0.0525 0.00506 0.000576 91.21 0.00123 Solid CP Pen Mdamp No friction -91.13 0.0357 0.00492 0.000392 78.88 0.00152 Solid Gen Pen No Friction -90.89 0.0630 0.00524 0.000693 106.07 0.00120 Solid Gen Pen No No friction -90.89 0.0336 0.00494 0.000370 91.32 0.00153

45

P-P material pair

Table 7: Results showing the effect of friction on the criteria for the P-P material pair

As the P-P results in Table 7 show, the use of friction:

- Decreases the height of the biggest bounce and the ratio of the height of bounce to the impact velocity. Looking at the numbers, friction has a bigger effect when the elements are solid.

- Decreases the duration of the biggest bounce, much more significantly for solid samples. Friction has a negligible effect on the duration of the biggest bounce for the Shell-CP-Pen combination.

- Increases the peak impact force.

- Has an insignificant effect on the duration of impact for all models except for the Shell-CP-Pen combination where friction causes a decrease in value.

The effect of the general contact definition is not noticeable for P-P as in St-St results. For P-P samples, the more noticeable results are those for the Solid-CP-Pen, which show the highest percent difference (between models with and without friction) for the height of the biggest bounce, the duration of the bounce, and the ratio of the height of bounce to the impact velocity.

Sample mesh

Contact definition

Contact

mechanics Damping Friction Impact velocity [mm/s] Percent Difference [%] Height of biggest bounce [mm] Percent Difference [%] Duration of biggest bounce [s] Percent Difference [%] Ratio of height of bounce to impact velocity Percent Difference [%] Peak impact force [N] Percent Difference [%] Duration of impact [s] Percent Difference [%] Shell CP Kin Mdamp Friction -75.59 0.0395 0.00529 0.000523 22.66 0.00171 Shell CP Kin Mdamp No friction -75.59 0.0452 0.00539 0.000598 20.82 0.00169 Shell CP Kin No Friction -75.40 0.0395 0.00526 0.000524 22.02 0.00171 Shell CP Kin No No friction -75.40 0.0450 0.00540 0.000597 20.63 0.00169 Shell CP Pen Mdamp Friction -75.59 0.0318 0.00519 0.000421 21.50 0.00175 Shell CP Pen Mdamp No friction -75.59 0.0332 0.00516 0.000439 20.53 0.00181 Shell CP Pen No Friction -75.40 0.0317 0.00517 0.000420 21.36 0.00175 Shell CP Pen No No friction -75.40 0.0329 0.00515 0.000436 20.43 0.00180 Solid CP Kin No Friction -75.41 0.0441 0.00540 0.000585 31.42 0.00163 Solid CP Kin No No friction -75.41 0.0557 0.00562 0.000739 28.39 0.00163 Solid CP Pen Mdamp Friction -75.59 0.0322 0.00518 0.000426 25.53 0.00168 Solid CP Pen Mdamp No friction -75.59 0.0461 0.00545 0.000610 24.04 0.00167 Solid CP Pen No Friction -75.41 0.0319 0.00517 0.000423 25.53 0.00169 Solid CP Pen No No friction -75.41 0.0459 0.00545 0.000609 23.98 0.00167 Solid CP Pen Sdamp Friction -75.43 0.0320 0.00518 0.000424 25.51 0.00168 Solid CP Pen Sdamp No friction -75.43 0.0459 0.00545 0.000609 23.95 0.00167 Solid Gen Pen No Friction -75.41 0.0444 0.00542 0.000589 30.59 0.00162 Solid Gen Pen No No friction -75.41 0.0567 0.00563 0.000752 27.41 0.00161

46

P-St material pair

Table 8: Results showing the effect of friction on the criteria for the P-St material pair

As the P-St results in Table 8 show, the use of friction:

- Decreases the height of the biggest bounce, the duration of the biggest bounce, and the ratio of the height of bounce to the impact velocity for all models except the Solid-CP-Kin and Solid-Gen-Pen configurations. For those, friction significantly increases the 3 above mentioned criteria. Comparing with the remaining models, it seems that the combination of solid with CP-Kin is what makes this difference. The only uncertainty is whether running Shell-Gen would give the same result, or if also for that, the solid meshing is what makes the difference.

- Increases the peak impact force.

- Increases the duration of impact for shell samples, and decreases it for

solid samples. For Solid-CP-Pen, however, the change is insignificant.

The effect of the general contact definition, and the kinematic contact mechanics is visible in this case and discussed for more visible results in section 10.3.

Sample mesh

Contact definition

Contact

mechanics Damping Friction Impact velocity [mm/s] Percent Difference [%] Height of biggest bounce [mm] Percent Difference [%] Duration of biggest bounce [s] Percent Difference [%] Ratio of height of bounce to impact velocity Percent Difference [%] Peak impact force [N] Percent Difference [%] Duration of impact [s] Percent Difference [%] Shell CP Kin Mdamp Friction -114.25 0.0514 0.00523 0.000450 58.56 0.00166 Shell CP Kin Mdamp No friction -114.25 0.0833 0.00554 0.000729 53.05 0.00149 Shell CP Kin No Friction -114.06 0.0524 0.00522 0.000459 54.96 0.00163 Shell CP Kin No No friction -114.06 0.0815 0.00551 0.000715 50.38 0.00151 Shell CP Pen Mdamp Friction -114.25 0.0647 0.00529 0.000566 47.64 0.00171 Shell CP Pen Mdamp No friction -114.25 0.0780 0.00547 0.000683 46.52 0.00167 Shell CP Pen No Friction -114.06 0.0630 0.00528 0.000552 47.42 0.00170 Shell CP Pen No No friction -114.06 0.0773 0.00547 0.000678 45.98 0.00167 Solid CP Kin No Friction -114.06 0.0725 0.00552 0.000636 70.27 0.00139 Solid CP Kin No No friction -114.06 0.0587 0.00534 0.000515 62.14 0.00155 Solid CP Pen Mdamp Friction -114.25 0.0696 0.00544 0.000609 53.51 0.00157 Solid CP Pen Mdamp No friction -114.25 0.0719 0.00551 0.000629 51.04 0.00158 Solid CP Pen No Friction -114.06 0.0696 0.00543 0.000610 54.01 0.00157 Solid CP Pen No No friction -114.06 0.0721 0.00549 0.000632 52.03 0.00158 Solid CP Pen Sdamp Friction -114.08 0.0698 0.00544 0.000612 53.92 0.00158 Solid CP Pen Sdamp No friction -114.08 0.0722 0.00547 0.000633 51.85 0.00158 Solid Gen Pen No Friction -114.06 0.0736 0.00550 0.000645 69.38 0.00135 Solid Gen Pen No No friction -114.06 0.0533 0.00522 0.000467 62.65 0.00156