Linköping University Department of Management and Engineering Master thesis 30 credits 2019|LIU-IEI-TEK-A–19/03390−SE
CAE modelling of cast aluminium in
automotive structures
Master Thesis in Division of Solid Mechanics
Sreehari Veditherakal Shreedhara&Subrat Raman Singh
Academic supervisor: Daniel Leidermark Industrial supervisor: Johan Jergeus
Examiner: Kjell Simonsson
CAE modelling of cast aluminium in automotive structures SHREEHARI VEDITHERAKAL SHREEDHARA
SUBRAT RAMAN SINGH
© SHREEHARI VEDITHERAKAL SHREEDHARA, SUBRAT RAMAN SINGH, 2019
This work is protected according to the copyright law (URL 1960:729)
Master thesis 2019
ISRN : LIU-IEI-TEK-A–19/03390-SE
Department of Management and Engineering Divison of Solid Mechanics
Linköping University SE-581 83 Linköping Sweden
Abstract
In automobile industry, there is a big push for the automotive car manufacturers to base engineering decisions on the results of Computer Aided Engineering (CAE) solutions, and to transform the prototyping and testing, from a costly iterative process, to a final verification and validation step. The variability in components material properties and environmen-tal conditions together with the lack of knowledge about the underlying physics of complex systems often make it impractical to make reliable predictions based on only deterministic CAE models. One such area is the CAE modelling of cast aluminium components. These cast aluminium components have gained a huge relevance in the automobile industries due to their commendable mechanical properties. The advantage of the cast aluminium alloys are being a well-established alloy system in manufacturing processes, their functional in-tegrity and relatively low weight. However, the presence of pores and micro-voids obtained during the manufacturing process constitutes a specific material behaviour and establishes a challenge in modelling of the cast materials. Furthermore, the low ductility of the mate-rial demands for the advanced numerical model to predict the failure.
The main focus of this master’s thesis work is to investigate modelling technique of a cast aluminium alloy component, a springtower, for a drop tower test and validate the predicted behaviour with the physical test results. Volvo Car Corporation currently uses a material model provided by MATFEM for cast aluminium parts which is explored in this thesis work, to validate the material model for a component level testing.
The methodology used to achieve this objective was to develop a boundary condition to perform component level tests in the drop tower and to correlate these with the obtained results found by using various modelling techniques in the explicit solver LS-DYNA. There-fore, precise and realistic modelling of the drop tower is crucial because the simulation results can be influenced by major design changes. A detailed finite element model for the spring tower has been developed from the observations made during the physical testing. The refined model showed good agreement with the existing model for the spring tower and observations from physical tests.
Keywords: Cast aluminium, Explicit FEA, LS-DYNA, Drop tower, Spring tower, FEM, CrachFEM, MFGenYld
Acknowledgements
We would like to sincerely thank Peter Nyström (Senior Manager, Principal Engineering, Volvo Car Corporation) for the thesis opportunity and for facilitating us with all the infras-tructure for the thesis work. We are very grateful to our supervisor at Volvo Cars Safety Centre, Johan Jergeus for his tremendous guidance and support during the thesis. His un-wavering supervision of our work and sharing his knowledge and experience made us fill up the expectations set during this thesis.
Our sincere thanks to Domenico Macri (Senior Simulation Engineer), Jan Hinder (Prin-cipal Engineer), Sandeep Shetty (Senior CAE Engineer) and Pulkit Sharma (CAE Engineer) for their invaluable technical assistance. We gratefully acknowledge the valuable feedback during this thesis work from Mathias Retzlaff, Mats Landervik, Per-Anders Eggertsen and Simon Rydberg in the safety center at Volvo Car Corporation. We would also like to thank Reino Frykberg, Andreas Bjurhult and Erik Axen in helping us conduct the physical test-ing. Finally, we would like to thank everyone at Volvo Car Corporation who have directly or indirectly supported our work and for making us feel pleasant all the time during our thesis work. Their motivation coupled with various ideas and suggestions has been crucial for the progress of the thesis.
We express our deep gratitude for our academic supervisor Daniel Leidermark and ex-aminer Kjell Simonsson at Division of Solid Mechanics, Linköping University for their valu-able feedbacks with academic writing and constant guidance throughout our work. Their insights and constructive criticism along with their valuable suggestions have made us improve a lot while performing this thesis work. Our heart-filled thanks to our parents, friends and well-wishers who were motivating and cheering us all along the journey.
Linköping, September 2019 Sreehari and Subrat
Nomenclature
Abbreviations and Acronyms
Abbreviation Meaning
BIW Body in white
VCC Volvo Car Corporation
CAE Computer Aided Engineering
GenYld General Yield
FE Finite Element
GISSMO General Incremental Stress State dependent damage MOdel
OEM Original Equipment Manufacturer
IGES Initial Graphical Exchange Specification
ELFORM Element Formulation
NIP Number of Integration Points
HGID Hourglass ID
Mathematical Symbols
Symbol Meaning
[ ] Rectangular Matrix or Square Matrix
{ } Column Vector . Time Derivative Latin Symbols Symbol Meaning A Area c Speed of Sound [C] Damping Matrix
[D] Degree of freedom of structure matrix
E Young’s Modulus
I1 First Invariant
L Length
[M] Mass Matrix
{R} Force vector on structure’s nodes
Si j Deviatoric stress tensor
t Time
v Velocity
m Mass
Greek Symbols
Symbol Meaning
β Stress State parameter
ε Normal Strain
η Stress triaxility
θ Shear stress parameter
εp Equivalent plastic strain
λ Plastic Multiplier
ν Poisson’s ratio
σ Stress
ω Shear fracture
Subscripts and superscripts
Abbreviation Meaning
i, j Row and column of the Matrix
p Plastic part cr Critical x, y, z Cartesian co-ordinates m Mean eq Equivalent vM von Mises
Contents
1 Introduction 1 1.1 Motivation . . . 1 1.2 Objective . . . 2 1.3 Delimitations . . . 2 1.4 Research question . . . 2 1.5 Resources . . . 3 1.6 Other Considerations . . . 3 2 Background 4 2.1 Cast aluminium . . . 4 2.2 Previous Work . . . 5 3 Theoretical background 6 3.1 Yield Criteria . . . 63.1.1 von Mises yield locus . . . 6
3.2 Strain-rate dependency . . . 8 3.3 Hardening rules . . . 8 3.3.1 Isotropic Hardening . . . 9 3.3.2 Kinematic Hardening . . . 9 3.4 Damage Mechanics . . . 9 3.5 Failure modes . . . 10
3.5.1 Ductile normal fracture . . . 10
3.5.2 Ductile shear fracture . . . 11
3.5.3 Instability . . . 11
4 Method 12 4.1 Physical test set-up development . . . 12
4.1.1 Finite Element Model . . . 12
4.1.2 Analysis of the developed model . . . 14
4.2 Physical test . . . 15
4.2.1 Drop tower test . . . 15
4.2.2 Drop test evaluation . . . 16
4.3 Improved boundary condition . . . 17
5 Results and Discussion 19 5.1 Results from Load case 1 . . . 19
5.1.2 With improved boundary condition . . . 21
5.2 Results from Load case 2 . . . 29
5.2.1 With improved boundary condition . . . 29
6 Conclusion and future work 32 Bibliography 33 A Appendix 35 A.1 Explicit analysis . . . 35
A.2 Termination time for various models . . . 36
A.3 Results from original setup . . . 37
A.4 Results from simulation of Load case 1 at 7.5 m/s . . . 38
A.5 Results from simulation of Load case 1 at 5 m/s . . . 40
1
|
Introduction
1.1
Motivation
With increasing demand of electric vehicles, and for those vehicles to achieve a higher driv-ing range there is a necessity for reduction of weight. For these reasons various new mate-rials are being used by the automotive industry. Matemate-rials such as carbon-fiber reinforced plastics are effective in reduction of the vehicles weight but are also expensive to produce and manufacture. On the other hand, materials such as aluminium has achieved lot of relevance, due to their low density and well-established manufacturing processes [1].
The combination of strength with ductility is the key performance parameter for auto-motive applications. However, the relationship between the material properties, strength, stiffness and weight of a component is very complex and can be strongly influenced by ge-ometry and the manufacturing process.
In Figure (1.1), various materials used in a current Volvo’s car body are shown. Green coloured parts are the aluminium components used in the body. Currently, cast aluminium components used in the body of Volvo cars are the spring towers. Apart from the body other components such as engine block, transmission casing and brackets are also manufactured
CHAPTER 1. INTRODUCTION
by castings of aluminium. As Volvo has announced that they will aim for 50 percent of their sales to be fully electric by 2025 [2], with an aim for long endurance which can be achieved substantially by reducing the weight without compromising the crashworthiness of the vehicle. This makes cast aluminium a good choice in the current period. Therefore, it is essential to study the behaviour of each cast aluminium component under certain load conditions.
1.2
Objective
The main purpose of this thesis work is to study, verify and improve the Computer aided engineering (CAE) modelling methods for cast aluminium structures in a car. This is done by investigating the current modelling and simulation techniques used at Volvo Cars in Gothenburg, Sweden. Commercially available finite element (FE) programs are used to simulate and analyze the load conditions to predict failure in components. For this study, one such component, a spring tower is analyzed. Volvo Cars Corporation (VCC) uses a ma-terial failure model, MFGenYld+CrachFEM developed by MATFEM, that is incorporated in a FE software. These material cards have parameters defined from numerous coupon tests, but they have not been verified exclusively for cast aluminium components at VCC. Thus, a verification of the material is to be done to ensure that the material failure card is applicable to cast aluminium components in FE simulation. Based on this, the following was the scope of the thesis:
• The initial step was to decide a test setup, rigid enough to hold the component during a physical drop tower test.
• The second part was to conduct a physical test on the component and to obtain the results from it.
• The final step was to correlate the physical results with CAE by evaluating various modelling methods.
1.3
Delimitations
The thesis work had to be completed in 20 weeks and to stay within the time frame, it was limited to:
• Analyzing and verification of only one cast aluminium component, a spring tower. • Verification only on a component level.
1.4
Research question
Based on the problem description and the purpose of this thesis the following research questions are answered during the analysis and conclusion:
• What is the best failure prediction method for cast aluminium based on the the cur-rent computational capabilities?
• Will shell and solid elements provide accurate results for relatively thin walled struc-tures ?
CHAPTER 1. INTRODUCTION
1.5
Resources
Following resources and computer tools were used to carry out the work: • LS-DYNA as finite element solver
• ANSA for pre-processing and modelling
• METApost for post-processing and reading and calculating results • MATLAB for plotting results
• Drop tower test rig at VCC
1.6
Other Considerations
No ethical or gender-related issues are aroused by the work. Nor does it have a direct con-nection to issues concerning the environment or sustainable development of society. How-ever, the work is on cast aluminium structures which reduces the vehicle’s weight hence producing more fuel economical vehicles.
2
|
Background
This chapter discusses more about cast aluminium alloys, their properties and general prob-lems while manufacturing them. Some previous work done their at VCC, regarding the spring tower testing is also mentioned here.
2.1
Cast aluminium
Traditional heavy steel components have been replaced by cast aluminium because of their various advantages. Even though they have certain advantages over conventional iron, they have some disadvantages as well, such as:
• Different regions of the same component show different properties • Variable imperfections in the part
• CAE methods and experience with cast aluminium is limited
The mechanical properties of these components are very much dependent on the micro-structure of the alloy. These properties can be modified by controlling the rate of cooling and also depend on the pressure maintained during the casting process. Cooling rate also influences the shape, size and spacing of the micro constituents. The casting method used depends on the amount of silicon percentage used in the component. Also, the rate of cooling is also very much governed by the amount of silicon used. Strength and ductility of a
CHAPTER 2. BACKGROUND
component will increase by adding silicon in the alloy. During the casting process, when the liquid solidifies the volume decreases and hence, the shrinkage results in the formation of porosity in the alloy. This is one of the reasons why aluminium can show varying properties in the same part. These pores are common in cast aluminium alloys and they can also lead to micro-crack initiation and propagation when subjected to external loading [3].
As mentioned earlier, use of these components should not affect the safety of the vehicle. Thus, it becomes very crucial to understand the behaviour of a car body and its components during a crash. A strut tower, also called a spring tower, placed in the car’s body as seen in Figure (2.1), is studied for its behaviour. The suspension and its mounts are fixed to this tower. For choice any cast aluminium component could be studied but based on its availability at VCC, this component was chosen. It was assumed that the modelling method chosen should work on all other cast aluminium components with appropriate boundary conditions.
2.2
Previous Work
In year 2012, drop testing was performed on the same cast aluminium component, but because of not having an optimal setup, desirable results were not obtained. Figure (2.2) shows the previous test setup. From the figure it can be seen that the spring tower and the ground were bounded by an adhesive and it was tightened using a bind belt. When the test was conducted using this boundary condition it was observed that during the impact the spring tower showed horizontal movement, as the fixture was not rigid and could not withstand the load. The results obtained from the test showed that there was no failure in the model as expected and the component showed only ductility.
3
|
Theoretical background
This chapter discusses the material and failure models typically used within commercially available finite element programs to predict the behaviour of a casting alloy. Due to low ductility of castings, it is essential to gain dependable analyses from crash simulations. However, modelling cast aluminium in FE environment has been a challenge for many years. Cast aluminium component shows isotropic behaviour, i.e. the material has the same properties in different directions.
The material and failure model MFGenYld+CrachFEM has been used in this thesis work provided from VCC and calibrated accordingly to suit the cast aluminium behaviour. It is an elasto-plastic and comprehensive failure model developed to work with various finite element software available [4]. This card uses two models, one for plasticity and the other for failure. MFGenYld (generalized yield model) combines a range of yield locus descriptions. In this work, a corrected von Mises yield criteria with an isotropic hardening rule is used and the failure model comprises failure modes such as ductile normal, ductile shear, and instability.
3.1
Yield Criteria
The yield criteria defines the limit of elastic behaviour throughout the loading history for a general three dimensional stress state [5]. Hence, it is pivotal in predicting whether or not a material will begin to yield under a given stress state, resulting in permanent deformation (shape change) or failure of the material.
3.1.1 von Mises yield locus
The von Mises equivalent stress describes a common yield locus for isotropic materials with symmetry under tension and compression, which is valid for a three-dimensional or genral stress state [4]. The von mises equivalent stress is described as:
σvM e (σi j) = r 1 2[(σx−σy) 2+ (σ y−σz)2+ (σz−σx)2+ 6(σ2x y+σ2yz+σ2zx)] (3.1) Here, σx,σy,σz are the normal stress components and σx y,σyz,σzx are the shear stress components in three dimensions. The above equation can also be written as :
σvM e (σi j) =
r 3
CHAPTER 3. THEORETICAL BACKGROUND
Here, Si j is the deviatoric stress, which is further defined as Si j=σi j−1
3σkkδi j (3.3)
where,σi j is the stress tensor,σkkis the trace of the stress tensor andδi j is the Kroneckers delta.
Corrected yield locus Base yield locus
Figure 3.1: Yield locus
In the Figure (3.1) the base von Mises and corrected yield locus are shown. For perfect plasticity the yield function for a material following von Mises yielding will be:
f (σi j) =σvMe (σi j) −σY= 0 (3.4) where, f (σi j) = 0 is the yield surface of the material in question andσY is the initial yield limit of the material. In modelling plasticity of materials an associated flow rule is com-monly assumed which is also considered in this work for which the plastic strain rate is perpendicular to the flow potential at the point of the current stress stateσi j. The base yield locus is sometimes further modified, where it is scaled in bi-axial tension and com-pression. In this work, the yield locus correction factor is defined as m, which can further be described as
m · f (σi j) =σvMe (3.5)
where, f (σi j) is the reference yield locus of the material and m = m(σi j). The correction factor depends on the principal in plane stresses represented by stress triaxility (η). Stress triaxility can be defined as :
η=q σ1+σ2+σ3
σ2
1+σ22+σ23−σ1σ2−σ1σ3−σ2σ3
CHAPTER 3. THEORETICAL BACKGROUND
The values forηchange under tension and compression, positive values denote tension and negative values compression [4].
Under tension, the correction is:
m = 1 −bT· b (Y L) T − 1 bT· b(Y L)T
· g(η) (3.7)
where, g(η) is an interpolation function, bT is a strain dependent bi-axial correction factor and bY LT is the ratio of uniaxial and bi-aixial tensile stress.
bT=σBT σvM e bY LT = σ vM e σ(Y L) BT (3.8)
In equations (3.8) σvMe is the von Mises equivalent stress, σBT is real stress in bi-axial tension as measured andσ(Y L)BT stress under bi-axial tension according to base yield locus. Under compression, the correction is:
m = n µ · 1 −bT· b (Y L) T − 1 bT· b(Y L)T · g(η) ¶ (3.9)
where, n =σT/ |σC | is a strain dependent asymmetry factor between tension and com-pression. Bi-axial corrections factors for compression is defined in the same manner as for tension. bT=|σBC| σvM e bY LT = σ vM e |σ(Y L)BC | (3.10)
where,σvMe is the von Mises equivalent stress,σBCis real stress in bi-axial compression as measured andσ(Y L)BC stress under bi-axial compression according to base yield locus [4].
3.2
Strain-rate dependency
The material behaviour and it’s properties are significantly dependent on the deformation rate. A strain rate dependent material will have high stress with fast loading rate. It is important on a large scale when simulating load cases at different velocities, as done in this study. Strain rate dependency is usually done by specifying various stress vs strain curves for different strain rates [4]. There are few factors which affect strain rate dependency, such as the composition of the aluminium alloy and the heat treatment the alloy has gone through.
3.3
Hardening rules
During plasticity, the yield surface changes as a result of the hardening that develops dur-ing the history of plastic deformation. The two basic models for hardendur-ing are isotropic and kinematic hardening. As mentioned the plasticity model MFGenYld uses isotropic harden-ing behaviour, which is explained in this section. Mixed hardenharden-ing models have also been developed which include both yield surface expansion and translation [5].
CHAPTER 3. THEORETICAL BACKGROUND
Figure 3.2: Isotropic hardening (left) and Kinematic hardening (right)
3.3.1 Isotropic Hardening
This model is based on the assumption that the yield surface expands uniformly in stress space as yielding occurs, as seen in Figure (3.2). The general form of the isotropic hardening yield criterion can be written as:
f (σi j,εp) =σvMe − k(εp) = 0 (3.11) where, f (σi j) is the yield function which defines the shape and orientation of the yield surface. The parameter k(εp) defines the size of the yield surface as a function of the plastic deformation history, equivalent plastic strain, and is initially equal toσY.
3.3.2 Kinematic Hardening
In kinematic hardening, it is assumed that the yield surface translates in the six dimension stress space but does not changes size or shape. It can be written as:
σvM
e (σi j−αi j) = r 3
2(Si j−αi j)(Si j−αi j) (3.12) where, the difference f (σi j−αi j) is called the overstress or reduced stress,αi jis a backstress tensor and Si j is deviatoric stress tensor.
3.4
Damage Mechanics
When considering material failure, damage is considered as important for the failure. Dam-age is characterized by nucleation, growth and coalescence of voids in the material until the load bearing area has eventually been reduced and material separation occurs. This reduc-tion in area leads to material softening. In various damage models, material softening is added to the constitutive relation either by porous plasticity, cf. Gurson (1977) or by con-tinuum damage mechanics cf. Lemaitre (1985). A method to consider damage evolution is based on a relationship between the damaged and undamaged area as shown in the Figure (3.3) [6]. The damage D is the measure for the reduction in area considering the voids or defects in the material. This damage variable can then be interpreted as the ratio between the damaged area, Ac and the initial or undamaged area Ar.
CHAPTER 3. THEORETICAL BACKGROUND
Figure 3.3: Undamaged area Ar, damaged area Acand normal direction of the surface n
From the equation (3.13) we can define D as: D = Ac
Atot =
Atot− Ar
Atot (3.14)
Ultimate fracture is expected when D becomes 1, the entire surface is damaged and there is no material left to maintain the connectivity. The reduction of the cross section has an effect on the stresses used to map the global response. So, the effective true stress for a damaged material would be as:
σtrue= (1 − D)σe f f ective (3.15) This is the basic damage theory and fundamentals for various failure models. The only differences is the damage accumulation. As discussed by Lian Xui [7], the simplest form of a damage accumulation is a linear function which means that damage evolution is propor-tional to the ratio of equivalent plastic strain rate to the equivalent plastic strain at failure. This is used for most the failures models such as Johnson and Cook.
˙ D =ε˙p
εf
(3.16) Here, ˙εp is the equivalent plastic strain rate and εf is the equivalent plastic strain at failure. Apart from these, other damage evolution laws are also available such as in the models like GISSMO and CrachFEM. Due to confidentiality the damage evolution model used in CrachFEM can not be described here.
3.5
Failure modes
Various failure mechanism are involved in fracture. Brittle fracture is without any prior plastic deformation while ductile occurs under plastic deformation. Fracture of glass or ceramic is typically of brittle type. The failure mechanisms observed in aluminium alloys are ductile fracture and localized necking. Localized necking is an instability phenomenon which emerges from local inhomogeneties [8]. Ductile fracture mechanism accounts for two types of failure, normal fracture and shear fracture [4].
3.5.1 Ductile normal fracture
Ductile normal fracture occurs when micro-voids grow and finally fuse under tensile loads. This surface is rough and lies perpendicular to the principal loading direction [4]. For a
CHAPTER 3. THEORETICAL BACKGROUND
Figure 3.4: Different failures modes
Plastic Instability(left), Ductile normal failure(center) and Ductile shear failure(right)
three dimensional problem a parameter defining the stress state is needed solely defined by the stress triaxility. With a combination of triaxility (η), see equation 3.6 and normalized principal stress, fracture can be described for stress state that have same triaxility but different fracture strains. This parameter relates the first principal stresses to the von Mises stress:
σ1= σ1
σvM e
(3.17) Tensile loads cause normal fracture and the normal fracture risk is accumulated only if the stress triaxility is positive, i.e.η> 0.
3.5.2 Ductile shear fracture
Shear band localization causes ductile shear fracture and can occur under various stress states. Unlike the stress state parameterσ1/σvMe in the ductile normal fracture model, the shear fracture model usesτmax/σvMe as a second parameter in the stress space. It is assumed that the fracture strain is a function of the shear stress parameterθwhich is defined as:
θ= σ vM e
τmax· (1 − kSF·η
) (3.18)
Here, kSF is a material dependent parameter describing the influence of stress triaxilityη on the shear fracture. The value of kSF is positive and under tension or compression maxi-mum shear stress is perpendicular to the sheet and shear fracture occurs under 45°over the sheet thickness [4].
3.5.3 Instability
The final criteria for failure is instability, seen in Figure (3.4). Equivalent plastic strain at onset of failure is a function of the deviatoric stress ratio,ζand the strain rate.
ε∗
M=ε∗M(ζ, ˙ε
p) (3.19)
where,ζ=σ2/σ1. Instability is due to localized necking, which is a failure under tension and occurs when the material becomes unstable and deformation localizes. The instability model is usually used with shells, with solids a very fine mesh can only account for local instabilities.
4
|
Method
In this chapter, the experimental test setup, it’s development and the finite element model is explained. The chapter also discusses about the co-relation work carried out to obtain desirable results.
4.1
Physical test set-up development
Developing a setup for the component’s testing was necessary to correctly correlate the results with FE analyses. As the component had to be tested in a drop-tower, the initial step was to develop boundary condition within the FE environment with the failure model. The component, the spring tower, was then physically tested with the here formulated boundary conditions and later correlated with the results obtained in the FE analysis. Thus, it was essential to obtain a boundary condition considering following factors:
• The boundary condition should not be difficult or impossible to construct in reality. • One or more load cases can be tested with very small or no change in boundary
con-dition.
• The test setup should be reasonably rigid during the impact for various load cases.
4.1.1 Finite Element Model
The CAD model of the spring tower was provided by Volvo at the beginning of the work. This model was imported to ANSA, a pre-processing tool, and a shell mesh was generated from the mid-surface. The finite element model had the spring tower and a cylindrical impactor. A quadrilateral element type was used having five through thickness integration points. The whole spring tower with impactor was modelled with 44390 shell elements with one of the commonly used element formulations, the Belytschko-Lin-Tsay formulation. This formulation is based on Reissner-Mindlin kinematic assumption theory with 5 degrees of freedom in local coordinate system and 6 degrees of freedom in global coordinate system [9]. Iterations were done with impactor striking the component from various sides of the spring tower. Contact forces were defined to make the parts interact with each other and the impactor was given an initial velocity. The impactor’s initial velocity and mass were iterated after each analysis. The purpose was to estimate the velocity and mass of the impactor, at which the component was damaged. As seen in the previous test the component showed only ductility and thus it became important during this study for the component to fail. Based on the simulations, two load cases were chosen as seen in Figure (4.1).
CHAPTER 4. METHOD
Load case 1 : For the first load case, with the impactor’s mass of 90kg and velocity at 6m/s, there was failure in the component. The impactor was placed along the x-direction as seen in the figure.
Load case 2 : The second load case had the same impactor mass but the velocity was reduced to 5m/s. Also, the impactor’s direction was rotated horizontally 90 degrees, i.e. aligned along the y-direction as seen in the figure.
Figure 4.1: Load case 1 and Load case 2
Figure 4.2: FE model spring tower, plates and bolts
Once the boundary condition was developed it had to be modulated in a way so that it fulfills the requirements. To do this, the three holes, shown in the Figure (4.2) were locked using M14 bolts. Two plates were also modelled as rigid bodies to hold the the spring tower firmly from either sides. There was a possibility to observe bending in the bolts during the impact which could effect the results by inducing cracks in the component. In such a case, cylindrical casing was constructed around the beam elements that would be used just to represent a bolt. This cylindrical casing was created in a way that it was similar to the physical dimension of the bolt and the contact between the bolt and spring tower could be assigned for proper simulation. In ANSA, these three M14 blots were modelled with pretension assigned to them. The bolt joints model consists of the shank, bolt head
CHAPTER 4. METHOD
and null beams. The bolt joints were modelled using two noded beam elements. The reason for the bolts joints to be modelled in the analysis was to verify if the bolts would withstand the impact without bending, during the physical test. From this it is possible to mimic the real test setup and get more accurate results that can be compared with the data from the physical test. Based on the time for testing and above mentioned requirements this test setup was considered appropriate.
4.1.2 Analysis of the developed model
The simulation model for the FE solver LS-DYNA was defined as a input deck with all the necessary files in it. This included control cards, material data and other part files. The bolt, plates and impactor were included as a separate part file. Explicit analysis (c.f. Appendix A.1) simulations were performed in LS-DYNA using 192 CPU cores at VCC. The total elapsed simulation time as seen in the Figure (4.3) are with MFGenYld + CrachFEM failure model. The total simulation time for spring tower was 3 minutes and 52 seconds.
CHAPTER 4. METHOD
After normal termination, the results stored in d3plot and binout files were read in METApost. D3plot stores the animation data from where deformations can be loaded while binout saves numerical data for plotting curves. The results at various states (time steps) were loaded and the plastic strain in the component were fringed. The component damaged near the impact points and the fringed contour can be seen in Figure (A.1). The elements failing were deleted here from the analysis, few elements were failing at the point of impact and some were failing because of bending caused by the flange.
4.2
Physical test
To determine and predict failure of each components in a vehicle body, it is necessary to study them individually. For this thesis work a physical test was done to correlate the experimental results with the simulations. It was important to get a failure in the compo-nent to understand the velocity at which the compocompo-nent is damaged and if the same can be achieved by CAE modelling techniques.
4.2.1 Drop tower test
There are numerous testing methods used by industries now to replicate an impact sce-nario. One of them is a drop tower test. It is generally used to test energy absorption by a component under a loading condition. An impactor with a certain mass is guided by columns of the tower. These guide rails have roller bearing attached to it for reducing the friction.
CHAPTER 4. METHOD
Force transducers and accelerometers are fitted on the impacting mass to measure the force and accelerations upon impact. Two accelerometers are also placed on the the rig and table to capture accelerations from it. The rig holds the specimen to be tested by some kind of fixture. This fixture should be rigid to hold the component on impact. Once the desired velocity is decided, the impactor is moved to the corresponding height and released from that distance to hit the component. The height from which the impactor is to be released, is approximated by equation (4.1) but due to friction and other external parameters the height must be adjusted. h = v 2 2g (4.1) K .E =1 2mv 2 (4.2)
Using the kinetic energy equation, the maximum impact energy is calculated.
4.2.2 Drop test evaluation
During the drop tower testing a few changes of the boundary condition were necessary, as mentioned below:
• M14 bolts used during the test had their head removed and were welded to the rigid plate.
• The impactor mass was changed to 88 kg from 90 kg.
• Before every test, a new height was calculated from equation (4.1) and it was almost impossible to achieve the same velocity at every impact. Such as for a 6 m/s the velocity in the drop test was 6.07 m/s.
• Due to time constraint the force transducer was not accordingly calibrated for the test. Hence, the results obtained from these transducers were wrong and for corre-lation those results were not considered. Only data from accelerometer was taken in account.
While conducting the drop tower test, 21 components were tested for both the load cases. Few tests with certain velocities were repeated twice or thrice to verify the consistency of obtained results. The impactor in the physical test had two accelerometers equipped in it and the results were obtained from the accelerometers as channel data. These data were imported in METApost and was filtered (with CFC-180) to remove noise from the curve. The acceleration was multiplied by the mass to obtain the force. Again, the filtered acceleration curve was integrated to obtain the velocity and from the obtained velocity, the displacement was plotted. Finally, a force vs displacement curve was plotted for each accelerometer. An average of all the force-displacement curves from these accelerometers were taken to compare it later with FE results. The energy curve was also calculated from this data. The absorbed energy in the test is the area under the graph in the force vs displacement curve, which was plotted by integrating the force-displacement curve.
CHAPTER 4. METHOD
4.3
Improved boundary condition
The FE model was modified to mimic the drop tower test setup. The whole setup of impactor with rollers, guide rail and mass was included in the FE model. In Figure (4.5) the whole setup is shown, the spring tower is meshed with solid elements here. The rollers had no constraints in any direction so that these rollers allow some movement after impact. The center nodes of all the rollers were connected with the nodes of the impactor as a nodal rigid body, which in turn was connected to a rectangular box modelled as a solid to assign mass to the whole setup. The rollers were assigned properties of steel and the guide rail was modelled as a rigid body. The impactor and ground plate were assigned mild steel property. The ground plate had elements constrained with the bottom most plate to lock the movement. The mass of impactor was also changed to 88 kg as it was altered during drop tower test. Finally, the position of the impact point in the FE model was modified as per the physical test. In FE simulations, the welding of bolts to the rigid plate was not taken in consideration. The important factor was for a pretension to act on the bolts holding the two plates together.
Figure 4.5: Modified setup for co-relation with solid meshed model of the spring tower
Correlation was also done by using solid and shell meshing for the spring tower. Various elements lengths were evaluated for both meshes. A few analysis were done with higher
CHAPTER 4. METHOD
order elements to incorporate bending. Iterations were done with different element formu-lations and friction coefficients to obtain better results. These simuformu-lations results are dis-cussed in Chapter 5. Simulation times with various elements are mentioned in Appendix, see Section A.2.
5
|
Results and Discussion
The drop tower tests were done at various velocities starting at 6 m/s for the first load case. There was a small crack in the component at 6 m/s. Few tests were also done at higher velocities to compare crack propagation in the physical test with CAE. After completion of the first load case, tests were done for the second load case with velocity starting at 5 m/s. All the results obtained from the drop tower tests were compared with that of FE and are presented in this section.
5.1
Results from Load case 1
5.1.1 With original boundary condition
Finite element simulations were fringed with plastic strain to check deformation in the component. Element deletion is identified as a crack or fracture in the component’s region.
Figure 5.1: Physical test compared to CAE model with rigid impactor at 6m/s
As seen in Figure (5.1), due to bending of the flange the component fractured. The area marked in red shows the critical location of the crack. This section tend to fold inwards due to impact, causing it to crack. With a rigid impactor setup, i.e., the initially developed FE model, shown in Figure (4.1), the initiated crack location was different compared to physical test, see Figure (5.2). The comparison was done with model meshed with shell elements of average length 3.4 mm using ELFORM02. It was considered that the setup revealed devia-tion in fracture predicdevia-tion due to rigid properties being assigned to the impactor. However, this deviation could also be because of other factors such as element types, formulations
CHAPTER 5. RESULTS AND DISCUSSION
and their sizes. For this reason it was necessary to study these factors.
(a) Crack in the component (b) Failure in CAE Shell model
Figure 5.2: Incorrect failure prediction in FE model with rigid impactor compared to physical test at 6m/s 0 10 20 30 40 50 60 70 80 90 Displacement (mm) -1 0 1 2 3 4 5 6 7 8 F o rc e (N ) ×104 Physical Simulation
Figure 5.3: Force vs displacement curve of physical and original FE model at 6m/s
The force-displacement curve of the impact at 6 m/s with original setup is plotted in Figure (5.3), where it can be seen that the correlation between the simulated and physical test curves also do not show such a good result. The physical test curve seen here is the average from two accelerometers used in the test. The peak force in the physical test is 75.46 kN but from the FE results it was 61.68 kN.
CHAPTER 5. RESULTS AND DISCUSSION
5.1.2 With improved boundary condition
To obtain better results with a shell mesh, analysis were done with the modified setup shown in Figure (4.5). As mentioned in section 4.3, the impactor was given elastic material properties from being rigid. Simulations were completed by changing various element for-mulations. In this case Element formulation (ELFORM) 16 was chosen, a fully integrated shell formulation. The simulation time was longer compared to the ELFORM02, as it has more integration points in the plane for each element. Results obtained from ELFORM16 were better compared to that from ELFORM02.
0 10 20 30 40 50 60 70 80 90 Displacement (mm) -1 0 1 2 3 4 5 6 7 8 F or ce (N ) ×104 Physical ELFORM16 3.4mm fs & fd 0.2 ELFORM16 2.5mm fs & fd 0.2 ELFORM16 3.4mm fs 0.25 & fd 0.2 ELFORM02 3.4mm fs & fd 0.2
Figure 5.4: Force-displacement curve of physical test with simulations pf shell model at 6 m/s
(a) 3.4 mm ELFORM16 (b) 2.5mm ELFORM16
CHAPTER 5. RESULTS AND DISCUSSION
Figure (5.4) shows few simulation curves to compare the FE results with the physical tests. It shows good correlation between the physical and simulated results. During corre-lation friction became an important parameter to be identified. Earlier simucorre-lations were done with static and dynamic frictional coefficient ( fs and fd) as 0.2, but as the surface of the impactor was rough, it’s value was increased to get better results. As friction was not measured experimentally, the value for static friction was estimated and changed to 0.25. Upon increasing the value for static friction better results were achieved. In the curve it can be seen that with ELFORM16 and 3.4 mm average element size with fs of 0.25 the maximum force is 76.18 kN, which is close enough to the value from experimental results. Also, the prediction of failure was better approximated. The FE model analysis had a failed element at the location of the crack, but the element was not deleted. In the physical test the crack was also very small and the correlation was considered good. As solid model was
0 10 20 30 40 50 60 70 80 90 Displacement (mm) -2 0 2 4 6 8 10 12 F or ce (N) ×104 Physical ELFORM16 3.4mm ELFORM16 2.5mm ELFORM16 1.0mm ELFORM16 2.5mm scatter ELFORM13 3.4mm ELFORM13 1.0mm
Figure 5.6: Force-displacement curve of physical test with simulations of solid model at 6m/s
also analyzed at the same velocity to compare the behaviour with shell elements. The force in the solid model was higher compared to that with the shell elements and the physical test. The solid model was analyzed with both first order and second order tetra elements. Simulations were first completed with first order tetra element using ELFORM13 and later higher order tetras were analyzed by ELFORM16. The number of integration points in the elements is important to obtain correct material behaviour during non-linear deformation. It is known that increased number of integration points results in longer simulation times but it was necessary in this task to study the results obtained from it. Also, these formula-tions are dependent on the element’s skewness, warping and angular deviaformula-tions and with a coarser mesh it is likely to obtain poor results. The correlation study was started with 3.4 mm element size and again to obtain better results, it was reduced to 2.5 mm and 1 mm. In contrast to the shell model, the solid model predicted failure more precisely.
CHAPTER 5. RESULTS AND DISCUSSION
(a) Crack in the component (b) Solid model 1 mm
(c) Solid model 2.5 mm (d) Solid model 3.4 mm
Figure 5.7: Failure in the FE model with different element size using ELFORM16 at 6 m/s
(a) Solid model 1mm (b) Solid model 3.4 mm
CHAPTER 5. RESULTS AND DISCUSSION
Figure (5.7) and Figure (5.8) shows failure predicted by solid model with various element sizes. ELFORM13 with a coarse mesh of 3.4 mm had poor results compared to ELFORM16. The crack prediction was also not accurate with 3.4 mm. Though after reducing element size the crack prediction was better but the force was still higher. The simulation time for the model with ELFORM13 had shorter simulation time compared to ELFORM16. These simulations were completed with selective mass scaling to further reduce the simulation time. Selective mass scaling is a method of adding few terms to the mass matrix and de-creasing highest eigenfrequencies of the system, while affecting the lower frequencies as little as possible [10]. In the Figure (5.6) the force vs deflection curve is plotted with re-sults obtained from various element sizes. Also, with ELFORM16 a stiff behaviour of the solid elements can be seen. The finer mesh seems to give better results than the coarse ones. The 1mm element size analysis gave peak load of 82.79 kN while the coarse mesh of 3.4 mm had 87.10 kN of maximum force. Both were higher compared to the physical i.e. 75.46 kN. Another simulation was completed using a high ductility scatter material card. The high ductility scatter MFGenYlD+CrachFEM card has its parameters determined from more ductile locations of a cast aluminium component. There were low ductility cards avail-able at VCC but as the behaviour with solid elements were already too stiff, this cards were not used in these simulations.
The velocity was increased to 9m/s for another impact and the component was damaged the way seen in Figure (5.1). Three impact tests were done at this velocity and a total of 6 force - deflection curve were obtained from two accelerometers fitted on the impactor. Again, an average of all the curves were taken and compared with CAE results.
0 20 40 60 80 100 120 140 Displacement (mm) -1 0 1 2 3 4 5 6 7 8 F or ce (N ) ×104 Physical ELFORM16 3.4mm fs 0.2 fd 0.2 ELFORM16 3.4mm fs 0.25 fd 0.2
CHAPTER 5. RESULTS AND DISCUSSION
Failure was at the same location as with 6 m/s, but the component fractured at this velocity. Elements were deleted near the location of impact. From the experimental results a maximum force of 78.69 kN was obtained. In the first simulation with dynamic and static friction values of 0.2, the maximum force was 74.96 kN. In this simulation the forces don’t drop gradually as in the case of the experimental curve. The first impact between the impactor and the component occurs at the displacement of 60 mm, as seen in the Figure (5.9). From the FE simulations it was seen that as the force builds up the first element in the shell model fails at the displacement of 65 mm. The peak force is reached at 124 mm after the impactor bounces back. The force was less compared to the experimental test, hence the changed friction values. After the change in fs, the first element fails when the impactor displaces 59 mm. At 100 mm displacement maximum force is achieved, i.e. 77.78 kN. Iterating higher friction values gives better curves but it is to be noted that some static and dynamic friction values also gives different location of fracture.
0 20 40 60 80 100 120 140 Displacement (mm) -2 0 2 4 6 8 10 F or ce (N ) ×104 Physical ELFORM16 2.5mm ELFORM16 3.4mm
Figure 5.10: Force-displacement curve with solid model at 9 m/s
However, the shell model showed better correlation with the experimental curves com-pared to the solid model. For the solid model the maximum force with average element size of 2.5 mm was 87.41 kN and with 3.4 mm element size it was 91.14 kN. Though from Figure (5.11) it can be seen that the solid model gave better prediction of the location of the failure. There is skinning of the material near the flanges as it folds, which is represented aptly with the solid model but with the shell model the location of the failure differs remarkably.
CHAPTER 5. RESULTS AND DISCUSSION
(a) Failure in the component after impact at 9 m/s
(b) Shell mesh 3.4 mm
(c) Solid mesh 3.4 mm (d) Solid mesh 2.5 mm
Figure 5.11: Failure in the physical and the FE model with different element size at 9 m/s
After the impact tests at 9 m/s, the velocity of impact was reduced to 7.5 m/s. The pur-pose was to check the behaviour of the component by reducing the velocity. The fracture behaviour was same as seen for 9 m/s. The shell model was analyzed with 3.4 mm and change in friction values. The behaviour was similar to the impact at 6m/s and 9m/s. It-erated static friction value gave better results with shell model. With solid model, higher values of the forces were obtained. The maximum force from the physical test was 77.55 kN and the best results obtained showed 77.20 kN after fsis 0.25. The location of the failure was predicted better by the solid model, as in the case of 9m/s. Results are included in the appendix section of the report, see Appendix (A.4). Now, the purpose of reducing the velocity even further was to identify the velocity at which there was no failure in the com-ponent. However, at this velocity, i.e. 5m/s a small crack was observed in the comcom-ponent. In the simulations of the shell model it was difficult to identify a crack in the component, see Figure (A.7). Solid model predicted failure better but with a finer mesh. A plotted force - deflection curve can be seen in Figure (A.5) and Figure (A.6) where the shell model with a static friction fs0.4 shows best results. Though this simulation predicts another location of the failure compared to the physical test. Again, the solid model gave a much better prediction of the failure location, see Appendix Figure (A.7)
Upon reducing velocity even further, there was no observable crack or fracture in the component at the velocity of 4 m/s. The force deflection curve, see Figure (5.12) and Figure (5.13), gave a good correlation with the experimental curve. The best result was achieved with a value of 0.4 for fswith ELFORM16 and an element size of 3.4 mm.
CHAPTER 5. RESULTS AND DISCUSSION 0 10 20 30 40 50 60 70 80 90 100 Displacement (mm) -1 0 1 2 3 4 5 6 7 F or ce (N ) ×104 Physical ELFORM16 3.4mm fs 0.4 & fd 0.2 ELFORM16 3.4mm fs & fd 0.2 ELFORM16 3.4mm fs 0.25 & fd 0.2
Figure 5.12: Force-displacement curve with the shell model at 4 m/s
0 10 20 30 40 50 60 70 80 90 100 Displacement (mm) -1 0 1 2 3 4 5 6 7 F or ce (N ) ×104 Physical ELFORM16 2.5mm ELFORM16 3.4mm
CHAPTER 5. RESULTS AND DISCUSSION
(a) Component after impact at 4 m/s
(b) Solid mesh 2.5 mm (c) Shell mesh 3.4 mm
CHAPTER 5. RESULTS AND DISCUSSION
5.2
Results from Load case 2
5.2.1 With improved boundary condition
After Load case 1 tests, the drop tower tests were performed with the second load case setup, discussed in Section 4.1.1 and shown in Figure (4.1). The first velocity of impact was chosen as 5m/s in this case. The component showed failure as predicted by FE simulations but the length of the crack was different.
Figure 5.15: Comparison of physical and FE simulation of load case 2 at 5 m/s
In the Figure (5.15) comparison is done between the physical impact and simulated shell model after improving the boundary condition. It is to be noted that failure criteria in CrachFEM are primarily used to predict the critical location in part. It cannot be reliably used to model crack propagation once the component had already failed [4]. This can also be seen in this load case, that failure is accurately predicted but not the length(size) of crack. Shell model with friction values fs and fd 0.2 showed best results. These results were obtained by using 3.4 mm element size and ELFORM16. As ELFORM16 showed best results with Load case 1, simulation with Load case 2 were analyzed using the same.
(a) Component after impact at 5mps (b) Shell mesh 3.4 mm with fs0.2 and fd0.2 Figure 5.16: Comparison of physical and FE model at 5 m/s
CHAPTER 5. RESULTS AND DISCUSSION 0 20 40 60 80 100 120 Displacement (mm) -1 0 1 2 3 4 5 6 F or ce (N ) ×104 Physical ELFORM16 3.4mm fs & fd 0.2 ELFORM16 3.4mm fs 0.25 & fd 0.2
Figure 5.17: Force-displacement curve with shell model at 5 m/s
The force-deflection curve plotted shows that, by changing the friction values for Load case 2 in the simulations, there was not much difference in the peak forces. From the curves in Figure (5.17) it can be seen that two peaks are obtained during the test. The first peak is at 33 kN, when the first impact is made. After the impact the forces drops and again builds up until the flange of the component fails. This behaviour can be because of the other flange being stiffer than the failed one. The physical test shows maximum force of 58.12 kN but the FE simulations shows 48.46 kN with fs, fdvalues as 0.2 and 47.31 kN with fs as 0.25, fd as 0.2. The correlation between the physical test and FE simulations were considered good regarding the failure prediction.
(a) Solid mesh element size 1 mm (b) Solid mesh element size 2.5 mm
CHAPTER 5. RESULTS AND DISCUSSION
Solid models with various element sizes were also analyzed for comparison. In this case both a solid model with 1mm element size and 2.5mm element size showed good results. The fracture with both the element sizes was aptly predicted but the crack size large in CAE compared to physical. The curves in Figure (5.19) show comparison of various element sizes with the physical test. Maximum force from the CAE results are 54.5 kN which is close enough to the maximum value obtained from the physical test, i.e. 58.12 kN.
0 20 40 60 80 100 120 Displacement (mm) -1 0 1 2 3 4 5 6 F or ce (N) ×104 Physical ELFORM16 1mm ELFORM16 2.5mm ELFORM16 2.5mm scatter
Figure 5.19: Force-displacement curve with solid model at 5 m/s
After the impact at 5m/s the velocity was further reduced to 4m/s. The results are included in the Appendix, see section A.6. From the results it can be seen that both the solid and the shell model showed higher forces compared to the physical test. Failure prediction was accurate with the shell model but again the solid model had larger crack size compared to the physical test.
6
|
Conclusion and future work
Results from Load case 1 and Load case 2 demonstrates how developing a right test setup plays a crucial role for correlating results. The provided material model (MFGenYld + CrachFEM) from MATFEM, for cast aluminium components shows good correlation be-tween the the physical test and the finite element results. Load case 1 was a better choice of test setup compared to the second load case. As with Load case 2, each impact on the component made the flange bend and later it was a tedious task to remove the component from the test rig for further tests. Also, the results presented in this report are specific to the component and it’s developed test setup. So, any other cast aluminium component from the vehicle’s body should also be analyzed in a similar manner to verify the material failure model.
Comparing the results of different solid and shell models in the performed simulations reveals both shell and solid elements were accurate enough to provide a good prediction of failure in the component. The simulations with solid elements were observed to pro-vide more accurate results but were strongly mesh dependent. These solid elements with a coarse mesh showed a stiff behaviour, hence overestimating the peak forces. The refine-ment of the solid tetrahedron mesh yields more improverefine-ment in results but the simulation time increases significantly. In order to compensate for it, this work has shown that com-putational speed-up can be obtained by applying the selective mass scaling in explicit finite element analyses with very little or almost no change in the results.
In an explicit dynamic analysis reduced integration element are favourable because of the low computation cost in contrast with a fully integrated element. Though, a fully integrated element formulation is more suitable for predicting failure, which had been ob-served through the comparison with the physical tests. As it is suggested that ELFORM16 is only suitable for moderate strain, alternatively ELFORM13 proves to be more effective for large strain but with finer mesh [11]. But findings from this work indicate that using ELFORM13 even with a fine mesh revealed poor results compared to ELFORM16. It is possible that this task had moderate strains and thus ELFORM16 showed better results. Therefore, additional studies are needed to conclude which element formulation would be well suited in these scenarios. Apart from this, frictional values used in this thesis work were estimated and before any further validation of cast aluminium component, frictional coefficients should be determined experimentally.
To conclude, using the current resources available, the proposed modelling technique for a cast aluminium component, a spring tower, shell mesh of average element size 3.4 mm or solid model with a fine mesh (2.5mm) proves to be quite accurate. Selective mass scaling is an attractive technique due to its simplicity and can be used in situations where only small parts of the model need to be mass scaled meanwhile reducing the computational time.
Bibliography
[1] Djukanovic G. Aluminium vs. steel in electric vehicles : the battle goes on. Aluminium insider; 2018. https://aluminiuminsider.com/ aluminium-vs-steel-in-electric-vehicles-the-battle-goes-on/.
[2] Lambert F. Volvo clarifies electrification plan, aims for 50 percent of sales to be fully electric by 2025. electrek; 2018. https://electrek.co/2018/04/25/ volvo-electrification-plan-fully-electric/.
[3] Dharwadkar N, Adivi KP. Modelling of Engine suspension components for crash sim-ulation. Division of Vehicle safety, Department of Applied mechanics, Chalmers Uni-versity of technology; 2014.
[4] Dell H, Gese H, Oberhofer G. MF GenYlD + CrachFEM User’s theory Manual. 2014;. [5] Stouffer DC, Dame LT. Inelastic deformation of metals. 1st ed. John Wiley and Sons;
1996.
[6] Björklund O. Ductile failure in High Strength Steel Sheets. Linköping University; 2014.
[7] Xue L. Damage accumulation and fracture initiation in uncracked ductile solids sub-ject to triaxial loading. 2006;.
[8] Brenner F, Buckley M, Gese H, Oberhofer G. Influence of discretisation on stiffness and failure prediction in crashworthiness simulation of automotive high pressure die cast components. 9th European LS-DYNA Conference. 2013 Feb;Available from: https: //www.dynalook.com/conferences/9th-european-ls-dyna-conference.
[9] Hallquist J. LS-DYNA theory manual. Livermore Software Technology Corporation, California; 2006.
[10] Olovsson L, Simonsson K, Unosson M. Selective mass scaling for explicit finite element analyses. 2005;.
[11] Erhart T. Review of solid element formulations in LS-DYNA. LS-DYNA forum 2011; 2011.
[12] Cook RD, Malkus DS, Plesha ME. Concepts and Applications of Finite Element Anal-ysis, 3rd Edition. 3rd ed. Wiley Student Edition; 2004.
A
|
Appendix
A.1
Explicit analysis
In LS-DYNA, nodal displacements are calculated by using the central difference method. The equation of motion can be represented as:
[M]{ ¨D}n+ [C]{ ˙D}n+ {Rint}n= {Rext}n (A.1) Here, [M] is the mass matrix, R is the load vector for internal and external loads and [C] is the damping matrix. { ¨D} and { ˙D} are the nodal acceleration and velocity vector. The index n represents the current time steps with known nodal values. The explicit algorithm can also be given in the general form:
{D}n+1= f ({D}n, { ˙D}n, { ¨D}n, {D}n−1, ....) (A.2) The above equation (A.2) combines with the equation of motion at time step n. Time step n can be approximated by the conventional central difference equations:
{ ¨D}n= 1
∆t2({D}n+1− 2{D}n+ {D}n−1) (A.3) Equation (A.3) can be obtained from Taylor series expansion of { ¨D}n+1 and { ¨D}n−1 about time∆tn: {D}n+1= {D}n+∆t{ ˙D}n+∆t 2 2 { ¨D}n+ ∆t3 6 { ... D}n+ ... (A.4) {D}n−1= {D}n−∆t{ ˙D}n+∆t 2 2 { ¨D}n− ∆t3 6 { ... D}n+ ... (A.5) Adding equations (A.4) and (A.5) we get equation (A.3), Now, if we substitute equation (A.3) in (A.1) we get: · 1 ∆t2M + 1 2∆tC ¸ {D}n+1= {Rext}n− {Rint}n+ 2 ∆t2[M]{D}n− · 1 ∆t2M − 1 2∆tC ¸ {D}n−1 (A.6) In equation (A.6), the values for n and n − 1 are known and the internal forces can be determined. Since, explicit analysis usually solves the lower order elements, mass matrix will be a diagonal matrix. So, for a non-linear problem it can be solved directly without iterations. This diagonal matrix is one of the important features which makes the explicit method efficient and practical. Equation A.6 is conditionally stable .i.e., calculations do not become unstable unless
∆t ≤ 2 ωmax=∆tcr or, since ω= 2πf =2π T , ∆t ≤ Tmin π (A.7)
APPENDIX A
Here,ωmaxand its period Tmincorresponds to the highest frequency. The critical time step for equation A.6 is independent of damping. Instability is caused by too large time steps and is recognized as an unbounded solution that may grow in each time step [12]. The time step for this method is limited by the time it takes for the the shock wave, that arises from the loading, to transmit through the smallest element in the mesh. In the equation A.8 Lminis the smallest distance between nodes in an element and c is the speed of the sound.
∆t =Lmin
c (A.8)
Speed of sound can be calculated by
c = s
E
ρ(1 −ν2) (A.9)
Here, E is the Young’s Modulus, ρ is the mass density andν is the poisson’s ratio. The solving method does not allow the shock wave to transmit a longer distance than Lmin during the time step and regulates it by slowing the speed of sound and increasing the density. Thus, to prevent solution from becoming unstable, the mass of the element is raised.
A.2
Termination time for various models
Mesh Type Mesh size (mm) Number of Elements Time Step(ms) Element Formulation Selective Mass Scaling CPU Time (minutes)
Solid 1.0 16738173 tet 0.54 16 Yes 669
Solid 1.0 16738173 tet 0.005 16 No 4237
Solid 2.5 1003840 tet 0.54 16 Yes 71
Solid 2.5 1003840 tet 0.01 16 No 1521
Solid 1.0 16737615 tet 0.54 13 Yes 261
Solid 1.0 16737615 tet 0.08 13 No 879 Solid 3.4 432863 tet 0.40 13 No 37 Shell 2.5 59768 quad 3596 tri 0.54 16 No 25 Shell 3.4 30119 quad 2701 tri 0.54 16 No 19
Note : All the above simulation presented here are for load case 1 with 960 CPU cores apart from Solid 3.4 mm ELFORM13, which was completed with 200 CPU cores.
APPENDIX A
A.3
Results from original setup
APPENDIX A
A.4
Results from simulation of Load case 1 at 7.5 m/s
0 20 40 60 80 100 120 Displacement (mm) -1 0 1 2 3 4 5 6 7 8 F or ce (N ) ×104 Physical ELFORM 3.4mm fs fd 0.2 ELFORM 3.4mm fs 0.25 fd 0.2
Figure A.2: Force-displacement curve with shell model at 7.5m/s
0 20 40 60 80 100 120 Displacement (mm) -2 0 2 4 6 8 10 F or ce (N ) ×104 Physical ELFORM16 2.5mm scatter ELFORM16 2.5mm ELFORM16 3.4mm
APPENDIX A
(a) Failure in the component after impact at 7.5 m/s
(b) Shell mesh 3.4 mm with fs0.25 (c) Shell mesh 3.4 mm with fs0.20
(d) Solid mesh 3.4 mm (e) Solid mesh 2.5 mm
APPENDIX A
A.5
Results from simulation of Load case 1 at 5 m/s
0 10 20 30 40 50 60 70 80 90 100 Displacement (mm) -1 0 1 2 3 4 5 6 7 8 F or ce (N ) ×104 Physical ELFORM16 2.5mm fs & fd 0.2 ELFORM16 3.4mm fs & fd 0.2 ELFORM16 3.4mm fs 0.25 & fd 0.2 ELFORM16 3.4mm fs 0.4 & fd 0.2
Figure A.5: Comparison of Force-displacement curve with shell model at 5 m/s
0 10 20 30 40 50 60 70 80 90 100 Displacement (mm) -1 0 1 2 3 4 5 6 7 8 F or ce (N ) ×104 Physical ELFORM16 2.5mm scatter ELFORM16 1mm ELFORM16 3.4mm
APPENDIX A
(a) Crack in component after impact at 5 m/s (b) Shell mesh with fs0.4
(c) Shell mesh 3.4 mm (d) Solid mesh 3.4 mm
(e) Solid mesh 2.5 mm (f) Solid mesh 1 mm
APPENDIX A
A.6
Results from simulation of Load case 2 at 4 m/s
0 10 20 30 40 50 60 70 80 90 100 Displacement (mm) -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 F or ce (N ) ×104 Physical ELFORM16 3.4mm fs & fd 0.2 ELFORM16 3.4mm fs 0.25 & fd 0.2
Figure A.8: Force-displacement curve with shell model at 4 m/s
0 10 20 30 40 50 60 70 80 90 100 Displacement (mm) -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 F or ce (N) ×104 Physical ELFORM16 1mm ELFORM16 2.5mm ELFORM16 2.5mm scatter
APPENDIX A
(a) Physical test results (b) Shell mesh element size 3.4 mm
(c) Solid mesh element size 1 mm (d) Solid mesh element size 2.5 mm
(e) Solid mesh element size 3.4 mm
