Strojniski vestnik - Journal of Mechanical Engineering 57(2011)4, 345-356 D01:10.5545/sv-jme.2009.074
Paper received: 24.06.2009 Paper accepted: 27.01.2011
Fatigue Failure Study of the Lower Suspension Vehicle Arm Using a Multiaxial Criterion of the Strain Energy Density
Abdelhamid Saoudi - Mohamed Bouazara* - Daniel Marceau Department of Applied Sciences, University of Quebec at Chicoutimi, Canada
The objective of this study is to evaluate the potential of light alloy mechanical part use in automobile industry by studying their fatigue life using various parameters such as effect of suspension dynamic, excitation type, geometry and mechanical part weight. The studied part is the lower suspension arm made from 7075-T6 aluminium alloy. The strain density energy approach enables us to compare two same order tensor: the multiaxial and uniaxial cases. The random displacement excitation is obtained analytically from the power spectral density PSD. The force excitation is obtained by a simple normalisation of spectrum displacement. To avoid the use of Newton-Raphson method during the partial fatigue life calculation step in all mesh elements, a Matlab interface to identify the critical elements is developed. The strain energy density (SENER) signal of the critical element is corrected to remove anomalies by WAFO Matlab interface algorithm. Rainflow cycles are extracted using Markov formulation in order to calculate the number of signal repetitions to failure, which is calculated from Miner law. ©2011 Journal of Mechanical Engineering. All rights reserved. Keywords: fatigue failure, vehicle, dynamic, suspension, aluminium
0 INTRODUCTION
Weight reduction not only improves the slip angle between tires and the road, reaction to turns and the stability of the vehicle, but it also makes driving more effective and safe over long distances and yields lower gasoline consumption [1] to [3]. In this framework, the research goals are to study the dynamic and vibratory effects on certain aluminium alloy parts, in particular fatigue life of the lower suspension arm. Fatigue causes cracking, which develops gradually under the action of random loading repetition. These random loadings can lead to rupture by fatigue during the application of stress levels. Failure occurs when a crack reaches a critical length lc and the stress intensity factor K [4] reaches a critical value Kc. The complexity of the fatigue problem led many researchers to approach the subject using several methods. Bazergui et al. [5] based their work on the curves of fatigue SN (stress S depending on the number of cycles to failure N) and two empirical models of fatigue to prevent rupture. They report rotational bending as a cheaper standardized test. To improve the ability of predicting fatigue life of a part, they practised on standard test-tubes before installing the operational part on apparatus. In practice, stress variation is often periodic, but it
is not always sinusoidal. The average value of the stress is not null because of the static contribution from the weight of the part and premature tightening. Cervello [6] analyzed and studied the design of railway wheels with weak noise. A numerical procedure was used for calculation of the loss factor. This procedure was checked on plates by means of experimental modal analyses. It allowed better treatment of smoothness and the choice of a commercially available arrangement with feasible technology. Elmarakbi et al. [7] studied the validity of multiaxial criterion of fatigue failure based on strain energy density equivalent to uniaxial case. They performed three-dimensional finite element analysis on SAE notch axis, which is used as a test component to evaluate the criterion of multiaxial damage caused by fatigue. They also equalized energy density of multiaxial case obtained from finite element commercial software. The analytical equivalent strain density energy in uniaxial cases is calculated and compared to numerical multiaxial one, since both of them have the same tensor order (order zero). March [8] based his study on the same criterion as Elmarakbi et al., density of the energy of cracking. March performed tensiontorsion fatigue tests and compared the results with the energy density criterion. The model aimed to
predict fatigue life of rubber particles using two approaches. One focused on crack initiation, giving the history of some parameters such as stress state and deformations. Other approaches were based on ideas of the breaking process and dealt with predicting the propagation of particular cracks, giving the history of released energy rate. Crack propagation approach was developed by Rivlin and Thomas using reference [8] who applied Griffith's criterion to rubber. The difficulty in applying the crack propagation approach to rubber is that it requires advance knowledge of the initial specimen and crack state which causes the final failure. Instead of part endurance prediction models that were based on Paris law, De-Andrés et al. [9] used element cohesive law. In this approach, the creation of new surfaces is the final result of a process of gradual loss of elasticity as separation increases.
In the present paper, prediction of an automotive part's fatigue life is suggested by developing a multiaxial elasto-plastic numerical model. Knowing mechanical properties of aluminium alloy 7075-T6, fatigue failure prediction will be simulated with commercial
finite element software Abaqus and a Matlab interface. The lower vehicle suspension arm of a quarter vehicle suspension system is studied in a critical case, where the part is embedded at suspension joint. The stress state gives values of shorter and safer fatigue life by underestimating their true value. Treatment and analysis of the numerical results include:
•	development of a dynamic vehicle suspension system model,
•	development of a road irregularity model,
•	development of a numerical model using commercial finite element software,
•	analysis of results of strain density energy time evolution of the multiaxial case in all mesh elements,
•	filtering of the results to extract the critical elements in the elasto-plastic case,
•	correction of critical element signal from anomaly of strain density energy,
•	extraction of rainflow cycles by Markov algorithm,
•	calculation of partial and total fatigue life of the part from the uniaxial criterion equivalent to the multiaxial case.
Fig. 1. Vehicle suspension control and the lower suspension arm
1 VEHICLE SUSPENSION MODEL
The vehicle part studied is the lower suspension arm of vehicle suspension system. It is necessary to evaluate the dynamic behaviour of suspension system and roadway profile model. Rahnejat [10] studied the dynamics of Macpherson suspension system in a quarter vehicle. He represented spring and shock absorber as a single element with stiffness constant. In the present study a simple model is developed as illustrated in Fig. 1.
In this model, partial stability of a vehicle is ensured by control suspension system. For an action chain this system has a tire of stiffness constant Kt and for a return chain in negative feedback a spring of stiffness Ks assembled with a shock absorber Cs. The excitation force F2 of road irregularity, is balanced through the tire by the F1 negative feedback. F1 negative feedback brings back the suspension arm to its place of equilibrium linearly. For cases of small disturbances and domination of elastic behaviour of material, it produces a state of low stress. Then, the following relation is deduced:
Li „
— = — ^ x,
L L
1 2
L
(1)
Eq. (1) gives the deflexion x2, compared to equilibrium position and caused by road profile, allowing us to deduce displacement x1. Unsprung mass, such as mass of spring, tire and shock absorber is negligible compared to the dynamic stress brought into play. However, the suspension parameters such as: the spring stiffness constant Ks =16000 N/m; Shock absorber Cs = 1000 Ns/m and tire stiffness constant Kt =160 000 N/m [11].
The forces applied to the lower suspension arm by a quarter vehicle are: F2 = Kt x2 exerted on the tire by x2 and F1 — Ksx1 + Csx1 negative feedback force exerted by spring and shock absorber.
Knowing the power spectral density (PSD) experimental values, x2 is deduced and consequently, F2, enabling x1 spectrum. Once x1 is known, complete determination of F1 values is done. The road profile spectrum is a function of the vehicle speed, which is obtained from PSD. Normalisation of the road profile spectrum is necessary to keep the same frequency band and to be able to transform it into a force excitation. The
factor, which must be multiplied by PSD values, depends on numerical values of suspended mass, tire stiffness constant and shock absorber coefficient, and the damping ratio of both the shock absorber and the unsprung mass. In this study, numerical value of this factor is estimated at 0.025. x1 is numerically given by:
. _ x1(t + At) - x1 (t) x, — 1	At
The excitation caused by the road profile has a random nature estimated from power spectral density PSD. The road profile is defined by function X(t) and by assuming the road surface is a random, stationary, Gaussian and centered process: all its statistical properties are invariable with any change in argument t, that distribution law of variable X(t) is a Gaussian law, and that the average ofX(t) for any tpertaining to [0, T], is null. The stationary random process and the Gaussian law X(t) can be regarded as a periodic function in time t, of amplitude a, circular frequency m and phase $ [11]. Thus, the random process X(t) can be written according to the following equation:
X(t) = acos(at -(
(2)
X(t) can be expressed by the sum of a series of harmonic functions called Fourier series where i = 1 to N:
X(t) = Yjai cos(® J )•
(3)
The phase fa is random and distributed in the interval [0, 2n]. So the average of Eq. (3) is given by:
n
e[x (t)] — XaE[cos(®,.t - j)] — 0, (4)
thus:
i*2ff 1 £[cos(^.] = f cos^I — )dtt = 0, (5) J0	2n
Consider the mean square of X(t):
r	,	N	M	N 1
E[X(t)2]— E[£a,cos(®,i- j)£ajCos(®/-j] —^-o^ . (6)
.	j—1	. 2
Eqs. (5) and (6) show that the process is stationary. In the case of a random process, stationary, Gaussian, the following equation must be satisfied by using the mean m and standard deviation 5:
2
s2 =J+C°SXX (m)dm- m2,	(7) as m = 0, Eq. (7) become:
n	n I
¿2 = 2^Sxx (œ,)Aœ = E[x2(t)] = £-tf ,	(8)
Eq. (9) is deduced:
ai =v4sxx )a(a) .
(9)
Finally, the process X(t) according to power spectral density is given by:
N
X(t) = 4SXX(^)Aw cos(^.t). (10)
For spatial frequency range spanning from 0.01 to 10 cycles/m, power spectral density Sxx(y) from experience can be represented by an exponential function [11] given by:
S (Y) =
S (Yo) s (Yo)
/ v«1
Yo j
/ \-«2 Y
Yo
if Y ^ Yo
if Y > Yo
(11)
with y0 = — (cycles/m), q1 = 3.14 ± 0.76 and q2 2n
= 2.11 ± 0.38. y is the spatial frequency linked to circular frequency w and the vehicle speed V by:
m = 2nVY.
(12)
From Eqs. (10) and (11), power spectral density and time history displacement can be plotted as shown in Figs. 2 and 3 [11]. The road spectrum thus obtained corresponds to a vertical random displacement. However, the present study takes into account direct excitation by a random dynamic force. Thus, it is necessary to transform the random displacement into a random force while keeping the same pulsation. Therefore, the minor road spectrum is normalised. The spectrum values are multiplied by a common factor estimated from the statistical data. For a minor road, the speed of the vehicle is lower than 75 km/h (21 m/s). For 75 km/h, a 960 kg vehicle has a momentum equal to 20160 kgm/s. Thus in 0.5 second, the maximum force excitation is estimated to be 40320 N. This force is distributed equally on
four quarters of the vehicle's suspension system. In the undamped cases, the vehicle weight opposes about 3750 N in the return chain to the maximum 4000 N random excitation value. The maximum value of the force of excitation through a tire of stiffness 160000 N/m is about 4000 N. Figs. 4 and 5 depict a comparison between damped and undamped spectrums. It also illustrates opposition phases between excitation F2 and force F1 of the return chain in vehicle control suspension system.
2 ELASTO-PLASTIC NUMERICAL MODEL
Aluminium alloys have an elasto-plastic behaviour with a plastic contribution of about 36% in certain cases (Table 1). Gbadebo [12] and [13] proposed an analytical solution for elasto-plastic case, where the matrix of total tensor deflection increments is the sum of elastic and plastic contributions. They extrapolated the nonlinear uniaxial behaviour of the material in multiaxial case to calculate strain energy density. In the plastic behaviour Gbadebo et al. [12] and [13] have developed the endochronic theory of plasticity. This approach is the same as the one used in the commercial finite element software Abaqus.
Table 1. Elasto-plastic behaviour of some aluminium alloys
Alloy	Code	Ultimate limit [MPa]	Yield stress [MPa]	Ratio [%]:
				
Aluminium alloy	2024-T351 2024-T4 Al 7075-T6 Al	455 476 578	379 303 469	16.70 36.34 18.85
This model takes into account two cases of nonlinearity: non-geometric linearity (variation of the geometrical configuration of the system in time) and material nonlinearity expressed by elasto-plastic coupling. The numerical diagram of integration depends on its stability, computing time and accuracy of the method. In this study implicit integration scheme is used [14] and [15]. The loading and boundary conditions specify the
	
	
	
	
	
	
	
Frequency [Hz]
Fig. 2. Power spectral density PSD as a function offrequency
2	3
Time [s]
Fig. 3. Displacement as a function of time
places of excitation and the degrees of freedom imposed on the system. The mesh takes into account the stress state and critical points: load location, and boundary conditions area. The analytical equivalent strain density energy in uniaxial cases is calculated, and compared to the numerical multiaxial one, since both of them have the same tensor order (order zero). The critical element is filtered through all the elements of the mesh, giving the maximum sum of positive variations. The critical element loading signal is corrected from anomaly in order to apply the Markov algorithm conditions. The rainflow cycles are extracted using Markov method. Due to the
non linearity of the Mansson-Coffin equation, partial fatigue life is calculated by the Newton method. The number of random cycle repetitions leading to part failure is determining at final stage using Miner law. To determine local stress state and fatigue criterion, Lachowicz [16] used multiaxial strain energy density criterion without specifying the calculation of density of plastic deformation energy. Pan et al. [17], Lee et al. [18] and Li et al. [19] compared experimental results of multiaxial fatigue with models of multiaxial criterion of strain energy density already existing as the criterion of critical plan which generalizes the Smith Watson Topper SWT uniaxial model.
Fig. 4. Road excitation F2 in undamped (- - -) and damped modes (—)
Fig. 5. Damped mode: road excitation F2 (—) and feedback chain F1 (- - -)
The local stress plane state, characterized by: o2, VP, and X = o2 / o1, that represents the two principal stresses, phase angle between o1 and local axis OX and the ratio of biaxiality [20], respectively. The stress state in a material can take the following aspects:
• uniaxial case: X = a2 / a1 = 0 and yp are constant in time. In this case an uniaxial model is used such as the SWT criterion;
•	proportional multiaxial case: varies in time and yp is constant. In this case the effective strain criterion is used;
•	nonproportional multiaxial case: both
X = a2 / 01 and yp vary in time. In this case the critical plane is used.
In all three cases the uniaxial criterion equivalent to the multiaxial case as suggested by Elmarakbi et al. [7] can be used. Indeed the uniaxial
-5.0E+03
-1.0E+U4
-1.5E+04
-1.0E+03
-2.0E+03
-3.0E+03
-4.0E+03
-5.0E+03
criterion of strain density energy equivalent to the multiaxial one, is used in all cases of X and fp, because it is a complete integral calculation in case of uniaxial analytical integration and all components of shear and distortion in multiaxial numerical case are calculated.
Since aluminium alloys present more significant plastic than elastic behaviour, the nonlinear profile is given by cyclic Ramberg-Osgood behaviour law as:
^ = — + (—) '
f
where K ' =—— is coefficient of endurance limit,
b f
n' = b/c cyclic hardness exponent, of coefficient of fatigue strength, c exponent of fatigue ductility, b fatigue strength exponent and ef coefficient of fatigue ductility. The strain energy density, relative to multiaxial case, can be expressed as:
i
(13)
a 2E n'+1 { K'
A resolution of Eq. (19) gives necessary constraint to produce same energy density as in the uniaxial case. If replacing the stress value in Eq. (16), gives the corresponding deformation, then fatigue life can be predicted by using the following Manson-Coffin equation:
f = E (2*(2N)c , (20)
where Ae = em
- en
is strain interval, o
f
coefficient of fatigue strength, N number of cycles to failure, E Young modulus, C exponent of fatigue ductility, b exponent of fatigue strength and ef coefficient of fatigue ductility.
3 ELASTO-PLASTIC NUMERICAL MODEL UNDAMPED RIGID CASE
Us = Ua or £Sjdey o^ . (14)
In the present study, the model of three-dimensional elastic finite element (FEM) analysis is used. The value of strain energy density obtained from the analysis in three-dimensional FEM, is equalized to the strain density energy of uniaxial case. The mathematical formulation of total uniaxial deformation energy per unit of volume, obtained by an exact integration is represented as [7]:
U a = Ua.
u=j;
oadea
We know that:
e„ =^ + n
thus:
E K '
ds„ 1 a„ n - = — + ■
daa E 1 n' k '
(15)
(16)
(17)
Eq. (15) becomes:
Ua =J«
1 I °an
E I
n k 'n
da„
(18)
and finally, the following relation is obtained:
According to Rahnejat [10] the suspension arm cannot support a negative feedback force of a quarter vehicle, which is superior to 3750 N, corresponding to the rigid case. This part is subjected to random excitation of the road described by PSD and to quarter weight of vehicle which opposes it. The right end of the part is linked to a suspension joint in order to allow a rotation around y-axis as illustrated by Fig. 6. The suspension arm is 2.5 cm thick, 50 cm long and 30 cm wide. As shown in Fig. 6, the concentrated forces are replaced by uniform surface pressures following:
P=f=F
(21)
Fig. 6. Loading and boundary conditions in a portion of the Al 7075-T6 alloy part
Abaqus-6.4 offers us tetrahedral elements. A controlled mesh is necessary to refine the
number of elements close to critical zones as illustrated in Fig. 7. Moreover, due to bending in a stress plane state, the mesh must be refined in thickness to represent state of stress plane better. Vertical partitions are created in the thickness. The linear elements of reduced integration tolerate distortion, therefore, a refined mesh of these elements is used in any simulation where the distortion levels can be very high. Thus, linear tetrahedral elements are chosen in an implicit integration stable scheme.
Fig. 7. Mesh of the lower suspension arm
The need for determining the critical element and its coordinates led to choose a strategy which allows isolation of the element having maximum positive variation. Contrary to the static case, the rough maximum value of strain energy density does not necessarily correspond to the critical element. This allows an application of the Newton-Raphson algorithm to only one element instead of applying it to all the mesh elements in the structure. In order to extract the number of cycles to failure, the nonlinear Mansson-Coffin equation is used. This filter generalizes the case where excitation is multipoint and shifted in time, giving a tangle of mesh material element signals. This filter is based on the following algorithm:
xc and zc are the critical element coordinates having:
tf
M(ic, j) = Mnc (Max(£ AU(m, n))),
where:
m = time parameter, k = simulation duration,
tf =^2=
n = spatial parameter to make ic and jc in the same column nc
AU(m,n) = U(l,n) - U(p,n), when U(l,n) > U(p,n).
Indeed, it is a variation of strain energy density which is involved in failure by fatigue in case of the dynamic stress and not absolute value of strain energy density as in the static case. Fig. 8, illustrates strain energy density profile (SENER) within the part. The SENER is critical in neighbourhoods of boundary conditions. Fig. 9, shows time evolution of strain energy density in centroid element 169. The critical element is located in the area where the part is embedded.
Fig. 8. Instantaneous strain energy density in the 25 mm specimen
k!
2!( k - 2)!
Fig. 9. Strain energy density time evolution of critical element no. 169
The critical element strain energy density time evolution signal has a random profile. To apply Miner law, rainflow cycle must be extracted.
m=0
Each rainflow cycle is described by direct passage of power spectral density PSD at each mesh element and was then solved by a rigorous theory based, on one hand, on the definition of a rainflow cycle suggested in 1987 by Rychlik [21] and on the other hand, on the theory of Markov chains. Indeed, a rainflow cycle, as illustrated in Fig. 10, can be mathematically characterized in the following way:
Let us consider o(t) where t e [0, T] and the stress maximum Mi of level K occur at time ti. The extents (m-, Mi) and (m, m+) are defined, where :
•	m- is the minimum of a(t) and Mi is between last passage to negative slope of a(t) and the maximum Mi. This minimum is on the left of Mi and occurs at time t-,
•	m+ is the minimum a(t) which is between Mi and the first passage to positive slope of a(t) by the level K. This minimum is on the right of M, and occurs at time t + .
Fig. 10. Mathematical characterization of rainflow cycle
If there is no passage of s(t) by the level K before or after time ti, then t- = 0 and t+ = T. The rainflow extracted at time ti is then defined either as the extent (mf, Mt), or (M , mf). This
minimum mfc is given by applying the condition:
max(mi , m+ ) f
m,
t; > 0 else
(22)
Fig. 11 shows the critical element SENER signal and rainflow cycles counting. When loadings are composed of various cycles of various amplitudes and various average values, it is necessary to measure total damage produced by these cycles by using the laws of damage calculation, which were developed from linear damage rule suggested by Palmgren in 1924 [18]. The mathematical formulation under which it is currently known was proposed by Miner in 1945 and it is expressed as Eq. (23).
1=1 Nt
(23)
Random stress history is described as a sequence of blocks of constant amplitude. Each block i is composed of N cycles of amplitudes. Partial fatigue life Ni corresponding to this stress amplitude is determined from the strain energy density approach. Failure is predicted when damage D is equal to 1. It is thus necessary to find the number of times that random loading is repeated, so that D is equal to unit. It is necessary to find the Bf number which one must multiply by D to reach rupture. Bf is calculated using Eq. (24).
mrfc =
Bf =\
' Ni
since,
BfD = Bf (V n) = 1.
f	, N,
(24)
Table 2 shows the metallurgical fatigue parameters of aluminium 7075-T6 alloy used in lower suspension arm.
Table 2. Metallurgical parameters of aluminium 7075-T6 alloy
Young modulus E [MPa]	n'	[MPa]	s'f	b	c
71000	0.106	1466	0.262	-0.143	-0.619
The relative cycle fatigue life part's is estimated as 8.86xlO11. Indeed critical element maximum Von Mises stress is about 31 MPa and corresponding maximum deformation is 4x10-4 These slightly low values can be further attenuated to make part safer by reinforcing it at critical points. Optimization of the part shape will make it possible to increase the number of repetitions needed to failure, make vehicle safer and further decrease vehicle weight in non critical places.
4 ELASTO-PLASTIC NUMERICAL MODEL, DAMPED AND RIGID CASES
Damping suspension largely decreases the applied force in the feedback chain middle part section. Vehicle weight is considered at 3750 N, not exceeding absolute value of F1max = 1200 N in the case of damping by a spring and shock absorber. This force is opposed in feedback chain by road random excitation deadened by tire stiffness Kt = 160000 N/m, and maximum value F2max = 4000 N. Feedback chain F2 is useful to stabilize the vehicle, but will not compensate F1 in terms of stress since difference between F1 and F2 is higher in case of damping than in rigid case. This is confirmed by a strain energy density (SENER) packing illustrated in Fig.12. The critical element filtered in both cases is element 169 located at the level of embedding as illustrated in Fig. 8. The model part and mesh are the same in both cases. Consequently, the fatigue life part in the damping case will be lower than in the rigid case. This shows that despite the comfort gain when using damper and spring, the fatigue life value of certain automotive parts such as a lower suspension arm of a vehicle decreases due to load order.
[R101]
6.00
4.00
2.00
a)
0.00
i
Jj
0.00
1.00
4.00
5.00
2.00 3.00
Time
Fig. 12. Critical element strain energy density time evolution; a) 25 mm rigid case,
b) 25 mm damping case
1
n
5 CONCLUSION
The principal objective of this research was to develop a hybrid model in order to study
the potential of an aluminium alloy as the new
material for a lower suspension arm of a vehicle. Analytical and numerical model were developed to simulate dynamic behaviour of a suspension system as well as the state of stress and strain energy density in the lower arm of vehicle suspension. The spectral aspects of fatigue and dynamic behaviour of a vehicle suspension system were studied. The multiaxial criterion of the uniaxial case strain energy density equivalent to the multiaxial one represents a rational approach. Moreover, strain energy density criterion, independent of average value of loading, is more practical than the Morrow model, which requires corrections due to effect of the average value of loading. The Morrow model needs some correction due to mean value, because it is a tensor of order 1, which depends on orientation.
In this study, the maximum value of force transmitted by tire and the feedback chain of the force were evaluated. Then a model of negative feedback and forward path of the vehicle's lower suspension arm, which is subjected to very important dynamic stresses, was established using the road irregularity model power spectral density.
To filter the critical element and to extract the fatigue life, Matlab interface was generated to locate automatically the critical elements without applying Newton-Raphson algorithm in every mesh element. This filter produces the case where excitation is multipoint and shifted in time giving a tangle of signal mesh elements.
The elastoplastic nonlinear case was described using the Ramberg-Osgood uniaxial relation binding stress to deformation. The commercial finite element software Abaqus version 6.4 extrapolates this relationship in the multiaxial case. The sensitivity of certain types of mesh elements makes the execution of the stable implicit integration scheme difficult. Ramberg-Osgood and Mansson-Coffin relation allow us to describe two areas: high-cycle fatigue area and low cycle fatigue area. In the case of damping then in the rigid case even despite the comfort gain when using a damper and a spring, the fatigue life value of certain automotive parts such as a lower
suspension arm of a vehicle decrease due to load order.
In the future, a complete model should be developed to study the various cases, most importantly the effect of the "Jounce bumper" rubber support of suspension.
6 REFERENCES
[1]	Bignonnet, A., (2001). A Global Approach for Vehicle Weight Reduction. Mechanics & Industries. 2, p. 173-180. (In French).
[2]	Miller, W.S., Zhuang, L., Bottema, J., Wittebrood, A.J., De Smet, P., Haszler, A., Vieregge, A. (2000). Recent development in aluminium alloys for the automotive industry. Materials Science and Engineering, A280 p. 37-49.
[3]	Saito, M., Iwatsuki, S., Yasunaga, K., Andoh, K. (2000). Development of aluminium body for the most fuel efficient vehicle. JSAE Review, vol. 21, p. 511-516.
[4]	Wilfred, K., Jean, P.M., Gérald, Z. (1991). Introduction of Material Sciences, 2nd edition Presse Polytechniques et Universitaires Romandes (EPFL Press). (In French)
[5]	Bazergui, A., Bui-Quo, T., Biron, A., McIntyre, G., Laberge, C. (1993). Résistance des matériaux. École Polytechnique de Montréal, 2nd ed., Montreal.
[6]	Cervello, S., Donzella, G., Pola, A. (2001). Analysis and design of a low-noise railway wheel. Proceedings of the Institution of mechanical Engineers. Part F, Journal of Rail, vol. 215, p. 79-92.
[7]	Elmarakbi, A., El-Hage, H., Bhattacharjee, S. (2002). Multiaxial fatigue crack initiation by strain energy density using finite element method. CSME Forum 2002, Kingston.
[8]	Mars, W.V. (2002). Cracking energy density as a predictor of fatigue life under multiaxial conditions. Rubber Chemistry and Technology, p. 1-17.
[9]	De-Andrés, A., J., Pérez, L., Ortiz, M. (1999). Elastoplastic finite element analysis of three-dimensional fatigue crack growth in aluminium shafts subjected to axial loading. International Journal of Solids and Structures, vol. 36, p. 2231-2258.
[10]	Rahnejat, H. (1998). Multi-Body Dynamics: vehicles, machines and mechanisms. Society of Automotive Engineers, Warrendale.
[11]	Bouazara, M. (1997). Étude et analyse de la suspension active et semi-active des véhicules routiers Ph.D. thesis, Laval University, Laval.
[12]	Gbadebo, M., Meera, O., Sing, N.K. (2006). A comparison between analytical models that approximate notch-root elastic-plastic strain-stress components in two-phase, particule-reinforced, metal matrix composites under multiaxial cyclic loading: Theory. International Journal of Fatigue, vol. 28, p. 910-917.
[13]	Gbadebo, M., Meera, O., Sing, N.K. (2006). A comparison between analytical models that approximate notch-root elastic-plastic strain-stress components in two-phase, particule-reinforced, metal matrix composites under multiaxial cyclic loading: Experiments. International Journal of Fatigue, vol. 28, p. 918-925.
[14]	Sun, J.S., Lee, K.H., Lee, H.P. (2000). Comparison of implicit and explicit finite element methods for dynamic problems. Journal of Materials Processing Technology, vol. 105, p. 110-118.
[15]	Rebelo, N., Nagategaal, J.C., Taylor, L.M. (1992). Comparison of implicit and explicit finite element methods in the simulation
forming processes. Numerical methods in industrial forming processes: Proceedings 4th international conference , p. 99-108.
[16]	Lachowicz, C.T. (2001). Calculation of the elastic-plastic strain energy density under cyclic and random loading. International Journal of Fatigue, vol. 23, p. 643-652.
[17]	Pan, W.F., Hung, C.Y., Chen, L.L. (1999). Fatigue life estimation under multiaxial loadings. International Journal of Fatigue, vol. 21, p. 3-10.
[18]	Lee, B.L., Kim, K.S., Nam, K.M. (2003). Fatigue analysis under variable amplitude loading using an energy parameter. International Journal of Fatigue, vol. 25, p. 621-631.
[19]	Li, B., Reis, L., De Freitas, M., (2006). Simulation of cyclic stress/strain evolutions for multiaxial fatigue life prediction.
International Journal of Fatigue, vol. 28, p. 451-458.
[20]	Haiba, M., Barton, D.C., Levesley, P.C. (2003). The development of an optimization algorithm based on fatigue life. International Journal of Fatigue, vol. 25, p. 299-310.
[21]	Pitoiset, X. (2001). Spectral Methods for Fatigue Analysis of Metallic Structures under Multiaxial Random Loading. Ph.D. thesis, Universite libre de Bruxelles, Brussels. (in French)