Nonlinear dynamic analysis of parabolic leaf springs using ANCF geometry and data acquisition

The parabolic leaf spring is widely used in modern vehicle suspension systems because it has many desirable features, such as weak interleaf friction and light weight. In this paper, the parabolic leaf spring is analyzed using automatic geometric data acquisition and the finite element (FE) absolute nodal coordinate formulation (ANCF). In order to account for manufacturing considerations in developing the virtual models, the relaxed uniform cubic B-spline is used to represent the leaf spring profile curves and geometry. Three-dimensional scanning techniques based on structured light and multiple images are explored in this study to automatically extract the leaf spring complex geometric data. A new procedure is proposed to develop the ANCF/FE mesh from the physical object. Both double-leaf uniform-thickness and double-leaf parabolic spring models are developed and analyzed using different ANCF elements. Using ANCF geometry, piecewise linearly tapered parabolic leaf spring models are constructed, accounting for the leaf pre-stress. The interleaf contact is enforced using a penalty approach and a smoothed Coulomb friction model. It is shown that the fully parameterized low-order beam and plate elements suffer from locking problems, while the thin plate element can lead to less accurate results. The use of the new strain split method (SSM) as a locking alleviation technique is also examined in this investigation. It is shown that while the current SSM implementation can be effective in solving the locking problem in the case of symmetric bending-dominant loading, it may not produce accurate results in the case of torsional loading. The comparative study performed demonstrates that the higher-order ANCF beam element is more suitable for developing the leaf spring models compared to other ANCF elements considered in this investigation. The numerical results obtained show that the friction effect in parabolic leaf springs is much weaker than that in the leaf spring with uniform-thickness.


Introduction
Leaf springs are commonly used in the suspension systems of cars, trucks, and railroad vehicles to connect the chassis or frame to the axles. They are designed with a stiffness rate sufficient to support the car body and transmit the load from the chassis to the axle, and long leaf springs can act as guide components in the longitudinal direction. As the suspension elastic elements, leaf springs serve as initial shock absorbers that reduce the impact and isolate the vehicle body from the imperfections of the road surface, thus improving the ride quality. The mechanical properties of the leaf spring can significantly influence the performance of the suspension and vehicle, including the ride comfort, the handling performance, and the vehicle stability. Therefore, it is important to study and evaluate the dynamic behavior of the leaf springs in order to improve the vehicle performance and ensure its safe operation.
The leaf spring was first applied to the horse-drawn carriage in full elliptic form, meaning two sets of leaf springs were assembled to approximate an ellipse [1]. However, the profile of leaf springs in modern vehicles is no longer elliptic. Additionally, the profile of the modern leaf spring is not a circular arc, despite the fact that leaf springs are usually designed based on a circular arc profile following the ideal uniform strength beam model. It should be noted that most modern leaf springs are not close to either the elliptic or the circular arc profile curves. Accurate description of the geometry of the leaf spring profile curve, however, is necessary in order to develop credible virtual prototyping mod-els that capture accurately the spring deformations and forces.
Because the profile curve of the leaf spring, as shown in Fig. 1, can be complex, it is advantageous to use three-dimensional geometric data acquisition techniques to extract the spring geometric information. To this end, automatic geometric data acquisition (AGDA) techniques are discussed in this paper. The three-dimensional scanning technology has been of particular interest due to its potential for rapid and accurate measurement of complex three-dimensional geometry. The technology has been applied in various engineering fields for uses such as surveying buildings [2], road surface condition evaluation [3], evaluation of production samples after manufacturing for quality control, as well as for reverse engineering applications [4]. Additionally, three-dimensional scanning could prove to be useful in industry for the development of more realistic FE meshes created from physical objects such as manufactured production samples. In doing so, reliance on geometric assumptions based on imported computer-aided design (CAD) system models or on data possibly provided by manufacturers of sub-systems or components could be reduced. Furthermore, eliminating additional steps by developing the analysis mesh directly by scanning the physical component has the potential to reduce the time required to develop complex flexible body models. Additionally, meshes created directly from production samples will provide an as-manufactured geometric data set, which does not rely on any geometric simplifications or assumptions made at the design/drawing state. This method also allows for obtaining the geometry of products which have no associated CAD files due to the age of the design, which could predate modern CAD software. Therefore, accounting for the complexity of the leaf spring geometry using an automatic procedure to extract the geometric information from the physical leaf spring can be advantageous in developing reliable virtual prototyping models. As shown in Fig. 2, several three-dimensional reconstruction techniques can be used [5]. Two methods based on the multiple images and structured light, which are classified as passive and active methods, are discussed in detail in this paper.

Background
Leaf springs have been the subject of many investigations focused on the study of their strength and durability as well as on their effect on the vehicle dynamics and stability. Because of the lack of accurate computational multi-body system (MBS) models, many researchers turned to experimental methods to study the behavior of leaf springs [6,7]. Another widely used approach is to model the leaf spring based on the beam theory, which was first proposed by Timoshenko [8] to obtain the equivalent spring stiffness coefficient. Later, the beam theory was further developed by SAE [9]. Common curvature or concentrated loading assumptions were typically used in the conventional leaf spring design process [10]. However, all the traditional meth-ods employ many assumptions aimed at simplifying the formulations and calculations, and in many cases the effect of the interleaf contact was neglected.
With the development of computational mechanics, the FE method is often used for developing leaf spring computer models, taking into account the effect of the interleaf contact. Zahavi [11] studied the effect of the contact forces on the leaf spring deformation using the FE method. Shokrieh and Rezaei [12] studied a composite leaf spring using the general-purpose FE software ANSYS with the goal of optimizing the width and thickness of the leaves. Duan et al. [13] developed an MBS dynamic model of a tandem suspension with a tapered leaf spring, which was created using the ADAMS/CHASSIS leaf spring module.

FE models
Commercial FE computer programs have been widely used in some recent leaf spring investigations. Many authors used the commercial FE software ANSYS to investigate parabolic leaf springs. The focus of some of these investigations has been on the mono-leaf parabolic leaf springs, which have no pre-stress as result of the assembly process [14][15][16]. A three-leaf parabolic spring was also the subject of an investigation by Kumar and Aggarwal [17]. However, the effect of the initial stress caused by the spring assembly was not considered. Considering the effect of the prestress can be a challenging problem, particularly when developing an approach that accurately accounts for the spring geometry. Therefore, in this paper, a double-leaf parabolic spring will be analyzed, considering the effect of the pre-stress resulting from the spring assembly process, which can be systematically represented virtually using ANCF elements based on continuum mechanics theory.

MBS models
Leaf spring models have also been developed based on the floating frame of reference (FFR) formulation [18,19]. The FFR formulation allows for using a reduced-order model by employing component mode synthesis techniques. While FFR is a widely used method in the field of MBS dynamics, existing FFR elements are not suited for accurately representing the complex leaf spring geometry. Another alternative is to use ANCF elements, which employ absolute nodal coordinates, for modeling leaf springs [18]. Instead of using the FFR infinitesimal nodal rotations as nodal coordinates, absolute nodal slopes are used to represent the rotation and deformation of the body when ANCF elements are used [20]. Additionally, ANCF elements have desirable features including the constant mass matrix and zero centrifugal and Coriolis inertia forces. Yu et al. [21] developed a leaf spring model using ANCF fully parameterized plate elements and the concept of the ANCF reference node. A doubleleaf uniform-thickness spring model was developed to demonstrate that ANCF elements can be systematically used to predict the spring nonlinear dynamic behavior.

Contributions of this study
This paper is focused on developing new FE models for the nonlinear dynamics and vibration analysis of parabolic leaf springs. The main contributions of this paper are: (1) demonstrate how automatic geometric data acquisition (AGDA) can be integrated with ANCF elements to develop accurate representation of the spring geometry; (2) demonstrate how the geometry mesh can be used as the analysis mesh without the need for any conversion that may lead to geometry distortion, thereby demonstrating the use of the integration of computer-aided design and analysis (I-CAD-A) in the simulation of leaf springs; (3) develop a procedure that accounts for the complex geometry and the effect of the pre-stresses resulting from the assembly of the parabolic leaf spring; (4) conduct a comparative study based on different ANCF leaf spring models in order to evaluate the performance of different ANCF elements in modeling leaf springs; (5) demonstrate the use of the proposed approach by simulating the scenario of axle articulation when the vehicle passes a unilateral bump; and (6) examine the use of the new strain split method (SSM) as a locking alleviation technique for leaf spring applications [22].

Organization of the paper
This paper is organized as follows: In Sect. 3, a brief description of the four ANCF elements used in this investigation to model the leaf spring is presented and their ability to capture complex geometries is demonstrated. In Sect. 4, two methods to extract geometric data automatically, the procedure to obtain the FE model from the extracted geometric data, and an alternative method to avoid scanning the assembled configuration are discussed. The discussion in Sect. 4 is necessary in order to explain how to develop the ANCF/FE geometry/analysis mesh from physical objects based on three-dimensional scanning techniques. Section 5 discusses the use of ANCF elements to represent the leaf spring geometry including the configuration of the spring, the profile curve representation, and the representation of the parabolic leaf spring geometry and prestress. In Sect. 6, the use of the ANCF reference node (ANCF-RN) to connect the leaf spring to the vehicle components and to model rigid elements is explained, and the fundamental differences between the ANCF-RN and the rigid body elements (RBE) used in the conventional FE literature are explained. In Sect. 7, the method used in this investigation for the simulation of the interleaf contact is described. In Sect. 8, the recently introduced strain split method (SSM) proposed as a locking alleviation technique is discussed. In Sect. 9, uniform-thickness, double-leaf spring models developed using different ANCF elements are compared. In Sect. 10, a parabolic leaf spring model is developed using the ANCF higher-order beam element in order to alleviate the locking problem. The effects of prestress and interleaf friction are shown, and the SSM performance in alleviating locking problems is evaluated using parabolic leaf spring models. In Sect. 11, conclusions drawn from the investigation are provided.

ANCF element geometry
Four ANCF elements including the three-dimensional low-order beam element with 24 nodal coordinates (LOBE24) [23], the high-order beam element with 42 nodal coordinates (HOBE42) [24,25], the thick plate element [26], and the thin plate element [27] are used in this paper to develop the leaf spring models. Two nodes and four nodes are employed by the beam and plate elements, respectively. The four elements can be classified into fully parameterized and gradient-deficient elements; the LOBE24, the HOBE42 and the thick plate element are fully parameterized elements. All the fully parameterized elements utilize all three gradient vectors as nodal coordinates and the corresponding elastic forces are derived based on the general continuum mechanics approach. In addition to the position and gradient vectors, three curvature vectors are used as vectors of nodal coordinates in the HOBE42. By contrast, the elastic forces of thin plate element, which is a gradient-deficient element since there is no gradient vector in the thickness direction, are developed based on the Kirchhoff plate theory that does not account for the transverse shear deformation [27]. Both LOBE24 and HOBE42 are based on cubic interpolation in the longitudinal parameter. Similar to the beam element, the interpolations of the plate elements in the two directions of the mid-surface are also cubic. While the interpolation in the thickness direction for the thick plate element is linear, the element has a full set of gradient vectors at the nodal points. For the thin plate element, no gradient vector is used in the thickness direction. It is important to note that in the case of ANCF elements, shell structures can be modeled using plate elements by an appropriate choice of gradient vector coordinates and that fully parameterized ANCF beam and plate elements retain the isoparametric property.

Low-order beam element (LOBE24)
This element has two nodes and each node has 12 coordinates that define the nodal position vector r and three position vector gradient vectors r x = ∂r/∂ x, r y = ∂r/∂ y, and r z = ∂r/∂z, where x, y, and z are the element parameters. The shape function matrix of LOBE24 used in this investigation is [23] where the shape functions s 1 k , k = 1, 2, . . . , 8, are defined as In these equations, l is the length of the element, and ξ = x/l, η = y/l, and ζ = z/l.

High-order beam element (HOBE42)
This high-order beam element has 42 nodal coordinates. In addition to the position vector and the full set of gradient vectors, three curvature vectors r yz = ∂ 2 r/∂ y∂z, r yy = ∂ 2 r/∂ y 2 , and r zz = ∂ 2 r/∂z 2 are used as nodal coordinate vectors. The HOBE42 shape function matrix is defined as [24,25] The shape functions s 2 k , k = 1, 2, . . . , 14, are defined as The dimensionless parameters that appear in this equation are the same as the ones used for the LOBE24 element.

Thick plate element
The ANCF thick plate element has four nodes, each of which has 12 coordinates that represent the nodal position vector r, and three position vector gradient vectors r x , r y , and r z , where x, y, and z are the element parameters. The element, therefore, has 48 coordinates.
The shape function matrix of this element is defined as [26] S 3 = s 3 1 I s 3 2 I . . . s 3 16 I The thick plate element shape functions s 3 k , k = 1, 2, . . . , 16, are defined as In these equations, a, b, and t are, respectively, the length, width, and thickness of the element, and ξ = x/a, η = y/b, and ζ = z/t.

Thin plate element
The ANCF thin plate element is a gradient-deficient element since only two parameters are used in the element interpolation. The element has four nodes, each of which has 9 coordinates that represent the nodal position vector r, and two position vector gradient vectors r x and r y , where x and y are the element parameters. The shape function matrix of thin plate element is defined as [27] S 4 = s 4 1 I s 4 2 I . . . s 4 12 I The plate element shape functions s 4 k , k = 1, 2, . . . , 12, are The element dimensions and dimensionless parameters that appear in these shape functions are the same as the ones defined for the thick plate element.

Element geometry
The ANCF elements described in this section will be used to develop leaf spring models in this investigation. Table 1 shows a comparison between these elements, which can be used to describe curved and tapered geometry that characterizes leaf spring designs. Table 2 shows some of the possible shapes that can be obtained using these elements. In this table, only one element is used to produce the shape based on the element nodal coordinates reported in the table. By using position vector gradients instead of rotations, no assumptions are made with regard to the magnitude of the rotations within the element. This allows complex shapes to be represented using a small number of ANCF elements.
To explicitly show the ability of ANCF elements to describe complex geometry, the single element model is used in the table to describe significantly curved and tapered configurations using the four types of the ANCF elements employed in this investigation.

Physical leaf spring geometric data
As previously discussed in this paper, the geometry of the parabolic leaf spring is complex and difficult to measure. However, because of the importance of accurate geometric representation for credible virtual durability analysis, the development of an efficient auto- Vector of Nodal Coordinates r, r x , r y , r z r, r x , r y , r z , r yz , r yy , r zz r, r x , r y , r z r, r x , r y matic procedure for extracting the geometric information from the physical leaf spring is necessary. The approach can be generalized to create the geometry and the FE analysis mesh from the physical product. In this section, methods that can be used to reconstruct the three-dimensional CAD/FE meshes from the physical leaf spring are discussed. The main steps used in these methods are schematically shown in Fig. 3.

Geometric data acquisition using structured light scanning
A simple and widely used type of three-dimensional scanning techniques is the structured light scanning. Structured light scanners project a two-dimensional banded pattern, such as narrow linear stripes, on the three-dimensional target surface in a specific frequency of light [28]. One or more cameras are then used to view the surface at an oblique angle. A bandpass filter can be used in the camera to only detect the frequency of light projected on the surface of the object [29]. When the pattern is projected on an irregular surface and viewed from the oblique angle, the linear banded pattern becomes distorted. The distortion of the projected pattern present in the viewing angle of the camera due to the contours of the surface is then used to determine the variation in surface geometry [30]. While a single scan using a few two-dimensional images can map the surface of a portion of an object, the means of determining the geometry relies on lineof-sight and two-dimensional images to acquire the vol-ume of the object. Thus, to obtain the surface geometry of a three-dimensional object, a number of these scans must be taken and integrated to form a data set which can describe the surface and volume of an object. It is important to note that the grid of data points for a given surface geometry is distributed uniformly on a plane normal to the direction of the scanner, so oblique angles will have lower point cloud density than surfaces normal to the scan direction. Thus, deep wells or features such as bolt holes, or in the case considered in this paper, the gap between two leaves of a leaf spring assembly will not be captured accurately.
In order to develop the full surface geometry, the data sets from each scan must be processed to form a single set. This data set includes a number of points obtained from the grid of surface points obtained during the scanning process and is referred to as a point cloud. The process of aligning three-dimensional point clouds into a single complete data set is referred to as registration [31]. Algorithms designed for integrating the multiple scans into a single set of points are often based on a pairwise registration scheme. In this scheme, pairs of point clouds are condensed into a single point cloud based on distinguishing geometric features common to both point clouds. For this method to work properly, distinguishing geometric features common to both point clouds are required to allow the registration to be completed [30].
A point in the cloud, referred to as a key point, is first selected on each of these geometric features. It is important that these points are located on a feature containing sufficient point density and uniqueness; that    is, it should not be part of any repeated geometry that would be indistinguishable from other features in the data set. Figure 4 shows an example of points on unique features of the leaf spring that could be used as key points. Once the key points have been estimated, they are used to define a larger set of points called a feature descriptor. The feature descriptor is composed of the key points as well as several adjacent points. Local position vectors can be used to define the location of the adjacent points in the cloud with respect to the key point. These sets of points and position vectors are used to determine the proper position and orientation of one point cloud with respect to another such that they can be combined into a single data set [29]. It is important to note that there is a certain amount of noise present in the data that will prevent exact matching of feature descriptors in the pair of data sets. Thus, the assembly of the data sets into a single cloud requires estimates to be made to determine potential matches between the point clouds. This step is referred to as the correspondence estimate. In this step, the similarity and position of the given feature are compared. If the correspondence is found to not be a good match, then the correspondence estimate is rejected. Once a correspondence estimate is successfully found, a transformation is determined to translate and rotate the bodies such that the common feature on each point cloud is coincident [29]. Figure 5 shows the process of registration. It is important to note that geometries which are smooth or have only repetitive patterns will have many possible ways in which the registration algorithm will attempt to assemble the two grids into a single point cloud, which can lead to improperly combined point clouds.

Geometric data acquisition using multiple images
The geometric data acquisition based on the multiple images technique is another important threedimensional reconstruction approach. The required tool, the camera, is relatively inexpensive and portable when compared with the device used for structured light scanning. Additionally, many modern cameras are high-resolution with several megapixels. The procedure to convert multi-view two-dimensional images to three-dimensional geometric data points consists of four steps: camera calibration, depth determination, registration, and material application. To accurately determine the depth, camera calibration is usually required [32]. This is followed by the determination of the depth, which is the three-dimensional component missing from the planar image. The registration that follows is performed to combine the obtained multiple depth maps to create a final mesh. As an optional step, the material application may need to be implemented to obtain the texture and color of the reconstructed object. The steps of the three-dimensional reconstruction approach based on the multiple images used in this study are outlined in Fig. 3. The reconstructed threedimensional mesh and the corresponding point cloud are shown in Fig. 3b, c.

FE mesh based on the acquired geometric data
Using either of the methods described above, a point cloud that describes a grid of points on the scanned surface can be obtained. Many problems can occur during the process of converting the point cloud to a surface, such as an inherent amount of noise present in the scanning process as a result of variation in the depth reading from the scanning device as well as the holes and missing data. Additionally, it should be noted that the point cloud describing the surface can have a large data set, which can be computationally expensive to process. Therefore, rather than processing the entire point cloud and converting the cloud to a surface with an existing algorithm, only some of the points are selected and used to interpolate the profile curve of the leaf spring based on the representation method of the spring profile curve introduced in Sect. 5.2. The procedure to obtain the FE mesh from acquired geometric data is also shown in Fig. 3.
In this paper, the multiple-image-based reconstruction technique is used to acquire the three-dimensional geometric data from a set of two-dimensional photos taken at different angles. A Nikon D200 digital single lens reflex camera is used, and the leaf spring surface is marked slightly with chalk to create black and white points, in order to produce a high-quality mesh from scanned surface. The obtained meshed surface is shown in Fig. 3b. Following the scanning, one can obtain the three-dimensional model and extract the point cloud, from which appropriate sample points are selected to interpolate the profile curve of the surface shown in Fig. 3c. However, to develop the FE mesh with beam or plate elements, the points on the mid-surface rather than the surface points are used. To determine the points on the mid-surface from the points on the two surfaces, the plane cross-section assumption is used. In other words, the cross section is assumed to remain plane after deformation. Because the dimension of the crosssection is uniform and the initial curvature is small, the mid-surface is assumed to coincide with the geometric central plane. Thus, the points in the middle of the fiber lines connecting two corresponding points on the upper and lower surfaces, which share the same material coordinates in the two directions of the mid-surface, can be assumed to be the points on the mid-surface. For convenience, points which are evenly distributed along the length of the leaf on the upper and lower surfaces are chosen as point pairs. In this implementation, the surface curve of the leaf is interpolated based on the relaxed uniform cubic B-spline, which will be introduced in the following section, using the selected points on the surface. Subsequently, the equally spaced points along the surfaces are determined and the corresponding point pairs are obtained, as shown in Fig. 3d. Then, according to the plane cross-section assumption, the points on the mid-surface are obtained at the centers of the fiber lines connecting the point pairs. Based on the obtained points on the mid-surface, the mid-surface curve is then interpolated to determine the gradients. The length of the fiber line connecting the two points on the surfaces is used as the thickness of the leaf spring at the point on the mid-surface. The interpolated curves including the surface and mid-surface curves are shown in Fig. 3e, and the final ANCF FE model is shown in Fig. 3f.

Pre-stressed reference configuration
Both the stress-free and pre-stressed assembled reference configurations are required in order to determine the pre-stress in the leaf spring. Acquiring the geometric data of the assembled leaf spring is much more difficult and time-consuming as compared to acquiring the data of the separate leaves due to the small varying gaps between the two leaves when they are assembled. These gaps lead to the problems of inaccurate point cloud meshes as previously mentioned. Instead of scanning the assembled leaf spring using the existing scanning technology, two alternate approaches can be used. The first is manual measurements of points on the assembled leaf spring, and the second is an iterative numerical procedure. Both methods have shown to produce similar results in this investigation.
The iterative numerical procedure can be used to obtain the pre-stressed reference configuration (assem- Fig. 6 Quasi-static simulation of the leaf spring assembly process. a Stress-free configuration. b Initial configuration e 0 before inner loop of static analysis. c Final equilibrium configuration e i after inner loop of static analysis bled spring) shown in Fig. 6c from the stress-free configuration of the leaves before assembly, shown in Fig. 6a, by simulating the quasi-static process of assembly numerically. As shown in Fig. 6b, the initial leaf spring assembly configuration is acquired by clamping the two leaves together at the center clamped segment. The initial penalty δ 0 at the two end contact elements can be obtained to start the iteration to determine the final equilibrium leaf spring assembly configuration as shown in Fig. 6c. The normal contact force F n j,i at contact point i of element j can be expressed in terms of the penetration δ j,i between the contact surfaces as F n j,i = K 0 δ j,i n j,i , where K 0 is the current penalty contact stiffness coefficient and n j,i is the corresponding normal vector. As shown in the flowchart of Fig. 7, the iterative algorithm consists of two loops: the outer and inner loops. The purpose of the outer loop is to adjust the current contact stiffness coefficient K 0 according to the obtained equilibrium configuration e i to ensure convergence to the desired value of the penetration δ tol . The inner loop is used to calculate the static solution under the current unbalanced force F u = F e − Q e , where F e = n e j=1 p i=1 S T F n j,i n j,i is the generalized force associated with the normal contact force calculated using the current stiffness coefficient K 0 and configuration e 0 , p is the number of contact points, n e is the number of contact elements, S is the associated shape function matrix, and Q e is the elastic force calculated using the current configuration e 0 . Convergence is assumed when the maximum penetration is smaller than a certain user-specified tolerance δ tol , that is, 0 < max (δ) < δ tol , and the obtained equilibrium configuration e i is considered as the final assembly configuration.

ANCF leaf spring geometry
While the leaf spring profile curve is often assumed a circular arc in the conventional design process, the actual profile curve is quite different from either a circular arc or an elliptic arc, as previously mentioned. Therefore, a simple assumed functional description does not capture correctly the spring profile geometry. In this section, the representation of the leaf spring geometry using ANCF elements is discussed. Based on the geometric features of the leaf spring, the relaxed uniform cubic B-spline is proposed to represent the spring profile curve using the scanning point cloud. A discussion on the conventional piecewise-tapered modeling approach used for parabolic leaf springs is also provided in this section. It is shown that, by using Fig. 7 Flowchart for analysis of quasi-static assembly process ANCF elements, the tapered geometry of the leaf spring can be systematically and accurately described.

Leaf spring configuration
The general leaf spring configuration is shown in Fig. 1. The leaves are clamped together by the center bolt, and the leaf spring assembly is fixed to the spring seat by a set of U-bolts, while the spring seat is rigidly attached to the axle. The leaf spring is connected to the chassis using the shackle at one end and a pin joint at the other end. The shackle can have a small rotation to allow for small longitudinal translation of the spring with respect to the chassis. Leaf springs are designed to have pre-stress in order to reduce the stress of the master leaf when the spring is mounted on the vehicle. This pre-stress can be achieved by changing the leaf curvature during the assembly process, as shown in Fig. 8a. Specifically, nips between the adjacent leaves, below the master leaf, are formed by successively reducing the radius of curvature of the leaves, as shown in Fig. 8a  [1]. Therefore, when the leaves are clamped together, as shown in Fig. 8b, the master leaf is subjected to a bending preload opposite in direction to that caused by the vehicle load, resulting in overall reduction of the spring bending stress during the vehicle operation.
To extend the life cycle of the leaf springs, they are usually designed based on the concept of the uniform strength beam. Depending on if the thickness or width of the leaf is kept constant in the ideal model, the designed leaf spring can be classified as a uniformthickness leaf spring or a parabolic leaf spring, respectively. The uniform-thickness leaf spring is shown in Fig. 9a. In the multi-leaf spring, the leaves with different curvatures are clamped together and the adjacent leaves contact along the full longitudinal length due to the assembly pre-stress. By contrast, the other method is to vary the thickness of the leaves, rather than the width, with respect to the length of the leaf in order to obtain uniform stress throughout the leaf spring. This spring type is called the parabolic leaf spring and is shown in Fig. 9b. In practice, the thickness of the center clamped segment, which is rigidly connected to the spring seat, is constant. Unlike the uniform-thickness multi-leaf spring, the parabolic leaf spring normally has only one or two leaves because even one single parabolic leaf can function as a uniform strength beam. The weight of the parabolic leaf spring is therefore less than that of the uniform-thickness multi-leaf spring for Fig. 8 The leaf spring before and after assembly. a Leaves before assembly. b Assembled leaf spring Fig. 9 Two types of leaf spring based on two uniform strength beam concepts. a Uniform-thickness multi-leaf spring. b Parabolic leaf spring the same stiffness. In practice, to increase the stiffness of the leaf spring, parabolic leaf springs with two leaves are commonly used, and a rubber or plastic spacer is often inserted between adjacent leaves at the ends to separate them and reduce the interleaf friction. As shown in Fig. 9b, there are gaps between leaves in the assembly and the leaves initially contact only at the two edges, unlike the uniform-thickness multi-leaf spring in which the leaves can have distributed contact along the length of the leaves. As a result, the friction between the parabolic spring leaves is much less than that in the uniform-thickness leaf spring, improving stability and ride comfort. Because the parabolic leaf spring has advantages of light weight and weak interleaf friction, it is widely used in modern suspension systems.

Relaxed uniform cubic B-spline
In previous publications, a polynomial has been used to interpolate the leaf spring profile curve. A piecewise cubic B-spline curve with a continuous second derivative is used in this investigation to describe the leaf spring profile curve. In order to represent accurately the spring geometry, the steps of the manufacturing and assembly process should be considered when constructing the profile curve. The leaf spring is formed by applying a distributed force on a heated semi-finished product along the longitudinal direction gradually with a finger press mold with no concentrated moments at the end points, and therefore, the curvature at the edges is zero. Consequently, the relaxed uniform cubic B-spline, which ensures zero second derivatives at the endpoints, can be used to represent the profile of the spring. The relaxed uniform cubic B-spline is constructed by connecting Bezier curves C i , i = 1, 2, . . . , n with C 2 continuity at the internal breakpoints and zero second derivatives at the first and last breakpoints [33]. Therefore, the relaxed uniform cubic B-spline P (u) for 0 ≤ u ≤ n can be defined using the given data points D i , i = 0, . . . , n. Based on the properties of the relaxed uniform cubic B-spline, the Bspline control points B i , i = 0, . . . , n can be obtained by solving the following system of linear equations: Equivalently, the equations can be written in a more general matrix form for the convenience of calculation as The B-spline is constructed with n Bezier curves with the ith Bezier curve defined by the four control points P 0 = B i−1 , With the four control points, the ith Bezier curve can be obtained as P (u) where u in this equation is the parameter of the Bezier curve, and 0 ≤ u ≤ 1 . Using the relaxed B-spline representation, the leaf spring profile curve can be systematically constructed using the scanned data points obtained using one of the procedures described in the preceding section.

ANCF parabolic leaf spring representation
The development of the geometry/analysis mesh of the parabolic leaf spring can be challenging due to its nonuniform-thickness, which is function of the square root of the longitudinal distance (arc length). While ANCF elements can be used to represent nonlinear geometry, linearly tapered elements are used in this investigation in order to provide simple and clear demonstration of how the ANCF position vector gradients can be used to obtain the desired shapes. Stretch and shrinkage of the element cross section can be conveniently realized by changing the norm of the gradient vectors [34,35]. For example, tapering an ANCF element in the z direction can be achieved by simply multiplying the gradient vector r z by scaling factors α and β at the first and second element node, respectively, as shown in Fig. 10.
For an arbitrary point with dimensionless coordinates in the longitudinal and lateral directions ξ and η, the position vectors of the corresponding points on the top and bottom surfaces defined, respectively, by r t , r b , and the cross-sectional thickness h ξ are where h 0 is the nominal thickness before altering the element shape, and ζ t and ζ b refer to the nondimensional thickness parameter whose definition varies with the type of element. It is clear that the thickness of the element changes linearly with respect to the length; thus, the parabolic leaf spring can be approximated with linearly tapered ANCF elements. Table 2 shows examples of geometries that can be obtained using a single element of the types of elements used in this study.

ANCF modeling of the spring pre-stress
In order to develop accurate representation of the leaf spring pre-stress, the three different configurations shown in Fig. 11 must be considered. These configurations are the straight configuration, the initial stressfree configuration, and the assembled pre-stress configuration. The initial stress-free configuration is a curved configuration as shown in Fig. 11. A line element dx in the straight configuration corresponds to a line element dr i in the initial configuration and to a line element dr in the current configuration. The relationship between dr and dr i is defined using the matrix of position vec- tor gradients J, the relationship between dr and dx is defined using the matrix of position vector gradients J e , and the relationship between dr i and dx is given by the matrix J i . One has the following relationships: The Green-Lagrange strain tensor is defined in terms of the position vector gradient matrices as where J e = ∂ (Se) /∂x, S is the shape function matrix, and e is the vector of nodal coordinates of the leaf spring in the current configuration. Two initial configurations of the leaf spring can be used depending on whether or not the pre-stress is considered in the analysis, which are the configurations before and after assembly, as shown in Fig. 8a, b. Denoting the vector of nodal coordinates of the leaf spring before and after the assembly as e u and e 0 , respectively, the corresponding matrices of position vector gradients are denoted as J u and J 0 . If the leaf spring assembly is used as the initially stress-free configuration, shown as the dashed line in Fig. 11, there will be no initial stresses. In order to account for the pre-stress caused by the assembly process, the configuration before assembly shown in Fig. 8a should be used as the initial stress-free configuration, the mapping relation is shown as the solid

Connection to vehicle components
As shown in Fig. 1, the leaf spring is connected to the chassis at one end using the shackle and at the other end using a pin joint. It is also rigidly connected to the axle at the spring seat by a set of U-bolts at the center section. In this investigation, it is assumed that the chassis, shackle, and leaf spring eye can be modeled as rigid bodies since they are relatively stiff compared to the flexible leaf spring. The chassis, axles, and leaf springs can all be considered as one ANCF/FE mesh in which the constraint equations and connectivity conditions, including the pin joint constraints, can be formulated using linear algebraic equations. This allows for eliminating the dependent variables at a pre-processing stage, thereby significantly reducing the number of coordinates and constraint forces that need to be determined during the dynamic simulation.
The development of such an ANCF geometry/ analysis mesh can be achieved using the concept of the ANCF reference node (ANCF-RN) which allows for including rigid components in the FE mesh. It is important to point out that the ANCF-RN concept is fundamentally different from the rigid body element (RBE) used in the commercial FE computer programs. The RBE approach is often used with elements that employ infinitesimal rotations as nodal coordinates and an incremental co-rotational procedure. The ANCF-RN approach does not impose any restriction on the amount of displacement, including rotation, and can be used with a non-incremental solution procedure. In the RBE approach, infinitesimal rotations are used, while in the ANCF-RN approach, gradient vectors are used. The constraint equations used in these two different approaches to formulate the rigidity conditions and other joint and connectivity conditions are fundamentally different.
In order to develop a geometry/analysis mesh that includes the chassis, axles, and leaf springs, the ANCF-RN approach is used to represent the rigid chassis, rigid axle, shackle end, and leaf spring eyes. In the case of a rigid body, six nonlinear ANCF-RN con- Fig. 12 Contact detection and distributions of contact points. a Contact detection. b Contact points inside the element surface. c Contact points at edges straint equations, |r x | = 1, r y = 1, |r z | = 1, r T x r y = 0, r T x r z = 0, and r T y r z = 0, are applied for each reference node during the dynamic simulation to ensure that the ANCF-RN position vector gradient matrix is an orthogonal matrix [36]. Since the constraints between the ANCF-RN and other leaf spring mesh nodes are linear in the ANCF coordinates, the dependent degrees of freedom associated with these linear constraint equations can be eliminated at a pre-processing stage before starting the dynamic simulation through the use of a velocity transformation matrix. The number of ANCF reference nodes can also be minimized by using the appropriate gradient constraint equations at the pre-processing stage; this will allow for the elimination of the largest number of dependent variables and constraint equations before the start of the simulation. The use of the new ANCF geometry approach is not only necessary for the development of new efficient and accurate leaf spring models, but is also necessary for the successful integration of computer-aided design and analysis (I-CAD-A) that will eventually eliminate reliance on three different software (CAD, FE, and MBS) which suffer from serious incompatibility problems.

Interleaf contact forces
In this study, the contact force between the spring leaves is formulated using a penalty approach. In this approach, the bottom surface of the upper leaf is considered as the master surface on which the contact points are pre-specified, and the upper surface of the adjacent lower leaf is considered as the slave surface. The minimal distance criterion is adopted to efficiently find the closest point on the upper surface of the adjacent lower leaf [37]. The contact detection method used in this study is shown in Fig. 12a. Let P be one of the specified points on the master surface, and Q be a potential contact point on the slave surface. After determining the nearest point Q, it is necessary to check the penetration δ by projecting the vector r P Q = r P −r Q along the unit normal vector n = r P,x × r P,y / r P,x × r P,y , where r P and r Q are the position vectors of points P and Q, r P,x and r P,y are the position gradient vectors defined by differentiating r P with respect to x and y, respectively. If δ = r Q − r P T n is smaller than a specified tolerance, points P and Q are assumed to be in contact.
It should be noted that the contact between leaves predominantly occurs at the edge of the shorter leaf [21]. This fact can be utilized in developing an efficient contact search algorithm. Depending on the type of the leaf spring, the distribution of the contact points may vary. For the uniform-thickness leaf spring, contact points on the element internal surface are employed as shown in Fig. 12b, while for the parabolic leaf springs, edge contact points are employed as shown in Fig. 12c. The normal force applied to the contact points based on a penalty approach is defined as F n = K δ + cδ n, where K is the penalty stiffness coefficient, c is the damping coefficient, δ is the penetration,δ is the contact point penetration rate, and n is the normal vector to the contact surface. Using the normal force vector, the tangential friction force F f can be calculated. In order to avoid discontinuity and improve the computational efficiency of the model, the following smoothed tangential friction force model is used [38]: where μ is the friction coefficient, t r is the unit vector along the relative tangential velocity v tr , and v tr = v r − (v r · n) n, v r is the relative velocity, v s is an assumed slip velocity, and γ = |v tr | /v s .

MBS equations and locking alleviation
The ANCF geometry/analysis meshes of the leaf spring, developed using different ANCF elements, can be simulated using general computational MBS algorithms designed for solving a system of differential/algebraic equations (DAE's). The MBS equations of motion can be derived using the principle of virtual work in dynamics. The number of constraint equations can be significantly reduced by using the ANCF-RN approach previously discussed in this paper. This approach allows for eliminating the system dependent variables at a pre-processing stage. The remaining nonlinear constraint equations, which cannot be eliminated at a pre-processing stage, are combined with the system differential equations using the technique of Lagrange multipliers. The strain split method (SSM) that was recently proposed for the ANCF fully parameterized beam and plate/shell (structural) elements is applied to the leaf spring example considered in this paper [22]. The basic idea of the strain split method is to use two different constitutive models with an additive split of the Green-Lagrange strain tensor, such that the high-order strain terms that arise from the ANCF coupled displacement interpolation functions are decoupled from the lowerorder terms. In case of the low-order ANCF beam element [23], the position field and the gradient vectors can be written as r = r c + yr y + zr z r x = r c x + yr yx + zr zx r y = r y , r z = r z (15) and the matrix of position vector gradients can be written and additively split into two parts as J = r c x + yr yx + zr zx r y r z = J c + J k = r c x r y r z + yr yx + zr zx 0 0 (16) where J c includes the terms associated with the beam centerline, and J k includes any remaining higher-order terms. The Green-Lagrange strain tensor can be written similarly into two parts as ε = ε c + ε k such that Using the two terms of the additive Green-Lagrange strain split, the second Piola-Kirchhoff stress can be written in Voigt form as and E is the Young's modulus, λ and μ are the Lamé's parameters, and k s2 and k s3 are the shear correction factors. Using the definition of the second Piola-Kirchhoff stress shown above and the Green-Lagrange strain, the virtual work of the elastic forces can be written in the general continuum mechanics approach as where σ is the second Piola-Kirchhoff stress tensor. The SSM details for handling initially curved geometry can be found in the work of Patel and Shabana [22]. The SSM implementation details discussed above are based on the procedure that follows the general continuum mechanics approach. An alternate implementation of the strain split method is based on handling the two stress matrices that arise from the strain split separately while formulating the vector of elastic forces. Starting from a potential function defined as the elastic forces can be formulated as Q k = ∂U /∂e, where e is the vector of nodal coordinates. Another implementation that is based on the continuum mechanics approach and a slightly modified constitutive interpretation is that of the element local frame. It was noted that in the case of structures with a small radius of curvature, SSM leads to slight overprediction of the element deformation. A local frame-based interpretation of the constitutive model yields accurate results in the case of structures with large curvature. Since the transformation of the constitutive coefficients is rather cumbersome, the strains can be transformed into the element local frame, the local stress can be evaluated using the SSM constitutive model and the global stress can be evaluated by using the element local frame to global frame transformation applied to the local stress matrix [39]. In this paper, since the curvature of the leaf spring is quite small, the deviation of the SSM solution is quite small. Accordingly, the SSM implementation for the leaf spring analysis does not include the local frame transformation.

FE comparative study
Before developing new parabolic leaf spring models, it is important to examine the behavior of different ANCF elements discussed in this investigation. To this end, a simple transient dynamic analysis, in which the leaf spring is subjected to the sudden application of a constant vertical force without damping, is performed. In order to focus on the behavior of the element, the pre-stress is neglected in this preliminary analysis. Furthermore, because constructing tapered leaf springs using the thin plate element is not straightforward, a uniform-thickness double-leaf spring with the same profile obtained from scanning the parabolic leaf spring, shown in Fig. 13b, is used in this prelimi- nary FE comparative study instead of the parabolic leaf spring. It should be noted that the tested leaf spring is not strictly a standard laminated leaf spring because the lengths of the two leaves are the same and laminated leaf springs generally include more than two leaves of varying length. The uniform-thickness double-leaf spring is employed for the convenience of comparison with the corresponding parabolic leaf spring.
Both beam and plate elements can be used to model the leaf springs [18], and therefore, two different ANCF beam elements and two different ANCF plate elements are tested in this numerical investigation. These elements are the thin plate element [27], the thick plate element [26], the three-dimensional low-order beam element with 24 nodal coordinates (LOBE24) [23], and the three-dimensional high-order beam element with 42 nodal coordinates (HOBE42) [24,25], which were discussed previously in this paper.
The leaf spring geometry shown in Fig. 13b is symmetric and has two full length leaves with uniformthickness t t = 0.02667 m and width b = 0.1016 m. The spring has the following dimensions, as indicated in Fig. 13a: A = B = 0.7357 m, C = 0.09429 m, and the pack thickness D = 0.05484 m. As previously discussed in this paper, the leaf spring profile curve is interpolated using the relaxed uniform cubic B-spline based on the geometric information directly extracted from the three-dimensional scanning data.
The results presented in this section are based on a 28-element mesh of the spring assembly, in which each leaf is modeled with 14 elements, and the end elements on the two leaves are assumed to be in contact. The Young's modulus of the leaves is assumed to be E = 2.06 × 10 11 Pa, the coefficient of friction is 0.6, and a force P = 31,121 N equal to half the weight of the chassis is applied on the leaf spring-chassis reference node. The vertical displacements of the chassis reference node predicted using different element models are compared.

ANCF beam element models
Both the LOBE24 and HOBE42 elements are used to analyze the leaf spring. The results obtained for the vertical displacement of the chassis reference node are compared and shown in Fig. 14. In this figure, the simulation results are obtained for two values of the Poisson ratio, ν = 0.27 and ν = 0, in order to examine the effect of locking on the two elements. As shown in Fig. 14a, the LOBE24 results do not agree with the HOBE42 results in the case of ν = 0.27. Specifically, the LOBE24 mesh exhibits smaller vertical displacement, which implies that the LOBE24 element model is stiffer. On the other hand, the results obtained using the two elements coincide when the Poisson ratio is zero, as shown in Fig. 14b. This is attributed to the fact that the LOBE24 element suffers from the Poisson locking, which is caused by the coupling of the longitudinal and transverse beam strains due to the Poisson effect and the varying orders of interpolation used in those directions [40]. The use of the HOBE42 element, on the other hand, eliminates the Poisson locking as the result of using the quadratic interpolation in the transverse directions [24,25].

Two ANCF plate element models
The leaf spring is also modeled using both thin and thick plate elements. A comparison between the results of the vertical displacement of the chassis reference node predicted using the two element models is shown in Fig. 15, where the simulation is performed with two different values of Poisson ratio, ν = 0.27 and ν = 0, in order to examine the effect of locking on the two elements. As expected, the thick plate model results in smaller vertical displacement compared to the thin plate model in the case of ν = 0.27, as shown in Fig. 15a. The thick plate model exhibits stiffer behavior because this element also suffers from the locking problem [41]. As shown in Fig. 15b, the vertical displacements obtained by the thin plate and thick plate elements match more closely but are still slightly different in the case of ν = 0. Two effects contribute to the slight differences in this case of zero Poisson ratio. First, the stress forces of the thin plate element are formulated using the Kirchhoff plate theory, which neglects the shear deformation in the thickness direction, while a general continuum mechanics approach is used to formulate the stress forces of the thick plate element. This contributes to the difference in the results despite the fact that the ratio between the thickness and the length of the leaf spring is small and the shear effect is not significant. Second, in the case of the gradient deficient element, contact points and forces are

Comparison between the beam and plate models
The results of the vertical displacement of the chassis reference node obtained using the beam and plate element models are also compared. As the thick plate element and the LOBE24 element suffer from the locking problem, the Poisson ratio is set to zero when comparing all the elements to eliminate the effect of locking. The results are shown in Fig. 16. While the results obtained using the thick plate, the LOBE24, and the HOBE42 elements agree well, the solution obtained using the thin plate model is slightly different for the reasons previously explained. Based on these results, it can be concluded that the thick plate and the LOBE24 elements suffer from the Poisson locking problem.
Additionally, the thin plate is not suitable for modeling the leaf spring contact. Therefore, based on the comparative study presented in this section, the high-order beam element, HOBE42, which does not suffer from the Poisson locking, will be used to obtain the reference solutions for the dynamic simulation results presented in the following section.

Parabolic leaf spring simulation results
Based on the simple transient analysis of the preceding section, the parabolic leaf spring is modeled using the high-order beam element and is approximated by piecewise linearly tapered elements. The geometry extracted from the scanning data is used to develop the parabolic leaf spring model, taking into account the effect of the pre-stress. The focus in this section will be on demon-strating the advantage of using the parabolic leaf spring in reducing the interleaf friction and on examining the use of the SSM locking alleviation technique when applied to leaf spring modeling.

Interleaf friction and pre-stress
The boundary and loading conditions of the leaf spring model used in the simulation are shown in Fig. 17.
In the test performed in this section, the configuration shown in Fig. 17 is used, and the vertical load is suddenly applied on the center clamped segment of the leaf spring. As shown in the figure, the leaf spring is pinned to two horizontal sliding hinged supports at the two eyes. The procedure to account for the effect of the pre-stress is previously described in this paper. The prestress in the longitudinal and thickness directions introduced by the assembly process is shown in Fig. 18. It is clear that desirable initial pre-stress, which is opposite to the stress caused by the vehicle load, is introduced in the master leaf. To explicitly show the advantage of the parabolic leaf spring in reducing the interleaf friction, a uniform-thickness double-leaf spring model with the same profile curve was also developed. The vertical displacement of the center clamped node of the leaf spring is shown in Fig. 19. It is clear that the amplitude of vibration decreases due to the interleaf friction for both the uniform-thickness and tapered leaf springs. After 10 periods of vibration, the amplitude decreases by about 30% and 61% of the initial amplitude for the tapered and uniform-thickness leaf springs, respectively. It is clear that the effect of friction in the tapered leaf spring is smaller than that in the uniform-thickness leaf spring. This is expected since the tapered leaves contact only near the two edges, whereas the uniform-thickness leaves may contact along the entire leaf. This contributes to decreasing noise and improving ride quality.

SSM locking alleviation
Three different leaf spring models are developed in this section to check the ability of the SSM locking alleviation technique when applied to leaf spring modeling. In these models, Young's modulus is selected to be E = 2.06×10 11 Pa and the Poisson ratio is ν = 0.3; the Poisson ratio was slightly increased in order to increase the locking effect. As in the case of the previous simu-  lations, the spring geometry is obtained from the scanning data. Both the single spring and dual-spring models are considered in this section. In the case of the dual-spring model, the two leaf springs are assumed to be connected by a rigid axle. The front axle suspension full load is used as the steady load. The load for a dual-double-leaf spring suspension is 14,000 lbs, and the corresponding force is Q = 62242.5 N. A static analysis of a single mono-leaf spring is first performed to check the SSM performance in static analysis and element convergence in the case of leaf spring modeling. Due to the symmetry, the half mono-leaf spring is used for the analysis presented in this section and one end of the leaf spring is assumed to be fully clamped, as shown in Fig. 20. The tip vertical load applied to the half mono-leaf spring is P = Q/8 = 7780.3 N. Different elements/methods are used to analyze the leaf spring, and the abbreviations of all the elements/methods are shown in Table 3; the solution obtained using the commercial FE computer program ANSYS is used as reference. The Timoshenko beam BEAM188 is used in ANSYS, and the results converge at 20 elements. The obtained corresponding tip vertical displacements f are shown in Table 4. As expected, it is found that the SSM solution converges with mesh refinement. From the view of practical engineering accuracy, the solution obtained with 12 ANCF beam elements can be considered as a converged solution. Accordingly, the 12 element mesh for the half leaf will be used in the following dynamic analysis. From the comparison between the GCM, SSM, HOBE42 and ANSYS, it is clear that the SSM method can successfully alleviate the locking phenomenon and leads to a solution that is in a good agreement with the commercial FE code ANSYS.

Dynamic analysis of the mono-leaf spring
The dynamic analysis of the mono-leaf spring is also conducted. The comparison of the GCM, SSM, and HOBE42 solutions for the tip vertical displacement is shown in Fig. 21. This comparison clearly shows that the SSM solution agrees well with the HOBE42 solution for the tip vertical displacement. By contrast, the GCM solution gives smaller values because of the locking phenomenon. Accordingly, it can be concluded that the SSM locking alleviation works well even for the dynamic analysis of the initially curved and tapered structure, which is consistent with the conclusion drawn from the static analysis.

Unilateral bump negotiation
Using the ANCF-RN concept, a simplified suspension assembly is built as shown in Fig. 22. The central ele- To clearly show the orientation of the reference nodes, the gradient vectors at the reference nodes r x , r y and r z are plotted as red, green, and blue arrows, respectively. Since the road in general is not perfectly flat, a bump on the road surface usually incurs larger impact load F d than the steady load F st during stable driving. Traffic calming measures are usually installed on the road such as the speed bump to slow the vehicle, especially in the residential areas. Therefore, it is important to perform dynamic simulation of the leaf spring corresponding to the vehicle driving on a rugged road or passing a bump. For example, when only one wheel passes over a bump, the loading condition of the leaf spring can be more severe compared to both wheels of the axle simultaneously passing over the bump. Therefore, the leaf spring is analyzed when only one wheel passes over the bump. There are no national standards on speed bump size and passing speed; however, according to some street design manuals [42,43], the length of the bump is shorter than 1 foot, the height is normally 3-4 in., and 20 miles per hour is a reasonable passing speed. While an actual speed bump can be smoothed with a gradual slope, for simplicity, it is assumed to be rectangular in this numerical investigation, with dimensions shown in Fig. 23, and the passing velocity is assumed to be 8.94 m/s. To determine the value of the impact force, the simple spring-damper quarter vehicle model shown in Fig. 23 is used. According to the static or damped dynamic analysis of the doubleleaf spring, the stiffness of the spring is obtained as K = 1.030 × 10 6 N/m and the damping coefficient is set as C = 5 × 10 4 N s/m. The mass of the quartercar body is m = 3172.4 kg. The quarter vehicle model is accelerated to 8.94 m/s to pass the bump. Then, the maximum amplitude of the force applied on the leaf spring F d = K δ l + Cδ l is obtained as the impact force, where δ l andδ l are the relative vertical position and velocity, respectively, between the car body and the wheel. The vertical position of the car body versus its horizontal position is shown in Fig. 24a, and the corresponding vertical force is shown in Fig. 24b. From  Fig. 24b, it is clear that the load on the leaf spring significantly increases from the steady value to an extremum when the vehicle passes over the bump before the sudden impact is damped. The obtained maximum impact force is F d = 84,166.6 N, which is more than twice the steady load F st = Q/2 = 31,121.25 N. The vertical upward impact load F d associated with the unilateral bump is applied on the corresponding central segment of the leaf spring and the other side is applied with F st to simulate the case of axle articulation when the vehicle passes a bump on one side of the vehicle. The vertical displacement of the center node of the impacted leaf spring is shown in Fig. 25. It can be noted that the largest deformation appears at around t = 0.086 s and the corresponding configuration of the leaf spring is shown in Fig. 26a, while the axle articulation is shown in Fig. 26b. Obvious rotation of the shackle can be seen in Fig. 26c. Due to the twisting of the leaf spring, different penalties at the contact points on the edges can be seen in Fig. 26d. It can be seen from Fig. 25 that while the SSM implementation can still alleviate the locking in the case of the dynamic analysis of axle articulation and produces larger displacement than GCM implementation, it cannot give a solution close to the HOBE42 solution. This is due to the large twisting distortion of the leaf spring associated with the axle articulation loading condition. Because of the linear interpolation in the beam cross-sectional direction, the low-order beam cannot capture the nonlinear distribution of the torsional deformation. Furthermore, the SSM implementation is not currently designed to tackle the locking associated with torsional deformation. Thus, neither GCM nor SSM approach can give the correct solution in the case of torsional deformation. This can also be proved by the fact that the SSM solution is correctly predicted in the case of symmetric steady loading, where no torsion can be observed. The distribution of strain ε x x obtained from the HOBE42 simulation is shown in Fig. 27. The obvious twisting effect can be seen through the strain distribution shown in this figure.

Conclusion
Leaf springs are widely used in vehicle suspension systems to connect the chassis and axles [44]. They serve to absorb the vibration and improve the ride quality and have been the subject of many numerical as well as experimental studies [7,[10][11][12][13][14][15][16][17][18][19]21,45]. Compared with the uniform-thickness leaf spring, the parabolic leaf spring is much lighter for the same stiffness and has the advantage of reduced interleaf friction, making it less noisy and more effective in improving the ride quality. In this investigation, the parabolic leaf spring was analyzed using ANCF elements in order to properly model the geometry and pre-stress [46]. Three-dimensional scanning techniques were used to extract the complex geometric information from the physical product and ANCF/FE models were created from the data acquired from the scanned geometry. An automated procedure for developing the ANCF geometry/analysis mesh from the physical object was proposed. The relaxed uniform cubic B-spline was first used to represent the leaf spring profile curve, and the parabolic leaf spring was approximated with piecewise linearly tapered ANCF elements. To select the appropriate ANCF element type for modeling the leaf spring, two ANCF beam elements and two ANCF plate elements were considered in this investigation to develop the uniform-thickness double-leaf spring, accounting for the interleaf contact. It was shown that the loworder fully parameterized beam element and the thick plate element suffer from the locking problem. Additionally, the thin plate element can lead to less accurate results because it neglects the shear and because the contact points and forces are distributed on the element mid-surface. Therefore, the high-order beam element was chosen in this investigation to perform the dynamic simulations of the parabolic leaf spring model developed with piecewise linearly tapered elements. The numerical results obtained demonstrated that the tapered leaf spring is more effective in reducing the interleaf friction as compared to the uniformthickness leaf spring. Furthermore, the use of an SSM implementation is examined as a locking alleviation technique. The results obtained, however, showed that while the current SSM implementation can be effective in reducing Poisson locking, it may not correct the beam element kinematic deficiencies when subjected to torsional loading. Future work will include the development of high-fidelity vehicle models that include the parabolic leaf springs based on the new ANCF geometry/analysis meshes proposed in this investigation as well as further developing the SSM technique in order to improve its performance in the case of torsional loading of ANCF beam elements.