Validation and optimization of the leaf spring multibody numerical model

In this paper, a multibody approach for leaf springs modelling that may be used for full vehicle simulations is presented. Two different configurations of the leaf spring are investigated. The first one reproduces an actual mounting configuration of a leaf spring suspension in a vehicle, including the shackle. In the second configuration, the shackle is eliminated by mounting the leaf spring on linear bearings that are able to move on a guide rail. This allows for better spring rate evaluation, without the influence of shackle’s length and position. Both configurations were modelled in commercially available software LMS Virtual.Lab Motion and tested experimentally. Results of investigations have proved that developed model correctly reproduces behaviour of the leaf spring for dynamic simulations, with better convergence achieved for configuration without shackle. Initial results for static tests have shown significant error in behaviour of the numerical model and experimental test rig. In order to minimize observed errors, the optimization process of the model parameters was performed in Noesis Optimus by means of single-objective optimization. It was proved that presented approach may be used for proper modelling of the leaf spring behaviour, and this type of models may be effectively used in full vehicle simulations.


Introduction
Computer-aided engineering methods are playing an important role in design process of many different products in modern engineering. Over the course of years, the demand on complexity of developed systems has increased, while time dedicated on the development is constantly reduced. This is especially visible in electronic and automotive industry. Increasing demands from clients to deliver a complete and reliable product in as short time as possible force developers to investigate new approaches which will not only reduce the development time of the new product, but would also allow to reduce costs associated with it. An example of such approach is 'virtual prototyping'. Virtual prototyping allows the investigation of many different designs of the new products by means of virtual prototypes, while the number of physical prototypes is limited only to most promising configurations. Such approach results in reduced time and costs of the development [10]. The creation of simulational model is a compromise between complexity of the model and performance. The model has to be complex enough, to provide accurate and valuable information about the performance of the system, while simultaneously be simple enough in order to minimize the computing time. These two objectives are usually contradicting themselves, and finding a compromise is not an easy task. In order to achieve satisfying results in a reasonable time, different modelling techniques with various degree of complexity have to be investigated.
Another factor that is very important especially in automotive and transport industries is the safety of the new designs. New types of the vehicles become faster, lighter and more comfortable, but they have to be also safer. One of the most important system that plays a crucial role in safety is the suspension [4]. One of the type of suspensions used in road or rail vehicles is based on leaf springs. This type of suspension is one of the oldest but is still commonly used in many types of vehicles. Although leaf springs are now rarely used in the passenger cars (i.e. Chevrolet Corvette), they are still widely used on sport utility vehicles (SUVs), four-wheel drive vehicles, all-terrain vehicles, load-carrying vehicles and trains [14,17]. Since leaf springs are nonlinear systems which dissipate energy through inter-leaf friction, they cause many problems in modelling [5,12]. Many of the CAE tools used for full vehicle modelling in automotive industry are based on multibody approach for vehicle dynamics' simulations [2,16], and methods used for modelling of the leaf spring suspension should be compatible with this approach.
There are several methods to model the behaviour of a leaf spring, each having some advantages and disadvantages. First attempts at modelling the leaf spring behaviour were based on incorporating a combination of springs, dampers and friction elements in order to create an equivalent system, which would replicate the behaviour of the leaf spring. These types of models are replicating each physical phenomena of the leaf spring individually, i.e. modelling the stiffness by the spring and hysteretic behaviour by the damper [5]. An example of such system could be two springs connected in series, with a friction model between them or three springs connected in one point, with two of them placed horizontal and one vertical to represent longitudinal and vertical force-displacement characteristic [2,5]. Next approach is based on discretization of the leaf spring by rigid elements, which are connected by revolution joints with associated torsional springs or dampers [2,5,8]. A typical example of such approach is so-called SAE X-link model. Such models are very simple, have only a few degrees of freedom and are relatively easy to implement in typical multibody software. However, they cannot represent the roll behaviour and lateral loads correctly, and also hysteresis from the inter-leaf friction is rarely included [5]. Most commonly, the SAE 3-link model is used, but there are also publications in which leaf springs were modelled by 5 or even 10 links [6,11] with the latter using spherical joints to connect different segments. Another approach that can be used for leaf spring modelling is finite elements method (FEM). This method is giving accurate prediction of the spring behaviour, allows for modelling of multileaves springs and inter-leaf friction in order to replicate the hysteretic damping. However, FEM models are usually computational heavy, and they are very hard to be incorporated in a multibody simulation of complete vehicle [5]. The last approach involves a flexible multibody model of the full leaf spring. It may be regarded as a combination of discrete and FEM model, as each leaf of the leaf spring is modelled by a series of segments connected by beam force elements with sectional properties of the leaf. This approach allows for incorporation of contact elements with friction into the model, in order to simulate the hysteretic properties of the leaf spring [5]. Beam models provide very good accuracy in comparison with discrete models, but they are also less computationally demanding than FEM models and can be implemented in multibody simulation of a complete vehicle.
The main objective of the study presented in the paper is validation of the leaf spring model built-in multibody environment LMS Virtual.Lab Motion and its optimization in order to tune up its properties to match with experimental results.
The paper is organized as follows. In Sects. 2 and 3, modelling method of the leaf spring in LMS Virtual.Lab Motion and model development are presented. Section 4 consists fundamental information about multibody dynamics equations used in model development. In Sects. 5 and 6, two test configurations are described and initial numerical and experimental results are presented. Section 7 consists of description of optimization procedure and numerical results obtained with optimized leaf spring model are discussed. Summary and conclusions drawn from the presented investigation are presented in Sect. 8.

Modelling method of the leaf spring implemented in LMS Virtual.Lab Motion
The popularity of the leaf spring suspension is due to its simple design, cheap manufacturing and simple assembly process. Despite this simplicity, it presents several challenges from modelling point of view. Namely, the leaves are in contact between each other, they slide relatively to each other with friction, there can be loss of contact in part of the leaf, the suspension undergoes large displacement, and additional elements can be present (such as clips, shackle and clamps). In the LMS Virtual.Lab Motion environment, it is possible to represent each single phenomenon by combining several basic elements: 1. The flexibility of each leaf can be modelled by discretizing the leaf and concentrating the flexibility at discrete locations along the leaf. 2. The contact between leaf can be modelled in different ways by exploiting the geometry of the leaf and a contact detection algorithm based on such geometry (CAD contact that is rather generic but also time consuming) or with more tailored algorithms that take advantage of the knowledge of the kinematics (such as sphere to extrude contact that assumes that one leaf is an extruded solid and the counterpart can be modelled as a sphere at the location of the contact. This method is more limited but guarantees low computational time). 3. Friction forces come together with the contact ones. 4. Clips can be modelled as rigid bodies with contact. 5. Connections to the frame, if compliant, are modelled by means of bushings.
In particular, to model the flexibility of each leaf, the element is discretized along its length in a number of user defined segments represented by rigid bodies (carrying the inertial properties) and connected by force elements that in the LMS Virtual.Lab Motion environment are named beam elements. These forces are characterized by a 6 × 6 constant stiffness and damping matrices connecting two axis systems. The relative linear and angular displacement (and velocity) between the two axis systems cause a force that is proportional to the stiffness and damping.
The element of the stiffness matrix is derived from the geometrical properties of the element and a simple linear beam theory [13]: where E is the Young's modulus, A is the cross section, l is the free length of the beam, G is the shear modulus and I x, y, z are the moments of inertia. The damping matrix can be defined in two ways, by means of a viscous damping (the user defines the damping ratio ξ ): or by means of structural damping (preferred choice in this case): For the contact between leaves, it has been chosen the sphere to extrude method (see Fig. 1). An equivalent sphere is placed in the location of the most probable point of contact on the lower leaf. The upper leaf is defined as the extruded surface that will contact with the sphere. The contact force is based on the penetration of the sphere into the extruded surface and its velocity (normal to the surface) and is computed according to a Hertzian formulation. The friction forces are computed using the velocity tangential to the surface. The parameters of the contact element are usually difficult to estimate or measure; therefore, in practice it is required to use them as parameters in a model updating procedure when experimental data are available.

Development of leaf spring model using multibody approach
For numerical and experimental tests, a front leaf spring from Suzuki Samurai, a small off-road vehicle, was chosen. The leaf springs installed in this car are typical semi-elliptical, non-tapered leaf springs consisting of four leaves made of steel. Simple design makes it an ideal candidate for model development.
First to develop the multibody model, the geometry of the assembled and disassembled leaf spring was acquired. The leaf spring was scanned using Gom Atos II optical three-dimensional scanner. Obtained geometry was exported as STEP/STL file to the CAD environment and is presented in Fig. 2. In CAD software, the   geometry measurements were taken including dimensions of the individual leaves, curvatures of each leaf and assembled leaf spring. Geometrical measurements showed some variations of leaves thickness throughout their length. Further investigation revealed that there were two sources of these imperfections-scan errors connected with the 3D scanner resolution and actual geometrical variation observed on leaf spring. Since leaf springs model is built by discrete solid segments connected through force elements, some simplification in regard to leaf thickness was made. An average value was used to ensure that thickness of all segments was constant for individual leaf model. Geometric properties of the leaf spring used in model are presented in Table 1.
A numerical model was created in LMS Virtual.Lab Motion using a dedicated leaf spring modelling interface: LMS Leaf Spring Tool (Fig. 3) that uses the SAE J1123 (Leaf Spring for Motor Vehicle Suspension Made to Metric Units) as reference for the definitions and nomenclature [15]. The model created with this interface is created based on three specified attachment points: location of the first leaf's connection to the chassis, the centre of the spring where all individual leaf are bracketed and attachment point of the shackle to the first leaf. Since LMS Leaf Spring Tool does not take leaf spring eyes into account, they were included in the model later, by connecting them to the earlier defined attachment points and replacing default Revolute Joints  with Bracket Joints. Next, rest of the test rig's geometry, described in the following section, was attached to the eyes through Revolute Joints, moving the rotation centre from first leaf ends to the its proper position in the eye (Fig. 4). Necessary changes to the attachment points coordinates were calculated using 3D scans and a CAD environment.

Multibody dynamics equations
A mechanism is modelled in LMS Virtual.Lab Motion by mean of a Newton-Euler formulation. A maximal set of coordinates including the translational coordinates x, y and z and the four Euler parameters e 0 , e 1 , e 2 and e 3 is used to describe the position and orientation of each body of the mechanism in space.
This type of formulation implies that a set of algebraic constraints needs to be imposed to remove the degrees of freedom that are linked by means of kinematic relationships (for instance, joints such as revolute, translational and spherical). This set can be expressed as the following holonomic constraint: The equations for the algebraic constraints are added to the Newton-Euler equations of motion, and a set of Lagrange multipliers, λ, is used. Moreover, in many practical applications, a mechanism is actuated by controls or hydraulic systems or other types of external active forces. The behaviour of such an element can be usually described by a set of first-order state equations in the state variable, χ . This variable could be the pressure in a hydraulic jack or current controlling a valve. When assembling all these state equations with the equations of dynamics and the constraints, the system is described by the following system of differential algebraic equations (DAE): where M is the mass matrix, Q a is a vector of applied forces and v is the vector of velocities (derivatives of 5). The formulation is fully described in the book of Haug [3].

Virtual and experimental tests
To identify properties of the leaf spring, two experimental setups (Fig. 5) were investigated during static and dynamic tests. First of the test configurations resembled the actual mounting configuration of the leaf spring as mounted in the vehicle (reversed) what is presented in Fig. 5a. The fixed end of the leaf spring was connected to the chassis, modelled as U-shape beam, by revolute joint. The second eye of the spring was connected to the shackle, which was then attached to the chassis. In both connections, revolute joints were used. However, the presence of the shackle has a strong influence on leaf spring characteristic, and it should not be used in situation when the characteristics of the spring alone are to be obtained [1]. Therefore, another test rig configuration was developed, where the shackle was eliminated and both spring eyes were mounted by revolute joints on linear bearings moving on a guide rail. This configuration allows for better spring rate evaluation, without the influence of shackle's length and position (Fig. 5b). Experiments on each of the test rig consist of static and dynamic tests. During static tests, the leaf spring was loaded axially by actuator until value registered by the force sensor reached 2 kN. Afterwards, the leaf spring was unloaded. Force and displacement during both stages of the test were measured and saved. The time of complete cycle was 120 s, loading and unloading lasting for 60 s each. Based on performed measurements, force-deflection static characteristics, for loading and unloading, were determined. In the second part of the experiment, the response of the leaf spring to dynamic operation was tested by means of cycling loading using sine wave signal with constant amplitude and frequency. Total number of loading cycles was set to 50, and test was repeated four times for each of the dynamic tests. Three different amplitudes of sine wave signal with 1 Hz frequency were used to excite the system. Also, one additional experiment with frequency 2.5 Hz was performed to allow better assessment of the acquired data.

Virtual test rig
Both of the test configurations were modelled and simulated in LMS Virtual.Lab Motion (Fig. 6).
As described in the Chapter 3, the numerical model of the leaf spring was used to create two different mounting configurations. The relative rotation of the spring eyes against supports and the shackle was modelled To replicate actual laboratory test rig configuration (Fig. 7), additional elements were added to the model: a plate rigidly attached to the leaf spring at its clamping area and an actuator that allows for Z-axis force excitation of the test rig. CAD contact element was used to model the contact between the actuator and plate. In created model, the friction forces between the plate and the actuator, as well as in the leaf spring bushings, were omitted. To excite motion of the virtual actuator, real displacement signal recorded during experimental test was used.

Experimental setup
In order to verify numerical results obtained from LMS VirtualLab Motion, experimental test rigs that reflect both of the configurations were built. Both of the test setups were carried out on Instron 8872 load frame equipped with 100-mm stroke hydraulic actuator, capable of generating forces up to 25 kN in both directions. Both of the test configurations are presented in Fig. 7. Proposed test setups consisted of a hydraulic actuator (1), pusher (2) used to load the spring, t-slot table (3), force transducer of 10 kN range (4), tested leaf spring (5), mounting U-shape beam (replaced in the second experiment by linear bearings and guide rail) ( (7). The purpose of the wooden board was to provide a flat surface for the actuator to act on.

Results of the numerical and experimental tests
During numerical and experimental tests, load-deflection characteristics for each of the configuration and frequency have been obtained. Because of the strong influence of the shackle, each of the tests was analysed separately. Figure 8 presents direct comparison of the load-deflection characteristics for configuration with and without shackle (left and right, respectively). Experimental results confirmed influence of the shackle on leaf spring characteristic: without shackle there is strong reduction in stiffness and amount of energy dissipated per cycle. Developed model correctly estimates an overall shape of experimentally identified force-displacement curve, including slightly nonlinear behaviour of the suspension when the shackle is present (Fig. 9). This nonlinearity is caused by changes in system configuration. Leaf spring characteristic is influenced by shackle's length, its angle, the camber of the spring and load applied to it. Since all of those parameters (except length) are changing during experiment with force applied to the leaf spring, the resultant spring characteristic is also changing as the test progresses [1]. Replacing the shackle with linear bearings has eliminated this effect on the spring characteristic. Damping capabilities of the spring are directly related to amount of energy dissipated during loading and unloading, which can be calculated by integration of both virtual and experimental hysteresis loop according to expression 8  where F force recorded during cycle, x displacement of the spring under load. Initial simulations have failed to correctly recreate hysteretic behaviour of the spring; dissipated energy per cycle is significantly lower in comparison with experimental data. Also, slight errors in stiffness estimation were present for both cases. A detailed comparison of achieved results is given in Table 2.
In the second set of tests, the influence of the dynamic load on the spring was estimated. In this case, the results achieved for all tested amplitudes and frequencies show much better estimation of the dissipated energy per cycle. Still, much higher errors can be observed for the model with the shackle. Errors in energy dissipation for most cases are on around 10 % level, as may be observed in Fig. 10.
It must be noted that friction in leaf spring bushings, shackle as well as friction between wooden board and the pusher was omitted during modelling phase. This may be one of the sources of better convergence of the initial results for second configuration (Fig. 11).
The difference in energy dissipation between static and dynamic simulation is caused by the model of friction forces present in contact elements used to simulate interaction between leaves. The evaluation of friction is speed dependent; thus, a very slow static test (over 2 min) gives improper inter-leaf friction evaluation leading to minimal energy dissipation. For dynamic loading with frequency of 1 Hz, the threshold value of the Transition Velocity parameter is exceeded allowing for proper friction evaluation. No speed dependency is observed above this value, as for 2.5 Hz energy dissipation is at the same level as for 1 Hz (Fig. 12).   Experimental test results have shown that leaf spring characteristics during dynamic tests are unaffected by changes in loading speed (Fig. 12a). Despite 2.5 times increase in loading frequency, hysteresis loop shape is preserved and amount of energy dissipated per cycle is virtually on the same level, as can be seen in Table 3. The same effect may be also observed in simulation results (Fig. 12b), with energy dissipation within 10 % error tolerance. Those small errors could be eliminated either by manual tuning of the friction coefficient, implementation of the friction between actuator and plate and by optimization procedure. In both dynamic and static simulations, there are present some errors in leaf spring stiffness evaluation. However, those errors are small, not exceeding 10 %, and may be caused by improper leaf thickness in the model. As stated in Sect. 3, there is a stochastic distribution in leaf thickness on the real leaf spring. Since the leaf thickness in Virtual.Lab model is constant for each segment of the leaf, an average value, calculated from 3D scan measurements, was assumed. This issue is further investigated in the following section of this paper, where correction factor for leaf thickness is considered as an input parameter during optimization procedure.
Despite those small errors in energy dissipation, it can be assumed that developed model is correctly reproducing leaf spring characteristic during dynamic operation under axial loading.
Achieved results for static simulation provide insufficient level of correlation with experimental data. The aim of the analysis presented in further chapters is the minimization of errors between simulation and experimental data. Improvement in model performance could be achieved by means of single-or multiobjective optimization, where parameters of the model would be optimized to better reflect experimental test results. The optimization will allow tuning up the two most probable sources of errors: leaves thickness and friction coefficient in contact elements.

Numerical model optimization procedure
Initial results of dynamic simulations show good correlation with conducted experiments. However, static simulations have shown significant differences between behaviour of the numerical models and experimental test rigs. Although a possible explanation has been given in previous sections of this paper, an attempt to improve the model's performance has been undertaken. To achieve this goal, a model updating procedure was carried out using Noesis Optimus optimization software. The model updating procedure (Fig. 13) for created leaf spring models was based on iterative parameter tuning by means of single-objective optimization. During each iteration, the experimental data were compared with results of the simulation to evaluate model performance by Optimus. Then, a new set of parameters, based

Optimization procedure in Noesis Optimus
The optimization procedure described in previous subsection was used to develop a workflow (Fig. 14) in Noesis Optimus that allowed for automatic model parameters selection, updating and saving of the model assembly, calling LMS Virtual.Lab Motion solver automatically for model calculation, results extraction from LMS Virtual.Lab Motion and post-processing in Mathwork's MATLAB.
The new model parameters calculated in Noesis Optimus were automatically introduced into LMS Virtual.Lab Motion model by CATIA macro written in VBScript. Subsequently, the macro was executed through system command line using LMS Virtual.Lab Motion built-in automation scripts. This allowed to avoid unnecessary opening of the GUI (Graphical User Interface), which would slow down the whole model updating process. When macro has finished changing the parameters of the model, an input file was generated and then passed to LMS Virtual.Lab Motion solver in a batch run. Solver was called again by automation scripts, so computation could be done in background without the GUI. Afterwards, results of the calculation were extracted by another VBScript macro and passed to MATLAB for post-processing. Finally, results of MAT-LAB's calculations were passed to Optimus, so they could be compared with experimental data. At this stage, objective function is evaluated to assess model performance, and new values of parameters for subsequent iteration are calculated.

Design variables
As it was mentioned in Chapter 4, design variables for optimization procedure were chosen based on model behaviour observed during initial static simulations. Achieved results for both models have shown vast differences in estimation of energy dissipated per cycle and small errors in stiffness in comparison with experimental data. Since inter-leaf friction is predominant source of energy dissipation in leaf springs, friction parameters in contact element need to be adjusted. Friction coefficient in Sphere-to-Extruded-Surface Contact elements was therefore selected as a first design variable. As mentioned in previous chapters, force, generated by model, is dependent on sectional properties of individual leaves. To alter spring stiffness, a change in individual leaf thickness can be made. Usually when model is based on particular object, such correction could make the model invalid, as it would not represent the real life counterpart of the model. However, during modelling phase, a stochastic distribution in thickness on entire leaf length was observed on 3D scans used to recreate the spring in multibody environment, so problems regarding stiffness evaluation in proposed leaf spring model may be solved by small changes in the model geometry, but in very limited range so it could still represent the actual spring that it was based on. Since an average value for each leaf was assumed in modelling phase, there is an uncertainty in regard of its correct value. Therefore, a correction factor of leaf thickness was assumed as a second design variable. To keep design variables amount small, the same correction factor was applied to all leaves. Furthermore, thickness correction was kept small (0.2 mm, which corresponds to around 3 % of the thickness) so the model would not deviate too much from its real counterpart.

Objective function and optimization algorithm
Optimization procedure involved minimizing the objective function described by expression 9 where n, number of output variables; t i , target value for ith output; tw i , target weight for ith output; R i , ith output value during iteration; S i , target scale, s i = |t i | + 0.001. Output variables selected for optimization procedure was maximum force F max registered during loading cycle and amount of energy dissipated per cycle L R . To speed up convergence rate, L R value was rounded up to second decimal place. Assumed target values for optimization procedure are listed in Table 4.
To minimize objective function 9, a global optimization algorithm, called Differential Evolution [9], was used. During the optimization procedure, algorithm is randomly generating a population of design vectors of size NP. Each vectorx (g) i , with size equal to the number of design variables, is a separate individual of generation g. For each individual, a fitness value is calculated based on objective function. In the next step, the algorithm is calculating a new vectorν i in a so-called mutation procedure. This is done by adding the weighted difference between two vectors from previous population to a third according to expression 10 [7]: where i, ith individual, i = 1 . . . N P; g, generation number; F, weight factor;x (g) Then, crossover betweenx (g) i (called parent) and mutantν i is introduced in order to create a trial vectort i . Each t i,n (n = 1, . . . , design variable amount) element of trial vector is selected according to formula 11 [7]: where CR, crossover rate; D, vector dimension; rand(0, 1), random number [0,1). L R = 8.9 (J) 0.6 L R = 6.05 (J) 0.6 Afterwards, an individual is selected between parentx (g) i and trial vectort i based on objective function value according to expression 12. Selected individual will become a part of the new generation [7] x g+1 i After the new generation is set, convergence condition are checked. If convergence is not met, then mutation, crossover and selection is repeated. Algorithm parameters values assumed during optimization procedure are given in Table 5.
During optimization procedure, 12 different simulations were conducted during single generation (algorithm iteration). Since total number of possible generations is 25, therefore maximum value of model calculations equals to 300. Between simulations model was automatically rebuilt, and changes were made to geometry of 86 bodies and parameters in 142 force elements.

Optimization results
As a results of the optimization procedure, new force-deflection curves for static simulation were acquired. A direct comparison between experimental data, initial simulation results and optimized model is shown in Fig. 15. Achieved values for design variables and output parameters are given in Tables 6 and 7

Discussions
Convergence criteria for suspension model with shackle were achieved after 196 cycles, while for model without shackle after 216 cycles. In both cases, achieved results show great improvement in leaf spring characteristics estimation, with errors within fractions of percentages (Tables 6, 7). New characteristics' shapes are almost the same as those acquired during experimental tests (Fig. 15 ). Appropriate energy dissipation was achieved only for high values of friction coefficient in contact elements. However, values this high do not have physical meaning, but still they allowed for almost error free match of hysteretic properties of the spring. It should be noted that dynamic simulation results discussed in Sect. 4 of this paper have shown that energy dissipation is connected to the loading speed. For starting friction coefficient values (μ = 0.1), during dynamic operation, much better results were achieved than for very slow (over 2 min) static test. The main source of this is the friction model used in Sphere-to-Extruded-Contact which is dependent on value of transition velocity parameter. For velocities below this value, the nominal friction force is multiplied by a hyperbolic tangent function to give a smooth transition through zero velocity. In case of velocity equal or close to zero, there is no or very little friction in modelled contact because of scaling imposed by this hyperbolic tangent function [7]. Decreasing the value of transition velocity parameter improves friction evaluation at lower speeds, but has drastic influence on computation efficiency, leading to excessive model calculation times. There are two possible solutions of this problem. First one is to increase the speed of the 'static' tests to the level where the friction evaluation is more realistic. Second one is to include into Sphere-to-Extruded-Contact element speed independent part of the friction model which would be responsible for proper modelling of the friction at very low speeds. It must be noted that such loading case as in proposed static test has no counterpart during operation of the real leaf spring in the car and that during dynamic loading of proposed model, friction evaluation was at much more plausible level. Described effects have influence on the results only for very slow loading of the spring, unlikely during dynamic simulation of complete vehicle. In both cases, thickness correction was applied, lowering the overall thickness of the spring assumed during modelling stage. Correction value is directly related with the second design variable which was friction coefficient. High values of friction coefficient assumed by optimization algorithm had a direct effect on leaf spring thickness. Therefore, leaf thickness correction should take into account not only initial stiffness error, but also overcome influence of the second design variable on the leaf spring stiffness. The resulting conclusion is the same as in previous paragraph.

Conclusions
In this paper, a numerical leaf spring model was developed using multi-body approach in commercially available software. Proposed model was created using LMS Virtual.Lab/Motion environment and tested with two mounting configurations. Simulations have proved, that developed model is correctly reproducing an overall shape of the leaf spring characteristics for dynamic simulations, with better convergence achieved for leaf spring model without shackle. Initial results for static tests have shown significant error in behaviour of the numerical model and experimental test rig. Those differences were minimized through iterative model tuning performed in Noesis Optimus by means of single-objective optimization. Leaf spring characteristics acquired after optimization process shows great improvement in evaluation of stiffness and energy dissipated per cycle, with errors within fractions of percentages. Presented in the paper results proved that implemented in LMS Virtual.Lab Motion tool for leaf spring modelling properly simulates its behaviour during axial loading in real road excitation conditions where frequency is usually above 1 Hz. Moreover, it provides fast and easy method for modelling of the leaf spring suspension in the software widely used for vehicle dynamics modelling, and it allows for significant reduction of the time required for full vehicle dynamics modelling.