Projects in Applied Mechanics 2017

Research report 2017:04 Projects in Applied Mechanics 2017 EDITED BY VALERY CHERNORAY AND HÅKAN JOHANSSON Department of Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Göteborg, Sweden 2017 Research report 2017:03 Projects in Applied Mechanics 2017 Edited by VALERY CHERNORAY AND HÅKAN J OHANSSON Department of Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Göteborg, Sweden, 2017 Projects in Applied Mechanics 2017 VALERY CHERNORAY AND HÅKAN JOHANSSON © T H E R E S P E C T I V E A U T H O R S , 2017 Research report 2017:04 ISSN 1652-8549 Department of Applied Mechanics Chalmers University of Technology SE-412 96 Göteborg Sweden Telephone +46 (0)31 772 1000 Preface The course TME131 Project in Applied Mechanics is a mandatory course within the Applied Mechanics Masters programme at Chalmers. The course was carried out during spring semester 2017. In total 7 projects were carried out, and 4 of these are collected here. Implementation of an engineering method for wind turbine aerodynamics Students: Zuher A BDEL M ALLAK Kári B REKASON Jakob E KERMO David E NGLÉN Oliver S JÖGREN Madhavan VASUDEVAN Supervisor: Dr. Hamidreza A BEDI TME131 Project in Applied Mechanics Group 3 D EPARTMENT OF A PPLIED M ECHANICS C HALMERS U NIVERSITY OF T ECHNOLOGY G ÖTEBORG , S WEDEN M AY 31, 2017 i Abstract Implementation of an engineering method for wind turbine aerodynamics was the objective of the project work. The Blade Element Momentum (BEM) method was implemented to calculate the aerodynamic loads and the results were validated with given values from the vortex method. The computational benefits of BEM over vortex method was realized and deviation of results at every stage were quantified. The blade of the wind turbine is discretized using cosine, half-cosine and equidistant mesh types. A convergence study was made for the three mesh types, where cosine and half cosine mesh gave the better convergence to consider for further studies on account of better resolution of higher gradient of forces near the tip of the blade. The classical BEM method was first implemented with no yaw misalignment involved and then compared with vortex method and BEM code of National Technical University of Athens (NTUA) for the NREL 5MW wind turbine. A series of yaw misalignments were introduced and the impact of sign and magnitude of the same on the power and thrust generated were studied and compared with the vortex method. The pitch angles were found crucial in higher velocities as for the yaw misalignment case, the pitch angles used for no yaw misalignment case were over-pitching the blades, producing less than the rated power. The normal and tangential forces were computed across the section of the blade for different velocities with and without yaw misalignment and were compared to the results from vortex method. In attempt to make BEM model replicate reality more, loss factors were associated to the calculations and it was striking to notice that the hub and the tip loss factor could not contribute commendably to help BEM results match the results of vortex method. The load dependence on the azimuth angle in yawed inflow condition opened the possibility of extending the current work to dynamic structural analysis and the impact of blade pitch angles were studied. iii Contents 1 2 3 4 5 Introduction 1.1 Background 1.2 Purpose . . . 1.3 Goals . . . . 1.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 2 2 Theory 2.1 The classical BEM theory . . 2.2 Corrections to classical BEM 2.2.1 Tip and hub losses . 2.2.2 Yawed inflow . . . . 2.3 Vortex method . . . . . . . . 2.4 Previous studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 7 7 8 8 9 Method 3.1 Blade discretization . . . 3.2 Classical BEM . . . . . . 3.3 BEM with yawed inflow 3.4 Power calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 10 11 12 12 . . . . . 13 13 14 18 19 22 . . . . . . . . Results and Discussions 4.1 Convergence study of mesh . . . . . . . . . . . . . 4.2 Comparison between BEM and the Vortex method 4.3 Loss correction factors . . . . . . . . . . . . . . . . 4.4 Effects of yawed inflow . . . . . . . . . . . . . . . . 4.5 Altering blade pitch . . . . . . . . . . . . . . . . . . Conclusion References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 25 v List of Figures 1 2 Aerofoil geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Overview of yawed inflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix x 2.1 2.2 Velocity triangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Forces acting on the blade element . . . . . . . . . . . . . . . . . . . . . . . 4 5 3.1 Overview of mesh types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 Power generated without yaw . . . . . . . . . . . . . . Thrust generated without yaw . . . . . . . . . . . . . Power for negative and positive yaw angles . . . . . . Thrust for negative and positive yaw angles . . . . . . Power generated with yaw . . . . . . . . . . . . . . . . Thrust generated with yaw . . . . . . . . . . . . . . . Normal and tangential force distribution without yaw Normal and tangential force distribution with yaw . . Loss factor comparison for power . . . . . . . . . . . . Zoom of loss factor comparison . . . . . . . . . . . . . Loss factor comparison for thrust . . . . . . . . . . . . P Pmax compared for different yaw angles . . . . . . . . Azimuthal variance in torque with yaw . . . . . . . . Azimuthal variance in thrust with yaw . . . . . . . . Effects of altered pitch angles . . . . . . . . . . . . . . 14 14 15 15 16 16 16 17 18 18 19 20 20 21 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables 3.1 Nominal operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.1 4.2 Power from BEM compared to vortex method . . . . . . . . . . . . . . . . Thrust from BEM compared to vortex method . . . . . . . . . . . . . . . . 14 15 vii Nomenclature BEM theory α Angle of attack β Inflow angle χ Skew angle γ Yaw angle λr Local tip speed ratio, λr = Ω Rotational speed of the rotor ω Rotational speed of the wake ψ Azimuth angle ρ Density of the local atmosphere σ′ Local solidity factor, σ ′ = θp Pitch angle θt Elemental twist angle a Axial induction factor a′ Tangential induction factor Ωr V Bc 2πr askew Axial induction factor corrected for yawed inflow B Number of blades c Chord length CD Drag coefficient CL Lift coefficient dD Elemental drag force dFθ Elemental lift/force tangential to blade disk dFx Elemental thrust/force normal to blade disk dL Elemental lift force dr Element span dT Elemental torque F Combined loss coefficient p Static pressure R Rotor radius r Local blade radius viii Rhub Hub radius V∞ Inflow velocity/Wind speed W Local/relative flow velocity ix Terms and definitions There are various terms used in this report that are common in the field of aerofoil studies. These terms will be explained in this section. Figure 1 shows a cross section of a wind turbine blade in motion. The chord line is the straight line that connects the leading and trailing edge. The flow angle at the blade is indicated by β. The angle of attack is the angle between the relative wind velocity and the chord line where it is defined as: α= π − β − (θt + θp ) 2 (1) To increase efficiency of the wind turbine the blade is twisted along its span and θt describes the local twist angle. For wind speeds higher than the turbine’s rated wind speed, in order to generate constant power, the local angle of attack must be reduced along the blade. This can be done by the mechanism of pitch angle regulation during the operation of wind turbines at such conditions. The pitch, θp is a global variable as the whole blade always has the same pitch while twist can be different depending on spanwise location. θ V∞ rΩ θt + θp W α β ch or d x F IGURE 1: An overview of terms and angles used for aerofoils. x The amount of lift and drag an aerofoil produces can be decomposed into the lift- and drag coefficients CL and CD . Within a certain Reynolds number range, CL and CD are interpreted as empirical constants that quantify the drag and the lift forces as functions of area and dynamic pressure for any given angle of attack. F IGURE 2: Overview of yawed inflow. An explanation of the yaw-, skewand azimuthal angle. An overhead view is shown on the left while the right side depicts a front view of the rotor disk. Figure 2 shows the case of yawed inflow and the yaw angle is marked as γ. After interacting with the wind turbine the outflow will have a slightly different angle than γ. This is the skew angle χ which is a function of the yaw angle γ. The azimuth angle ψ defines the position of the blade, where the start position is located on the horizontal line. F J 1 1 Introduction In this section the background for this project is presented as well as the purpose and the main goals of this project. Finally the limitations are listed. 1.1 Background As a consequence of growing concerns for climate change today, more and more emphasis is put on renewable energy sources when generating power for modern society. One of the more developed and used techniques belonging to this subgroup of energy sources is wind power, of which the usage has increased exponentially over the last decade. This trend of increase in popularity demonstrates the need for analysis tools that allow for optimized designs of wind turbine blades. Modern analysis tools are all based on computational efforts. In the domain of fluid mechanics, the most popular method is Computational Fluid Dynamics, in which one solves the governing equations of conservation of quantities of interest iteratively on a discretized domain. The proven disadvantage of this method is that it is very expensive in terms of computational times and power for industrial applications. Easier and more efficient implementation is achieved with the vortex method, which is based on incompressible, inviscid and irrotational flow theory. Using the latter method, computational duration can be decreased from a few months to an hour. In the domain of fluid mechanics, simplifications and assumptions are often useful when aiming to construct an easily and efficiently implemented model that display a sufficient level of accuracy. By introducing a number of simplifying assumptions describing the behaviour of the flow, the Blade Element Momentum (BEM) method can be derived. This is a highly efficient method, in which computational times and efforts can be more or less eliminated. The basic BEM method does not take into the account the possibility of the inflow being yawed, thus not aligned with the rotor axis. As this often is the case despite active control on modern wind turbines, the BEM method needs to be extended to take this effect into account. 1.2 Purpose The purpose of this project is to validate BEM analysis results for a given wind turbine, with respect to corresponding analysis results achieved using the vortex method. Furthermore, the BEM method is to be extended to also take into account yawed inflow, to provide a fast and cheap alternative to the more extensive aforementioned methods. Chapter 1. Introduction 2 1.3 Goals The goal of this project is thereby divided into three subtasks on order to meet the purpose stated above. These are the following; • Implement the classical BEM method in MATLAB in order to predict loads such as power and thrust as well as distribution of tangential and axial forces along the blades • Extending the classical BEM code to take into account yaw misalignment • Comparing the results to vortex method for validation and benchmarking 1.4 Limitations This project will be restricted to only analyzing the NREL offshore 5-MW baseline wind turbine blade. 3 2 Theory In the following chapter the theory for the BEM method is presented along with the corrections made to extend this method. Lastly, the vortex method and previous studies is addressed. 2.1 The classical BEM theory Before proceeding to the following derivations, it is noted that the procedure is the same as in [4]. The blade element momentum method is applied through two main procedures, one where a momentum balance is applied on a rotating annular stream tube passing through a turbine rotor plane and the second where the aerodynamic loads are computed on the aerofoils using the lift and drag coefficients associated with the different elements of the blade. These coefficients for different aerofoils are available as a function of attack angle from the provided inputs of the project. Usually this data is obtained from wind tunnel testing. Interest is put in the second formulation whichhttps://www.sharelatex.com/project/ signifies the calculation of the thrust and power output from the turbine. The most simplistic situation considered in this project is steady-state flow completely aligned with the rotor axis. The velocity profile is considered uniform, thus the viscous effects of the ground are neglected. An introduction of four different stations is made. The first station is far upstream of the rotor plane, the second and third just before and after the rotor plane, respectively, and the fourth is far downstream of the rotor plane. Assuming inviscid flow as well as p1 = p4 and V2 = V3 , Bernoulli’s equation can be applied to represent the momentum balance which yields; 1 (2.1) p2 − p3 = ρ(V12 − V42 ) 2 Furthermore, when considering the forces, the assumption is made that the radial span can be divided into annular stream tubes of width dr and area as seen by inflow velocity dA. The axial force in terms of flow properties is defined as; dFx = (p2 − p3 )dA (2.2) The axial induction factor a is defined as a= V1 − V2 V1 (2.3) Moreover, the following expressions are easily derived; V2 = V1 (1 − a) , V4 = V1 (1 − 2a) (2.4) 4 Chapter 2. Theory Combining these expressions yield the following expression for the incremental axial force obtained from momentum conservation; 1 dFx = ρV12 [4a(1 − a)]2πrdr 2 (2.5) As the flow enters the rotor plane, a rotation described by the angular velocity ω is induced by the movement of the blades which rotate at the angular velocity Ω. From a consideration of physics, the following expression for the incremental torque dT is obtained; dT = ρV2 ωr2 2πrdr (2.6) The tangential induction factor a′ is defined as a′ = ω 2Ω (2.7) Combining these expressions yield the equation dT = 4a′ (1 − a)ρV Ωr3 πdr (2.8) To summarize, the conservation of momentum has resulted in expressions for incremental axial and tangential forces. When proceeding to the aerodynamic analysis a number of assumptions are made to simplify the model. These are that there are no aerodynamic forces interacting between blade elements and that the forces are considered only as functions of the lift and drag coefficients. x W V(1-a) r r 2 blade rotation r r wake rotation F IGURE 2.1: Velocity triangle on a blade element. The flow angle β is the angle between rotor axis and relative flow. The initial inflow has no rotation but as the flow passes over the rotor disk a rotational component of velocity is induced. To take into account the generated angular momentum of the wake the tangential velocity over the airfoil is assumed to be the average of the inlet and the outlet of the rotor disk. Thus the tangential velocity becomes Ωr + 21 ωr where ω is the angular velocity of the wake. Together with the tangential induction factor the component of the relative velocity W in Figure 2.1 becomes: Ωr + ωr = Ωr(1 + a′ ) 2 (2.9) 5 Chapter 2. Theory From the definition of the axial induction factor it follows that; tan β = λr (1 + a′ ) Ωr(1 + a′ ) = V (1 − a) (1 − a) (2.10) where V represents the incoming velocity. Also from observing Figure 2.1 we have; W = V (1 − a) cos β (2.11) The next step in the procedure is to consider the aerodynamic forces acting on the blade element. In accordance with common practice, the lift and drag forces are defined as perpendicular and parallel to the incoming flow, respectively. This is illustrated in Figure 2.2. Note that β is the angle between the axial inflow velocity and the velocity relative to the airfoil, Θt and Θp are the twist and pitch angles, respectively, and α is the attack angle. θt + θ p x V(1-a) α r r F L D Fx F IGURE 2.2: Lift and drag forces acting on the blade element. To be able to relate to the expressions for axial and tangential forces obtained previously from consideration of momentum conservation, the lift and drag force per span are translated in these directions to form the incremental tangential and axial force components. It follows from observing Figure 2.2; dFθ = dL cos β − dD sin β , dFx = dL sin β + dD cos β (2.12) Where dL and dD are the incremental lift and drag forces respectively. To proceed with eliminating these unknown variables from the expressions above, equations for these are introduced as follows; 1 (2.13) dL = CL ρW 2 c dr 2 1 dD = CD ρW 2 c dr (2.14) 2 Thus after elimination, the following expressions for tangential and axial forces are obtained, where B represents the number of blades; 1 dFx = B W 2 (CL sin β + CD cos β)c dr 2 (2.15) 6 Chapter 2. Theory 1 dFθ = B W 2 (CL cos β − CD sin β)c dr (2.16) 2 Upon multiplying the tangential force by the radius, the torque dT is obtained as; 1 dT = B W 2 (CL cos β − CD sin β)cr dr 2 (2.17) Introducing the local solidity factor σ ′ and utilizing Equations 2.10 and 2.11, the following expressions are derived for the axial force and torque; V 2 (1 − a)2 (CL sin β + CD cos β)rdr cos2 β (2.18) V 2 (1 − a)2 (CL cos β − CD sin β)r2 dr cos2 β (2.19) dFx = σ ′ πρ dT = σ ′ πρ In order to take into account the losses that occur at the blades, a loss correction factor, F, is introduced into the expressions for axial force and torque generated from momentum theory, Equations 2.5 and 2.6, to obtain the following equations; dFx = F ρV12 [4a(1 − a)]πrdr (2.20) dT = F 4a′ (1 − a)ρV Ωr3 πdr (2.21) The correction factor F is further described in section 2.2.1. To summarize the undergone procedure so far, the objective has been to generate two sets of equations for the incremental axial force and torque. The first set was obtained from a consideration of momentum conservation, whereby these quantities are expressed in terms of axial and tangential induction factors. The second set was generated from expressing the aerodynamic forces acting on the blades in terms of the lift and drag coefficients CL and CD which, for the sake of repetition, are assumed to be known functions of the flow attack angle from the input data. As these alternative expressions are equated, Equations 2.20 and 2.18 as well as Equations 2.21 and 2.19, the characteristic, coupled axial and tangential induction equations arise; σ ′ [CL sin β + CD cos β] a = 1−a 4F cos2 β (2.22) a′ σ ′ [CL cos β − CD sin β] = 1−a 4F λr cos2 β (2.23) As the iteration process is concluded and converged values of induction factors and thus also the pertinent velocities are achieved, the next step in the analysis is calculating the power and thrust output. These are the two desired outputs of the whole process, and will thus also be of use when analyzing numerical aspects of the procedure such as mesh convergence behaviour. When proceeding to analyze power output, it is stated that the power contribution from each individual annulus is; dP = ΩdT (2.24) 7 Chapter 2. Theory Subsequently, the total power is found by means of integration across the radial span, using the above formulation; P = Z R dP dr = rh Z R ΩdT dr (2.25) rh Here, rh is the hub radius. The power coefficient is defined as the ratio between the produced power and the available power in the wind, thus it can be interpreted as a form of efficiency. RR P r ΩdT dr (2.26) = 1h 2 3 CP = Pw 2 ρπR V By combining this expression with equation 2.19, it is possible to directly develop an integral expression for the power coefficient by means of algebra; Z P 8 CD CP = = 2 tan β]dλr (2.27) λQλ3r a′ (1 − a)[1 − Pw λ λh CL 2.2 Corrections to classical BEM The plain BEM derived in Chapter 2.1 can only account for wind speeds that are uniform and normal to the rotor disk. It also neglects the losses present in the turbulent region close to the wind turbine nacelle as well as at the blade tips. In order to correct this, various alterations are implemented. These will be further discussed in this section. 2.2.1 Tip and hub losses One way of correcting the classic BEM theory is by considering the aerodynamic effects of losses for the hub and tip part of the rotor blades. The main purpose of the loss correction factors is to account for vortices generated by the blades. The influence of vortices on the induced velocity field is important since this effects the power output, especially at the tip of the blades. This correction will therefore lead to a more realistic approximation of the wind turbine performance. [6] In this project, a loss correction factor according to Prandtl’s formula have been used. This factor varies between 0 to 1 and represents the decrease of forces at the blades. Prandtl’s loss factor corrects the assumption of an infinite number of blades for the classic BEM method to a rotor with finite number of blades. The correction factor F, defined in Equation 2.28, is integrated in the equations for axial and tangential induction factors, Equations 2.22 and 2.23. In order to take into account the losses that occur at the tip of the blades, Equation 2.29 is inserted into Equation 2.28. [3] [6] F = ftip 2 cos−1 e−f ; π B R−r = 2 r cos β (2.28) (2.29) Just as for the tip of the blades, the hub section can be corrected with a loss factor considering the impact by vortices on the induced velocity field. The formula is similar to the 8 Chapter 2. Theory one used for tip loss. To account for the losses at the hub, Equation 2.30 is inserted into Equation 2.28. [7] fhub = B r − Rhub 2 r cos β (2.30) To account for both tip and hub losses the factors are simply multiplied to a total correction factor, F. [7] F = Ftip ∗ Fhub (2.31) 2.2.2 Yawed inflow In yaw conditions a correction to the induction factor a is introduced, askew . Since flow conditions with yaw varies dependant on blade position, this variable is calculated for each position.[7] It is important to emphasize that askew is not calculated by iterations like the previous induction factor in Equation 2.22. Rather, Equation 2.22 and 2.23 are used to find converged induction factors for the component of wind normal to the disk plane. Now a is replaced with askew : 15π r χ askew = a [ 1 + tan cos ψ ] (2.32) 32 R 2 where χ is the skew angle calculated from the yaw angle γ: χ = (0.6a + 1)γ (2.33) Realizing that askew will take different values for a specific element depending on the blade position, askew will actually have a impact on the flow angle. Looking at Figure 2.1 it is apparent that the vector V∞ is now changing with blade position so the flow angle β is best solved for in 3D. For each blade position we know in which direction rΩ and V∞ are oriented. Projecting V∞ onto the cut-plane of the blade gives us the local velocity V∞,loc . The angle between rΩ and V∞,loc is the flow angle and can be calculated with: − → −−−→ ! rΩ · V∞,loc π β = − cos−1 − (2.34) → −−−→ 2 |rΩ||V∞,loc | −· ) represents vector quantities. Where (→ 2.3 Vortex method The traditionally used BEM method is a simple and fast method but does not cover all aerodynamic effects and is only useful for limited flow conditions. The BEM method can be extended with various corrections but the assumptions are often based on simplified physics. Using a vortex method results in a more accurate and physical prediction of the aerodynamic forces on a wind turbine with a wider flow condition range. Chapter 2. Theory 9 In a vortex method, the rotor blades are replaced by a lifting line or surface. The vortex flow generates finite lift and is modelled as circular streamlines around a fixed point. The vortical form of a wake, caused by the rotating blades, is modelled either by ring elements or horseshoe formed trailing vortices. The main objective of the vortex method is to account for the effects of wake rotation and a finite number of blades when predicting wind turbine performance.[1] There are several different types of vortex methods. The free wake vortex method is the most complex vortex method since it builds on few confining assumptions. This method assumes inviscid, incompressible and irrotational flow where the geometry of wake vorticity is based on a finite number of helical vortex elements whose position depends on the local velocity field. [2] [8] In the prescribed wake model, which is another type of vortex method, the blades are modelled either by a lifting line or a lifting surface. The vortices generated from the blades create a helical vortex sheet with a constant pitch angle in the axial direction behind the rotor. Unlike the free wake vortex model the prescribed wake method does not allow for wake expansion. [1] [2] 2.4 Previous studies The use of BEM in order to assess the power output of wind turbines is established and a widely used method in the industry. There are even commercial as well as open source codes available that utilize it. One of these open-source codes is the FAST code offered by the National Renewable Energy Laboratory (NREL). It makes use of unsteady BEM as a foundation but then further expands into the realms of aeroelastic calculations and control theory. This goes to show that BEM is in fact a sound model, although it has its limitations. Previous studies have been made to compare BEM with more complex models such as the vortex method. These models have even incorporated alterations to account for different yaw angles or even shear inflow. Using BEM with these alterations generally works well but as soon as conditions become extreme, it results in relatively poorer predictions and eventually a divergence of the iterative processes involved. One aspect that has been studied is that BEM, relying on iterative solutions, becomes unstable at higher tip speed ratios. Compared to the vortex method, BEM should provide relatively lower power outputs and predict well the reduced output for yawed inflow.[8] 10 3 Method The output for wind velocities between 5 and 25 m/s was calculated using corresponding pitch and rotational velocities provided by the NREL 5MW normal operational data, which is illustrated in Table 3.1. For wind velocities less than 12 m/s, there are no pitch angles introduced. Yaw angles were varied from negative 20◦ to positive 20◦ for all wind speeds using different correction factors. TABLE 3.1: Operating conditions. U [m/s] 12 13 14 15 16 17 18 θ p [◦ ] 4 6.65 8.70 10.46 12.07 13.55 14.93 U [m/s] 19 20 21 22 23 24 25 θ p [◦ ] 16.24 17.49 18.7 19.87 21.01 22.1 23.17 U [m/s] 5 6 7 8 9 10 11-25 rot.speed [rad/s] 0.627 0.753 0.878 1.003 1.129 1.255 1.267 3.1 Blade discretization The very first step is to discretize the domain and assign values of given data, that is needed to perform the calculation, to each element. The discretization of the span of the blade will be carried out using the cosine-rule, half cosine rule and equidistant spacing. Discretizing the domain with the cosine rule makes the element size near the hub and the tip become significantly smaller than the rest of blade. The half cosine rule only refines the mesh at the tip. In this way the large aerodynamic gradients near the ends are captured in a better way compared to an equidistant mesh with an equal number of element. The difference between the three mesh types for 25 number of elements are illustrated in Figure 3.1. As mentioned before this study will only encompass analysis of one particular type of turbine blade, more specifically the blade on the NREL 5MW baseline wind turbine. All dimensions of the blade and needed data were provided for this project. Using these dimensions a discretized blade was formed through linear interpolation of chord length and twist angle as well as aerofoil type. Tables of lift and drag was provided for all angles of attack for the numerous different aerofoil types. [9] 11 Chapter 3. Method Full Cosine rule y [m] 5 0 -5 0 10 20 0 10 20 0 10 20 y [m] 5 30 40 x [m] Half Cosine Rule 50 60 70 30 50 60 70 50 60 70 0 -5 y [m] 5 40 x [m] Equidispaced 0 -5 30 40 x [m] F IGURE 3.1: Different mesh types for 25 number of elements. 3.2 Classical BEM The following list illustrates the steps in the iteration process to calculate both the a and a′ . Each iteration process occurs for a certain wind velocity, rotational speed and each element of a blade. ′ 1. Set a and a equal to zero. 2. By Equation 2.10 determine the angle of flow, β. 3. Equation 1 gives the angle of attack, α. 4. One can interpolate a CL and CD value from data based on the α calculated in the previous step. ′ 5. By Equations 2.22 and 2.23 determine a and a . new 6. If | a−a anew | and | step 2. ′ ′ ′ a −anew ′ anew | ≤ 10−5 continue with step 7 otherwise start over from 7. When a and a are known for all elements, continue with the power calculation. 12 Chapter 3. Method 3.3 BEM with yawed inflow Since the flow angle β now changes with blade position this will in turn affect the angle of attack and finally the lift and the drag forces of each blade position. The steps for implementing yaw are listed below: 1. Follow the steps for no yaw condition to get a and a′ for all elements using only the wind component normal to the rotor disk. 2. Solve Equation 2.32 with Equation 2.33 where ψ is the current position of a blade. 3. Use askew to calculate the local flow angle β to get new CL and CD . 4. Continue to power calculations using local β, CL and CD . Note that askew should be used instead of a in the power calculations. 3.4 Power calculations From the process of calculating a and a′ the local flow angle β as well as the resulting lift and drag coefficients have been established. This is also true for yaw conditions. The total power and thrust generated are constant for classical BEM but taking yaw into consideration a time-averaged output will have to be computed. The time averaged output can be computed from the average tangential and normal forces of each blade element travelling one full revolution in the disk. The time averaged (denoted (·)) net force along the tangential and normal direction of any individual element i, is: dF θ,i = dLcosβ − dDsinβ , dF x,i = dLsinβ + dDcosβ (3.1) Where the instantaneous lift and drag force of each element in any blade position is given by Equations 2.13 and 2.14. Once the contribution of each element is known they can be used to sum up the total power and thrust generated by B blades: P =B N X i=1 ri dF θ,i Ω , T =B N X dF x,i (3.2) i=1 Where the summation is across all elements. These are the common calculations involved in the load assessment on individual blades, whether it be with or without yaw misalignment. 13 4 Results and Discussions The methods listed in Chapter 3 provided results which are contained within the following five sections. To make the presentation and discussion of the results more easily read, the discussions of the results have been contained alongside the results. 4.1 Convergence study of mesh Figure 4.1 illustrates the total power and thrust calculated with respect to number of elements respectively, for three different mesh types. As can be seen both the power and thrust generation display a highly fluctuating behaviour before settling around the converged value visualized as black lines. A plausible reason to this is that the values of chord length, twist as well as lift and drag coefficients are interpolated from the provided discretized data which cause the solution to be highly mesh dependent for coarser meshes. However, the power and thrust pertaining to the full cosine mesh seem to stabilize faster around the converged value, only requiring about 20 elements to display a satisfying level of accuracy. A solution based on equidistant mesh is apparently more sensitive as it settles around the converged value at around 30 elements. Thus the benefit of using the full cosine mesh is that a higher rate of convergence will be achieved than that of the other two mesh types. Convergence of Power, V = 15 [m/s] 5.6 Full Cosine Half Cosine Equidistant 5.5 Convergence of Thrust, V = 15 [m/s] 440 Full Cosine Half Cosine Equidistant 430 Thrust [kN] Power [MW] 5.4 5.3 5.2 420 410 5.1 400 5 4.9 390 5 10 15 20 25 30 Number of elements 35 40 5 10 15 20 25 30 Number of elements F IGURE 4.1: Mesh convergence for total power and thrust with incoming velocity of V = 15 m/s and zero yaw 35 40 14 Chapter 4. Results and Discussions 4.2 Comparison between BEM and the Vortex method A comparison between the results of the Vortex method, NREL turbine results and the results from the BEM method has been made for both the power and thrust generation of the wind turbine. These results are illustrated in the Figures 4.2 and 4.3. All results show the same trend. The power increases with increasing wind speeds but upon reaching 11 m/s the power generated is constant while the thrust decreases. This is because of the added pitch which aims to keep the power generated constant by reducing angles of attack. These reduced angles of attack will then generate less drag and as a consequence, less thrust. One can see that the difference between the turbine results and the BEM method results is small both for the power and thrust. When comparing the results from the vortex method with the results from the BEM method in Figure 4.2, the former provides almost equal or slightly higher values of power generation. However, as seen in Figure 4.3, the vortex method always gives lower values of thrust generation as compared to the BEM method. 0 ° yaw 6 Our code NREL turbine results Vortex Method 700 5 600 Our code NREL turbine results Vortex Method 4 Thrust [kN] Power [MW] 0 ° yaw 800 3 500 400 2 300 1 200 0 100 5 10 15 20 25 5 10 Wind velocity [m/s] 15 20 25 Wind velocity [m/s] F IGURE 4.2: Power generated without yaw. F IGURE 4.3: Thrust generated without yaw. TABLE 4.1: Power output obtained from the vortex and BEM methods. Method/Yaw angle vortex Method [MW] BEM [MW] Difference [%] 0◦ 5.27 4.95 6.1 5◦ 5.22 4.9 6.1 10◦ 5.06 4.74 6.3 15◦ 4.80 4.49 6.5 20◦ 4.47 4.15 7.2 In Table 4.1 the power output relating to both the BEM and vortex method when the wind velocity is at 11 m/s is illustrated. One can see that the difference increases as the yaw angle increases and that the vortex method predicts lower power output than the BEM method. The thrust at the same wind velocity for both the BEM and vortex method is illustrated in Table 4.2. It is observed that the difference between the two methods is the smallest for no yaw condition and increases with increasing yaw angle. Interesting to note is that the vortex method seem to predict higher thrust than the BEM method for this particular wind speed contrary to the power characteristics. 15 Chapter 4. Results and Discussions TABLE 4.2: Thrust output obtained from the vortex and BEM methods. Method/Yaw angle vortex Method [kN] BEM [kN] Difference [%] 0◦ 688 706 2.5 5◦ 684 702 2.6 10◦ 670 691 3.1 15◦ 647 672 3.9 20◦ 615 646 5.0 In Figures 4.4 and 4.5 the influence of the sign of the yaw angle is visualized. The comparison is made for low magnitude yaw of 5◦ and high magnitude of 20◦ . The power and thrust generated by the positive and negative counterparts are calculated to be exactly equal in the case of the BEM method and the vortex method predicts a small deviation between the two. This could be accredited to the sensitivity of the two different methods, and the vortex method proves to be more susceptible to changes in the direction of yaw angle introduced. 6 Yaw sign- Thrust comparison Yaw sign- Power comparison 800 700 5 600 4 Thrust [KN] Power [MW] yaw angle =-5 yaw angle = 5 yaw angle =-20 yaw angle =20 vortex-yaw angle =-5 vortex-yaw angle = -20 vortex-yaw angle = 5 vortex-yaw angle = 20 3 yaw angle =-5 yaw angle = 5 yaw angle =-20 yaw angle =20 vortex-yaw angle =-5 vortex-yaw angle = 20 vortex-yaw angle = 5 vortex-yaw angle = 20 2 1 10 15 20 Wind velocity [m/s] F IGURE 4.4: Power with ±yaw. 400 300 200 0 5 500 25 100 5 10 15 20 25 Wind velocity [m/s] F IGURE 4.5: Thrust with ±yaw. Figures 4.6 and 4.7 illustrate a comparison between different yaw angles of 0◦ , 5◦ , 10◦ , 15◦ , 20◦ with the BEM code and the vortex method. There is significant reduction in the power and the thrust output in the both the methods but the difference in the higher velocities is more striking in the BEM code. It is to be noted that the pitch angles introduced to maintain constant power output in the case of no-yaw was adopted for the cases with the yaw angles introduced as well. The vortex method still produces a constant power output which the BEM does not. 16 Chapter 4. Results and Discussions Power comparison-yaw angles Thrust comparison-yaw angles 6 800 700 5 600 4 Thrust [KN] Power [MW] yaw angle = 0 yaw angle = 5 yaw angle = 10 yaw angle =15 yaw angle = 20 vortex- yaw angle =5 vortex- yaw angle =10 vortex- yaw angle =15 vortex- yaw angle =20 3 yaw angle = 0 yaw angle = 5 yaw angle = 10 yaw angle =15 yaw angle = 20 vortex- yaw angle =5 vortex- yaw angle =10 vortex- yaw angle =15 vortex- yaw angle =20 2 1 10 15 400 300 200 100 0 5 500 20 5 25 10 15 20 25 Wind velocity [m/s] Wind velocity [m/s] F IGURE 4.7: Thrust comparison with yaw. F IGURE 4.6: Power comparison with yaw. In Figure 4.8 the normal and tangential force per length for each segment over the blade is plotted. The operational conditions are; an incoming velocity of 11 m/s and no yaw angle. Important to point out is that the vortex method data used in this comparison does not extend as close to the hub as the data obtained with the BEM method does. However up until about x/R ≈ 0.17 the contribution from the elements to the overall power is very little as one can see in Figure 4.8. This is due to that the cylindrical airfoils near the hub will give no lift at all but only drag. The negative tangential force seen in Figure 4.8 to the right for the blue line, is a result of the dominating drag force projected down to the rotor plane acting in opposite direction to the rotation of the turbine blade. Normal force distribution, zeros yaw 8 BEM Vortex Method 0.6 F T [kN/m] 6 F N [kN/m] Tangential force distribution, zero yaw 0.8 4 2 0.4 0.2 0 0 BEM Vortex Method -0.2 0 0.2 0.4 0.6 x/R 0.8 1 0 0.2 0.4 0.6 x/R F IGURE 4.8: Normal and tangential force distribution for velocity V = 11 m/s and yaw angle γ = 0◦ . 0.8 1 17 Chapter 4. Results and Discussions For the normal force the BEM corresponds well with the vortex method in the middle of the blade and close to the tip. Closer to the hub the solution starts to diverge from the vortex method data, especially for the tangential force. In the area of largest normal and tangential forces respectively, the BEM method seem to underestimate the value consistently. This underestimation further increases as the yaw angle increases which can be seen in Figure 4.9. In the range of around 0.22 < x/R < 0.4 the BEM seem to give higher values than the vortex Method for both normal and tangential force distribution. This could imply that the hub correction factor for the BEM method utilized here fails to represent the true loss close to the hub. More on that in section 4.3 Tangential force distribution, yaw = 20 ° 0.8 F T [kN/m] 0.6 0.4 0.2 0 BEM Vortex Method -0.2 0 0.2 0.4 0.6 0.8 1 x/R F IGURE 4.9: Tangential force distribution for velocity V = 11 m/s and yaw angle γ = 20◦ . 18 Chapter 4. Results and Discussions 4.3 Loss correction factors As one can see in Figure 4.10, where the wind velocity is plotted against generated power for the classic BEM method (no yaw), the impact of the loss correction factors results in a reduced power output. For the case with only the hub correction factor considered, the power reduction for wind velocities below 11 m/s is negligible compared to when no correction factor is used. At low wind velocities, the tip loss correction results in a slightly decreased power output, but when reaching velocities around 10 m/s there is a significant decrease compared to the case with no loss. For wind velocities above 12 m/s, where the pitch angle is introduced, the hub correction deviates from the case without losses and decreases the generated power, while the tip loss correction has a upgoing tendency with increasing wind velocity. This is shown in Figure 4.11. For the highest wind velocity of 25 m/s the tip correction factor will cause an even higher power output than the case with no loss corrections. The reason for this is that at this high velocity the pitch angle actually turns the blade causing the tip to experience a negative angle of attack. So by considering losses at the tips these negative contributions of lift will be damped resulting in a higher generated power than with no tip losses. When implementing loss correction factors for both the tip and hub this results in a power curve that follows the tip correction, but levels out for higher velocities as it is influenced by the power decreasing hub factor. Loss factor comparison 6 Loss factor comparison 5.5 5 3 Power [MW] Power [MW] 5.4 No Loss Hub Loss Tip Loss Hub & Tip Loss 4 2 5.3 5.2 No Loss Hub Loss Tip Loss Hub & Tip Loss 5.1 1 5 4.9 0 5 10 15 20 Wind velocity [m/s] F IGURE 4.10: Loss factor comparison for power 25 12 14 16 18 20 Wind velocity [m/s] F IGURE 4.11: Zoom of loss factor comparison. 22 24 19 Chapter 4. Results and Discussions Figure 4.12 illustrates how the thrust depends on loss factors with respect to the wind velocity. Just as for the power comparison the hub loss correction will have a marginal impact on the thrust except for higher wind velocities where it has a decreasing tendency. The tip loss correction factor overall results in a small thrust reduction, but slowly starts to increase towards the case without loss factor for velocities above 12 m/s due to the increasing pitch angle. Loss factor comparison 800 No Loss Hub Loss Tip Loss Hub & Tip Loss 700 Thrust [kN] 600 500 400 300 200 100 5 10 15 20 25 Wind velocity [m/s] F IGURE 4.12: Loss factor comparison for thrust. 4.4 Effects of yawed inflow P In the Figure 4.13 the normalized power, Pmax has been plotted with respect of the tip speed ratio, λ for a number of yaw angles. The power for each yaw angle is normalized with the Betz limit which has the following Equation 4.1.[5] Pmax = 16 ρπR2 3 v 27 2 (4.1) The Betz limit is the maximum power that a wind turbine can extract from the incoming wind. In Figure 4.13 the wind velocity is 11 m/s, which is right before the pitch kicks in. One can see that the power output decreases with an increased yaw angle. One can also notice that for a tipspeed ratio from one to about four, the normalized power for all the different yaw angels are equal. 20 Chapter 4. Results and Discussions P/P max for different yaw angles and tip speed ratios 0.7 0.6 0.5 P/P max 5° yaw 0.4 -5° yaw 10 ° yaw 0.3 -10 ° yaw 15 ° yaw 0.2 -15 ° yaw 20 ° yaw 0.1 -20 ° yaw 0 1 2 3 4 5 6 7 8 9 10 r F IGURE 4.13: P Pmax for different yaw angles vs. tip speed ratios. Figure 4.14 shows the instantaneous torque generated as a function of location in the rotor plane. Without yaw, there is no variation in the torque generated and so it is constant throughout the full revolution of a blade. However most of the torque is generated across the central region of the blade with the regions close to tip and hub becoming almost negligible in comparison. Adding yaw to the system greatly alters the torque generated as now the wind blows 15◦ from the left. Most of the torque is now generated on the left side of the rotor plane and there is a large variation in torque generated as a function of blade position. Tip and hub regions are relatively unaffected as their torque generation is still negligible in a global scope. F IGURE 4.14: Comparison of azimuthal variance in torque generation. All values have been normalized with the maximum local torque generated without yaw. Comparison is with 11 m/s and a yaw angle of 15◦ . Chapter 4. Results and Discussions 21 For thrust there is a similar relation as with the torque and this is shown in Figure 4.15. Just as expected there is no azimuthal variance in thrust generation without yaw. With the yaw introduced, there is a greater generation of thrust on the windward side of the rotor disk. Even visible is the greater thrust on the bottom half of the rotor disk, caused by the rotor blades rotating against the flow as opposed to with the flow at the top of the disk. Figure 4.15 also demonstrates the possible extension of BEM into the realms of structurally dynamic calculations as it shows the distributed load of the blade and how it varies with time. F IGURE 4.15: Comparison of azimuthal variance in thrust generation. All values have been normalized with the maximum local thrust generated without yaw. Comparison is with 11 m/s and a yaw angle of 15◦ . Chapter 4. Results and Discussions 22 4.5 Altering blade pitch Altering the pitch effectively alters the power output and thrust. Figure 4.16 demonstrates this as it shows the power and thrust for different pitch angles at an incoming wind velocity of 23 m/s. The black lines represent values for the nominal rotational speed of the NREL turbine. Also marked in magenta is the normal pitch angle associated with the wind speed. Altering the pitch angle from its normal operating angle we see that increasing the pitch lowers the power and decreases the thrust while decreasing the pitch has the opposite effect. It might seem intuitive that the power produced should always F IGURE 4.16: The net effect of altering pitch is lower power output and thrust. Nominal rotational speed i 1.267 rad/s and normal operating pitch angle is marked in magenta. increase when reducing the pitch. This is not true however since the lift generated for an increased angle of attack will depend on whether the blade is operated before or after its stall. When increasing the pitch with a stalled blade, there will be a greater generation of lift and vice versa. This stall point is indicated on the power plot of Figure 4.16. To the left of this point, the blade will generate more lift for an increased pitch angle. To the right of this point the blade will generate less lift for an increased pitch angle. Altering the rotational speed of the rotor also has an effect on power generated. It’s effect depends on whether the blades are completely stalled or not. If the blades are well below their stall angle, increasing rotational velocity will decrease the power output and decrease thrust. This is because the increased rotation will reduce the angle of attack. Decreasing the rotational velocity at this point will have opposite effects for the same reason. Should the blades be well above their stalling point the effects are switched. This happens since increasing the angle of attack beyond the stall point reduces the lift generated. Since the blade are twisted and composed of different aerofoil sections, different sections can stall at different times. The point where lines meet in Figure 4.16 is characterized by the fact that approximately equal contributions of blade span is stalled vs not stalled, so altering the rotational speed has minimal effects on both power and thrust. 23 5 Conclusion Code for implementing the BEM method was developed using MATLAB, the power and thrust characteristics were studied for a velocity range of 5 m/s to 25 m/s. Three different types of meshes were employed; namely the half cosine, full cosine and equidistant mesh. The results obtained from the parametric study carried out in this project gave a qualitative insight into the convergence characteristics of the BEM method. The convergence of power and thrust generation with an increasing number of elements were studied for all the mesh types at a constant fixed wind velocity. Strong fluctuations in convergence were observed and accredited to the discontinuous distribution of blade parameters along the span of the blade. This distribution is a direct consequence of the discretization of the blade, where different element values of parameters were applied in different discrete regions of the span. Observing the results pertaining to the convergence rate relating to the different mesh types, it is concluded that the half cosine and full cosine provides a higher rate of convergence for the same number of elements than the equidistant mesh. The reason to this behaviour is that larger gradients of power and thrust generation as well as different parameters at the tip of the blade are resolved to a greater extent when a finer mesh in this region is applied. Given that the computational time of the BEM method is negligible, a low sensitivity of the convergence rate to the mesh type and number of elements was observed. However, if a more extensive method was used and the mesh was enriched, the trend would be more noticeable. The characteristic power and thrust curves obtained from half cosine mesh type with 30 elements were compared to the performance of 5 MW NREL wind turbine. Equal contributions by the blade in power and thrust generation was a major take away in the consideration of classical BEM method. The accuracy of BEM is a crucial property to consider, while acknowledging the computational benefits it offers as compared to the alternative methods. The required computational time using the vortex method for one combination of velocity and yaw angle is roughly one hour. In this project, with the specific parameters studied, the BEM method required less than five seconds to analyze the cases of one yaw angle combined with 21 different wind velocities, respectively. Hence, the choice between the vortex method and the BEM method is a concern of choosing between the accuracy of the solution and the computational efforts of the analysis. The yaw angles were introduced as an attempt of extending the BEM method to take into account deviations of wind velocity directions. Comparisons were made on the sensitivity of the sign of the yaw angles and the BEM method proved to be unaffected by the direction of the yaw angle, unlike the vortex method which developed a marginal difference. The magnitude of the yaw angles were altered and it was observed that the power and the thrust generation of the blades decreased as the yaw angles increased. A similar trend was observed with the vortex method results as well and the maximum deviation between the results obtained with the two methods was calculated to be 7.2 % at a wind velocity of 11 m/s. At a velocity of 23 m/s the absolute deviation between the results of the BEM and vortex methods were large. This could be associated to the rigidity of the nature of BEM and vortex method codes to generate constant power for different yaw conditions with the same pitch angle used for that of now yaw condition. The vortex method claims high rigidity in this regard. In efforts to increase the fidelity of Chapter 5. Conclusion 24 the implemented BEM method, physical phenomena like hub and tip losses were introduced. It was observed that the impact of accounting for these losses is not noticeable in the correction to benchmark vortex method results. In the study of the power coefficients using different tip speed ratios one could observe that every yaw condition has an optimal tip speed ratio where maximum power is delivered. These ratios were not the same for any given yaw condition but rather the ratio decreased for an increased yaw misalignment. This iterates that the BEM method is more appropriate for smaller tip speed ratios and no yaw misalignment model. As a characteristic of the yaw condition where the blades develop unequal power, a distribution of one blade over one complete revolution was studied and distinguished with a no yaw counter-part. The unsteady loads on the blade through its revolution could be a point of departure for structural analysis regarding flutter and fatigue. However, the assumption of time-averaged behaviour in this project limits the scope of the analysis. The results obtained could serve as a point of departure for future work relating to the study of aeroelastic properties. For future work it would also be interesting to implement corrections for shear inflow. Furthermore, a comparison of this BEM model with the more complex CFD model would be an interesting study. As previously mentioned, the code assumes time-averaged behaviour so dynamic aspects of the operation are neglected. Should the code be made instantaneous in time, these aspects could be utilized in flutter analysis and/or the regulation of pitch. It would therefore be interesting to expand this study in that direction. 25 References [1] [2] [3] [4] [5] [6] [7] [8] [9] Hamidreza Abedi. Aerodynamics Loads On Rotor Blades. 2011. Hamidreza Abedi. Development of Vortex Filament Method for Wind Power Aerodynamics. 2016. Martin O. L. Hansen. Aerodynamics of Wind Turbines second edition. 2008. Grant Ingram. Wind Turbine Blade Analysis using the Blade Element Momentum Method, version 1.1. 2011. Anthony L. Rogers James F. Manwell Jon G. McGowan. Wind Energy Explained: Theory, Design and Application, second edition. 2009. Emrah Kulunk. Aerodynamics of Wind Turbines. 2011. A.C Hansen P.J Moriarty. AeroDyn Theory Manual. 2005. J. Gordon Leishman Sandeep Gupta. Comparison of Momentum and Vortex Methods for the Aerodynamic Analysis of Wind Turbines. 2005. J Jonkman. S Butterfield. V Musial. G Scott. Definition of a 5-MA Reference Wind Turbine for Offshore Developement. 2009. Investigation of Wind-Induced Noise to Optimize Masts Project in Applied Mechanics J ONATHAN FAHLBECK J OHAN F ORSGREN S ANKAR R AJU N ARAYANASAMY M ARTIN OTTOSSON DAVID W INQVIST Department of Applied Mechanics C HALMERS U NIVERSITY OF T ECHNOLOGY Gothenburg, Sweden 2017 Project in Applied Mechanics 2017 Investigation of Wind-Induced Noise to Optimize Masts Jonathan Fahlbeck Johan Forsgren Sankar Raju Narayanasamy Martin Ottosson David Winqvist Department of Applied Mechanics Division of Fluid Mechanics Chalmers University of Technology Gothenburg, Sweden 2017 Investigation of Wind-Induced Noise to Optimize Masts Jonathan Fahlbeck, Johan Forsgren, Sankar Raju Narayanasamy, Martin Ottosson, David Winqvist © Jonathan Fahlbeck, Johan Forsgren, Sankar Raju Narayanasamy, Martin Ottosson, David Winqvist, 2017. Supervisors: Dr. Hua-Dong Yao and Prof. Lars Davidsson, Applied Mechanics, Chalmers Examiners: Dr. Håkan Johansson and Dr. Valery Chernoray, Applied Mechanics, Chalmers Project in Applied Mechanics 2017 Department of Appied Mechanics Division of Fluid Mechanics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: An isosurface of the λ2 criterion with the velocity magnitude vizualised along the surface, for mast geometry (b). Gothenburg, Sweden 2017 iv Abstract In modern society the awareness of disturbing noise has increased. Thus it is of importance to find the cause of noise generation and prevent the designs which produces them. In this report the wind-induced whistling noise from two sailboat mast profiles have been evaluated. It has previously been observed that one of the profiles generates a whistling noise while the other one does not. The predetermined conditions were a free stream velocity of 20 m/s and a yaw angle of 33◦ . The commercial software STAR-CCM+ has been used along with aero-acoustic features. The simulations showed that no whistling noise could be identified. Although several interesting features regarding the flow around the mast geometries have been discovered, no root cause of the whistling noise could be completely established. In order to fully understand the noise generation, additional work needs to be performed on the subject. Thus it is recommended to examine more yaw angles and velocities of the incoming flow. In this manner it is the authors believe that it would be possible to finally find the answer of what is causing the whistling noise. Keywords: Whistling, Noise, STAR-CCM+, Mast, FW-H, Aero-Acoustics, Vortex Shedding v Acknowledgements First of all we would like to thank our supervisor Hua-Dong Yao. Without his continuous support and advice, the overall understanding and progress within the aero-acoustic field, would not have been attained. His happy spirit helped us through many dark hours of computational errors. The computer resources have been of utter most importance to secure results. A special thanks is thus directed to Chalmers University of Technology. Finally we would like to thank Seldén Mast. Without their initial idea and interest, this project would not have been possible. William Holt deserves extra acknowledgement for his commitment to the progress of the project. Jonathan Fahlbeck, Johan Forsgren, Sankar Raju Narayanasamy, Martin Ottosson, David Winqvist, Gothenburg, May 2017 vii Contents List of Figures xi List of Tables xi Nomenclature xii 1 Introduction 1.1 Background . . . 1.2 Case Description 1.3 Project Goal . . . 1.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Theory 2.1 Acoustic Theory . . . . . . . . . . . . . . . . . . 2.2 Computational Theory . . . . . . . . . . . . . . 2.2.1 Basics of Fluid Dynamics . . . . . . . . 2.2.2 Basics of Computational Fluid Dynamics 2.2.3 Ffowcs Williams-Hawking formulation . 2.2.4 DNS, LES, RANS, DES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 2 2 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 5 5 6 6 7 3 Method and Computational Settings 3.1 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Acoustic Domain Setup . . . . . . . . . . . . . . . 3.1.2 Mesh generation . . . . . . . . . . . . . . . . . . 3.1.3 FW-H Surface and Receivers . . . . . . . . . . . . 3.1.4 Mesh Dependency Study . . . . . . . . . . . . . . 3.2 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Boundary Conditions and Material Properties . . 3.2.2 RANS Simulation . . . . . . . . . . . . . . . . . . 3.2.3 DES Simulation . . . . . . . . . . . . . . . . . . . 3.2.4 Stopping Criteria . . . . . . . . . . . . . . . . . . 3.3 Post-processing . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Fast Fourier Transform and Sound Pressure Level 3.3.2 Overall Sound Pressure Level . . . . . . . . . . . 3.3.3 λ2 Criterion . . . . . . . . . . . . . . . . . . . . . 3.3.4 Contour Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 8 9 11 11 12 12 12 13 13 13 14 14 14 14 4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . 15 ix Contents 4.1 4.2 4.3 4.4 4.5 4.6 Velocity . . . . . . . . . . . . . . Sound Pressure Level . . . . . . . Overall Sound Pressure Level . . λ2 Criterion and Vorticity. . . . . Additional Results . . . . . . . . General Discussion of the Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 16 17 18 20 21 5 Conclusion 24 Bibliography 25 x List of Figures 1.1 Mast profiles, (a) has a whistling noise, (b) does not have a whistling noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 3.1 3.2 Domain and the mast profile (b). . . . . . . . . . . . . . . . . . . . . 9 View of a directed volume mesh with 5.6 million cells, for mast (b). . 10 4.1 4.2 Velocity contour plots. . . . . . . . . . . . . . . . . . . . . . . . . . A-weighted SPL. Data recorded in two receivers. See Figure 4.4 for the location. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . PSD of static pressure in the probes located close to the tips. Comparison between the two configurations. . . . . . . . . . . . . . . . . Overall Sound Pressure Level around the two different mast configurations. Note that 0◦ refers to downstream of the mast. . . . . . . . Isosurface of λ2 = −250000 1/s2 , for mast configuration (a) (top) and for mast configuration (b) (bottom). . . . . . . . . . . . . . . . . . . Vorticity contour plots. . . . . . . . . . . . . . . . . . . . . . . . . . Mast configuration (a). . . . . . . . . . . . . . . . . . . . . . . . . . Vorticity near the mast. . . . . . . . . . . . . . . . . . . . . . . . . 4.3 4.4 4.5 4.6 4.7 4.8 . 15 . 16 . 17 . 18 . . . . 18 19 20 20 List of Tables 3.1 3.2 Domain design parameters. . . . . . . . . . . . . . . . . . . . . . . . . 9 Material data for dry air at 15 °C and atmospheric pressure [12]. . . . 12 xi Nomenclature AA ASZ CAA CFD CFL CTH DES DNS FFT FW-H LES OASPL PSD RANS SPL Aero-Acoustics Acoustic Suppression Zone Computational Aero-Acoustics Computational Fluid Dynamics Courant–Friedrichs–Lewy Chalmers University of Technology Detached Eddy Simulation Direct Numerical Simulation Fast Fourier Transform Ffowcs Williams-Hawkings Large Eddy Simulation Overall Sound Pressure Level Power Spectral Density Reynolds-Averaged Navier-Stokes Sound Pressure Level 1 Introduction 1.1 Background In today’s society, people are much concerned about noise pollution. Hence industries strive to make products less noisy in order to comprehend with the awareness and sensitiveness of their customers. The aim of this project is to investigate the unwanted whistling noises created by sailboat masts when the boat is moored and the sail is stowed inside the mast. An annoying whistling sound may arise when the flow velocity and the yaw angle are large enough in certain boat mast designs. The company Seldén Mast AB desires to know the root cause of this whistling noise in order to enable them to design the masts in accordance to prevent it. Seldén Mast provided Chalmers University of Technology (CTH) with two mast segments, of which one has been experienced to produce a higher noise in comparison to the other. In Figure 1.1 the two mast profiles are illustrated. Since the observations of the whistling sound is based on experiences rather than experimental data the conditions of when this behaviour occurs is not fully known. It has been estimated to a free stream velocity of 20 m/s and a yaw angle from 20◦ to 45◦ . Figure 1.1: Mast profiles, (a) has a whistling noise, (b) does not have a whistling noise. The aim of this collaboration between CTH and Seldén Mast is to achieve a winwin situation. By using computational fluid dynamics tools, the group at CTH aimed to help Seldén Mast in their understanding on the whistling noise generation. Thereby enabling the company to reduce unwanted sounds generated by the masts and making necessary changes in the design. The proposed analysis of the mast designs will hopefully help the company in making them more competitive. In 1 1. Introduction addition to the assistance provided to Seldén Mast, the group at CTH intended to improve their skills in project work as well as to enhance their knowledge in fluid mechanics and aero-acoustics. 1.2 Case Description The two mast configurations in Figure 1.1 are investigated regarding noise generation. Mast (a) has a cross sectional dimension of 148 mm × 266 mm and mast (b) has a cross sectional dimension of 270 mm × 140 mm. As noise is generated by turbulence, 3D computations are performed. The height of the two profiles are set to 250 mm, and thus the height of the 3D domain. The flow velocity is 20 m/s and the yaw angle 33°. The fluid is dry air at 15 °C and atmospheric pressure. As 3D computations are challenging, instability issues such as numerical errors, convergence problem and flow field development are expected. Mesh generation is also an area where challenges might be encountered. 1.3 Project Goal Use acoustic wave modeling to explain the effects of the tip1 geometry in a mast regarding noise generation. The project should result in recommendations to consider while designing a mast to prevent noise generation. 1.4 Limitations In order for the project to be feasible within the time frame a number of boundaries and limitations are applied. The project was to be carried out from March to May 2017. The limits are determined by either absolute factors, such as the amount of time available, or in order to assure that an overall quality of the work is attained. Since unsteady 3D computations was performed, which are both time and computationally demanding, time available and computational resources are of importance. Hence it is hard to establish a fine balance between number of cells and computational time. The field of aero-acoustics was new for the project group. In the narrow time frame a large amount of the time available had to be used in order to fully understand the underlying problems in the field of aero-acoustics. Thus the task is limited to investigate noise caused by turbulence. Other types of sounds are not to be investigated. All the simulations are performed in the commercial software STAR-CCM+. The software is also used in the analysis part of the project in combination with MATLAB. 1 The tip geometry refers to the tips at the opening of the mast, see Figure 1.1. 2 2 Theory 2.1 Acoustic Theory Acoustics refers to the science of sound and originates from the Greek word for hearing [1]. Sound is of high significance in our every day life when it comes to communication, awareness of our surroundings and orientation. In recent years the awareness of sound and noise has increased. One needs to note that not all sound is bad. We still want to hear a click when we turn the key in the door, music or sounds alerted by warning devices. The sound which today is unwanted is the one that is "out of place" and disturbing [2]. It is known that sharp bangs or long time exposure to loud noises can damage our hearing. In fact noise can also cause other health issues or even death. For example, flies that are exposed to sound levels around 160 dB die after a short period of time. To get this noise level into perspective one can compare this with the noise produced by a refrigerator ~50 dB, a chain saw ~100 dB or a shotgun ~140 dB. Humans who spend much time in areas with loud noise usually have a higher blood pressure and heart rate. Noise has also been proven to affect sleeping patterns and just simply annoys people. All these facts have made companies and whole industries more aware of the effects of noise. Today acoustics is an important aspect in product development [3]. In order to understand how noise travels, one can explain it by looking at what happens when you throw a rock into a lake. When the rock hits the water surface, the rock will be slowed down. The kinetic energy of the rock then is transferred into the water and causes ripples to form on the water surface. The ripples are gradually transported outwards and heating the water slightly as they pass and fade away. In the same way the noise travels from, for instance, a sudden clap of your hands. Energy will spread from the clap in a series of sound waves which consists of regions with increased pressure. The fluid particles will move closer together for a short period of time as the wave passes. For a louder clap, the fluid particles move tighter together and a higher pressure difference is created. Like the ripples on the water, the sound waves of the clap will die away as the energy eventually is transferred into heat. Since the losses are small, sound waves can travel a very long distance before they fade away [2]. There are different ways by which noise or sound can be created. For instance, the one mentioned above which was as a "clap". Others are frictional noise or vibrations 3 2. Theory noise from structures. Since the cause of these sounds are fairly straight forward, the industries ability to dampen these sounds has come pretty far. In recent years the removal of the "out of place" sounds has therefore been focusing on other sources of sound such as aero-acoustics. An example of the importance of aero-acoustics is when you are in a car driving at high velocities. If you open the window the aeroacoustic noise will be the dominating one. Aero-acoustic noise generation can in a simple manner be explained by the flow interacting with geometrical irregularities of a car, an airplane, etc. The interaction creates unsteady turbulent flows which are often detached and thus in turn generates noise. The noise generation can be divided into two phenomena: impulsive noise and turbulent noise [3][4]. Impulsive noise is generated by the movement of surfaces or a surface in a nonuniform flow condition. The non-stationary load on the body causes pressure fluctuations to occur, which are generated as sound. This noise can be estimated in a fairly simple manner from aerodynamic simulations. Turbulent noise, on the other hand is harder to predict but it is quite common and predominantly exists. Since the turbulence is stochastic by nature, it has a broad frequency spectrum. The turbulent noise therefore creates a broadband noise consisting of many frequencies. Turbulent energy is most efficiently transferred into acoustic energy in the presence of sharp edges. The sharp edge forces two flows of different velocities to have a sudden blend. An example of this is when a flow crosses a bluff body. Unless the body is aerodynamically shaped the fluid will separate from the body and cause a wake to form behind the body. The wake will consist of fluid with low velocity. The pressure difference between the still fluid and the fluid crossing the body will cause them to blend. This causes strong local equalizing flow which in turn results in peak pressures, i.e impulsive noise. Behind a bluff body the pressure will characteristically alternate to the different sides of the wake in order to comprehend with the velocity of the moving flow on both sides of the wake. This effect is a Von Kármán vortex street which refers to the fluid twirls of alternating directions forming behind an object. Vortex shedding is the phenomenon of every fluid twirl "coming loose" from the body [3][5]. As stated in previous paragraph, vorticity is generated at boundaries by the relative velocity of two surfaces, such as fluid and wall. When a certain threshold is reached the two-dimensional wakes formed by vorticity, shed from the bluff body surface. This results in a transition into the third dimension. The feature explains the fact that the vorticity can be stored as a vortex street, leading to conservation of vorticity in a system. Since the vorticity is mainly created just downstream of a bluff body the vorticity further downstream is merely a response of the rest of the fluid, trying to adjust to the instability that the bluff body initialized. Vortices tend to merge downstream of the bluff body, creating larger sections of rotating fluid with lower velocity. As a result one could state that the vorticity generated immediately downstream of the bluff body is of higher importance [6]. Since turbulent broadband noise always exists when turbulence is present, one can state that Aero-Acoustic (AA) noise consists of a broadband noise. Sometimes it is 4 2. Theory accompanied by a narrow-band impulsive noises. As previously stated, though the impulsive noise is fairly easy to predict, in order to predict the turbulent noise one has to estimate the turbulence. Hence, in order to fully estimate the AA noise one has to use CFD-tools. This is a fairly recent methodology that has been evolved with the utilization of computers. This field of science is called Computational AeroAcoustics (CAA) [3][4]. There are two main methodologies in CAA to compute an acoustic field. The first one being "direct method", which is considered to be the most exact and would be the equivalent to DNS in the CFD field. In the direct method the governing equations for the flow and acoustic field are solved over the entire domain, from the aerodynamic effective area to a far-field observer. The fact that the domain of interest is very large, makes the direct method extremely expensive when it comes to the number of cells and time steps. The more common way to perform AA computations is to use hybrid methods. The sound generation in the aerodynamic effective area is decoupled from the transport of the sound to the far-field. In more simple terms one can say that one method is used for the sound generation and another method is used for the transport process. There are several methods available, both for sound generation and the transport process [3]. Even though the possibility to fully estimate the aero-acoustic noises is fairly new due to improvement in technology (computers) today, a lot of the ground work was performed in the 1960s. Many of the equations still used today are based on the acoustical transport techniques from half a century ago. The most frequently mentioned ones are the Lighthill analogy and the Ffowcs Williams–Hawkings (FWH) equation. These are used as transport methods while performing calculations with the hybrid techniques [3]. The FW-H equation is explained in the section 2.2.3. 2.2 Computational Theory In order of fully understanding CAA some basic principles need to be addressed. As stated in section 2.1, CAA is a mixture of acoustic theory and CFD. In this chapter the theory behind these aspects will be explained. The aim is that this chapter will help the reader in the understanding of the final results. 2.2.1 Basics of Fluid Dynamics The principle of fluid dynamics is that it describes a pattern of flows. In fluid mechanics it is common to investigate a flow within a fixed control volume. The flow is governed by the continuity (eqn 2.1), Navier-Stokes (eqn 2.2) and energy (eqn 2.3) equations [7]. dρ ∂vi +ρ =0 (2.1) dt ∂xi 5 2. Theory " dvi ∂p ∂ ∂vj 2 ∂vk ∂vi ρ =− + + − δij µ dt ∂xi ∂xj ∂xj ∂xi 3 ∂xk ρ !# + ρfi du ∂vi ∂qi = σij − dt ∂xj ∂xi (2.2) (2.3) Here ρ is density, v velocity, x length, t time, p pressure, u internal energy, q conductive heat flux, f body force and σ the stress tensor. 2.2.2 Basics of Computational Fluid Dynamics In CFD the variables of a flow field are solved for inside a domain of interest using a computational software. The domain is discretized into a number of finite volumes or cells in which the governing equations are solved in an iterative process. Depending on the type of flow, different numerical schemes are used to discretize the differential equations. 2.2.3 Ffowcs Williams-Hawking formulation This formulation is a more general form of Lighthill’s equation [8][9], where moving walls are allowed to be present inside the domain. The formulation is based on the assumption that a volume B(t) exists which is enclosed by the surface S(t). The surface needs to be sufficiently smooth to allow the definition of a smooth function x, t) such that h(x    > 0 if x ∈ B(t) x, t) = = 0 if x ∈ S(t) h(x    < 0 if outside B(t). One can note that the Heaviside function of this function, H(h), will be zero inside the volume B(t) and unity outside B(t). If the mass conservation equation and momentum equation is multiplied by H(h), use Lighthill’s procedure for acoustic variable, p′ = p − p0 , apply Green’s theorem and use the free-space Green’s function, equation 2.4 can be obtained [3]. # # " " ∂ Z fH ∂ 2 Z (ρvi vj − σij )H x, t) = dVy − dVy p (x ∂xi ∂xj 3 4πr ∂xi 3 4πr τ =te τ =te ′ R 2 + ∂ ∂t2 Z " R3 R (p ′ /c20 − ρ )H 4πr ′ # dVy (2.4) τ =te # # " " ∂ Z p′ ni − σij nj ρ0b · n ∂ Z dS − dS + ∂t 4πr(1 − Mr ) τ =te ∂xi 4πr(1 − Mr ) τ =te S(te ) S(te ) x − y ||, f σij is the viscous stress, n is the normal to the surface, H = H(h), r = ||x is the body force, c0 is the speed of sound in the fluid surrounding the listener, b is x − y )/rc0 and te = t − r/c0 . the velocity of the moving surface, Mr = b· b·(x 6 2. Theory The two surface integrals in equation 2.4 are later used to set up the FW-H surface integral as equation 2.5. See section 3.1.3 for further information. # # " " ∂ Z ∂ Z p′ ni − σij nj ρ0b · n x, t) = p (x dS − dS ∂t 4πr(1 − Mr ) τ =te ∂xi 4πr(1 − Mr ) τ =te ′ 2.2.4 (2.5) S(te ) S(te ) DNS, LES, RANS, DES The major models that are utilized for solving turbulence in CAA or CFD as a whole are Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), Reynolds Averaged Navier-Stokes Simulation (RANS) and Detached Eddy Simulation (DES) [7]. DNS solves all the turbulence scales, making it highly expensive. Thus it is not suggested to be performed for objects with larger dimension with the computational facilities currently available. LES resolves all large turbulent scales, i.e, it is volume averaged. The principal idea behind LES is to reduce the computational cost by modelling or exclude the smallest scales, which are the most computationally expensive to resolve. This is achieved via averaging across time and space, which effectively removes small-scale information from the numerical solution. RANS models all turbulence. It is time-averaged and different turbulence models are applied. Thus making it faster than other methods to compute. The most common turbulence models are the k − ε or k − ω. DES is a combination of advantages of both RANS and LES. In DES, the scales close to the wall are solved using RANS (or unsteady-RANS), within the boundary layer. Further away from walls are the scales resolved using LES. Thus DES behaves as a hybrid RANS-LES model. 7 3 Method and Computational Settings While performing the analysis of the two mast designs a methodology divided into three stages was used. The different stages were Pre-processing, Simulation and Post-processing. Most of the work was performed in the commercial software STARCCM+. The work was initially performed on mast profile (a) (see Figure 1.1) due the simplicity of the geometry. All through the process, thorough documentation of the choices made were performed. This was later used when repeating the set up for mast (b). In this fashion a precise comparison between the two designs could be established. 3.1 Pre-processing When performing analysis using numerical methods, the quality of the final result is highly dependent of the work done prior to pushing the "run button". A well planned domain using a suitable mesh design to account for the important features is therefore crucial. In this section the different parts of the Pre-process will be explained. 3.1.1 Acoustic Domain Setup The outer domain design was based on the mast mean diameter Dm . In Figure 3.1 and in Table 3.1 the design features are visualized. The principle of the design in Figure 3.1 is that the white area (outer domain) is formed by a coarser mesh while the gray area (inner domain) has a finer mesh. This feature was applied in order to capture the small-scale flow structures downstream of the mast. The bold line which is dividing the areas is a FW-H surface, see section 3.1.3. A RANS simulation was performed and gave a good estimation of the size and location of the turbulent flow field behind the mast. This simulation was used to define the gray shaded refinement zone and the FW-H surface. In general is it necessary to use a very large domain when numerically calculating acoustics. This is because the pressure fluctuations that noise is built from are small and therefore sensitive to disturbance. Any presence of reflected noise could affect the result. The following methods were used to reduce the risk of reflecting noise and numerical errors: • To reduce the reflections of sound with high frequency an increased cell size near the boundary was used. This makes the mesh unable to capture the 8 3. Method and Computational Settings sounds since approximately 20 cells per wavelength is needed to capture the sound correctly [10]. This design also has the advantage that it reduces the computational effort. • Regarding the sound with low frequency an acoustic Suppression Zone (ASZ) was used around the boundary. ASZ dampens the acoustic waves before they reach the outer boundaries to prevent them from reflect back into the domain. The zone can therefore be thought of as a sponge absorbing the sound. The zone size, l, could be estimated from the speed of sound, c0 , and predicted minimum frequency, according to equation 3.1. c0 = 340 m/s, fminimum = 500 Hz ⇒ l= c0 fminimum = 0.68 m (3.1) Table 3.1: Domain design parameters. Parameter Used size Dm (D + d)/2 R 15Dm L 25Dm r 3Dm w 13Dm D 0.27 m d 0.14 m (a) Domain seen from above (top) and from (b) Illustration of mast profile (b), the side (bottom). with geometrical lengths. Figure 3.1: Domain and the mast profile (b). 3.1.2 Mesh generation In order to find a suitable mesh for performing the simulations, a couple of parameters needed to be taken into account. Firstly, the mesh should not exceed more than 10 million cells. This was to ensure that the computational resources provided were used effectively. Secondly, the mesh needed to be fine enough in the areas where small-scale vortex shedding occurs. This is where sound generation is expected. The main areas of interest were the region around the tips (Figure 3.2c) and the region behind the mast (Figure 3.2b). The mesh in these regions was controlled by a wake 9 3. Method and Computational Settings refinement operation. Here one could specify the direction of the anticipated flow, to which the geometry will be exposed. The length, angle and the cell growth rate of the zone were controlled. The mesh can be seen in Figure 3.2. (a) Top view of the whole domain. (b) Zoomed on the wake region. (c) Zoomed on the tips. (d) Cut plane around the mast. Figure 3.2: View of a directed volume mesh with 5.6 million cells, for mast (b). To properly resolve the acoustics, the cell size must be small enough. One good practise in estimating the cell size is to use the fact that at least 20 cells per wavelength is preferable [10]. With the speed of sound, estimated to be the standard value at sea level, c0 = 340 m/s, combined with an approximate value for what the frequency of the noise might be, fwhistling ≈ 1000 Hz, the cell size can then be calculated according to equation 3.2. ∆CAA = c0 20fwhistling = 0.017 m (3.2) This cell size was desirable to transport the sound correctly in the inner domain. In section 3.1.1 a couple of ways to minimize the risk of reflected noise are listed. These strategies were included into the mesh. In addition to the mesh design in regards where to have large/small cells, the actual mesh design needed to be accounted for. In 3D simulations, as in this case, three common mesh types are tetrahedral, hexahedral and polyhedral. Tetrahedral is the least complex with four vertices, four faces and six edges. The other ones are more 10 3. Method and Computational Settings complex and are built from several tetrahedral cells. This results in the tetrahedral mesh being the fastest to generate, but lacks in accuracy when solving the actual problem, since the error decreases as the number of cell faces increases. In acoustic simulations, a large number of iterations is in general required in order to attain a correct solution. Therefore, even though the mesh takes longer time to generate, the polyhedral mesh was the most suitable candidate, due to a higher accuracy at a shorter computational time. The whole mesh generation process was performed in an iterative manner where different mesh strategies were tested. Methods that were tested and later discarded were for example mesh designs using volumetric control for the entire domain, another one using surface controls. The finally chosen mesh generation method was a directed mesh. The advantage of this strategy is that the mesh is easy to control regarding number of cells. The mesh is also faster to compute in comparison to the other used methods. The fact that this thesis aims to compare two different geometries was an argument for the choice of an easily controlled mesh and thus the directed mesh. 3.1.3 FW-H Surface and Receivers In order to capture the noise generated by the mast, a permeable FW-H surface was applied within the computational domain. The data from CFD computations were then extracted along the FW-H surface as suitable source terms in the acoustic equations [11], see section 2.2.3 and equation 2.5. The FW-H surface was applied as an internal interface between inner and outer region within the fluid domain, see bold line in Figure 3.1. Once AA had been enabled in the used software, the interface was treated as a permeable surfaces within the "FW-H surface" menu. A number of FW-H receivers (36) was adopted outside of the fluid domain to record the acoustic data. The receivers were installed in a circular pattern, with a radius of 10 m, around the mast. Since the acoustic transport equations are decoupled from the flow equations it is possible to place the receivers outside of the domain. 3.1.4 Mesh Dependency Study A mesh dependency study is planned in order to establish mesh independent results. It is desirable to find a mesh giving accurate results without wasting computational resources. The study was planned so that the general design of the mesh were unchanged and the only difference was the number of cells. Three different meshes were created for mast configuration (a) (see Figure 1.1) with 2.2, 4.4 and 7 million cells. Probes were inserted into the domain to measure velocity and pressure to later be able to compare the results. The plan was then to choose one mesh and apply this 11 3. Method and Computational Settings set up also for mast configuration (b). Unfortunately problems were encountered regarding finding a converged solution. A lot of time was spent on fixing this and when a converged solution was obtained the time was too short to analyze all the meshes. Therefore, a fully completed mesh dependency study was not performed. For the final results, a 4.4 million cell mesh and a 5.6 million cell mesh were used for mast configuration (a) and mast configuration (b), respectively. The difference in cell number between the mast configurations is due to varying complexity in the geometry. 3.2 Simulations This section consists of different settings and options used for the simulations in STAR-CCM+. The settings and preferences used in all simulations are based upon [10], discussions with supervisor and the project group. 3.2.1 Boundary Conditions and Material Properties The names of the domain boundaries can be found in Figure 3.1. The Inlet boundary was set as a velocity inlet with a laminar incoming flow.Sides and Outlet were set as pressure outlet. Atmospheric pressure at sea level was prescribed for the pressure outlets. The Upper and the Lower boundaries were set to periodic boundary condition where the flow propagates from one boundary to the other. The fluid is dry air at 15 ◦ C and atmospheric pressure (101325 Pa). Other material properties under consideration are presented in Table 3.2. The ideal gas law is a viable assumption for the cases investigated. It is a fair approximation, since the regular air at non-extreme temperatures and pressure are used in these simulations. Table 3.2: Material data for dry air at 15 °C and atmospheric pressure [12]. T [K] 288 h i kg ρ m 3 1.226 h µ 10−6 Nm·s2 17.96 i h W κ 10−3 m·K 25.24 i h J Cp kg·K 1006 i Pr 0.7159 T temperature, ρ density, µ kinematic viscosity, κ thermal conductivity, Cp specific heat capacity, P r Prandtl number. 3.2.2 RANS Simulation Initially a RANS simulation was performed to understand the turbulent flow field around the mast. This simulation was performed to get an initial solution, from which a more comprehensive DES simulations could be performed, see section 3.2.3. The RANS simulation was performed as a steady state simulation using segregated flow solver but without any AA features. The SST k − ω model was used together with all y + wall treatment. During the RANS, first order schemes were used so that a solution for the flow field could be attained rapidly. 12 3. Method and Computational Settings 3.2.3 DES Simulation DES and AA tools were utilized in order to find and determine the source of the noise. The initial solution was obtained from the steady state RANS simulation. In contrast to the RANS the DES was carried out with an implicit unsteady scheme where a low CFL number was desirable. A Segregated flow solver was applied, just as for the RANS. The utilized model was SST k − ω with all y + wall treatment enabled. In the DES, second order schemes were used to achieve a higher accuracy. The convection was treated with a Hybrid-BCD scheme. During DES, the number of iterations per time step was used as an important factor to govern the over all convergence. It is critical to attain convergence within every time step and then consequently convergence in time as well. The number of iterations per time step was set to 10 due to promising convergence while monitoring the residuals. To develop the correct unsteady behaviour in the flow field the fluid needs to pass through the domain several times. In practice this meant that two different options could be used. Either the simulation needed to go on for a long time, which would be expensive regarding computational power. Another alternative is to start at a large time step and develop the field and then gradually decreasing it until a sufficiently low time step is reached. The second option was used during this project. The time step was first set to 10 ms and then decreased in several steps to a final time step of 0.1 ms. The final time step gave a maximum CFL number of around 10. This was considered the lowest possible time step due to the computational power available and was therefore good enough for this project. The enabled acoustic options were "FW-H unsteady", along with an Acoustic Suppression Zone (ASZ). The FW-H was used with the "On-the-fly" model option which provides noise prediction in the receivers in parallel with the CAA. The Inlet, Sides and Outlet boundaries had the ASZ set to 0.68 m. At the other boundaries no ASZ was applied. The ASZ was configured to suppress low frequencies since high frequency is dampened by increasing cell size near the boundaries, see Section 3.1.1. 3.2.4 Stopping Criteria The stopping criteria were based on residual plots, monitor plots and receivers in the far field. It was of importance to let the simulation run for sufficiently long time for the receivers to capture enough noise data. 3.3 Post-processing In the Post Processing stage, the focus is to analyze the gathered data. This section explains the methods of comparison when analyzing the wind induced noise level for the different masts. 13 3. Method and Computational Settings 3.3.1 Fast Fourier Transform and Sound Pressure Level For analyzing the aero-acoustics, a Fast Fourier Transform (FFT) algorithm was used to process the signals recorded by the receivers. The receivers record the pressure fluctuations over time and the FFT algorithm transforms the time dependent signals from the time domain into the frequency domain. After the FFT the Sound Pressure Level (SPL) was plotted as function of frequency to distinguish tonal peaks in the signal. This operation was efficiently performed using an FFT tool in STARCCM+. The SPL was calculated using the A-weighting function. This means that low frequency noise is filtered out in the similar way as in the human ear. In order to analyze if the noise recorded in the receivers was generated from vortex shedding produced by the tips of the masts, probes were introduced in these areas to monitor the static pressure. An FFT algorithm computation was then performed on the pressure-time history from the probes. This procedure was similar to that of the data obtained from the receiver. The difference is that instead of the SPL, the Power Spectral Density (PSD) was computed. PSD describes the signal power per unit frequency. To confirm if the sound originates from this region, the pressure spectra from the receivers and that of the probes should be of similar pattern. If a pressure peak is present in the SPL from the receivers, a similar peak should appear in the PSD from the probe data. This would imply that the noise is generated in the region of the probe. 3.3.2 Overall Sound Pressure Level In order to visualize the direction in which the highest level of noise is produced, the Overall Sound Pressure Level (OASPL) was computed in every receiver location. The OASPL is based on the root mean square value of the pressure fluctuations. MATLAB was utilized for these calculations. 3.3.3 λ2 Criterion The λ2 criterion is a technique that was utilized to visualize the vortices in the turbulent flow using isosurfaces. An isosurface is a surface in space where a certain variable is constant, in this case λ2 . This method utilizes the eigen value of the strain rate tensor, S and vorticity tensor, Ω. λ2 is then the median eigenvalue from S 2 + Ω2 , i.e, λ1 ≤ λ2 ≤ λ3 [13] [14]. Hence it was of interest to investigate λ2 , as the vorticity produced in the wake region is known to have a strong correlation to the generation of noise. 3.3.4 Contour Plots Velocity plots were utilized to visualize the flow field. It was important to analyze the vorticity contour plots, as they depict the vortex shedding which is correlated to noise generation. 14 4 Results and Discussion This chapter presents the results and the discussion of those. To simplify the understanding of noise generation, the flow field is discussed first. Gradually the various physical aspects of noise generation are presented and analyzed. The results are based upon the 4.4 million mesh for mast (a) and a 5.6 million mesh for mast (b). (a) Mast configuration (a). Top view (top) and side view (bottom). (b) Mast configuration (b). Top view (top) and side view (bottom). Figure 4.1: Velocity contour plots. 15 4. Results and Discussion 4.1 Velocity To get a better understanding of the flow field around the masts, the velocity contour plots are shown in Figure 4.1. As expected, the flow fields look similar for both configurations. A wake pattern has been formed downstream of the mast and shows a Von Kármán vortex street. The vortex street is caused by the unsteady flow separation. From Figure 4.1 one could see that the vortices in the wake are created from the bluff body and not from the tip geometry. Mast configuration (b), which was not expected to create a high level of noise, also seems to create a larger unsteadiness than mast configuration (a). This indicates that mast configuration (b) should produce a higher level of noise than mast configuration (a), according to the theory in section 2.1. 4.2 Sound Pressure Level As stated in section 2.1, the noise consists of a turbulent and an impulsive part. The turbulent noise is often referred to as broadband noise and consists of noise in several frequencies. This project has been aimed to investigate the source of the tonal noise which originates from the impulsive noise. In Figure 4.2 an FFT has been performed on the pressure-time data from the receivers and is presented as an A-weighted SPL. The spectra for both mast configurations show a broadband noise without any obvious protruding peaks, indicating that no tonal noise can be observed. The graphs show that the SPL for configuration (b) is higher than for mast geometry (a). Although, the sound pressure levels are considered to be relatively low. (a) Receiver at 90°. (b) Receiver at 270°. Figure 4.2: A-weighted SPL. Data recorded in two receivers. See Figure 4.4 for the location. In order to evaluate the origin of the noise, probes were deployed in the domain for data comparison with the receivers. The location of the probes can be seen in Figure 4.3a. In this case the receivers do not indicate any tonal noise generation. In order to verify if this result is correct, it is possible to look at the pressure fluctuations 16 4. Results and Discussion in the probes, see Figures 4.3b, 4.3c and 4.3d. Here it can be seen that no obvious spectral peaks occur. Hence it is possible to state that the recordings in the receivers are reliable and no tonal noise is created in any of the mast configurations. (a) Probe locations. (b) Probe 1. (c) Probe 2. (d) Probe 3. Figure 4.3: PSD of static pressure in the probes located close to the tips. Comparison between the two configurations. 4.3 Overall Sound Pressure Level The Overall Sound Pressure Level for both mast configurations can be found in Figure 4.4. The Figure is oriented as such that the inflow comes from the left hand side (180°). The Figure shows the OASPL in each of the 36 receivers deployed around the mast. The shape of the noise profile is similar for both the cases. The highest noise is found diagonally downstream of the mast. The lowest noise is found upstream of the mast, which is reasonable as the mast is blocking the sound which is created in the wake. It is also evident that the sound level is higher for mast configuration (b). 17 4. Results and Discussion Figure 4.4: Overall Sound Pressure Level around the two different mast configurations. Note that 0◦ refers to downstream of the mast. 4.4 λ2 Criterion and Vorticity. The λ2 criterion is visualized for the two mast configurations in Figure 4.5. As can been seen mast (b), contradictory to real world experience, generate an intenser vortex shedding than (a) near the tips. Thus it could be claimed that mast configuration (b) will generate a larger noise at this yaw angle. Also in this visualization it can be seen that vortices seems to be created from the mast body and not from the tips. This gives a hint that, for this case, the tips do not generate vortex shedding and thus hardly a tonal noise, according to the theory. Figure 4.5: Isosurface of λ2 = −250000 1/s2 , for mast configuration (a) (top) and for mast configuration (b) (bottom). 18 4. Results and Discussion The vorticity contours, see Figure 4.6, show that the mast (b) generates stronger vorticity than configuration (a) near the tips. As vorticity is a good indicator of noise, one could claim that configuration (a) produces less noise compared to mast (b). Also, it could be found that the generation of vorticity is largely governed by the bluff body, rather than the tips at this angle of attack. (a) Mast configuration (a). Top view (top) and side view (bottom). (b) Mast configuration (b). Top view (top) and side view (bottom). Figure 4.6: Vorticity contour plots. 19 4. Results and Discussion By both the visualization of the vorticity in Figure 4.6 and the λ2 -criterion in Figure 4.5 one can see that major vorticity occurs near the mast for the (b) case and not in (a). This could be some indication of that vortex interaction with the tips of configuration (b) will not create any tonal noise, since the interaction occurs but the mast still not produces any tonal noise according to Figure 4.2. For configuration (a) nothing about the correlation between vortex interaction with the tips and tonal noise can be stated, since the vortices are completely separated from the tips. 4.5 Additional Results While the project progressed, several different simulations were performed. Due to the low rate of acoustic simulation experience, many of these simulations were in purpose of practice. Consequently the results from the simulations have been of varying quality. During this period of trial and error a few interesting findings have been made. The results from using an initial coarser mesh with 2 million cells, showed signs of producing more noise than the final result. It is of importance to note that due to the mesh being more coarse, these results cannot fully be trusted. Nevertheless a few interesting details can be observed. In Figures 4.7-4.8 a comparison has been made between an initial result for a coarse mesh and the final result. The sound pressure level for the two meshes can found in Figure 4.7. In Figures 4.8a and 4.8b the vorticity is plotted close to the tips for both meshes respectively. As previously stated in section 2.1, noise and vorticity does in theory have a strong correlation. It can be seen in Figures 4.7-4.8 that the higher sound pressure level coincides with a higher degree of vorticity around the tips of the mast. (a) Coarse mesh. (b) Final mesh. Figure 4.8: Vorticity near the mast. Figure 4.7: Mast configuration (a). 20 4. Results and Discussion 4.6 General Discussion of the Results To summarize the results, it is found that mast configuration (b) generates a higher noise level than (a). This is in contrast to what was expected from real life observations. According to Seldén Mast the noise from mast (a) should be louder than the noise from mast (b). No significant tonal noise has been recorded during the simulations. The only noise that has been captured is broadband noise. The captured noise level has also been very low, approximately 20 dB. It is of interest to evaluate why the sound levels found in the simulations do not coincide with the expectations. A topic of interest is the angle of attack which is evaluated. This angle is set to be 33°. As stated previously, it does not look like the vorticity reaches the tips of the mast (a). From the observations made prior to the project it was expected that the noise was occurring when the flow was coming in at an angle of attack of around 20°- 45°. If another angle of attack had been examined, the results might have differed. A source of error could be that in real life the flow is never laminar. There is always turbulence and fluctuations in the wind. In the presence of turbulent flow, the separation of the flow occurs at a later stage from the surface of the bluff body. This would perhaps convey the fact that the vortices would be created closer to the tips of the mast. Hence the tips of the mast would be exposed to larger pressure fluctuations and noise would be generated. The delayed separation in combination with the fluctuating wind direction would imply that the noise generation in real life is less sensitive to the angle of attack in regard with the simulated case. It is of interest to discuss why there is a difference in result between the initial solution with 2 million cells and the final result for configuration (a). It would have been of great interest to perform a full mesh dependency study. Due to lack of time, this could not be done and hence both results are interesting to discuss. The mesh was in both cases performed in a very similar fashion and the only difference is the number of cells, which were 2 million and 4 million. The difference in cell size is therefore two times in volume, but regarding the length of one cell this results in a difference of 21/3 which is hard to visualize graphically. The first possible reason for the difference in result between the two meshes is that the coarser mesh is simply not fine enough and hence the faulty results are found. Another one could be that due to the size of the mesh, the coarser one ran for a longer time. Hence the flow could be picked up at a later time step. There is a possibility that the final mesh for configuration (a) has not yet reached the time where the vortices occur close behind the mast. Perhaps the most likely reason is that due to the perfect case of laminar flow that hits the mast, the finer mesh resolves this very well and a large wake is created that never gets close to the tips. The coarser mesh on the other hand does not solve the flow as well and there is an error that causes the vorticity to occur close to the mast. Interesting to note here is that to some extent the reason for the vortices acting close behind the mast is of less importance. By using the result from the coarser mesh we can discuss what would happen if the vortices actually 21 4. Results and Discussion did effect the tips. In section 4.5 it was shown that the initial results with a coarse mesh gave a higher SPL than the final results. Here the correlation between SPL and vorticity could be studied. The fact that more noise is created from the mast when there is a higher level of vorticity around the tips is of interest. Once again it is of great importance to remember that these results are not fully verified. But apart from the coarser mesh, there is nothing intending that these results are untrustworthy. These results would imply that the theory as well as the statement from Seldén are verified. That mast configuration (a) is indeed creating a larger noise than configuration (b). This can be seen if comparing Figure 4.7 and Figure 4.2. Presuming that this is correct it is of interest to look at Figure 4.6b and compare that to Figure 4.8a. It seems that in both cases the noise is generated by the vorticity created by the tips. It is then interesting to note that configuration (a) seems to generate more noise. Hence one could state that the design of mast configuration (b) helps to reduce the noise level. From Section 2.1 it is know that sharp edges in general generate noise. This is due to the sudden blend of fluids at different pressure levels. If the Figures 4.6b and 4.8a are compared it is possible to see that there is vorticity around the tips in both cases. Note that the tips that are referred to in configuration (b) are the two enclosing the smaller hole for an extra sail. By comparing the sound level of these two cases, it can be seen that configuration (a) generates a higher noise. The main difference in the tip geometry is that in configuration (a) the distance between the tips is much larger than in configuration (b). From this one could claim that the amount of space after the edge would effect the capacity of noise that can be created. Another interesting aspect to investigate is to regard the area inside the mast, where the two tips create an opening, as a resonance box. In music instruments, such as guitars, the strings create the vibration and with the help of the resonance box the tune is created. The empty space within the mast could possibly work as a resonance box. This would imply that it would be of interest to have as little empty space as possible within the mast in order to reduce the creation of noise. An example of noise that is created in a similar fashion is when you blow into a bottle. The more liquid that the bottle contains, the higher frequency the noise has. This is an effect of the amount of open space in the bottle. A full bottle barely produces any noise. This would therefore imply that it is of interest to make a sail, which when it is folded inside the mast, will fill out the entire cavity all through the mast. Since numerical methods have been used in this project there is always a chance that there are some errors in the simulation set up. Especially so since the time frame was short in regards to the prior knowledge within aero-acoustics of the group. Unsteady 3D computations are demanding and require large resources, therefore the computations were time consuming. It would for instance have been desirable to use a larger domain than the one used in these simulations. Especially in the span-wise direction, where the domain had to be limited in order to attain a number of cells 22 4. Results and Discussion that was low enough. The mesh itself may have been too coarse to fully resolve the correct aero-acoustic features. As a part of the simulations there was an aim to conduct an extensive mesh dependency study. This was partly abandoned due to difficulties in generating a mesh with sufficiently low number of cells. Thus the mesh generation process became very time consuming and valuable time was taken away from the simulations and post processing. In the start-up of the simulations large convergence problems occurred. This was mainly due to difficulties when setting up the physics for capturing the acoustics. Various boundary conditions and properties were therefore tested in order to acquire a converged solution. It should be emphasized that none of the group members had any previous experience in the aero-acoustic field. 23 5 Conclusion The project goal was to explain the effects of the tip geometry in a mast regarding noise generation. It was concluded that the two mast geometries indeed do generate noise. However, from the case examined it cannot be verified that any of the two configurations creates a tonal noise. The results show that the vortex shedding was mainly caused by the bluff mast body itself and not the tips. Aero-acoustic theory suggests that the whistling should occur due to vortex shedding caused by sharp edges. Since the vorticity is low in the tips-region, no major whistling should be generated. It is the believe of the authors that other yaw angles might generate a whistling noise caused by the tips. But as for as now, one could not determine the root cause of the whistling noise. During the project the experienced tonal noise could not be reproduced. No recommendations can therefore be made with certainty to Seldén Mast based on the results. However some interesting aspects regarding the flow around the mast have been identified and analyzed. From this study alone, it cannot be concluded that mast (a) in general generates a higher tonal noise than mast (b). In order to fully understand the noise generation features it is thus recommended to conduct a more comprehensive investigation. A profound study should include several yaw angles and varying free stream velocities. A more extensive mesh dependency study would also be of interest. This is to ensure that the results are reliable enough. 24 Bibliography [1] Junker Miranda, U. (2009) Bonniers Uppslagsbok. ©Bonnier Fakta AB. [2] Goldsmith, . (2012).DISCORD: The story of noise Oxford University Press. ©Mike Goldsmith 2012. [3] Albrecht Wagner, C. Hütti, T. Sagaut, P (2007).Large Eddy Simulation for Acoustics. ©Cambridge University Press 2007. [4] Anselmet, P. Mattei, P.O. (2016).Acoustics, Aeroacoustics and vibrations ©ISTE Ltd 2016. [5] D. Olivari. (2012), Von Kármán vortex shedding, Retrieved from https://www.encyclopediaofmath.org [6] Powell A.,Why Do Vortices Generate Sound?, ASME. J. Mech. Des. 1995;117(B):252-260. doi:10.1115/1.2836464. [7] Lars Davidson, (2017) Fluid mechanics, turbulent flow and turbulence modeling. Unpublished. [8] Lighthill, M. J. (1952). On sound generated aerodynamically. I. general theory. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 211(1107), 564-587. [9] Kaltenbacher, M., Escobar, M., Becker, S., and Ali, I. (2010;2009;). Numerical simulation of flow-induced noise using LES/SAS and Lighthill’s acoustic analogy. International Journal for Numerical Methods in Fluids, 63(9), 1103-1122. doi:10.1002/fld.2123 [10] Axel Kierkegaard (2016). Best Practices for Direct Noise Calculations. [11] Hua-Dong Yao Lars Davidson Lars-Erik Eriksson Shia-Hui Peng Olof Grundestam Peter E. Eliasson (2014), Surface integral analogy approaches for predicting noise from 3D high-lift low-noise wings, pp. 326-338. [12] F.J. McQuillan, J.R. Culham and M.M. Yovanovich (1984). PROPERTIES OF DRY AIR AT ONE ATMOSPHERE. Microelectronics Heat Transfer Lab University of Waterloo Waterloo, Ontario [13] Dong, Y., Yan, Y. & Liu, C. (2016), New visualization method for vortex structure in turbulence by lambda2 and vortex filaments, APPLIED MATHEMATICAL MODELLING, vol. 40, no. 1, pp. 500-509. [14] J. Jeong, F. Hussain (1995),On the identification of a vortex, J. Fluid Mech. 285 69–94. 25 Fatigue Life Improvement of Welded Joints Final report JAKOB ALM EFREM EFREMOV NIKLAS FLÖE TOBIAS MATTSSON SIRI RYDMAN WILLIAM STÅHLBERG Supervisor: Moyra McDill Client: Lennart Josefsson May 31, 2017 Abstract The use of high frequency mechanical impact (HFMI) in the post-weld treatment of certain structural steels is the subject of a benchmark exercise in a specialist committee, Materials and Fabrication Technology, for the International Ship and Offshore Structures Congress (ISSC) 2018 conference. As part of the Chalmers University of Technology contribution to the benchmark exercise, a numerical investigation of the effects of cyclic fatigue loading on the stability of compressive residual stress fields induced by HFMI after welding, is performed. The benchmark exercise specifies the use of S355, a structural steel with a nominal yield strength of 355 MPa, and a particular geometry, that of a stiffener, fillet welded to a membrane-loaded plate. This geometry, provided as a finite element model by the benchmark exercise committee, is known to be sensitive to fatigue as it has a high stress concentration at the weld toe. Temperaturedependent material properties as well as mechanical properties for the simulation of HFMI and the cyclic analysis are also specified. A progressive finite element analysis (FEA) is performed, in which the welding process is first simulated. The resulting residual stress field is then imported into a sequential mechanical FEA of the HFMI treatment and the subsequent cyclic loading. Both the welding process and the HFMI treatment are simulated using simplified approaches. The welding process is simulated by applying a thermal cycle, while the HFMI treatment is simulated by a quasi-static contact analysis. Two different cyclic load cases are considered; the first load case uses a series of constant-amplitude cycles preceded by peak over-and-under loads, while in the second, the peak over-and-under loads follow the constant-amplitude cycles. Throughout the course of the progressive FEA, a number of deficiencies and concerns with the benchmark exercise are identified. In particular, inconsistencies in the S355 material models used in the analysis of the weld versus those specified for the HFMI and fatigue stages are found to be of concern. Further, the benchmark exercise does not yet specify the exact nature of the cyclic loading that is to be applied. The results of the progressive FEA show that the simplified approach to modelling both the welding process and HFMI treatment gives results that correlate well with the experimental and FEA data available in the literature. The results of the cyclic loading demonstrate that the beneficial effects of the HFMI treatment in terms of compressive residual stress in the vicinity of the weld toe are almost eliminated if the test specimen is subjected to unexpected peak loads during the cyclic loading. It is found that the stability of the HFMI-induced residual stress field after constantamplitude cyclic loading cannot be determined. The study of the effects of constant-amplitude loads on the compressive residual stress field from the HFMI is recommended once the benchmark exercise has fully defined the nature of the fatigue loading that is to be used. Because of the inconsistencies in the S355 material models, it is also recommended for these to be clearly specified in order to make results between benchmark participants more comparable. i Contents 1 Introduction 1.1 Goal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 2 Model 2.1 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Benchmark exercise issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 3 3 Material properties 3.1 Thermomechanical weld simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Simulation of HFMI treatment and cyclic loading . . . . . . . . . . . . . . . . . . . 3.3 Benchmark exercise issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4 5 6 4 Progressive FEA 4.1 Thermomechanical weld simulation . . . . . . . . . . . . . . 4.1.1 Boundary conditions . . . . . . . . . . . . . . . . . . 4.1.2 Numerical results . . . . . . . . . . . . . . . . . . . . 4.1.3 Benchmark issues and error sources . . . . . . . . . 4.2 Simulation of HFMI treatment . . . . . . . . . . . . . . . . 4.2.1 Test specimen boundary conditions . . . . . . . . . . 4.2.2 Indentation tool geometry and boundary conditions 4.2.3 Contact properties . . . . . . . . . . . . . . . . . . . 4.2.4 Numerical results . . . . . . . . . . . . . . . . . . . . 4.2.5 Benchmark exercise issues . . . . . . . . . . . . . . . 4.2.6 Error sources . . . . . . . . . . . . . . . . . . . . . . 4.3 Simulation of cyclic loading . . . . . . . . . . . . . . . . . . 4.3.1 Boundary conditions . . . . . . . . . . . . . . . . . . 4.3.2 Load cases . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Numerical results . . . . . . . . . . . . . . . . . . . . 4.3.4 Benchmark exercise issues . . . . . . . . . . . . . . . 4.3.5 Error sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 8 8 9 10 10 11 11 12 12 13 14 14 14 15 15 18 18 5 Discussion 19 6 Conclusions and contributions 6.1 Recommendations and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 20 References 21 Nomenclature γ1 , γ2 Dynamic recovery terms in Chaboche model µ Friction coefficient ν Poisson’s ratio σ0 Initial yield stress in Chaboche model σA Load amplitude C1 , C2 Kinematic hardening parameters in Chaboche model E Young’s modulus h Film coefficient R Cyclic load ratio, ratio is minimum stress to maximum stress. AW As welded C3D6 Continuum 3D element with 6 nodes C3D8 Continuum 3D element with 8 nodes CA Constant amplitude FEA Finite element analysis HFMI High Frequency Mechanical Impact ISSC International Ship and Offshore Structures Congress L Longitudinal OU Over and under R3D4 Rigid 3D element with 4 nodes Room temperature 20°C T Transverse TT Through-the-thickness 1 Introduction Welding, a common and effective method of joining steel structures, is used extensively in the vehicle and shipbuilding industries. Metallurgical and other changes associated with welding often make the weld area prone to damage and crack development from cyclic loading, i.e. fatigue [1]. Of particular concern are the high tensile residual stresses often found in and around the welded area. There are a number of treatments that can be used to counteract this, for example the use of low temperature transformation electrodes, post-weld stress relief and high frequency mechanical impact (HFMI) [1]. In an HFMI treatment of a weld, a cylindrical pin with a radius of approximately 2 mm is used to continously impact the weld zone. These impacts, which are performed at frequencies of at least 90 Hz [2], cause the area around the weld to deform plastically. This plastic deformation results in compressive surface stresses that improve the weld toe geometry and prevent crack initiation, which in turn increases the fatigue resistance of the weld [2]. The use of HFMI in the post treatment of welds is the subject of a benchmark exercise in a specialist committee, Materials and Fabrication Technology, for the International Ship and Offshore Structures Congress (ISSC) 2018 conference. The benchmark exercise specifies a particular geometry, that of a stiffener, fillet welded to a membrane-loaded plate, as shown in Figure 1. This geometry is known to be sensitive to fatigue as it has a high stress concentration at the weld toe. The benchmark exercise combines nonlinear finite element analysis (FEA) and experimental studies. The current project, undertaken as part of TME131, will contribute to the Chalmers University of Technology submission to the benchmark exercise undertaken by the Materials and Fabrication Technology specialist committee through a nonlinear FEA of a benchmark model for the particular geometry. Figure 1: Test specimen used in the benchmark exercise [mm]. [3] It has been shown that HFMI treatment on the welded area works better in higher strength steels [2]. In this project, an S355 low-alloy steel (yield strength of 355 MPa) was considered, which is weaker than the materials that are usually used, such as S700. The benchmark test focused on determining the stability of the introduced compressive stresses, i.e., how well the stress field near the welded area maintains its compressive state after cyclic loading. Lower strength steels have the advantage of being less expensive, which means their use is preferable when possible. A progressive nonlinear FEA was conducted using ABAQUS [4]. In order to test the stability of the compressive stresses introduced by the HFMI treatment, the analysis was divided into two stages: an initial thermal cycle to simulate the thermal stresses from welding and a coupled simulation of the HFMI treatment to introduce a compressive stress field followed by constant-amplitude load cycles with some over-and-under loading. 1 1.1 Goal The goal of the project was to determine the stability of the compressive stress field induced by HFMI in the welded joint to contribute to the benchmark exercise in a specialist committee, Materials and Fabrication Technology, for the International Ship and Offshore Structures Congress (ISSC) 2018 conference. Specifically, the intention was to examine the effects on the stress field after being subjected to fatigue loading. HFMI treatment has been shown to be effective in improving the fatigue life in S700, a high-strength steel [5] and the aim of the project was to examine if the same is true for S355, a low-strength steel. 2 2 Model The model that was used during the course of this project can be seen in Figure 2. It is a half model of the geometry shown in Figure 1 which has been developed specifically for the benchmark exercise and is used by all benchmark participants in order to ensure consistent results throughout the exercise [6]. The mesh was previously used by [3] and [7] in similar studies. In this project a quarter-symmetric model would have been sufficient, but the half-symmetric model was used to be consistent with the benchmark exercise. Figure 2: The meshed model that is used by all benchmark participants. The coordinate system that is used throughout the report is shown in red, where L is the longitudinal direction, T is the transverse direction, and TT is the through-the-thickness direction. 2.1 Geometry The provided model is half-symmetric about the T-TT plane. It consists of 22,464 linear hexahedral elements (C3D8) with 2 × 2 × 2 Gauss points and 504 linear wedge (C3D6) elements with 2 Gauss points. The mesh is refined in the critical area of the weld toe where stress concentrations are known to occur [6]. It was assumed that the mesh had been tested for mesh sensitivity by the benchmark participants. 2.2 Benchmark exercise issues During the course of the project, a number of issues were identified in the benchmark exercise mesh. The following issues were found: • 191 out of 22,968 elements are distorted with poor aspect ratios and distorted geometry. In particular, there are 6 hexahedral elements and 10 linear wedge elements with an aspect ratio greater than 10. This can lead to unrealistic results such as shear locking in hexahedral elements if bending is present [8]. • The mesh is not completely symmetric about the L-TT plane due to possible discretisation errors. • Although the mesh is refined in critical areas it is still very coarse, restricting the ability to fully capture the complex stress fields in the vicinity of the weld. These issues were taken into consideration and worked around, but no changes were made to the actual model, as the nature of the benchmark requires the same model be used by all the participants. The possible effects are discussed in terms of the validity of the results. 3 3 Material properties In the benchmark exercise the stiffener, plate and the fillet weld were modelled as S3551 , a low-alloy structural steel with an expected minimum yield strength of 355 MPa for thicknesses up to 16 mm [10] 2 . In describing the benchmark exercise, the study by [3] uses a yield strength of 315 MPa in the thermomechanical model with isotropic hardening. However, a yield strength of 435 MPa is specified for the HFMI and fatigue simulation, using the Chaboche mixed hardening model [3]. While the benchmark exercise specifies the use of S355, the mechanical properties used in [3], [7], [11], and [12] are by no means consistent. The steel used in [11] is normalized, utilising properties from the SYSWELD3 [13] database. Furthermore, [7] considers phase transformations which can increase the yield limit in the welding simulation, together with a yield strength of 435 MPa for the HFMI and fatigue analysis. The study by [12] used different yield strengths of 339 MPa and 435 MPa in the Chaboche model. This led to some concerns when comparing various results, whether FEA or experiment, in the available literature. 3.1 Thermomechanical weld simulation The thermomechanical weld simulation used the same isotropic temperature-dependent material properties as in the benchmark exercise described by [3], having an initial yield strength of σy = 315 MPa at room temperature. The properties are shown in Figure 3. The choice of isotropic hardening was suggested by [5] to give accurate results when simulating the welding process. Using the same properties as [3] enabled direct comparison of results. Since [3] only provides yield stresses at zero plastic strain, yield stresses at higher strain rates were calculated using linear interpolation from the initial yield stress. The tangent modulus was set to 0.5% of the Young’s modulus based on [14]. Both the Young’s modulus and tangent modulus are temperature-dependent. Figure 3: Temperature-dependent material properties from Kim [3] for S355. The stress-strain behaviour was obtained by cyclic loading of a plane-stress plate in ABAQUS with prescribed displacement, by averaging the stress and strain over the entire region. The resulting stress-strain behaviour is shown in Figure 4 where the isotropic hardening can be observed. This hardening is what causes the resulting tensile stress field after cooling from the welding process. The weld expands as it heats, causing it to yield plastically. As it cools it contracts and hardens and as a consequence a tensile stress field is developed. 1 S355: C%: 0.23 max, Mn%: 1.60 max, Pn%: 0.05 max, Sn%: 0.05 max, Sin%: 0.05 max [9] higher thickness (80-100 mm) a yield strength of 315 MPa can be seen [10]. 3 SYSWELD is an FEA tool for the simulation of phenomena such as heat treatment and welding. 2 In 4 500 400 Stress [MPa] 300 200 100 0 -100 -200 -300 -400 -500 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Strain [%] Figure 4: Cyclic behaviour at 20°C of the isotropic hardening model used in the welding simulation using material properties from [3] and [14]. The isotropic hardening behaviour of the model can be seen. 3.2 Simulation of HFMI treatment and cyclic loading Unlike the thermomechanical weld simulation, the mechanical analysis was not dependent on temperature and used mechanical properties at room temperature. The Young’s modulus and Poisson’s ratio were set to E = 210 GPa and ν = 0.3, respectively. When modelling the hardening behaviour, the Chaboche combined isotropic-kinematic hardening model was implemented [15]. This material model was also used in the benchmark exercise by [3]. The Chaboche model has been shown to perform well in both implicit displacement-driven HFMI simulations and elastoplastic fatigue simulations [16]. As shown in Figure 5a, two sets of material parameters were compared against experimental data in [17] and [18] in order to find the parameter set with the best hardening behaviour. The comparison was made by evaluating the model behaviour for a single plane-stress element under monotonic tensile loading, where the loading was applied as a prescribed uniaxial displacement. Table 1 shows the two parameter sets, where [3] uses the parameters described in [19]. Table 1: Parameter sets for Chaboche kinematic hardening model. Steel grade S355 [19] S355 [12] σ0 [MPa] 435 339 C1 [MPa] 8971.8 33582 γ1 218.65 508.9 C2 [MPa] 12654.88 10061 γ2 106.98 94.24 Furthermore, the cyclic behaviour of the two parameter sets was evaluated in ABAQUS, shown in Figure 5b. It can be noted that the parameter set in [1] is more nonlinear in the strain hardening region. Also, it can be seen that the model captures the Bauschinger effect, i.e an increase in yield strength in one direction occurring at the expense of the yield strength in the other. 5 800 700 600 600 400 Stress [MPa] True stress [MPa] 800 500 400 300 0 200 0 100 = 435 MPa 200 0 -200 -400 = 339 MPa 0 -600 S355 Exp (de Jesus) S355 Exp (Forni) 0 0 = 435 MPa = 339 MPa -800 0 2 4 6 8 -3 10 -2 -1 0 1 2 3 Strain [%] True strain [%] (a) Tension test for single element compared against experimental data from de Jesus [17] and Forni [18]. (b) Cyclic behaviour for Chaboche model with parameter sets from Table 1 (single element). Figure 5: Stress-strain behaviour during monotonic loading and cyclic loading for S355 using the two parameter sets. It was found that the parameter set in [19] with a yield limit of 435 MPa has a closer agreement to the experimental data with regards to yield limit and hardening during monotonic loading and as such was chosen as the parameter set used in the simulation of HFMI-treatment and simulation of cyclic loading. Furthermore, this parameter set has been used previously by [3] and [7] in similar studies. 3.3 Benchmark exercise issues A possible source of error throughout the simulation is the fact that the welding process made use of a simpler material model with a lower yield strength, while the HFMI treatment and cyclic fatigue simulations used the Chaboche model with a higher initial yield strength. This was not considered to be an issue in the current study, since the purpose of the thermomechanical weld simulation was to obtain a representative stress field of the welding process, as verified by comparing to available experimental and FEA data. For the sake of consistency, the benchmark exercise should have provided material data and material models to be used by all participants. 6 4 Progressive FEA As with earlier research related to the benchmark, the focus was on determining the longitudinal stress, i.e., in the direction of the cyclic loading, along two specific paths A and B, starting from the weld toe. The paths, shown in Figure 6, are used in the benchmark exercise by [3] and in related research, e.g. [11]. Path B was also used by [7] in experimental X-ray measurements of residual stresses after welding, HFMI treatment and cyclic loading. These are the basis for experimental comparison. (a) Path A starts from the surface of the weld toe and continues in the longitudinal direction up to 25 mm from the weld toe. (b) Path B starts from the surface of the weld toe and goes through the entire thickness of 5 mm. Figure 6: Paths along which the longitudinal stress is of particular interest in the benchmark exercise. A complete progressive FEA simulation was performed by separating the analysis into two main stages: a thermomechanical weld simulation and a coupled simulation of the HFMI treatment and the cyclic loading. The thermomechanical weld simulation consisted of a forward-coupled transient thermomechanical FEA with the thermal strains driving the stress. The coupled simulation of HFMI treatment and subsequent cyclic loading was performed as a quasi-static mechanical FEA with the deformed shape, residual stresses and material hardening induced by the simulated HFMI treatment. The purpose of the thermomechanical weld simulation was to generate a tensile residual stress field around the weld representative of a welding process. This stress field was then imported as an initial condition into the coupled HFMI treatment and cyclic loading simulation. Other variables, such as the deformed geometry, were not imported as only the tensile stress field was deemed to have a representative effect on the subsequent stages. The purpose of the simulated HFMI treatment was to produce a representative compressive residual stress field at the weld toe, and carry forward the improved weld toe topology and material hardening induced by the HFMI treatment. All stages of the analysis used the full Newton-Rhapson iterative scheme to solve the nonlinear FE models [8]. Post-processing of results from the progressive FEA was performed using MATLAB [20]. The run times were approximately 30 minutes for the thermomechanical weld simulations, 15 minutes for simulation of HFMI treatment and 45 minutes for the computation of cyclic loading.4 4 The simulations were run with ABAQUS 6.14 on a Windows 7 machine with Intel Core i5-3470S CPU @ 2.90GHz and 16 GB memory. 7 4.1 Thermomechanical weld simulation The weld simulation consisted of a forward-coupled thermomechanical analysis where the entire test specimen had an initial temperature of 20°C. The nodal line inside the weld area where the stiffener meets the plate, shown in Figure 7a, was subjected to an incremental temperature increase up to 1300°C 5 and held there for 8 seconds. The nodal line was then quickly cooled to room temperature during the course of 45 seconds, shown in Figure 7b. This step was done as a transient analysis with a minimum and maximum increment size of 0.01 seconds and 60 seconds respectively, which resulted in approximately 80 increments of which most were located at the peak temperature. The specimen was then cooled down to 20°C in a steady-state step. The analysis considered both material and geometrical nonlinearities due to large deformations. 1400 Temperature [C] 1200 1000 800 600 400 200 0 0 5 10 15 20 25 30 35 40 45 50 55 60 Time [s] (a) (b) Figure 7: (a) The prescribed nodes at the bottom of the weld representing the nodal line stretching along the entire weld, in red. (b) Thermal cycle based on Kim [3]. 4.1.1 Boundary conditions The boundary conditions implemented for the welding simulation were such that a node at the bottom of the quarter-plane of symmetry was prescribed in the T-direction, see Figure 8a, and the node at the centre bottom of the free end of the plate was prescribed in the TT-direction which can be seen in Figure 8b. The entire plane of symmetry was prescribed in the L-direction as seen in Figure 8c. The heat flux was prescribed on the entire surface of the test specimen, with a convection coefficient of h = 1000 × 10−6 W/°C mm2 , except at the plane of symmetry, which was insulated. Initially, a film coefficient used by [14] of h = 20 × 10−6 W/°C mm2 was used, but a higher coefficient gave residual stress values closer to similar studies by [3] and [16]. (a) Node fixed in the T- (b) Node fixed in the TT- (c) Nodes fixed in the Ldirection. direction. direction. Figure 8: Prescribed displacements for the thermomechanical weld simulation. 5 This is just below the austenite and δ-ferrite phase transformation temperature. 8 4.1.2 Numerical results Figure 9 shows the residual von Mises stress field after cooling to room temperature shown side-byside alongside compared with [3]. Although, [3] did not provide a color legend, it can be observed that the stress fields are similar with respect to the stress variation around the surface. It can be seen that the same areas, circled in yellow, exhibit high stresses in both the simulation and the results by [3]. These areas do not directly affect the weld toe, but these similarities in the overall model further validate the simplification. In the simulation there were additional areas of high stress in the flange. In the stiffener, a diagonal stress gradient is observed where the lowest values are found at the top, circled in red. Figure 9: Comparison of the von Mises stress distribution after cooling to room temperature with the result to the left and compared to Kim [3], without legend, to the right [Pa]. High stresses are circled in yellow and low stresses are circled in red. In the longitudinal direction, Figure 10, high stresses are found in the circled area around the stiffener, but more specifically also around the weld toe. These tensile stresses at the weld toe are the source of the crack development due to fatigue loading. It can also be seen that there is a high compressive stress in the flange, as opposed to [3], seen in blue. It should be noted that the legend was scaled with respect to the surface stresses, with stresses above these levels located inside the model. Both Figures 9 and 10 are representative of the expected outcome of the thermomechanical analysis, with high tensile residual stresses in the vicinity of the weld. Figure 10: Comparison of the longitudinal stress distribution after cooling to room temperature with the result to the left and compared to Kim [3], without legend, to the right [Pa]. High stresses are circled in yellow. Figures 11a and 11b show the longitudinal stress along paths A and B after welding. Along Path A, where the stress is compared to numerical results from [3] and [11], it can be observed that the largest stress is located at the weld toe and decreases with increased distance from the weld. The results are similar to, [3] but differences are seen in comparison to [11], which has a lower and as 9 such a different stress behaviour. It is important to note that [11] used the standard SYSWELD library and did not state the yield stress used. However, the minimum expected 355 MPa, if used by [11], is higher than the 315 MPa from [3] and this reduces the ability to make a direct comparison of the results. Nevertheless, the as-welded values found along Path B are compared to numerical results from [3] and to experimental X-ray measurements in [7]. The study by [7] used a larger number of data points compared to the FE simulation and the FE can not therefore capture the full behaviour of the longitudinal stress along Path B, especially at the weld toe. However, it is in the same range as [7], especially after about 0.2 mm in depth and matches [3] closely. 500 400 400 AW Kim AW Stoschka AW FEA 300 Longitudinal residual stress [MPa] Longitudinal residual stress [MPa] 350 250 200 150 100 50 0 300 200 100 0 -100 -200 AW Kim AW Exp. Leitner AW FAE -300 -400 -50 -500 0 5 10 15 20 25 0 Distance from weld toe (Path A) [mm] 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Depth from weld toe (Path B) [mm] (a) Longitudinal as-welded residual stress versus distance from the weld along Path A compared to results from Kim [3] and Stoschka [11]. (b) Longitudinal as-welded residual stress versus distance from the weld along Path B compared to results from Kim [3] and Leitner [7]. Figure 11: Resulting residual stresses after welding along studied paths from the weld toe. The purpose of the thermomechanical welding simulation was to generate a representative tensile residual stress field. The results are similar to those presented by [3] in the benchmark exercise, and with the exception of the zone closest to the weld toe, a reasonable match to the FEA in [11] and the experimental results of [7] was obtained. Thus, even with a simplified thermomechanical analysis, i.e., of a simultaneous thermal cycle around the whole weld, the full residual stress field was determined to be suitable for use as the start condition in the HFMI analysis. 4.1.3 Benchmark issues and error sources There were significant simplifications made to the benchmark exercise welding process. The process was simulated by heating the entire weld simultaneously, as opposed the implementation of a heat source going around the weld joint. The simulations also did not reach temperatures at which phase transformation takes place. In addition, due to the use of forward-coupling, there was no contribution from heating during plastic deformation, although they would be expected to be small. Finally, the weld geometry provided by the benchmark is symmetric, which is unlikely in the case in a real welded joint. 4.2 Simulation of HFMI treatment After the thermomechanical weld simulation, the residual stress fields were imported into the FE model of the HFMI treatment, which consisted of an implicit quasi-static contact analysis with an indentation tool (in four parts) that pressed simultaneously and deforms the entire weld toe area6 . The purpose of the simulation of the HFMI treatment was to introduce compressive residual stresses in the vicinity of the weld toe. 6 Several different indentation tools were developed and tested. The four-part tool was selected because it enabled an indentation angle of 75° as in [3]. 10 The solver setup consisted of an initial step that defined the imported residual stress field and a contact step with a total step time of 8 seconds. The initial increment size was set to 0.01 seconds, with a minimum increment size of 10−6 seconds and a maximum of 0.1 seconds. The indentor was moved into contact with the weld for the first five seconds of the step and in the last three, it was returned to its original position. The numerical analysis included both material nonlinearity through the Chaboche hardening model as well as geometrical nonlinearity due to large deformations. Heating due to plastic deformation during the HFMI treatment was not considered. 4.2.1 Test specimen boundary conditions Essential boundary conditions on the test specimen were chosen to eliminate rigid body modes and to avoid bending during the indentation. Figure 12 shows the location of the nodes with certain displacements prescribed to zero. (a) Nodes fixed in the L- (b) Node fixed in the T- (c) Nodes fixed in the TTdirection. direction. direction. Figure 12: Prescribed displacements for the HFMI simulation. 4.2.2 Indentation tool geometry and boundary conditions The geometry of the four-part indentation tool was developed in CATIA V5 [21] and then exported to ABAQUS, where it was modelled with rigid 3D 4-node elements (R3D4) with an average element size of 0.5 mm. This mesh size matched the smallest elements in the contact zone of the slave surface. Figure 13a shows the geometry and positioning of the four-part indentation tool while Figure 13b shows the profile of an indentation tool, with a tip radius of 2 mm as in [3] and [7]. Figure 13c shows the path along which the four-part indentation tool indent the weld. In particular, note that the 75° indentation angle used in [3] was implemented. (a) The four-part indentation tool used for simulation of the HFMI treatment. (b) Profile of the indentation tool [mm]. (c) Illustration of the indentation path. Figure 13: Positioning, profile geometry and indentation path of the four-part indentation tool. The four-part indentation tool had prescribed axial displacements of 0.15 mm on reference points defined on the mid-plane of the indentation tools. This displacement was chosen as a mean value 11 of the indentation depths used by [3] and [7]. Local cylindrical coordinate systems were defined for the two bends along the weld, in order to achieve radial displacement for the curved indentation tools. To ensure that the indentation tools only moved in the prescribed direction, all other degrees of freedom on the reference points were set to zero. 4.2.3 Contact properties The contact interaction between the four-part indentation tool and the test specimen was defined as surface-to-surface contact with the test specimen as the slave surface and the four-part indentation tool as the master surface with no interaction between the four parts themselves. The tangential behaviour was governed by a friction penalty with a friction coefficient µ = 0.15 as in [3] while the normal-to-surface behaviour was set as hard contact on default setting in ABAQUS. 4.2.4 Numerical results The longitudinal residual stress field after the HFMI treatment simulation in the vicinity of the weld can be seen in Figure 14a. The compressive residual stresses are localized to the area where the weld meets the plate surface, which is the expected outcome of the HFMI treatment. Figure 14b shows the resulting plastic deformation after the HFMI treatment simulation in a cross section at the center of the weld toe. The indentation results in a smoother corner, effectively reducing the stress concentration factor at the weld toe. (a) Overview of the resulting longitudinal stress field after simulated HFMI treatment in the vicinity of the weld [Pa]. (b) Plastic deformation at the cross section of the weld toe after simulated HFMI treatment [mm]. Figure 14: Longitudinal residual stress fields after applied HFMI treatment. Figures 15a and 15b shows the longitudinal residual stress along Path A and Path B after welding and after HFMI treatment. Along Path A, it can be seen that the tensile stress field from the thermomechanical weld simulation was completely eliminated to a distance of approximately 10 mm from the weld toe. Note that [3] did not measure the residual stress along path A after HFMI treatment. Along Path B, the results are compared to numerical results from [3]. It should be noted that [3] uses a HFMI indentation depth of 0.1 mm, compared to the 0.15 mm indentation depth used in this analysis. It can be seen that the surface residual stress after HFMI treatment is approximately the same, however the compressive residual stress field from the current study is limited to a depth of around 2.25 mm. 12 500 500 AW FEA (Kim) AW FEA HFMI-treated FEA 400 Longitudinal residual stress [MPa] Longitudinal residual stress [MPa] 400 300 200 100 0 -100 -200 -300 -400 -500 300 200 100 0 -100 -200 AW FEA (Kim) AW FEA HFMI-treated FEA HFMI-treated FEA (Kim) -300 -400 0 5 10 15 20 -500 25 0 1 Distance from weld toe [mm] 2 3 4 5 Depth from weld toe surface [mm] (a) Longitudinal residual stress versus distance from the weld along Path A compared to results from Kim [3]. (b) Longitudinal residual stress versus distance from surface along Path B compared to results from Kim [3]. Figure 15: Resulting residual stresses after HFMI treatment along studied paths from the weld toe. Figure 16 shows the longitudinal residual stress along Path B with a comparison to experimental X-ray measurements made in [7] which had more data points into the depth of the test specimen compared to the FEA in the current study. It can also be seen that [7] measured negative values at the surface after welding. This means that the residual stress field along Path B is more complex than what the FEA in the current study was able to capture. Overall, the longitudinal stress field after HFMI treatment behaves impressively similar to the experimental data. 500 Longitudinal residual stress [MPa] 400 300 200 100 0 -100 -200 AW FEA AW Exp. (Leitner) HFMI-treated FEA HFMI-treated Exp. (Leitner) -300 -400 -500 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Depth from weld toe surface [mm] Figure 16: Longitudinal residual stress versus distance from surface along Path B compared to experimental X-ray measurements by Leitner [7]. The purpose of the simulated HFMI treatment was to generate a representative local compressive residual stress field in the vicinity of the weld toe for the simulation of cyclic loading. Although the simulated HFMI treatment was simplified, the results agree remarkably well with experimental data and other published FE results. The results could therefore be considered representative of the HFMI treatment process. 4.2.5 Benchmark exercise issues The mesh of the slave surface in the contact zone is too coarse, which can lead to penetration errors if a node-to-node contact interaction is used. With surface-to-surface contact, as used in the current study, penetration errors are in general not an issue. However, using surface-to-surface contact is more computationally expensive compared to using node-to-node contact [8]. Refining the mesh in the contact zone could allow the use of node-to-node contact which would decrease the computational time. 13 During the process of finding a suitable indentation tool, a few distorted elements were discovered inside the stiffener. The residual stress values of these particular elements were unrealistically high for some of the tested tools. Even though this was not experienced in the case of the four-part tool that was ultimately used, it should still be considered as an error source since these particular elements were known to behave badly. 4.2.6 Error sources During the contact, the effective von Mises stress in the contact zone was approximately of 1 GPa. The hardening behaviour of S355 with the Chaboche model using the determined parameter set, as shown in Figure 5b, showed hardening up to around 600 MPa before saturating. This raises questions about the applicability of the Chaboche model during the contact, since it is likely that the material would fail before reaching the effective stresses seen during the contact. This problem is also seen in [7], where the von Mises stresses during contact reach approximately 1.1 GPa. The four-part indentation tool was approximated as a rigid body and does not deform. Although a real HFMI pin is stronger than the material it impacts, there will be some elastic deformation present in the HFMI pin. The deformations are small, but the approximation used in this analysis still introduces an error in the solution. Furthermore the contact between the indentors and the benchmark model is defined as hard contact, which means there is no limit to the amount of contact pressure during the contact, and that the two objects separate when the contact pressure is equal to zero. In reality, there might be some cohesiveness between the HFMI pin and the test specimen. Therefore, the hard contact condition gave rise to another error source in the solution. 4.3 Simulation of cyclic loading Originally the plan was to run the simulation of cyclic loading as a separate FEA. Instead, the simulation of cyclic loading was coupled with the HFMI treatment in order to incorporate the material hardening and improved weld toe topology along with the compressive residual stress field generated by the HFMI treatment. 4.3.1 Boundary conditions For the cyclic loading analysis, the plane of symmetry was prescribed to zero displacement in the L-direction due to the symmetric boundary condition, as shown in Figure 17a. To eliminate rigid body modes and bending in both tensile and compressive loading, two nodes as shown in Figure 17b were prescribed to zero displacement in the T- and TT-directions and the node on the bottom on the free end was prescribed to zero displacement in the TT-direction, as shown in Figure 17c. The nodes in Figure 17b and 17c are located on the quarter-symmetry plane. (a) Nodes fixed in the Ldirection. (b) Node fixed in the T- and TT-directions. (c) Node fixed in the TTdirection. Figure 17: Prescribed displacements for the fatigue simulation. 14 4.3.2 Load cases To investigate the stability of the compressive residual stress field induced by HFMI treatment, a cyclic load was applied in the L-direction on the entire area of the free end. Unfortunately, the specific cyclic loading cases required for the benchmark exercise have not yet been determined by the specialist committee (ISSC). Instead, the load cycles were based on the load cases in [2] with input from the client [6]. The load cases consist of 20 constant-amplitude (CA) cycles and one over-and-under load (OU) cycle. Both the CA and OU cycles had a time period of two seconds. The OU cycle was set as the first cycle in the first load case and as the last cycle in the second load case. The stress ratio for both load cycles was R = −1. The load amplitudes, denoted σA , that were used can be seen in Table 2 and the shape of the cycles can be seen in Figure 18. The cyclic loading was simulated in two steps, where one step was used for the OU cycle and one for the CA cycle. The duration of the OU step was 2 seconds while the duration of the CA step was 40 seconds. To capture the cyclic behaviour and to obtain even step increments, the initial and maximum increment size was set to 0.1 seconds. The minimum number of increments was 422. The analysis involved the same nonlinearity as the simulation of HFMI treatment, i.e. material nonlinearity through the Chaboche model and geometrical nonlinearity due to large deformations. Table 2: Load amplitudes and stress ratios for the CA and OU for both cases. CA 37.5 -1 Load amplitude, σA [MPa] Stress ratio, R 125 125 Applied load, OU-first Applied load, OU-last 100 100 75 75 50 50 Applied load [MPa] Applied load [MPa] OU 106.5 -1 25 0 -25 -50 25 0 -25 -50 -75 -75 -100 -100 -125 -125 0 4 8 12 16 20 24 28 32 36 40 0 Time [s] 4 8 12 16 20 24 28 32 36 40 Time [s] (a) Load case with the OU load cycle as the first cycle. (b) Load case with the OU load cycle as the last cycle. Figure 18: Load cases showing amplitudes and OU loads. 4.3.3 Numerical results Figure 19 shows the longitudinal residual stress along Paths A and B after cyclic loading compared to the longitudinal residual stress after the simulated HFMI treatment. Along both paths, the longitudinal residual stress after the cyclic loading is nearly identical for both load cases. Placing the OU load at the start or the end of the load cycle makes only a small difference in the longitudinal residual stress. It is interesting to note that on the surface along Path A, the compressive residual stress produced by the simulated HFMI treatment, was completely eliminated. However, along Path B, some of the compressive residual stress field was retained. Along Path A the stress shifts from compressive to tensile twice in the first nodes before stabilising to approximately 20 MPa further away from the weld toe, as shown in Figure 19a. The cyclic loading reduced both tensile and compressive stresses along Path B, as such, the shift from com- 15 pressive to tensile stress remains at the same depth, approximately 2.25 mm, as shown in Figure 19b. Furthermore, the overall stress field is more uniform as compared to before cyclic loading. 400 Longitudinal residual stress [MPa] Longitudinal residual stress [MPa] 200 100 0 -100 -200 -300 OU-last, FEA OU-first, FEA HFMI-treated FEA -400 300 200 100 0 -100 -200 OU-last, FEA OU-first, FEA HFMI-treated FEA -300 -400 0 5 10 15 20 25 0 1 Distance from weld toe [mm] 2 3 4 5 Depth from weld toe surface [mm] (a) Longitudinal residual stress along Path A. (b) Longitudinal residual stress along Path B. Figure 19: Longitudinal residual stress after cyclic loading along studied paths from the weld toe. Figure 20 shows the changes in the longitudinal stress at the weld toe after each cycle for the first node in Path A. After the simulated HFMI treatment, at 0 load cycles, the longitudinal residual stress in the first node on surface is -309 MPa. The longitudinal residual stress after the last load cycle for the OU-first case is -23 MPa and -27 MPa for the OU-last case. These changes are compared in Figure 20 with experimental data from [7]. The shape of the stress curve for the experimental data should be compared with the OU-first load case, since [7] used a CA that is more comparable to the OU load used in the current study. It can be seen that the curves follow a similar behaviour for the first load cycle. The experimental data showed that much of the positive effect of the HFMI treatment is kept after 5 load cycles, while the FEA data showed that most of the positive effect is eliminated. This was expected, since the OU load in the current study was higher than the load amplitude used in the experimental study. The results from the OU-last case were however unexpected, since the CA is in the elastic range of the material. Therefore, the longitudinal residual stresses should not have changed, as there was no plastic deformation of the test specimen. 50 Longitudinal stress [MPa] 0 -50 -100 -150 -200 -250 -300 CA, exp (Leitner) OU-last, FEA OU-first, FEA -350 -400 0 3 6 9 12 15 18 21 Cycles Figure 20: Longitudinal stress in the node at the weld toe versus number of load cycles, compared with experimental data from Leitner [7]. Figure 21 shows the longitudinal stress-strain response for the two load cases. Since the cyclic loading simulation is coupled with the simulation of HFMI treatment, the longitudinal strain does 16 300 300 200 200 Longitudinal stress [MPa] Longitudinal stress [MPa] not start from 0% but 4.34%. In Figure 21a the OU load is applied first in the load cycle which can be seen in the figure since the material yields first in tension and then in compression, before it goes into linear elastic response from the CA cycles. When the OU load is placed last in the load cycle, as in Figure 21b, the material yields a small amount in each CA cycle which results in a hardening effect that pushes the curve upwards. The extremes in the curves are due to the OU cycle giving a large stress-strain response. 100 0 -100 -200 -300 OU-first, FEA Start Stop -400 -500 4.25 4.3 4.35 4.4 4.45 4.5 4.55 100 0 -100 -200 -300 OU-last, FEA Start Stop -400 -500 4.25 4.6 Longitudinal strain [%] 4.3 4.35 4.4 4.45 4.5 4.55 4.6 Longitudinal strain [%] (a) OU load in the beginning of the cycle. (b) OU load in the end of the cycle. Figure 21: Stress-strain response in the node at the weld toe for the two different load cycles. Figure 21b shows that the elastic range was shorter than what was expected for the Chaboche model, which should be at least ±435 MPa, if not higher, due to local hardening effects. And indication of the Bauschinger effect can also be observed, both in the compressive and tensile load peaks in the CA load interval, even though the longitudinal stress interval in the studied node is only around 200 MPa. This indicates that stress components in the other directions had larger magnitude than expected, and may have contributed to the yielding. When the material yields in tension on the OU-first case the strain decreases, as shown in Figure 21a. One concerning anomaly seen in Figure 21a is that the Young’s modulus appears to have changed since the slope of the stress-strain curves is not the same throughout the OU load. Although S355 steel is known to be weaker than materials previously tested in other studies, it was not expected that it would start to yield at such a low load amplitude. In the OU-last case, the loss of longitudinal compressive residual stresses was almost linear for CA loading with an amplitude of σA = 37.5 MPa, which is significantly below the yield limit. Since the width of the plate reduces from 80 mm on the free end to 50 mm at the weld toe, the stress amplitude is increased by a factor 80/50 = 1.6, resulting in a nominal stress amplitude of 1.6 · σA = 60 MPa and as such the material was not expected to yield. The unexpected results may be due to peaks in stress concentration in the node closest to the weld toe, but may also be due to numerical errors as the node is on the edge of the element and in a sharp corner. To resolve this, it would have been necessary to investigate the implementation of the material model and which boundary conditions should be used to properly simulate cyclic loading on plates with welded stiffeners. However time did not permit a detailed study of these issues. 17 4.3.4 Benchmark exercise issues Unfortunately, there were no load cases defined in ISSC benchmark exercise, which made it difficult to find a reasonable load level. In [2], different loads were applied to the edge of the plate, while [7] presented loads as a nominal stress range at the weld toe. This issue needs to be targeted, since this study and the work by [22] shows that the HFMI-induced stress field is sensitive to peak loads. In addition, another benefit would be that other results in this benchmark exercise could be more easily compared if different loads for different materials were better defined. 4.3.5 Error sources The mesh used in the benchmark is rather coarse which may reduce the ability to capture rapidly changing gradients [23]. The use of linear hexahedral (C3D8), linear wedge (C3D6) elements and full integration, which were used in all of the simulations, can cause problems in the form of shear locking if any unexpected bending occurs. In FEA, geometries are implemented as discrete connected points, which can lead to sharp corners and edges. This is something that is not normally exhibited in physical reality. Corners with a zero radius result in non-realistic stresses and plasticity occurring in these locations [23]. The results in the simulation are computed at the Gauss points in each respective element of the mesh, while values at other positions have to be extrapolated. When extrapolation is done at an edge of an element on the edge of a structure, where stress concentrations are likely to occur and are of high interest, it can lead to significant errors in the results when the stress field changes rapidly if the mesh density is too low [23]. This may prove to be a problem when evaluating data from a single node on a sharp part of the model, as was done in the simulation of cyclic loading when plotting the results from the node closest to the weld toe in Figure 20. 18 5 Discussion The progressive FEA consisted of a thermomechanical analysis of the welding process, followed by a quasi-static coupled analysis of HFMI and cyclic loading. The purpose of the progressive FEA was to determine if the compressive residual stress field induced by HFMI treatment on the low-strength steel S355 was stable under cyclic loading. In order to evaluate this, two load cases with 20 load cycles were considered. The evaluation of the stability was made by analysing the longitudinal residual stresses along two specific paths as a part of the benchmark exercise, Path A and Path B. The purpose of the thermomechanical welding simulation was to generate a representative tensile residual stress field compared to available experimental and FEA data. A simplified thermal cycle was used, where the weld was heated and subsequently cooled in order to generate thermal strains and a residual stress field in the test specimen. The residual stress results of the thermomechanical welding simulation corresponded well to available experimental and FEA data, and therefore it was deemed to be suitable for use in the HFMI simulation. Furthermore, it was shown that a simplified method for simulating the welding process provided representative results. When comparing to similar studies with more advanced simulation methodology, the results showed that it may not be necessary to simulate the entire welding process in order to obtain representative residual stress fields. Using the tensile residual stress field from the thermomechanical weld simulation as an initial condition, the simulation of HFMI treatment resulted in a representative compressive residual stress field in the vicinity of the weld toe. Furthermore, the compressive residual stress field was achieved with much simpler methods than used by similar simulations of HFMI treatment. It was found that the residual stress field after simulated HFMI treatment was in good agreement with both experimental data and other FEA data. Other effects of the simulated HFMI treatment were the local hardening of the weld toe, as well as improved weld toe topology, e.g. a smoother transition from weld to plate, resulting in a decreased stress concentration factor. The compressive residual stress field from the HFMI simulation was used as initial condition for the simulation of cyclic loading. In order to incorporate local hardening effects and weld toe topology improvements, the cyclic loading simulation was added as a subsequent step to the HFMI simulation. Two load cases were considered: the first load case was a constant-amplitude load cycle with a OU load as the first cycle while the second load case had the OU load as the last cycle. It was found that the compressive residual stress field was unstable regardless of the OU load occurring in the first or last cycle. The residual stress field after 20 load cycles was still compressive, however most of the positive effect from the HFMI treatment had been eliminated. The behaviour of the residual stress field under constant-amplitude loading was unexpected, as the material appeared to yield despite using load amplitudes much lower than the yield limit of the material. Because of these unexpected results, it was not possible to make a definitive statement on the stability of the compressive residual stress field subjected only to constant-amplitude cyclic loading. A possible reason for the unexpected results from the cyclic loading simulation is the lack of sufficient refinement of the computational mesh in critical regions such as the weld toe. Furthermore, it is not known whether the unexpected yield behaviour during cyclic loading was due to the Chaboche material model or other factors, such as ill-defined essential boundary conditions. Some bending was observed during the cyclic loading, whereas the test specimen should have been subjected to pure membrane loading with no bending. For this reason, it is important to evaluate which boundary conditions should be used in the benchmark exercise, as none were provided. 19 6 Conclusions and contributions The main conclusions of the current study relates to the goal of the project, i.e. determining the stability of the residual stress field induced by HFMI treatment after cyclic loading. • The compressive residual stress field induced by the HFMI treatment is very sensitive to cyclic loading with unexpected peak loads. Most of the positive HFMI effect is eliminatedn if the load cycle begins or ends with a peak OU load. The significance of this result is that the industrial use of the HFMI treatment methods must take into consideration unexpected peak loads that can completely remove the beneficial effects of the HFMI treatment when using S355. • It is not possible to draw any conclusions on the stability of the compressive residual stress field after constant-amplitude cyclic loading for S355 because of the unexpected nature of the results. There are also secondary conclusions, not related to the goal of the project, that contribute to the benchmark exercise. In particular, the methods used in the thermomechanical welding simulation and the simulation of HFMI treatment are of interest to the benchmark exercise as they are much simpler compared to methods used in similar studies. • Results from the thermomechanical weld simulation and the simulated HFMI treatment show that a simplified approach can give representative results compared to experimental and FEA data in the available literature. Comparing to similar studies with more advanced methods, such as in [3] and [7], it can be seen that the results are very similar. Utilising the simplified techniques from this project could allow for more efficient simulations as a part of the benchmark exercise. 6.1 Recommendations and future work There are some issues with the benchmark model which may cause unrealistic results that should be resolved before proceeding with further benchmark studies. As seen when comparing experimental data from [7] along Path B, the mesh in the TT-direction at the weld toe is too coarse to be able to represent the complex stress behaviour near the surface. It is also recommended that the material properties to be used in future studies are better defined. The current study and previous studies have used different material parameters in the different stages, which can lead to different results. For example, it was shown in [3] that there is a significant difference between isotropic hardening, kinematic hardening and strain-rate dependent hardening for the residual stress field after HFMI treatment. The mesh on the weld surface is too coarse for a proper contact analysis, therefore there exists a possibility of penetration errors, given that node-to-node contact is used. It is therefore recommended to use surface-to-surface contact if the current mesh is to be used in further studies in the benchmark exercise. As for normal-to-surface interaction behaviour, it is recommended that future work investigate if choosing another type of contact would benefit the results. One previous HFMI study by [19] used a finite stiffness of 250 kN/m instead of hard contact. Using this instead of hard contacts could possibly lead to more representative results and faster convergence. 20 References [1] L. Suominen, M. Khurshid, and J. Parantainen. Residual stresses in welded components following post-weld treatment methods. Procedia Engineering, 66:181–191, 2013. https://doi.org/10.1016/ j.proeng.2013.12.073. [2] E. Mikkola, H. Remes, and G. Marquis. A finite element study on residual stress stability and fatigue damage in high-frequency mechanical impact (HFMI)-treated welded joint. International Journal of Fatigue, 94:16–29, 2016. http://dx.dol.org/10.1016/j.ijfatigue.2016.09.009. [3] M. Hyun Kim. Numerical simulation of high frequency mechanical impact (HFMI) treatment in fillet welded joints. Meeting Presentation to the Specialist Committee on Materials and Fabrication Technology (ISSC 2018), Oslo, January 31, 2017. [4] Dassault Systèmes. Abaqus/CAE, Version 6.14, Accessed March 31, 2017. https://www.3ds.com. [5] M. Leitner, M. Khurshid, and Z. Barsoum. Stability of HFMI-treatment induced residual stresses during fatigue. Document no. XIII-2633-16, International Institute of Welding, Paris, 2016, 9 pgs. [6] L. Josefsson. Personal communication. Gothenburg, March 23, 2017. [7] M. Leitner, M. Khurshid, and Z. Barsoum. Stability of high frequency mechanical impact (HFMI) post-treatment induced residual stress states under cyclic loading of welded steel joints. Engineering structures, 143:589–602, 2017. http://dx.doi.org/10.1016/j.engstruct.2016.04.013. [8] Dassault Systèmes. Abaqus Theory Manual, Version 6.12, Accessed May 11, 2017. https://www.3ds. com. [9] N. Gilbert. Structural Steel - S235, S275, S355 Chemical Composition, Mechanical Properties and Common Applications. Online, Accessed May 16, 2017. http://www.azom.com/article.aspx? ArticleID=6022. [10] MEADinfo. Material Properties of S355 Steel - An Overview. Online, Accessed May 14, 2017. http://www.meadinfo.org/2015/08/s355-steel-properties.html. [11] M. Stoschka, M.J. Ottersböck, and M. Leitner. Integration of phase-dependent work-hardening into transient weld simulation. Civil-Comp Press, 39, 2014. http://dx.doi.org/10.1016/j.engstruct. 2016.04.013. [12] M. Khurshid, M. Leitner, Z. Barsoum, and C. Schneider. Residual stress state induced by high frequency mechanical impact treatment in different steel grades - Numerical and experimental study. International Journal of Mechanical Sciences, 123:34–42, 2017. http://dx.doi.org/10.1016/j. ijmecsci.2017.01.027. [13] ESI Group. SYSWELD Toolbox, Version 16, 2014. [14] A. A. Bhattia, Z. Barsouma, H. Murakawac, and I. Barsouma. Influence of thermo-mechanical material properties of different steel grades on welding residual stresses and angular distortion. Materials & Design (1980-2015), 65:878–889, 2015. [15] J.L. Chaboche. A review of some plasticity and viscoplasticity constitutive theories. International Journal of Plasticity, 24:1642–1693, 2008. http://dx.doi.org/10.1016/j.ijplas.2008.03.009. [16] M. Leitner, D. Simunek, S.F. Shah, and M. Stoschka. Numerical fatigue assessment of welded and HFMI-treated joints by notch stress/strain and fracture mechanical approaches. Advances in Engineering Software, 2016. http://dx.doi.org/10.1016/j.advengsoft.2016.01.022. [17] A.M.P. de Jesus, R. Matos, B.F.C. Fontoura, C. Rebelo, L.S. da Silva, and M. Veljkovic. A comparison of the fatigue behaviour between S355 and S690 steel grades. Journal of Constructional Steel Research, 79:140–150, 2012. http://dx.doi.org/10.1016/j.jcsr.2012.07.021. [18] D. Forni, B. Chiaia, and E. Cadoni. Strain rate behaviour in tension of S355 steel: Base for progressive collapse analysis. Engineering structures, 119:164–173, 2016. http://dx.doi.org/10.1016/j. engstruct.2016.04.013. [19] J. Foehrenbach, V. Hardenacke, and M. Farajian. High frequency mechanical impact treatment (HFMI) for the fatigue improvement: numerical and experimental investigations to describe the condition in the surface layer. Weld World, 60:749–755, 2016. http://dx.doi.org/10.1007/ s40194-016-0338-4. [20] The MathWorks, Inc. MATLAB, Version/2016b, Accessed May 11, 2017. https://mathworks.com/ products/matlab.html. [21] Dassault Systèmes. CATIA V5, Accessed May 11, 2017. https://www.3ds.com. [22] E. Mikkola and H. Remes. Allowable stresses in high-frequency mechanical impact (HFMI)treated joints subjected to variable amplitude loading. Welding World 61: 125, January 2017. http://doi:10.1007/s40194-016-0400-2. [23] D.J. Grieve. Errors arising in FEA. Online, Accessed May 13, 2017. http://www.tech.plym.ac.uk/ sme/mech335/feaerrors.htm. 21 Acknowledgements We would like to thank our supervisor, Moyra McDill, Prof. Em. (Carleton University), for her invaluable direction and help during the course of the project. We would also like to thank Lennart Josefsson, Prof. (Chalmers University of Technology), for his excellent guidance and technical knowledge. Human finger response to high frequency vibrations Wave propagation in human tissue from transient vibrations Report in Applied mechanics RON GEORGE NIKHIL MASINETE ANDREAS NÄKNE ALEXANDER OLSSON LINDA PIPKORN RISHAB RANGARAJAN Department of Applied Mechanics Project in Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2017 REPORT IN APPLIED MECHANICS Human finger response to high frequency vibrations Wave propagation in human tissue from transient vibrations RON GEORGE NIKHIL MASINETE ANDREAS NÄKNE ALEXANDER OLSSON LINDA PIPKORN RISHAB RANGARAJAN Department of Applied Mechanics Project in Applied Mechanics CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2017 Human finger response to high frequency vibrations Wave propagation in human tissue from transient vibrations RON GEORGE NIKHIL MASINETE ANDREAS NÄKNE ALEXANDER OLSSON LINDA PIPKORN RISHAB RANGARAJAN c RON GEORGE, NIKHIL MASINETE, ANDREAS NÄKNE, ALEXANDER OLSSON, LINDA PIPKORN, RISHAB RANGARAJAN, 2017 Department of Applied Mechanics Project in Applied Mechanics Chalmers University of Technology SE-412 96 Gothenburg Sweden Telephone: +46 (0)31-772 1000 Cover: The deformation and strain in the finger model. Chalmers Reproservice Gothenburg, Sweden 2017 Human finger response to high frequency vibrations Wave propagation in human tissue from transient vibrations Report in Applied mechanics RON GEORGE NIKHIL MASINETE ANDREAS NÄKNE ALEXANDER OLSSON LINDA PIPKORN RISHAB RANGARAJAN Department of Applied Mechanics Project in Applied Mechanics Chalmers University of Technology Abstract ISO 5349 claims that vibrating tools with frequencies in the range of 8-1250 Hz must be regulated. However, frequencies higher than 1250 Hz can also cause damage to nerves and blood vessels. There is thus a discrepancy between the regulated safe work and caused damage. Consequently, every year approximately 350 000 workers in Sweden are affected by hand-arm vibration syndrome. The research institute, Swerea IVF, is proposing a change of the ISO-standard with the aim of proving that high-frequency vibrations are causing damage to the fingers. As a part of the aim to change the ISO-standard, this project was investigating how a change of different load cases as well as skin layer geometry and material stiffness was affecting the strain and strain rate. This was done by performing 60 simulations with a 2D cross-section finite element model of a finger in LS Dyna. The evaluation of the strain and strain rate for the different analysis resulted in the following conclusions; acceleration frequency was having a greater impact than the acceleration amplitude, the shear modulus of stratum corneum was more important than the shear modulus of epidermis and soft tissue, no trend was found regarding the thickness of the skin layers, a change of bulk modulus and decay constant was showing no significant effect and finally the angle of vibration mattered for both distribution and magnitude. The project showed that high-frequency vibrations were resulting in maximum strain and strain rates near the fingerprint where nerves are present. In order to show that this is actually resulting in HAVS, damaging magnitudes of strain and strain rates is needed. That was recommended as future work. Keywords: HAVS, vibrant work, high-frequency, ISO 5349 i Nomenclature Abbrevations HAVS Hand-arm vibration syndrome TTS Temporal Threshold Shift Explanations Epidermis The outer skin layer Stratum corneum The outermost part of the epidermis, containing dead skin cells Soft tissue Everything else in the finger model which is not stratum corneum and epidermis ISO 5349 The international standard which limits usage of tools over a certain time interval node A, B, C, D Nodes selected within the finger model for our analysis Transient A sudden change of energy in a system resulting in a short-lived ”spike” response Variables β Decay constant θ Direction of resultant acceleration with respect to the finger base ǫ Strain ǫ̇ Strain rate ǫx Strain in x-direction ǫy Strain in y-direction γ In-plane shear strain. amax Amplitude of the applied vibration amax,x Amplitude of the acceleration parallel to the finger base. amax,y Amplitude of the acceleration perpendicular to the finger base. T Time period of the vibration t Simulation time Ev Voight Shear Modulus Er Reuss Shear Modulus Evrh Voight Reuss Hill Shear Modulus vfi Volume fraction of each tissue G(t) effective shear modulus G0 effective shear modulus G∞ effective shear modulus ii Contents Abstract i Nomenclature ii Contents iii 1 Introduction 1.1 Background . . . . . . . . . . . . 1.2 Project description . . . . . . . . 1.2.1 Project goal . . . . . . . . . . . 1.2.2 Provided material . . . . . . . 1.2.3 Assumptions and delimitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 2 2 2 2 2 Theory 2.1 ISO 5349 . . . . . . . . 2.2 HAVS . . . . . . . . . . 2.2.1 Vascular disorders . . 2.2.2 Neurological disorders 2.3 Finger model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 3 3 3 4 3 Methodology 3.1 Literature study . . . . . . . . 3.2 Pre-processing . . . . . . . . . 3.2.1 Applied loads . . . . . . . . . 3.2.2 Material parameters . . . . . 3.2.3 Geometry . . . . . . . . . . . 3.3 Simulations . . . . . . . . . . . 3.3.1 Main simulation scheme . . . 3.3.2 Extended simulation scheme 3.4 Post-processing . . . . . . . . . 3.4.1 Data collection . . . . . . . . 3.4.2 Data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 5 5 5 5 6 7 7 7 8 8 9 . . . . . . . . . . rate distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 . . 11 . . 11 . 13 . 13 . 13 . 15 . 17 . 17 . 18 . 18 . 19 . 20 4 Results and discussion 4.1 Literature study results 4.2 Overall strain and strain 4.2.1 Node D . . . . . . . . 4.3 Main simulation . . . . 4.3.1 Load . . . . . . . . . . 4.3.2 Shear modulus . . . . 4.4 Extended simulation . . 4.4.1 Geometry . . . . . . . 4.4.2 Bulk modulus . . . . . 4.4.3 Decay constant, β . . 4.4.4 Biaxial vibration . . . 4.5 General discussion . . . . . . . . . . . . . . . . . . 5 Conclusion and recommendations 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Recommendations for future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 22 22 A Simulation scheme with data 23 iii 1 Introduction This project is part of the course TME131 - Project in Applied Mechanics given at the master’s program in Applied Mechanics at Chalmers University of Technology during spring 2017. The project was proposed by the research institute Swerea IVF and is a part of their work towards a change of the ISO 5349 standard to reduce the number of work-related hand-arm vibration injuries. 1.1 Background Many industrial workers spend significant time on handling hand-held power tools. The vibration from these tools, which often are of transient nature, can cause injury to the user. This issue has a few names like white finger syndrome and Raynaud syndrome, but will in this report hereby be referred to as Hand-Arm Vibration Syndrome (HAVS). It is one of the most commonly occurring occupational injuries [1]. HAVS is a neurodegenerative disorder and causes a lot of physical discomfort [2]. The finger, on a macro level, is broadly classified into three main layers of tissue. The epidermis contains the stratum corneum (a layer of dead epidermal cells) which is of interest to the extent of this project. The next two layers are the dermis and hypodermis, which will be referred to as soft tissue for the rest of the report. The neuron-rich area, which is of interest when studying a neurodegenerative issue such as HAVS, is understood to be located at the interface between the soft tissue and the epidermis and this area will, therefore, be of greatest focus in this project [3]. In Figure 1.1.1 the schematic overview of the finger can be seen and how it has been translated into the simplified 2D cross section FE model provided by Swerea IVF. (a) Human finger tissue layers [4]. (b) FE model of the finger. Figure 1.1.1: Model of the skin and how it is translated into the FE model in LS Dyna. The International Standards Organization (ISO) regulates the daily exposure limits to the vibrations from these tools. The regulation, ISO 5349, is based on a subjective assessment that averages out the frequency, states that frequencies of the range 8-1250 Hz are dangerous and frequencies higher than 1250 Hz are not. The suggestion that high-frequency vibrations are not causing damage is contradictory to the fact that neural symptoms have been found among dentists (as one example), using high-speed drills (with high vibrational frequencies) as part of their daily work. This validates the imminent need for a more effective and accurate method for estimating risk when handling high-frequency tools [2], [5]. To minimize the risk of this occupational injury, the next step in the research is to try to investigate how these high-frequency vibrations cause damage. The most pertinent question is the source of damage; is it the high acceleration, the daily exposure time or the frequency content [2]? Swerea IVF has launched a project aiming at proving the fact that high-frequency vibrations actually cause damage. And consequently, that ISO 5349 can not be used to set limits for exposure to high-frequency vibrations and needs to be updated to 1 create guidelines for tool manufactures and users, and through this help more than 350 000 workers, only in Sweden, who work with vibrating tools for more than 2 hours each day [1]. 1.2 Project description As part of the work towards changing the ISO 5349 to also regulate work with high-frequency tools that Swerea IVF is working at, this eight week project aimed at investigating the strain and strain rate in a FE model of a 2D cross-section of a finger created by Swerea IVF. The outcome will then be to guide future work to the areas that seems to be the most critical and that will need further investigation. 1.2.1 Project goal The goal of the project was to assess how the strain and the strain rate in the finger are affected by high-frequency transient vibrations and different material and geometry parameters by performing finite element simulations. At the end of the project, a table containing all calculated values were presented and should be regarded as the ultimate deliverable. Sub goals The project more specifically aims at providing answers to the following questions: 1. 2. 3. 4. 5. 6. 7. 8. Where in the finger does the maximum strain and strain rate occur? How does a change of acceleration amplitude affect the strain and strain rate? How does a change of acceleration frequency affect the strain and strain rate? How does the thickness of the different skin layers affect the strain and strain rate? How does the shear modulus of each skin layers affect the strain and strain rate? How does a change in bulk modulus affect the strain and strain rate? How does a change of the decay constant, β, affect the strain and strain rate? How does a change in load angle, θ, affect the strain and strain rate? 1.2.2 Provided material The following material has been provided to this project and has thus been made, investigated or produced earlier by others: • A 2D finger cross-section FE model, containing: – – – – The geometry Mesh Contact conditions Material model • Mechanical material data (for the model) • Suggestion of study on strain and strain rate • Load cases 1.2.3 Assumptions and delimitations Since this particular project had a very limited time frame and is part of a bigger project at Swerea IVF, assumptions and delimitations need to be specified. • Assumptions – The provided material is regarded as good, e.g. the finger model, the mesh, the boundary conditions et cetera is not to be evaluated. – The skin thickness is coupled, e.g. when the epidermis is thin the stratum corneum is also thin. • Delimitations – The high computationally demanding simulations limits the simulations to fractions of milliseconds. – Only 2D studies will be performed. – Only one lower and one upper bound value for the material parameters will be investigated. 2 2 Theory For the extent of this project, it is important to understand the theory or background of ISO 5349 and the human physiological factors. The ISO 5349 was studied to deduce the issues with the standard and how the analysis should be proceeded to avoid the similar discrepancies. The human physiological factors of HAVS will help determine what would be the possible issue and area of concern where failure initiates. 2.1 ISO 5349 As stated earlier, there exists a standard that regulates the daily exposure to vibrational forces. Low-frequency vibrations, that mainly affect joints and ligaments, are fairly well assessed by the standard. However, when the frequency increases there the set exposure time increases, meaning the labourer can use the tool for a longer time. It is believed this high-frequency vibrations which are left unchecked by the ISO 5349 causes HAVS. The ISO 5349 standard is adequate for frequencies of the range 8-1250 Hz. For frequencies above 1250 Hz, the damage is localized to the tissues in the distal phalanges. This proves that ISO 5349 has to be revised for higher frequencies [5]. The standard uses a method of frequency weighting. The acceleration signal is filtered and an average (RMS) is calculated. This method works fine for continuous signals but not as good for transient vibrations (which is common in a lot of tools) [3]. This method is viable for lower frequencies. However, for the higher frequencies, this is not an appropriate method. The mathematical model is valid on a macro level stating that the damage to the bones and ligaments are not present. While considering a micro level, specifically at the tissue level, this method of analysis fails. In reality, neuro-degenerative disorders develop leading to necrosis of the neurons [6]. Furthermore, the standard also extensively discusses the necessity to consider multi-axial vibration with acceleration in x, y and z-direction due to the different relative orientation of tool and fingers resulting in different types of strains. Vibration transducer (accelerometer) is used to measure large peak acceleration. 2.2 HAVS The hand arm vibration syndrome is mainly due to chronic work with vibrational tools. Although it is not clear how HAVS occurs the differential diagnosis is easily identified, consisting of neurological and vascular disorders [6]. The most tangible way in discerning risk of damage is a pain. This is a more immediate method of determining damage subjectively. An objective method is by actual analysis of the tissue at the structural or in some cases the cellular level. In the case of HAVS, the frequency is very high and the amplitude is relatively low. Nociception is best explained by the gate control theory [7]. This states that there is a range of which pain is interpreted by the brain. Anything below and above that is not transmitted by the neurons [8]. 2.2.1 Vascular disorders The superficial blood vessels experience intense vasoconstriction. There is a clear demarcation between the affected and normal parts. The affected sections of the finger are characterized with a white or dusky yellow appearance. Typically, the circumference of the finger is numb and will feel dead on touch during an episode [6]. 2.2.2 Neurological disorders Here the sensorineural receptors are affected. The disorder is identified by a tingling sensation or numbness of the finger tips. This is due to a reduced tactile sensitivity. Upon chronic exposure to the vibrations, one can develop a significant loss in manual dexterity [6]. 3 2.3 Finger model The computational FE model of the finger was done by Swerea IVF. The finger is modelled as a 2D cross section including three skin layers; soft tissue as the main inner part, epidermis as the middle layer and the stratum corneum as the outer layer, see Figure 1.1.1. A previous project conducted an experiment to determine the shape of the fingerprint. This was done by pressing a finger into a mould, obtaining the shape of the fingerprint. This mould was then scanned to give a computer model of the fingerprint geometry. Any gaps in the finger geometry were filled computationally. From the fingerprint, the geometry for the rest of the finger was extrapolated for the three tissue layers [9]. The FE model also includes a bar that represents the tool which effects the finger [9]. This bar is modelled to move vertically upwards in order to simulate the behaviour of one vibration pulse. The condition between the tool and finger is a simple contact condition. There are no fixed points or similar, only contact condition that makes the surfaces to not merge. The model is using plane strain condition. The model uses a mixed mesh i.e different element types and element sizes. The mesh is created with a tool in LS Dyna with some manual refinements from Swerea IVF. The material model is of viscoelastic type and the material parameters that need to be specified for each skin layer are: bulk modulus, shear modulus and a decay coefficient β. 4 3 Methodology This chapter describes the method used in order to answer the questions specified in the project goal. The method, in general, consisted of a literature study and investigation of the existing FE model to gain the knowlperiphery needed in order to define the project, pre-processing, simulations and post-processing. LS Dyna was used as a solver, which was proposed by Swerea IVF since the existing model was created there. The pre-processing included the setup of the FE model; specifying loads, material- and geometry properties. The post-processing included extracting the strain and strain rate and analyzing the data. The pre- and post-processing were executed in LS-PrePost. 3.1 Literature study Since there is still a lot of uncertainty about what in the finger that actually causes damage when subjected to high-frequency transient vibrations it was of interest to see how other fields, working with human tissue, are assessing the damage. The literature study therefore aimed at providing an understanding of these fields, to give a foundation for discussion to the possible future implementation of these assessments into the ISO 5349 standard. Furthermore, the literature study also aimed at providing knowlperiphery of similar work with vibrating tools and some mechanical data for human tissue that had not already been provided to the project. 3.2 Pre-processing The pre-processing focused on defining the inputs to LS Dyna in the three different areas; applied loads, material parameters and geometry. This is essential for the simulations to be run. 3.2.1 Applied loads The finger in the model was exposed to four different acceleration load pulses (corresponding to a transient), by prescribing an acceleration to the tool. These were decided in a discussion with the supervisor since Swerea IVF has experimental results for the same load pulses. In LS Dyna the load cases were defined as a discrete number of x-y-pairs describing the sinusoidal acceleration followed by a straight line with the constant value of zero. The four pulses as they were inserted into LS Dyna can be seen in Figure 3.2.1. The curves have different amplitudes (|amax |) and frequencies (T = time period) according to the list below, and are all affecting the model perpendicular (uniaxial) to its contact area with the tool: Load case 1: |amax | = 10 km/s2 , T = 0.01 ms Load case 3: |amax | = 10 km/s2 , T = 0.1 ms 3.2.2 Load case 2: |amax | = 100 km/s2 , T = 0.01 ms Load case 4: |amax | = 100 km/s2 , T = 0.1 ms Material parameters The mechanical properties of the finger layer vary with various aspects such as age, consumption of tobacco and alcohol, chemical exposure and environmental factors. As a person ages, there is thinning of the epidermal layer. This thinning is because of a reduced moisture content of the tissue and leads to an increased sensitivity to vibrations [10]. Therefore the material parameters in a finger are highly dependent on the individual and there are no parameters that fit every finger. Therefore to have an extensive study, parameters such as the thickness, shear modulus, the bulk modulus and decay constant have been varied. The parameters for the mechanical properties of interest can be seen in Table 3.2.1 and Table 3.3.1. These values were provided to the project by Swerea IVF. The bulk modulus and the decay constant β of the three layers are 2190 MPa and 215 s−1 respectively. The shear modulus was varying over a range for all three layers. For soft tissue, it varies from 0.0017 to 0.017 MPa. The minimum calculated shear modulus for the epidermis was 0.021 MPa and maximum calculated shear modulus was 0.2 MPa. The values varies between 0.6 and 10 MPa for the stratum corneum. 5 Load case 1 10 Load case 2 100 Prescribed motion on tool Acceleration [km/s2] Acceleration [km/s2] Prescribed motion on tool 5 0 -5 -10 50 0 -50 -100 0 0.01 0.02 0.03 0.04 0.05 0.06 0 Time [ms] Load case 3 10 0.01 0.02 0.03 100 0.05 0.06 Prescribed motion on tool Acceleration [km/s2] Prescribed motion on tool Acceleration [km/s2] 0.04 Time [ms] Load case 4 5 0 -5 -10 50 0 -50 -100 0 0.05 0.1 0.15 0.2 0.25 0.3 0 Time [ms] 0.05 0.1 0.15 0.2 0.25 0.3 Time [ms] Figure 3.2.1: The four different load pulses the finger model was exposed to. The material parameters were changed using the MAT option in the keyword manager of LS Dyna. To input, the given shear modulus, the short (G0 ) and long-time shear modulus (G∞ ) had to be considered. The short time shear modulus is given in Table 3.2.1 and the long-time shear modulus was calculated using a factor of multiplication obtained from earlier simulations done by Swerea IVF. Table 3.2.1: Engineering data for human tissue [9], [11]. Skin layer Soft tissue Epidermis Stratum Corneum 3.2.3 Shear modulus [MPa] 0.0017 − 0.017 0.021 − 0.2 0.6 − 10 Bulk Modulus [MPa] 2190 2190 2190 β [s−1 ] 215 215 215 Geometry The geometry or specifically the thickness of the different skin layers of the finger also varies between individuals. In order to analyze how the strain and strain rate were affected by this, the analysis was performed on three different finger skin thickness. The thickness of stratum corneum and epidermis for the three different finger models are given in Table 3.2.1. The area of the soft tissue was kept constant and the thickness of the other two layers was varied. The area of the soft tissue is fixed at 307.29 mm2 with an effective diameter of 19.78 mm. Table 3.2.2: Finger Geometry for different models [9], [11]. Geometry Type Thin Skin Intermediate Skin Thick Skin Tissue Layer Epidermis Stratum Corneum Epidermis Stratum Corneum Epidermis Stratum Corneum 6 Thickness [µm] 80 25 145 32.5 210 40 Area [mm2 ] 5.06455 1.39954 5.81096 1.48592 6.55406 1.57438 3.3 Simulations For the simulations, which were performed on the three different geometries, the different load cases were combined with different material parameters in order to assess how the different combinations affect strain and strain rate. In the finger model the following parameters had to be investigated: ∗ Soft tissue – Bulk modulus • Skin thickness – Thin skin thickness – Intermediate skin thickness – Thick skin thickness • Material dampening – Decay constant, β • Material stiffness • Vibration Parameters – Shear modulus – Amplitude of acceleration, |amax | – Time period of pulse, T – Uniaxial and biaxial ∗ Stratum corneum ∗ Epidermis The simulation was planned in two different stages, the first stage which thoroughly evaluates how the different shear modulus for different skin layers affect strain and strain rate and the second stage which briefly evaluates how a change of the remaining parameters affect the strain and strain rate. This division was based on the fact that some material data and geometry data was given from the start of the project whereas some were obtained at later stages. The first simulation stage is hereafter called the main simulation and the second simulation stage is called the extended simulation. All in all a total of 60 unique simulations were run. 3.3.1 Main simulation scheme The main simulation (simulation number 1 − 32 in Table A.0.1 in Appendix A) consists only of the intermediate skin thickness subjected only to uniaxial loads and aimed to provide ground for evaluating the importance of the different skin layers, shear modulus, the acceleration amplitude and frequency. With four load cases and one upper and one lower bound for the three different skin layers the total amount of unique simulations for the main simulation were 4 ∗ 23 = 32. 3.3.2 Extended simulation scheme The extended simulation (simulation number 33 − 60 in Table A.0.1 in Appendix A) used four of the simulations from the main simulation (number 8, 16, 24 and 32), all using the same material parameters but different load cases, as comparison basis. These four simulation cases were then modified in different ways, according to Table 3.3.1 and the three different geometries, to enable evaluation of the remaining pertinent model parameters. Thus, four simulations were run for the thick skin and four simulations for the thin skin using otherwise the exact same settings as for the intermediate skin. Another twelve simulations were run with all settings the same, but this time with the load applied in a biaxial way. For the β and the bulk modulus the load case 4 was chosen, since it was assumed to give the highest strain and strain rates, and then the parameters were changed. Four simulations for each parameter were run. Biaxial vibration As ISO 5349 discussed extensively effects of multi-axial vibrations, simulations with different acceleration directions were investigated to evaluate the combined effect of compression and shear wave propagation in the finger. Similar investigations were conducted in seismology (study of earthquake) to access P-waves (Compression) and S-waves (Shear) [12]. 7 Since this project used a 2D-model only biaxial load cases were investigated. With the intention to keep the resultant acceleration at the same amplitude as for the main simulation, the x and y components of the biaxial acceleration can be found by q |a(t)| = a2x (t) + a2y (t) (3.3.1) ax (t) = a(t) cos θ ay (t) = a(t) sin θ where θ is the angle defining the resultant direction of the biaxial acceleration, measured clockwise from the parallel direction of the finger to tool contact surface. Thus meaning that θ = 00 corresponds to a case parallel to the finger and with θ = 900 is the same loading direction as the main simulation. However, θ is different from incident angle analyzed in [12]. With four cases for the thick skin, four for the thin skin, four for the beta constant, four for the bulk modulus and finally twelve for the biaxial vibration, the total amount of unique simulations for the extended simulation were 4 ∗ 4 + 12 = 28. The values that were tested can be seen in Table 3.3.1. Table 3.3.1: Varying parameters and its values in extended simulations. θ 150 450 750 900 3.4 β [s−1 ] 2.15 21.5 2150 21500 Bulk Modulus [GPa] 1095 1971 2409 3285 Post-processing When all simulations were performed the data (strain and strain rate) were collected and analyzed in order to be able to answer the research questions. This section will provide information about how this was done. 3.4.1 Data collection After the simulations have been run the overall behaviour of the finger was assessed by observations. This was done in order to get an overview of how the strain and strain rate distributes over the finger and to see where the highest values appear. This knowlperiphery in combination with the suggestion from supervisor could then be used in order to select a number of nodes from where the maximum strain and strain rate could be collected. These nodes were selected for the final study: (A) 119498 (B) 457634 (C) 445538 (D) 415203 Figure 3.4.1: Location of selected nodes. Here it was assumed that the finger could be seen as symmetric and thus only half of the finger needed to be analyzed in order to assess how the entire cross-section was affected. Three of these nodes were selected in the soft tissue close to the epidermis while one of them is further in the finger. This is to see how the strain and strain rate varies along the fingerprint and also how it decays through the finger. Another reason for the selected area of interest is, as mentioned earlier, that many nerves and blood vessels are located in that 8 area. The strain and strain rate was selected in LS Dyna as ”Mean IPT effective strain” and ”Effective strain (v-m)-strainrate”. These values were selected since the effective strain and strain rate gives a good picture of where the maximum values are occurring. It was also mainly of interest to assess how the strain and strain rate changes and not the specific values, thus the specific extraction type is not important. 3.4.2 Data analysis When the strain and strain rate had been collected the data were analyzed. The analysis was divided into three parts; an overall assessment of the distribution of the strain and strain rate, the data analysis based on the main simulations and the analysis based on data from the extended simulations. The overall strain and strain rate distribution When all the simulations had been run the results were assessed by checking the graphical representations in LS Dyna as well as the graphs of strain and strain rate for each node. In some cases, the simulation time had to be extended in order for the actual maximum strain or strain rate to occur. When a satisfying result was obtained for all simulations the overall distribution of the strain and strain rate was assessed by observations of the graphical representations. It was observed where the maximum strain and strain rate occurred and also how the model behaved, i.e. did the finger leave the tool or not, did the finger deform and so on. Main simulation scheme The purpose of this analysis was to assess how a change of the shear modulus for the different skin layers as well as different acceleration amplitudes and frequencies affect the strain and the strain rate in the finger. The strain and strain rate for each simulation was plotted to get a basic understanding of the data to perform an in-depth analysis of each parameter. The obtained results were ordered based on the loading cases. Shear modulus analysis The 32 simulations were ordered based on similar changes in the shear modulus of individual skin layers. The next step was to understand the change in the strain based on the change of the shear modulus. The results were analyzed for all three nodes for a pair of simulations. Acceleration amplitude and frequency The analysis of how the acceleration amplitude and frequency affect the strain and strain rate was done in two ways. The first to understand how each node behaves with each load case. Second, the simulations were sorted in an ascending order of strain in each node. The same procedure was repeated for all load cases. To show how the load changes for a given material stiffness the Voight, Ev , and Reuss, Er , modulus was obtained using the area to give the volume fraction, vfi . This gives the upper and lower limits of the stiffness. Using the upper and lower limit, the Voight-Reuss-Hill modulus, Evrh , (mean of Voight and Reuss) was calculated [13]. The simulations were then ordered based on the Voight-Reuss-Hill modulus. From this, it can be understood how the strain changes for the varying load cases. The following equations were used: Ev = i X vfi · Ei (3.4.1) 0 Er−1 = i X vfi 0 Evrh = Ei Ev + Er 2 (3.4.2) (3.4.3) Extended simulations The purpose of this analysis was to assess how a change of skin thickness between the skin layers as well as a change in bulk modulus, β and θ affect the strain and the strain rate in the finger. 9 Geometry The three thicknesses were tested with the four load cases (simulations 33-40). The data was analyzed by first divide it into groups having the same load case and then into another set of groups using the same node. Firstly one could see whether the different nodes had equal behaviour and secondly see if the behaviour were the same in each node for the different load cases. The data analysis focused on comparing the relative difference between the different cases, rather than comparing their absolute values. This meant that parameter importance for strain and strain rate could be identified. Bulk modulus and β The intermediate thickness was subjected to four different β and bulk modulus. As there was only one load case, in this scenario, the data was sorted directly with the same nodes. The behaviour at each node with varying β and bulk modulus was observed and plotted. It is important to mention that the shear modulus of each layer is dependent on time and β. With those values the effective shear modulus G(t) is given as [14]: G(t) = G∞ + (G0 − G∞ )e−βt (3.4.4) In the above equation G0 is the short time shear modulus and G∞ is the long-time shear modulus. Further the strain and strain rates at these nodes were weight averaged as explained below. Biaxial vibration Nodes A, B and C were analyzed for maximum strain and strain rate for different θ values. The load cases were studied separately were strain and strain rate were plotted against θ to see how the change of resultant vibrational direction affects the values of strain and strain rate. Weighted mean strain and strain rate Furthermore, a study of the weighted mean strain was performed for the studies of bulk modulus, β and biaxial vibration. The objective was to determine how the fingerprint (area of interest) experiences strain and strain rate on a more global level. Since the other studies, regarding the nodes, only were looking at very local effects a global study could show how the entire area of interest is affected on a mean basis. This was done by weighing the area each node represents by comparing their relative length to the entire length (rather area compared to area, however since the height is the same for all cases and irrelevant the relative difference will be the same). The relative area was found, in the 2D cross-sectional FE model, by identifying the coordinates of the points at the center of the fingerprint in the horizontal direction as well as the corresponding point at the fingerprint periphery. The total horizontal distance between these two points was then calculated. The fingerprint was then divided from start (center of the finger) to the end (periphery of the finger) by splitting the length between the new (center and periphery) and existing (A, B, C) nodes. These three sub-lengths were then divided by the entire length and an area representation was found. Figure 3.4.2 illustrates how the area representation has been performed. In Table 3.4.1 the results from this area representation weighting can be found. Figure 3.4.2: Node area representation Table 3.4.1: Weights corresponding node area representation. Node Node A Node B Node C Area weight 0.19 0.26 0.55 W eighted M ean Strain = P W eighted M ean Strain Rate = P 10 (Straini ) · (Area W eighti ) P (Area W eighti ) (Strain Ratei ) · (Area W eighti ) P (Area W eighti ) (3.4.5) (3.4.6) 4 Results and discussion In this section, the results of the literature study, the simulations and data analysis will be presented and discussed. The overall result for the total of 60 simulations can be found in Appendix A in Table A.0.1. The strain and strain rates for each node A-D can also be found in the table. 4.1 Literature study results A few journals related to the effect of vibration on soft tissue were studied. From [15], the finger and palm were subjected to a similar loading condition. The angular frequency of the finger was also considered as an output parameter since it could be determined and analyzed physically. Another study was conducted [16] where an axial load was applied on brain tissue at different rates. The resulting displacement at the ends of this tissue was calculated through an attached micrometer. A mathematical model for energy function was also validated from the results obtained. However, only a few parameters could be analyzed through the experimental analysis. Further, a journal [17] conducted a computational analysis of the finger tip subjected to periodic loading. A FE model of the finger, similar to the one used in this project, was subjected to different frequencies was analyzed. The temporal threshold shifts (TTS) was used as a laboratory approach to investigate the exposure-response relationship [18]–[20]. From a biomechanical perspective, the effect of the vibration exposure on TTS is associated with the dynamic distribution of strain in the fingertip induced by vibration. This journal served as a motivation to consider effective strain as a parameter for this analysis [17]. Finally, in the seismological research [12], 3D shear waves (S-waves) of different incident angles were assessed to a greater extent using FEA. The impact of this research also motivated to evaluate the effect of shear and compression in the finger. 4.2 Overall strain and strain rate distribution Figure 4.2.1 show how the strain propagated within the finger at different times. This figure represents simulation case 32 (one of the worst cases regarding the deformation of the finger) and was therefore not representative for all simulations. At time t = 0, the finger was resting on the plate (no load). At the time instant when the acceleration was prescribed to the tool it started to move in the positive y-direction. Due to the viscoelastic behaviour of the finger, the region near the edge of the finger came in contact with the tool. This behaviour can be observed in Figure 4.2.1c. At t = 0.0545 ms, there was maximum contact area between the finger and the tool. By this time the tool reached its maximum velocity and from the next instant it retarded. Due to the inertia of motion, the finger continued to move in the positive y-direction. Hence, at t = 0.0696 ms the finger was no longer in contact with the tool. A pattern observed in all simulations was that the strain initiated uniformly along the whole fingerprint and later on propagated towards the edges of the fingerprint. The critical regions highlighted in red were generally found in the soft tissue close to the epidermis in the area of interest. The strain did not seem to decay over time, but instead, the location of the maximum strain seemed to change location (concentration of red elements changed location). When considering strain rate, a pattern similar to the strain can be observed, see Figure 4.2.2 (the same simulation case as for the strain). However, the strain rate decreased much faster and the absolute maximum occurred earlier and closer to the fingerprint than for the strain. The maximum seemed to occur at the boundary between the epidermis and the soft tissue, i.e. in the area where a lot of nerves are. By observing the graphical representations after the simulations were done it could be seen that some of the load cases (especially load case 1, simulations 1-8) resulted in a nearly undeformed finger model. The model did not lift from the tool. For some load cases, the finger got highly deformed and lost contact with the tool (especially load case 4, simulations 24 - 32). This can be seen in Figure 4.2.3. This behaviour of the model was questionable, especially the total compression when it seemed like elements merged together. The strain and strain rates from these simulations may not be accurate. This could have caused numerical errors and thereby affected the results. This needed to be considered when evaluating the results. 11 (a) Time = 0.0212 ms. (b) Time = 0.0424 ms. (c) Time = 0.0788 ms. (d) Time = 0.133 ms. Figure 4.2.1: Strain in the finger at four different times, for simulation 32. The scale is from blue = 0 to red ≥ 0.70. (a) Time = 0.0212 ms. (b) Time = 0.0424 ms. (c) Time = 0.0788 ms. (d) Time = 0.133 ms. Figure 4.2.2: Strain rate in the finger at four different times, for simulation 32. The scale is from blue = 0 to red ≥ 40000. (a) Simulation 5. (b) Simulation 29. Figure 4.2.3: Comparison of the deformation for different simulations at time = 0.0364 ms It also seemed as numerical errors were present at some locations along the fingerprint since high values of strain and strain rate were occurring at random nodes, as can be seen in Figure 4.2.4. These were all located along the boundaries between two of the skin layers. The fingerprint geometry included sharp corners which might be an explanation. Another reason for these high local strains might be the rapid change of material stiffness between the two layers. In a real finger, the change would probably be more gradual. 12 Figure 4.2.4: Singular point. The distribution of the strain and strain rates also seems to match up with the assumption of that the most critical area seems to be along the fingerprint. Since HAVS is a result of degenerated nerves these high strains and strain rates in this area which was assumed to have a high concentration of nerves might be a reason for the issue. It is however of importance to remember that this is a simplified model of the human finger and if other parts of the finger would be included the overall strain state might change. One issue that should be discussed is that the data collection was done on a subjective basis, where the graphs for the strain and strain rate for each node was observed to see if the maximum value had occurred. Since this was assessed by observing when a graph reached a maximum point and then decreased it might be possible that local maximums had been found (if the graphs would start to increase again). However, it is assumed that this was not the case. For some of the extended simulations, the model deformed a lot and left the tool, hence the validity of the model for longer simulation times were unclear. One factor that would possibly be needed to take into consideration for longer times would be the fact that a human might use muscles to press the finger against the tool to prevent it from leaving the tool. Therefore the method used was assumed to be valid for the scope of this project. 4.2.1 Node D The results for node D were similar for all simulations. The strain increased linearly but never reached a maximum within the simulation time. Therefore the maximum values found in the result table, see Appendix A, and were thus not global maximum values. Some simulation times were extended in order to see if a global maximum could be obtained. These simulations, however, took 2 − 4 hours to run and it was therefore not possible to redo all simulations in order to obtain the maximum strain and strain rate of node D. It was also, as mentioned earlier, not clear about how valid the model was for longer simulation times. Therefore a decision was made to not include values for node D in the analysis. The result can, however, tell that node D reacted slower to the load, which also was expected since it was located further within the finger. The same delay also occurred for the strain rate but was more erratic. It was observed that all the values for node D (strain and strain rate) were much lower than for the other selected nodes. This showed that the strain and strain rate decrease further into the finger. 4.3 Main simulation The result of the data analysis of the main simulation will be presented in this section. The analysis mainly focuses on the strain and strain rate for varying load cases and shear modulus. 4.3.1 Load The increase in acceleration frequency resulted in a greater increase in strain rather than an increase in acceleration amplitude, as can be observed in Figure 4.3.1. For load case 1 and 2 there was a change in acceleration amplitude but the frequency was kept constant. For load case 2 and 4, however, the frequency was changed but the amplitude was kept the same. In the same figure, it can be concluded that node A seemed to be subjected to the highest strain. This might be because node A was located in the peripheral region of the fingerprint. Node A experienced the 13 highest strain probably since there existed shear as well as normal strains, unlike at node C, in the center of the finger, which probably only experienced normal strains. Figure 4.3.1: Strain variation across simulations. A more detailed view of the nodal behaviour due to the different load cases can be seen in Figure 4.3.2. It was seen that node B followed either node A or C depending on the frequency. For load cases 1 and 2 the frequency of vibration was high. For these cases, node B seemed to experience strains similar to node A. This indicated that there were a lot more factors that had to be further studied in order to draw a conclusion. For load cases 3 and 4, node B followed the same trend as node C. Figure 4.3.2: Increase in Strain for each load case. The strain rate followed a similar trend but is in the figure not ordered increasingly as in the case of the strain, see Figure 4.3.3. This was because the absolute maximum value of the strain rate was not reached for the same simulation case as for the strain. 14 Figure 4.3.3: Strain rate variation across simulations. According to ISO 5349, low-frequency vibrations cause greater damage to the human body compared to high-frequency vibrations. The strain rate analysis showed that high-frequency vibrations caused assumed high strain rates. Comparing the strain rates from load case 2 and 3 in Figure 4.3.4, the common range was 1000 − 3500s−1 . This indicates that high-frequency vibrations may affect the human tissue and possibly result in HAVS. Figure 4.3.4: Strain rate variation for a given load case. 4.3.2 Shear modulus The next analysis, shown in Figure 4.3.5, was to study how the load affected the strain but in an isolated method i.e. when there was a constant shear modulus. To study this, the effective shear modulus was determined using the Voight-Reuss-Hill average. It was found that the strain more or less increased exponentially from load case 1 to load case 4. 15 Figure 4.3.5: Variation of strain (y-axis) for the different load cases for a constant shear modulus (x-axis). In Figure 4.3.2 it can be seen which simulation resulted in the highest and lowest strain for each load case. These simulations were for the high value of the shear modulus in stratum corneum but low values for the other two skin layer’s shear modulus. By intuition, it might be thought that the lowest material stiffness would result in the highest strain which was not seen here. An explanation could be that stratum corneum (which was the outermost layer) when very stiff behaved like a rigid body and only translated the full load into the other skin parts resulting in higher strains. The result of including volume fractions of each skin layer can be seen in Tables 4.3.1 and 4.3.2. A comparative analysis of the data validates the initial finding of how a change in the soft tissue affects the strain to a greater extent than the stratum corneum. However, when considering the volume fraction of each tissue layer, the strain is more dependent on the change in shear modulus of the stratum corneum rather than the soft tissue. As an example take simulation pairs based on loading cases and the shear modulus of the soft tissue and stratum corneum being constant. When the soft tissue layer changes from 0.0017 to 0.017, the strain is increased to approximately 1000 % for every simulation. However, taking the same set of pairs with the epidermal shear modulus shifting from 0.017 to 0.0017, the strain drops down with 10 %. Table 4.3.1: Increasing shear modulus trend. Tissue Soft tissue Stratum corneum Epidermis Shear modulus increase [%] 1000 1667 952 Average change in strain throughout the finger [%] 61 181 89 Table 4.3.2: Decreasing shear modulus trend. Tissue Soft tissue Stratum corneum Epidermis Shear modulus decrease [%] 10 6 11 Average change in strain throughout the finger [%] 187 92 128 16 4.4 Extended simulation The result of the data analysis of the extended simulation will be presented in this section. 4.4.1 Geometry When looking at how the thicknesses of the skin layers affect the strain and strain rate at the selected nodes, one can see in Figure 4.4.1 and 4.4.2 that it is hard to find a pattern of the behaviour for different load cases. For example, when looking at the strain in node A, in Figure 4.4.1a, one get different behaviour for different load cases. It is however only in node C one can find some kind of pattern where one can see that a thinner skin gives higher strain. Another interesting observation was that the intermediate thickness seemed to give the lowest strain. For node A it also seemed to be a higher strain for thin skin, but lowest for thick skin. For node B it is the other way around. The different behaviour in node A, B and C can depend on many things. For instance, the different locations of the nodes will, of course, play a big part in the measurements. However, the differences in the trends are hard to explain. Between all the graphs it is hard to find a pattern relating the skin thickness to the measured strain and strain rate. Whether this is a result of bad data due to previously mentioned error sources or actually is a valid result describing how a real finger would respond is hard to tell. With too little information and knowledge to properly assess this issue, it could be assumed that the results may be at least somewhat faulty. One other aspect that was not specifically tested in this project is how the material parameters affect the result based on the skin thickness because then the parameters have a bigger area to act on. For example, a change of shear modulus in a thick skin probably affects more than the same change in a thin skin. But with thicker skin one also have more skin to work with which maybe can the compensate. Comparison of strain behavior at node b Comparison of strain behavior at node a Comparison of strain behavior at node c Strain Strain Load case 1 Load case 2 Load case 3 Load case 4 Strain Load case 1 Load case 2 Load case 3 Load case 4 Load case 1 Load case 2 Load case 3 Load case 4 Thick Intermediate Thin Thick Intermediate (a) Thin Thick Intermediate (b) Thin (c) Figure 4.4.1: The curves are moved together to easy make an comparison of the behaviour of the strain in nodes A, B and C with varying thickness for the four different load cases. Comparison of strain rate behavior at node a Comparison of strain rate behavior at node b Thick Comparison of strain rate behavior at node c Intermediate (a) Thin Thick Strain rate Load case 1 Load case 2 Load case 3 Load case 4 Strain rate Load case 1 Load case 2 Load case 3 Load case 4 Strain rate Load case 1 Load case 2 Load case 3 Load case 4 Intermediate (b) Thin Thick Intermediate Thin (c) Figure 4.4.2: The curves are moved together to easy make an comparison of the behaviour of the strain rate in nodes A, B and C with varying thickness for the four different load cases. 17 4.4.2 Bulk modulus It could be observed from figure 4.4.3 that the node A experiences a much higher strain than the other two nodes. As mentioned above in Section 4.3.1 the reason is that the two waves propagated into the finger originate in a region where the node A lies. This induces both compressive and shear stresses acting on this node which in turn causes higher strains. The bottom portion of the finger undergoes a concave expansion. The extent of this expansion is lesser at the center and more near the edges. As the node B is nearer to the edge it experiences more strain than node C. From intuition it could be understood that as the bulk modulus increases the finger becomes more stiff, which resulted in reduced strains. The same behaviour was observed in the analysis. For all the three nodes and the weighted average, the strain graph was decreasing with increasing bulk modulus. A similar pattern could be observed with the strain rates. Figure 4.4.3 depicts that as the bulk modulus increases the strain and strain rates decrease. Since the bulk modulus is a measure of how easy it was to compress a volume these results make sense. Variation of strain with bulk modulus Node A Node B Node C Weighted Variation of strain rate with bulk modulus 3.5 Strain rate [s -1 x104] 0.5 Strain 4 0.4 0.3 3 Node A Node B Node C Weighted 2.5 2 1.5 0.2 1000 1500 2000 2500 3000 1 1000 3500 1500 Bulk modulus [GPa] 2000 2500 3000 3500 Bulk modulus [GPa] (a) Strain. (b) Strain Rate. Figure 4.4.3: Variation of strain and strain rates with respect to the bulk modulus. 4.4.3 Decay constant, β In the case of the strain when the behaviour was analyzed with respect to decay constant, it follows a definite pattern. From Figure 4.4.4a it could be understood that, as the β values increases, the strain graph falls and reaches an optimum value. Once this optimum value is reached the graph stabilizes. Unlike the strain graph, the strain rate at node A behaves in an opposite manner. The strain rate increases and then stabilizes. The same pattern was observed for node B, node C and weighted average graph. Variation of strain with β Strain 0.5 Strain rate [s -1 x104] Node A Node B Node C Weighted 0.6 Variation of strain rate with β 3.5 0.4 0.3 Node A Node B Node C Weighted 3 2.5 2 1.5 0.2 1 0 2 4 6 Log scale of β 8 10 0 (a) Strain. 2 4 6 Log scale of β 8 (b) Strain Rate. Figure 4.4.4: Variation of strain and strain rates with respect to the decay constant. 18 10 4.4.4 Biaxial vibration From Figure 4.4.5, when θ = 45o and θ = 75o , higher weighted-mean strain and strain rates were obtained for most of the load cases. However, when θ = 15o and θ = 90o , strain and strain rates were relatively low. Moreover, the amplitude of the wave seems to be crucial for higher strain rates in load case 4. Variation of strain with θ , load case 1 x10 -3 2.5 Strain 600 Node A Node B Node C Weighted 2 1.5 1 0.5 40 60 θ [Degrees] 300 200 80 20 Variation of strain with θ , load case 2 0.02 0.015 3000 2000 0.005 0 20 40 60 θ [Degrees] 80 20 Variation of strain with θ , load case 3 3000 Node A Node B Node C Weighted 0.03 0.02 40 60 θ [Degrees] 80 Variation of strain rate with θ , load case 3 Node A Node B Node C Weighted 2500 Strain rate [s -1 ] 0.04 Strain Node A Node B Node C Weighted 1000 0.01 2000 1500 1000 0.01 500 20 40 60 θ [Degrees] 80 20 Variation of strain with θ , load case 4 3.5 Node A Node B Node C Weighted 0.4 0.3 40 60 θ [Degrees] 80 Variation of strain rate with θ , load case 4 Node A Node B Node C Weighted 3 Strain rate [s -1 x104] 0.5 Strain 80 4000 Strain rate [s -1 ] 0.025 40 60 θ [Degrees] Variation of strain rate with θ , load case 2 Node A Node B Node C Weighted 0.03 Strain 400 0 20 0.05 Node A Node B Node C Weighted 100 0 0.035 Variation of strain rate with θ , load case 1 500 Strain rate [s -1 ] 3 2.5 2 1.5 1 0.2 0.5 20 40 60 θ [Degrees] 80 20 40 60 θ [Degrees] 80 Figure 4.4.5: Variation of strain and strain rates with θ in all load cases. As shown in Figure 4.4.6, when decreasing θ from 90o to 45o , dominant compression waves caused more compressive strain (ǫy ) and shear strain (γ) at the periphery of the finger base. However, decreasing θ below 19 45o , dominant the shear waves increased the shear strain (γ) and the tensile strain (ǫx ) at the center of the finger base. At θ = 45o , damage could be highest, since there was significant amount shear strain (γ) at the periphery and center of the finger due to the equal amplitude of both shear and compression waves propagated into the soft tissue. As mentioned in table 3.2.2, soft tissue has least shear modulus compared to other layers making it susceptible to interlayer shear damage to the high concentration of nerves and blood vessels at soft tissue-epidermis boundary. Furthermore, the strain rate followed a similar trend as the strain. (a) θ = 15o Strain (b) θ = 15o Strain rate (c) θ = 45o Strain (d) θ = 45o Strain rate (e) θ = 75o Strain (f) θ = 75o Strain rate (g) θ = 90o Strain (h) θ = 90o Strain rate Figure 4.4.6: Contours of strain and strain rates with varying θ at load case 4. 4.5 General discussion Many assumptions have to be made in order to proceed with this project. The finger model is a simple one, where many parts of a finger are not considered. For the purpose of this project, which is in a pretty early stage in considering a full-scale finger model, it is probably enough to consider a 2D model of the finger. Since the main focus of the project is to evaluate which of the load cases and parameters, that could be changed in the model, is affecting strain and strain rate the most. Thus keeping a constant geometry, boundary conditions and mesh is important. It is, however, important to understand the simplicity of this model. One major part that is not considered in the model is the bone. However, since this project is focusing on the strain and strain rate at the interface of the soft tissue and the epidermis, computationally there is no need for the finger model to have a bone. However, the resonance frequency of the finger model will be altered in presence of the bone. If parameter such as pressure wave propagation within the finger is to be considered, it will be pertinent to include the bone. Another complication of working with the human material is the great variation among individuals. For the model the fingerprint is based on experiments on different individuals and the average is used which is one way used in order to represent the difference between individuals. The range of the shear modulus for the different skin layers also has to be seen as a way of representing different individuals. Stratum corneum, the outer layer, has the greatest difference which also makes sense since that is the skin layer that can be affected by how the hands and fingers are used. Also, other factors can have an effect on the geometry and properties of the skin. As an example on ageing, the epidermis becomes marginally thin. This is because the moisture content of the cells reduces. The dermis (the layer of soft tissue rich in nerves), reduces in thickness by 20 %. This leads to a decreased tactile sensitivity. As a result, it seems that there cannot be a clear demarcation between the symptoms of HAVS and the normal ageing process. The project has provided results that hopefully can be used as a guidance for where future work should focus on. The range of shear modulus for stratum corneum is still the largest compared to the soft tissue and the epidermis. This range needs to be further scrutinized with different combinations of the thickness of the skin layers. The problem that has been seen in this project is however that there are many factors affecting so it might be that for some loads (since vibrational tools might have components of both different 20 frequencies and amplitudes) some skin layer compositions might be better than others. This project, comprising of 60 simulations with different material, geometrical and vibration parameters, has also proved that propagation of high-frequency vibrations with high amplitude results in high strain and strain rates in the finger model. However, the extent of physical damage due to these high strain and strain rates will be assessed in the rat-tail experiment conducted by the supervisor in the near future. 21 5 Conclusion and recommendations In this section, the conclusions will be stated and the recommendations for future work will be brought up. 5.1 Conclusions The maximum strain and strain rate seemed to occur along the fingerprint in the area where it was assumed to be a high concentration of nerves. The strain and strain rate were affected greater by the increase in acceleration frequency than acceleration amplitude. The shear modulus of stratum corneum affected the result the most if consideration was taken to volume fraction, otherwise, the soft tissue was affecting the most. The thickness of the different skin layers was affecting differently dependent on load cases and location in the finger (no trend could be seen). The bulk modulus and the decay constant had not a great impact on the strain and strain rate. An increase in both bulk modulus and decay constant was, however, resulting in a decrease of strain and strain rate. The orientation of the vibrational force affected both magnitude and distribution of strain and strain rates along the finger base. 5.2 Recommendations for future work There is still a long way to go before an actual proof of that high-frequency vibrating tools are damaging the human fingers. It is still unclear what is actually causing the damage and if it would be the strain and the strain rate, as the assumption was for this project, there would still be a need of having actual values for when those parameters are causing damage. This is a recommendation for future work. In addition, the finger model could be improved in different ways to obtain more accurate results. The sharp corners between the skin layers should be smoother, boundary conditions (contact between tool and finger, friction etc), as some examples, should be looked up and investigated. The finger model is as mentioned earlier still far from a real human finger. How this simplification is affecting the result should also be investigated. 22 Geometry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Thick Thick Thick Thick Thin Thin Thin Thin Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Intermediate Time Amplitude Time period Beta Bulk Modulus Epidermis Epidermis Strat Corn Strat Corn Soft tissue Soft tissue [ms] [MPa] G0 [MPa] GI [MPa] G0 [MPa] GI [MPa] G0 [MPa] GI [MPa] [ms] [km/s^2] [1/s] 1 0.5 0.5 0.5 0.06 0.06 0.5 0.06 0.1 0.1 1 1 0.06 0.06 0.5 0.06 0.15 0.15 0.15 0.15 0.2 0.15 0.2 0.15 0.15 0.15 0.65 0.15 1 0.15 0.2 0.15 0.06 0.06 0.15 0.15 0.06 0.06 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.06 0.06 0.15 0.15 0.06 0.06 0.15 0.15 0.06 0.06 0.15 0.15 10 10 10 10 10 10 10 10 100 100 100 100 100 100 100 100 10 10 10 10 10 10 10 10 100 100 100 100 100 100 100 100 10 100 10 100 10 100 10 100 100 100 100 100 100 100 100 100 10 (θ = 45) 100 (θ = 45) 10 (θ = 45) 100 (θ = 45) 10 (θ = 15) 100 (θ = 15) 10 (θ = 15) 100 (θ = 15) 10 (θ = 75) 100 (θ = 75) 10 (θ = 75) 100 (θ = 75) 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.01 0.01 0.1 0.1 0.01 0.01 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.01 0.01 0.1 0.1 0.01 0.01 0.1 0.1 0.01 0.01 0.1 0.1 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 215 2.15 21.5 2150 21500 215 215 215 215 215 215 215 215 215 215 215 215 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 1095 1971 2409 3285 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 2190 0.021 0.021 0.021 0.021 0.2 0.2 0.2 0.2 0.021 0.021 0.021 0.021 0.2 0.2 0.2 0.2 0.021 0.021 0.021 0.021 0.2 0.2 0.2 0.2 0.021 0.021 0.021 0.021 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.0200004 0.0200004 0.0200004 0.0200004 0.19048 0.19048 0.19048 0.19048 0.0200004 0.0200004 0.0200004 0.0200004 0.19048 0.19048 0.19048 0.19048 0.0200004 0.0200004 0.0200004 0.0200004 0.19048 0.19048 0.19048 0.19048 0.0200004 0.0200004 0.0200004 0.0200004 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.19048 0.6 0.6 10 10 0.6 0.6 10 10 0.6 0.6 10 10 0.6 0.6 10 10 0.6 0.6 10 10 0.6 0.6 10 10 0.6 0.6 10 10 0.6 0.6 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 0.58068 0.58068 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 9.678 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.0017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.017 0.0005732 0.005732 0.0005732 0.005732 0.00057324 0.0057324 0.00057324 0.0057324 0.00057324 0.005732 0.00057324 0.0057324 0.00057324 0.0057324 0.0005732 0.0057324 0.00057324 0.0057324 0.00057324 0.0057324 0.0005732 0.0057324 0.0005732 0.0057324 0.00057324 0.0057324 0.00057324 0.0057324 0.005732 0.0057324 0.005732 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 0.0057324 Strain A Strain rate A 0.0020842 0.0012392 0.0054875 0.0038125 0.0011814 0.0010868 0.0044728 0.0025409 0.027813 0.017007 0.064007 0.046302 0.01636 0.015425 0.055059 0.032741 0.12648 0.089838 0.086139 0.063565 0.10458 0.060212 0.084958 0.040275 0.47278 0.40414 1.107 0.41787 1.079 0.44915 0.85975 0.51485 0.0029902 0.031665 0.034178 0.5075 0.003074 0.032097 0.043674 0.54696 0.54526 0.51712 0.51294 0.50666 0.64158 0.51976 0.51353 0.50227 0.0013594 0.02064 0.025447 0.40921 9.811 548.0000 12.697 17.587 2.348 30.658 37.281 52.587 41.465 32.333 128.75 115.64 65.04 63.939 143.91 189.13 572.1 464.65 1648.6 1443 1014.8 1003.6 1927.2 2419.2 3613 3434.6 2687.9 3083.1 3817.4 3613.1 2429.7 2537.2 18907 26247 22690 25474 28784 45393 25806 32959 206.97 2135.6 1901.6 33936 217.54 2280.9 3084 28007 37819 33595 32431 30728 28571 33499 32836 31627 98.39 1476.1 1319.8 23517 10.677 370.27 555.03 8236.4 19.172 2246.4 2334.6 31777 Strain B Strain rate B Strain C Strain rate C 0.0013053 0.00055697 0.0045296 0.0010621 0.00087941 0.00083891 0.0037552 0.00096261 0.010384 0.0083416 0.059974 0.01448 0.012546 0.011245 0.050041 0.015388 0.12139 0.063298 0.075269 0.040292 0.081866 0.055452 0.052589 0.037962 0.78665 0.4936 0.82306 0.3949 0.63962 0.42127 0.47578 0.36311 0.0015124 0.015144 0.038368 0.38197 0.0012103 0.014403 0.032921 0.33598 0.3901 0.36724 0.3582 0.35005 0.45287 0.34981 0.36405 0.37437 0.0011025 0.021518 0.045043 0.4731 13.375 19.561 27.642 31.978 97.972 17.926 43.952 23.504 37.093 33.714 103.54 83.122 271.45 286.77 188.37 438.01 690.13 672.8 971.94 810.18 2676 2252.4 3285.6 3092.2 4153 3622.6 2906.6 2047.6 4905.2 6386.8 2480.5 2407.7 32957 30441 35859 23165 40124 43123 25606 21459 618.84 2652.8 3125.2 29384 833.87 3569.7 1691.1 20133 29835 22412 20958 19067 24937 22542 21439 20712 523.54 3653.9 2329.4 28905 23.262 3793.8 1692.4 15826 369.19 2885.1 2369.3 27026 0.0013053 0.00055697 0.0035549 0.0017315 0.00082412 0.00077923 0.0022792 0.00075657 0.0078198 0.0070025 0.043627 0.021413 0.0120603 0.011631 0.026595 0.0076552 0.058431 0.035442 0.037847 0.020394 0.038368 0.030514 0.020474 0.016443 0.37965 0.26854 0.45459 0.21475 0.49312 0.213 0.20007 0.16385 0.0007161 0.0082139 0.018006 0.16892 0.0010353 0.010987 0.019202 0.21789 0.18284 0.16673 0.16105 0.15475 0.19597 0.15821 0.16432 0.16732 0.0014515 0.019065 0.03184 0.32773 0.00049305 0.019941 0.027991 0.30131 0.001087 0.01018 0.023525 0.43354 39.043 42.107 67.477 101.64 155.95 128.6 96.163 213.39 910.42 773.85 648.31 547.22 1772.9 1743 1342.4 2147.6 2607.5 2426.2 1285 785.84 4161.6 4226.8 1457.7 1198.9 19531 17759 445538 9726 10733 26726 13400 10868 267.61 1580.1 1051.6 9762.4 314.96 2808.9 981.99 13427 13377 11420 10926 10172 13308 10650 10766 11553 365.82 3922.1 1707.8 16464 77.198 4134.6 1730.2 15 667.0 272.81 2662.9 1438.2 13 807.0 Strain D Strain rate D 0.00052167 6.2011 25.566 0.00009056 0.00010942 8.2727 9.8889 0.00034115 0.0017136 0.0019477 0.004165 0.00435 0.0011566 0.0012144 0.056175 0.0043048 0.01384 0.015035 0.011421 0.01398 0.022126 0.020704 0.014533 0.020244 0.093726 0.098914 29.408 42.124 58.28 89.67 102.16 82.946 82.235 340.35 100.24 188.23 236.07 112.59 434.89 191.15 725.59 218.4 561.67 1334 1859.9 0.13586 4621.3 5769.8 6140.9 1931.8 5708.6 18.787 105.99 619.87 6188.2 25.529 102.82 504.39 5814.6 6170.6 5773.2 5703.2 5428.4 2991.4 5863.3 5492.2 3903.6 34.716 79.065 1076.8 9655.7 4.7558 44.363 838.47 8677.2 172.7 96.433 815.54 7931.4 0.21263 0.12479 0.19577 0.00040189 0.0046881 0.02093 0.20646 0.00036386 0.0039139 0.018467 0.21789 0.20238 0.19448 0.19924 0.19466 0.094304 0.20138 0.18981 0.1526 0.0002049 0.0030069 0.038159 0.3685 0.000041309 0.001055 0.028335 0.33185 0.00034867 0.0041573 0.028836 0.27953 23 Table A.0.1: All simulation cases with corresponding data. The main simulation consists of an analysis where the load cases and the shear modulus for the three different skin layers have been varied based on provided mechanical material data for human tissue. The extended simulation analyze the change of vibrational direction, β, bulk modulus and skin thickness. For the extended simulation, the only parameter change that has real world support is the skin thickness. All other parameter changes are done by percentage changes, that may or may not have a physical representation. In the Table A.0.1 all data for all simulation cases can be found. The green part reassembles the simulations for the main simulation, whereas the blue part reassembles the simulations for the extended simulation that are based on the darker green simulation cases from the main simulation. Simulation scheme with data A # References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] H. Lindell, Project proposals in applied mechanics, Wave propagation in human tissue from transient vibrations, 2017. S. G. Raju, O. Rogness, M. Persson, J. Bain, and D. Riley, “Vibration from a riveting hammer causes severe nerve damage in the rat tail model”, Muscle & Nerve, vol. 44, no. 5, pp. 795–804, 2011, issn: 1097-4598. doi: 10.1002/mus.22206. [Online]. Available: http://dx.doi.org/10.1002/mus.22206. Discussions with supervisor hans lindell at swerea ivf, Mar. 28 2017 & Apr. 28 2017. Modified by cropping the full picture and using only the important part., Jun. 1, 2017. [Online]. Available: https://en.wikipedia.org/wiki/Human_skin#/media/File:Skin_layers.svg. Svenska standard ss-iso 5349-1. Swedish Standards Institute, 2001. R. McCaig, Abc of occupational and environmental medicine, D. Patel and D. Snashall, Eds. John Wiley & Sons, Ltd., Publication, 2014, isbn: 978-1-4443-3817-1. R. Melzack and P. D.Wall, “Pain mechanisms”, Science, vol. 150, no. 3699, pp. 971–979, Nov. 1965. [Online]. Available: http://www.paincenter.pitt.edu/sites/default/files/Melzack%20and% 20Wall%20Gate%20Control%20Theory%20_.pdf. M. Hansen and J. Cheng, Understanding pain: What you need to know to take control: What you need to know to take control, R. D. U. Alan D Kaye M.D., Ed. ABC-CLIO,LLC, 2011, isbn: 978-0-313-39603-8. H. Berg, M. Demant, J. Edwijn, and D. Gustafsson, An investigation of the shear modulus for stratum corneum and epidermis, 2016. S. V. Saxon, M. J. Etten, and E. A. Perkins, Physical change & aging. Springer Publishing Company, 2015, isbn: 978-0-8261-9864-8. A. Pena, M. Arronte, E. De Posada, P. Luis, and T. Flores, “Non-invasive optical method for epidermal thickness estimation”, Online Journal of Biological Sciences, vol. 14, no. 2, pp. 163–166, 2014, issn: 1608-4217. doi: 10.3844/ojbssp.2014.163.166. [Online]. Available: http://www.thescipub.com/ ojbs.toc. J.-q. Huang, X.-l. Dua, M. Zhao, and X. Zhaoa, “Impact of incident angles of earthquake shear (s) waves on 3-d non-linear seismic responses of long lined tunnels”, Engineering Geology, vol. 222, pp. 168–185, 2017. [Online]. Available: https://doi.org/10.1016/j.enggeo.2017.03.017. G. Mavko. (). Effective medium theories, [Online]. Available: https : / / pangea . stanford . edu / courses/gp262/Notes/9.EffectiveMediumTheories.pdf. Ls dyna keyword users manual-volume 1. Livermore Software Technology Corporation(LSTC), 2007. K. Hamoudaa, S. Rakhejaa, and P. Marcotte, “Fingers vibration transmission performance of vibration reducing gloves”, International Journal of Industrial Ergonomics, vol. 15, no. 5, pp. 795–804, 2016, issn: 0169-8141. doi: 10.1002/mus.22206. [Online]. Available: http://dx.doi.org/10.1016/j. ergon.2016.11.012. K. Miller and K. Chinzei, “Mechanical properties of brain tissue in tension”, Journal of Biomechanics, vol. 35, no. 4, pp. 483–490, 2002. doi: 10.1002. [Online]. Available: https://doi.org/10.1016/S00219290(01)00234-2. J. Wu, K. Krajnak, D. Welcome, and R. Dong, “Analysis of the dynamic strains in a fingertip exposed to vibrations: Correlation to the mechanical stimuli on mechanoreceptors”, Journal of Biomechanics, vol. 39, no. 13, pp. 2445–2456, 2006. [Online]. Available: https://doi.org/10.1016/j.jbiomech. 2005.07.027. N.Harada and M.J.Griffin, “Factors influencing vibration sense thresholds used to assess occupational exposures to hand transmitted vibration.”, British Journal of Industrial Medicine, vol. 48, no. 3, pp. 185–192, 1991. [Online]. Available: http://www.jstor.org/stable/27727152. S.Maeda and M. J. Griffin, “Temporary threshold shifts in fingertip vibratory sensation from handtransmitted vibration and repetitive shock”, British Journal of Industrial Medicine, vol. 50, no. 4, pp. 360–367, 1993. [Online]. Available: http://www.jstor.org/stable/27727616. S. Kihlberg, M. Attebrant, G. Gemne, and A. Kjellberg, “Acute effects of vibration from a chipping hammer and a grinder on the hand-arm system”, Occupational and Environmental Medicine, vol. 50, no. 11, pp. 731–737, 1995. [Online]. Available: http://www.jstor.org/stable/27727616. 24