An Integrated Approach for Instability Analysis of Lattice Brake System Using Contact Pressure Sensitivity

This paper presents an attempt to evaluate the sensitivity of lattice parameters on the contact pressure instability of newly lattice brake disc. The instability characteristics are investigated through theoretical representation, modal analysis experiments, and nonlinear finite element thermo-mechanical analysis. Lattice properties are defined concerning the lattice truss angle geometry in the unit cell and periodicity of the lattice cell on the lattice plate. Motion dynamics of lattice plate concern the principal coordinates on the rotating disc are presented to use in the braking system. The modal behaviour of vanned and lattice brake disc/pad systems are defined through experimental modal analysis at free-free boundary conditions, and results are used as inputs of nonlinear finite element models as it goes through a partial simulation of the SAE J2521 drag braking noise test. Subsequently, the dynamic instability analysis of the brake disc is detailed by using the complex eigenvalue extraction technique concerning the contact pressure and lattice parameters effect. The sensitivity analysis of the brake instability respected to the mass fraction factor and relative density of the lattice structure is presented by using the average standard deviations of the contact pressure force. The likelihood of instability occurrence is quantified by definition of a single indicator derived from the system eigenvalues. The analysis indicates that the higher relative density with lower mass fraction factor of lattice structure is led to a higher temperature at the disc and pad surface. Mutually, the higher mass fraction factor with lower relative density is led to the lower temperature. The maximum contact pressure is observed in the model with less mass fraction factor and more uniform pressure distributions are observed at higher values of the mass fraction factor. The instability analysis points out that high instability frequencies are predicted at lower mass fraction factor and higher relative density.


I. INTRODUCTION
The braking system is one of the main parts of automotive vehicles with the responsibility for safety. Instability analysis of brake disc has been the focus of a lot of researches as the main subject of safety duties and comfort of riding in the braking system for many years [1]. Automotive vehicle manufacturers consider the warranty claims and subsequently, the profits loss which are related to the instability and noise The associate editor coordinating the review of this manuscript and approving it for publication was Hamid Mohammad-Sedighi . characteristic generated by the braking system whether it works correctly and in reliable condition [2], [3]. An effort to eliminate the brake noise is highlighted as it causes anxiety and discomfort for the vehicle occupants and pedestrians. Therefore a variety of experimental studies, numerical and finite element modelling techniques have been developed to cope with the problem. Many researchers have investigated the dynamic instability of the brakes concerning friction characteristics and show that it can be produced by the contact surface roughness and the nature of frictional mechanisms. They have mentioned that frictional forces varying by time and act as the leading role for brake noise [3]- [7]. Also, a significant number of studies have recorded that overheating of the brake system and deterioration of friction coefficient can cause brake noise [8].
On the other hand, the thermo-elastic expansion is developed by the frictional heat generation during braking and sliding with friction material on brake disc which is affected the contact pressure distribution [9]. On the thermo-mechanical point of view, the local high-temperature areas are generated by the high speed of sliding and inappropriate ventilation properties which are known as hot spots and the system can be unstable. The non-uniformity of contact pressure may cause vibrational frequencies at a low and high level which are known as hot roughness or hot judder or brake squeal [10], [11]. Therefore, braking system components are needed to have sufficient cooling properties to improve the braking performance which is more significant in high-performance passenger vehicles [12]- [16].
Manufacturing technologies improvement has enabled us to design, fabricate and use complicated structures in different kinds of science. Therefore, periodic cellular metals (PCMs) have considerably developed and allow us to use these structures in automotive engineering concerning lightweight properties and the capability to absorb energy, create unilateral fluid flows and improve heat dissipation. PCMs have open cell topologies with sandwich panel structure which enables them to heat dissipation through the thin truss ligaments and cell areas [17]. This capability is significant that can be used to absorb the contact energy through the sandwich structure's faces and transfer it to cell cores to improve the thermal transport. This is conducted through highly porous structures with 20% or less of their interior volume which has occupied by metal [18], [19]. On the other hand, PCMs with a series of three-dimensional truss ligaments that are arranged in a distinct pattern and repeating in all directions with a symmetrical geometric shape are named lattices. Contiguous media of lattice structures are made them a high potential structure to thermal transportation by coming up with a high thermal conductivity area and transferring heat to a cooling fluid or coolant flow through the lattice media [20]- [23]. The main significant advantage of lattice structures is to set up structures with considerable stiffness, strength, and load-bearing by using as little material as possible [24].
From this point of view, new designs to improve the effectiveness of thermal properties in passenger vehicles with the ventilated braking system are conceivable. Accordingly, a lattice design with high porosity and thermal capabilities have used as the core of brake disc to develop the ventilation and thermal capabilities of the braking system [25], [26]. The brake disc with lattice structure has shown more efficiency in absorbing heat and transferring it through its core cells [25]- [28]. The instability studies have indicated that it has a lower propensity to instability and also, at the low range of instability frequencies have shown that is more reliable [27]- [30].
In this research, a lattice brake disc with I-type design is modelled and the main parameters of I-type lattice are presented. The dynamics of motion concerning the principal coordinates and excitation caused by the distributed axial load is introduced to indicate the modal response of the structure and define the propensity of the lattice brake disc to instability. Referring to the model, validation is done by the impact hammer test to identify the experimental modal characteristics of the brake disc and compared with FE analysis at a free-free State condition. Subsequently, the validated model is used to nonlinear thermo-mechanical stability analysis conducting by the complex eigenvalue method to evaluate the contact pressure behaviour of new brake disc with lattice design. Sensitivity analysis of independent lattice parameters concerning the thermal and contact pressure properties is performed and the likelihood of instability occurrence is quantified by the definition of a single indicator derived from the system eigenvalues.

II. BRAKE INSTABILITY LITERATURE
Brake instability mechanisms have been investigated by many researchers about the different aspects of the problem. An important and essential study of disc brake vibration has done to examine the influence of the geometry on the modes of vibration of the brake rotor [31]. Moreover, the modecoupling type of instability has been used to investigate the prediction of unstable vibration modes using the complex eigenvalue analysis (CEA) [32]- [34]. The closeness of squeal natural frequencies and mode shapes of a rotating brake disc with the natural frequencies and mode shapes of a stationary rotor have studied by Nishlwakl et al. [35] and Monteil et al. [36]. Instability analysis of the brake system has been carried out through the non-linear contact pressure analysis by several studies [37]- [42]. Brake disc instability concerning different frequency ranges has been studied by Dunlap et al. [43].
The complex eigenvalue analysis (CEA) adoption with finite element analysis and linearized stability analysis of brake squeal have been done by some researchers [44]- [47] Lindberg et al. [48] studied the effect of vehicle speed and brake pressure on the roughness noise inside the vehicle Hetzler and Willner [49] and Choi et al. [50] developed the stability analysis with respect to the contact properties. Also, the transient dynamic analysis of brake system has been used to investigate the instability of the system concerning the divergent vibration of time response, by converting the time domain information into the frequency domain and use of the Fast Fourier Transform (FFT) technique [51], [52].
On the lattice point of view, dynamic analysis of periodic structures is necessary to define the natural vibration analysis of the structure consist of substructures with unit cell geometry that repeated symmetrically. Hence, linear differential equations of motion concerning the dynamic matrices of symmetric substructures can be described as a set [53], [54]. Also, eigenfunctions of the resulting periodic unit transfer matrix can be used to obtain the frequency response of VOLUME 8, 2020 the complete structure without increasing the analysis variables. Considerable studies on the forced response due to a rotating force on rotationally periodic structure have been done by Wildheim [55], [56] Noor and Mikulas [57] investigated the status of continuum modelling for large repetitive lattice structures, characterization of the continuum model, and application of continuum models to stability problems. Anderson and Williams [58] studied on vibration analysis of arbitrary lattice structures having repetitive geometry in any combination of coordinate directions. Karamoozian et al. investigated the dynamics and vibrational behaviour of lattice brake disc and have evaluated the instability behaviour of the lattice brake disc through the complex eigenvalue extraction technique [27]- [29].

III. LATTICE BRAKE DISC MODEL
Truss topology in lattice design has the main effect on the lattice loads supporting properties. Therefore a new design of lattice which is named I-type has introduced. The lattice design is strong enough to support the mechanical loads, also has proper ventilation and heat transferring properties which have been validated by Karamoozian et al. [27]- [29], [59] to use in the braking system. The shape of the unit cell and important parameters of the lattice model are presented in Figure 1.
The unit cell of I-type lattice structure is used to create the brake disc structure as shown in Figure 1(a) and define the micro-macro relations of the lattice structure. The Lattice beam with I-shaped periodicity cell is shown in Figure 1(b), to show the case of the structure in existence of rotating elements, where H is the height of the I-type lattice and t 0 is the thickness of the up and down plates of the I-type lattice sandwich structure. It should be noted that the coordinate systems which are at an angle with respect to the principal coordinate systems are used. Subsequently, the rotated elements will show monoclinic-like behavior in spite of the sandwich structure made of the orthotropic material, if the principal material coordinates not related to them.
The main lattice parameters consist of relative density and truss mass fraction factor are presented in Table 1 with respect to the geometry, truss length, truss angles and crosssectional areas in the unit cell. Relative density (ρ) indicates the quotient of truss metal volume to the unit cell volume and the truss mass fraction factor (η) specifies the lattice core elements to the stiffness contribution of the lattice structure concerning performance requirements.
It should be noted that l is the length of truss, ω is the angle of an inclination concerning to the truss member and the base of the unit cell, w is the width of a rectangular cross-section, b is the width of brazed node area.

IV. CONTACT DYNAMICS
It is considered that a disc with periodic lattice design rotating and lined on an elastic body as can be seen in Figure 2. An oscillatory and distributed motion moves across the structure and externally excites it. Also, the vibration of an elastic FIGURE 1. Schematically: (a) the I-type lattice unit cell (b) periodicity cell of the unit cell in principal directions [27] (c) the lattice brake disc [28]. body in neighbouring can be caused to oscillatory motion which moving in the lining surface and the contact area. Also, it is considered that the lattice disc and the back-plate of the brake pad interact in terms of the zero nodal circle, and interaction with the friction material of the pad is instant. The relative sliding of the elastic body on the surface of lining concerning friction material is caused to the distributed frictional traction and the distributed axial load µp can be obtained by transforming the traction surface onto the rotating lattice disc and area of the lining [60], [61]. The friction coefficient in the contact area between the lining surface and  [28].

FIGURE 2.
Schematics of the coordinate on the rotating disc and beam model considered in the present study, beam subjected to a distributed axial load due to friction concerning to free body diagram of the beam element, G = 0, G = 0 [28]. the rotating elastic body is considered as µ and the normal pressure in the contact area is considered as P. It should be noted that concerning to the interface tribology complexity different models of contact friction can be used [62].
To evaluate the instability situations on the contact problem concerning to the rotating lattice disc, the following assumptions are applied: (1) the contact between the rotating elastic body and lining surface is conformal, so the interfaces have no contact loss; (2) the lining area material properties and friction coefficient are constant during the moving; (3) the coupling of motions in axial and transverse directions is ignored. It should be mentioned that concerning the natural character of friction, the distributed axial load can be considered as a non-conservative and slope-dependent follower-type force [60], [61]. The schematics of the freebody diagram of the rotating lattice disc respected to the distributed axial load is presented in Figure 2, where θ is considered as the spatial variable,t as a temporal variable, and r e as the transverse displacement of the beam [60], [61].

A. CONTACT STIFFNESS
Theory of periodicity for the beams and rods is used to define the stiffness of the lattice disc structure concerning the construction of lattice plate periodicity cell by beams and rods. The mechanical properties of lattice plates have an interest in studies of several researchers regarding the theoretical and computational methods using the inverse approximation of framework modelled continuum structures [63]. The homogenization method by incorporating the classical strength theory is used for expanding the periodic lattice plates theory [64]. The homogenization method extracts the kinematic and equilibrium theories concerning to relations of the joint nodes in a finite-dimensional cell problem and principal coordinates [65]. Therefore, in this study, the stiffness properties of the lattice plate with structure and periodicity cell ( Figure 1) are evaluated through the method which has been investigated by Noor [63], Kalamkarov and Kolpakov [65]. The I-type lattice plate bending stiffness concerning the mass fraction factor and homogenized theory can be found in [27] as: where young's modulus of the up and down plates is E 0 , plates Poisson's ratio is ν 0 , thicknesses of plates is t 0 , the lattice truss thickness is t and the lattice truss young's modulus is E, respectively. ω is the inclination angle between the truss member and the base of the unit cell and the height of the lattice plate is H . (Figure 1 describes the geometric parameters). It should be noted that the lattice plate has symmetric geometry, so the coupling stiffness is neglected because it has zero value. Therefore, contact problem can be modeled through the validated method used by Ripin and Bin [42], Kang and Tan [60], [61] to define the stiffness of the contact interface by assuming the unit cells as the spring elements as shown in Figure 3.

B. PRESSURE LOAD DISTRIBUTION
It is needed to define the distributed load in the polar system. Based on the assumptions, the instantaneous distributed normal pressure P can be written as the sum of the static (P s ) and dynamic (P d ) normal pressures [60], [61]: Considering G and K as the shear and transverse moduli of the lining, respectively, P can be shown with the relative displacement of the lining and its elastic properties, as: where the neighbouring elastic body displacement in the transverse direction is indicated by r e . The lining deformation VOLUME 8, 2020 related to the static normal pressure (P s ) is represented by and the deformation affected just by the loads of local pressure on the condition of G = 0, as illustrated in Figure 2. Nevertheless, the deformation profile is usually considered for most elastic materials as represented in Figure 2. Noted that the lining shear modulus has a considerable effect on the instability behavior of the model and should be considered [61]. On the other hand, with consideration of polar coordinates, the general equation of a circle with a center at (r 0 , γ 1 ) and radius A, is: where A = r 0 , γ 1 is the curve assessment angle and has a relation γ 1 = Arctan (s), the curve slope is defined by m and the equation can have a solution for r in the general case. It should be noted that many mechanisms involving contact between two elastic bodies can be used to set up a coupled system, therefore, it is concentrated to investigate the parametric instability of rotating lattice disc in nonconservative condition with respect to the rotating distributed load as the excitation source. On other hands, assuming that the elastic body vibrates in terms of a sinusoidal traveling wave, the displacement r e in the transverse direction can be explained as: where the travelling wave component is shown by the sine term, the amplitude is a constant value and is described by a 0 . The angular velocity and the wave number are represented by the and n, respectively. The transverse vibration of the moving body at frequency ω f is represented by the cosine term. In the same way, by replacing the explained equations, the distributed axial load is given: (10) It should be noted that the vibration interaction on the contact area is created by a rotating disc and the pad affects by the source of the excitation applied to the disc or pad. In the present study, it is assumed that the excitation is caused by the dynamic force which is employed to initial excitation on the disc on rotating coordinates and the relations with the brake instability and other cases will look into future studies. Also, it is assumed that the vibration in the lining surface transmitted between the rotating disc and pad. The proportional vibration creation in one has resulted in the creation of displacement on the other one. It should be mentioned that the pad backplate is affected by the static pressure and concerning to the rotating disc it can be caused oscillation in the transverse direction, so the excitation considered as a moving distributed contact load f nm (r, θ, t). Therefore, θ is considered as the principle rotating coordinate and rotates with the disc as shown in Figure 2. Subsequently, a dynamic distributed load excites the rotating disc and the response of the disc in terms of mode pairs is evaluated.

C. DYNAMICS OF MOTION
In this section, the necessary basic motion equations of rotationally periodic lattice structure are introduced and are applied to a rotationally symmetric structure. Concerning the principal coordinates, equations of motion are [28]: where M nm is the constant for normalizing the corresponding eigenfunction, and the natural frequencies and the generalized forces are ω nm and Q nm (t), respectively. The normalization Rω nm (R)/M nm = 1 is used for simplification. On the other hand, consider a disc rotating at angular velocity and an axial distributed load f nm (r, θ, t) is acting dynamically on the disc at the radius R. Concerning to the symmetry of disc and rotation axes, the forcing function can be written with respect to the principal coordinates rotating with disc: where δ[θ − n ( t + α)] is the Dirac-delta function. It represents a function which is zero for all θ except = ( t + α), that it represents a unit response. Hence the force distribution is given by: The force on the plate can be defined by an analytical method where the forcing function can be expanded concerning to the principal coordinates rotating with the disc and integrated with r and θ directions: The generalized forces in terms of the nm mode are: where directions on the polar coordinates are r and θ, the mode shapes are represented by ϕ nm (r, θ) and are characterized by setting the eigenfunctions at n nodal diameters and m nodal circles (ϕ nm (r, θ) = ϕ cnm (r, θ) + ϕ snm (r, θ)) . The ϕ nm (r, θ) is given by ϕ cnm (r, θ) = cos nθ and ϕ snm (r, θ) = sin nθ. Hence, the generalized force has the nm mode terms of sinus and cosinus for the disc [28], [29]: The Eqs. (15) and (16) represents the generalized forces for the modes. Subsequently, the convolution integral is used to define the generalized responses (or normal responses) at the very lightly-damped system. In order to define the response, the convolution integral should be calculated by substituting Eqs. (10) and (16) to the convolution integral. The terms of natural frequencies ω nm are existed with the equations of the transient response and the investigation is concerned with the steady-state problem, which is considering only the terms are relevant to the excitation.
Knowing the generalized coordinates for the nm modes, the response is calculated using the mode-summation formula X dnm (r, θ, t) = ∞ n=0 ∞ m=0 q dnm (t) ϕ dnm (r, θ). Therefore, the response modes of the disc in terms of rotating coordinates is obtained: For another term: And using the solution at steady-state condition results in the disc response as: (20) where the m disc nm is related to the disc and is the the modal mass or generalized mass of nm modes. Eq. (20) describes the response of rotating disc at the nmth mode which is excited by the distributed axial load. It is a complex vibration at the frequency (n ) and consists of 'fixed vibration' terms. For the disc, it can be shown that ξ 1 = ξ 2 , so the response is a pure backward traveling wave. The response can be simplified to: Substituting the θ by (θ 0 + t) in Eq. (21), is resulted in the response equations of the rotating disc in the stationary coordinate θ 0 , after simplification have: From the stationary observer perspective, it could be understood that there are two standing waves which may be deformed the disc rim with a cosine shape. On the other hand, paying attention to the Eq. (22) represents that the standing wave is independent to the angular velocity and VOLUME 8, 2020 can be found at all rotating speeds. This kind of excitation can be formed on both sides of the brake disc and is due to performing a non-uniform pressure distribution in the contact area.
Hence, the general equations are determined for the modes of a disc rotating past the distributed axial load to simulate the travelling waves. Referring to Eq. (20), the generalized forces have sinus and cosinus terms with a harmonic force and frequency (n ) that are indicated the forcing function of a normal mode. Concerning to the nodal points positions of the nm modes, two equal harmonic forces are used to define the responses. Noted that the resonance possibility is occurred when the force is constant in time (ω f = 0 ) and the possibility occurs for n = ω 1n , which is indicated that the natural frequency is for backward rotating mode (backward is determined in relation to the direction of body rotation).
It could be derived that for any excitation frequency, there is a rotation speed at which resonance occurs. In other words, two apparent resonance frequencies occur at each rotation speed which one is below and the other one is above the actual natural frequency. Generally, the occurrence of all resonance cases at the same time is impossible but changing the forcing frequency, the angular velocity or both of them will result in each of the resonance cases.

V. EXPERIMENTAL TESTS AND SIMULATIONS
Concerning the rotationally coordinate, the excitation of a moving body at a frequency ω f on the rotating disc can be resulted in two frequency components (see Eq. (21)) which are indicate the modal behavior of rotationally disc. Therefore the response frequencies at an excitation frequency near to a natural frequency are analyzed and are separated from any adjacent modes in real conditions. At the present study, there are close nm modes in pairs and the mode numbers 2 to 8 are considered. It should be noted that it can be applied to the other diametral modes. Experimental modal analysis (EMA) is established to verify the accuracy of the analytical and FE models through the experimental determination of natural modes and frequencies at the free-free condition. The stability analysis is applied to define the natural frequencies and unstable modes. The modal characteristics of the brake rotor/pad are defined through the experimental impact hammer test. The brake disc parts natural frequencies information are evaluated by using the transfer function method. Therefore, The Fast Fourier Transform (FFT) is used to determine the functions of frequency responses. Figures 3 and 4 are presented the real time FFT plots for the vanned type brake disc/ pad system. The first peak frequency amplitude for the brake disc is turned up at 920 Hz ( Figure 3) and for the pad is occurred at 1153 Hz ( Figure 3). The complete procedure and details of experimental modal analysis (EMA) can be found in a previous study by Karamoozian et al [28], [29].

A. ANALYTICAL AND FE FREE-FREE MODAL ANALYSIS
The modal characteristics of the brake disc/pad are analyzed through the finite element technique at free-free conditions  to define the inputs of the FE model for the lattice type. Brake components with homogeneous material properties characteristics are applied as detailed in Table 2. The accuracy of FEA modal analysis to define the natural frequencies is enhanced through conducting the mesh density tuning on the components. The brick elements with 8-noded (C3D8) are applied to mesh the structure. The natural frequencies of the system are extracted by the eigenvalue technique. Noted that the quantity of model degrees of freedom should be equal to the quantity of natural frequencies of the finite element model [66], [67]. Accordingly, the Lanczos method is preformed, concerning the degrees of freedom in the system. Table 3 and Tables 4 and 5 are present the analysis outcomes of vanned and lattice types brake components respected to the 15 kHz extraction and nodal diameters from one to 8. It can be observed that all nodes and anti-nodes are spaced equally because of the axisymmetric of discs. Noted that the bolt holes of the top hat section demonstrated an effect on the diametral mode without any significant effect on the rubbing (contact) surface. A comparison is made between the EMA natural frequencies and the natural frequencies are obtained by the FE free-free modal analysis. The investigation demonstrated a good correlation between the experimental and finite element simulation results. Therefore, the FE model is validated for executing instability analysis.

VI. NONLINEAR THERMO-MECHANICAL FE MODEL
Comprehensive numerical modelling of braking instability has required the influence of thermal effects under different operating conditions that include pressure, velocity and   coefficient of friction. Temperature is changed due to frictional heat generation caused by the axial and radial deformation of the rotor along with pad expansion [40]- [47], [68]. The contact area between the rotor and pad surface is affected by changes in the shape of the rotor or pad and subsequently is influenced the load distribution at the friction interface.
Therefore the modelling of the braking instability should be undertaken using a thermo-mechanical model which can be applied to define the effect of thermal deformations when the pad and disc are in contact [69].
The ABAQUS software is used to present a non-linear finite element model for brake disc and pad at the contact interface which is represented the vanned type model at first and subsequently is used to analysis of the lattice brake model. The vanned type brake disc model has been used by a number of researchers to investigate the contact pressure distribution [31], [42]- [46], [41], [70], [71]. In the present research, the vanned type brake disc model is used to investigate the performance of the newly developed methodology before implementation within the more complex lattice-type brake disc model. Same identical control strategy properties concerning to contact control, material properties, structural constraints and thermal constraints are used in FE models.
In common with accepted practice by Liles [44], the stability analysis of the brake system is determined by the process of complex eigenvalue analysis (CEA). The sign of the real part is indicated stability and the imaginary part has defined the frequency of vibration of the unstable mode [38]- [72]. The mechanism is responsible for any unstable behaviour lies in the asymmetry of the stiffness matrix caused by frictional coupling at the disc-pad interface [37]. The initial thermomechanical contact simulation is carried out using the validated vanned type brake disc model. The three main areas are addressed within the context of the contact problem through the thermal and contact pressure analyses. This phenomenon contributes to the nature of the engagement that occurs when the pad is brought into contact with the disc. The thermomechanical FE model of the braking system concerning lattice brake disc parameters is established. Instability analysis is done concerning the effect of lattice parameters on the strength, thermal behaviour and contact pressure which are led to instability propensity.
The thermal analysis has described the temperatures predicted during the brake application whereas the contact analysis results are focused on the thermal contact pressure distribution. The results are predicted the transient thermal and structural performance of the system, including temperature gradients, contact pressure distribution of the system as it undergoes the SAE J2521 drag braking procedure. The results of the thermo-mechanical analysis are provided with the temperature-dependent transient contact pressure distribution which is an essential element of the brake instability simulation. Subsequently, the instability analysis of the braking system is investigated through the complex eigenvalue method on the proper frequency range between 1000 Hz to 15 kHz. The non-linear thermo-mechanical analysis is used to derive the response of the system concerning the transient contact pressure distribution at particular times of braking period. The unstable modes are defined by evaluation of the several unstable mode shapes and the associated frequency of them.
The overview of the non-linear thermo-mechanical contact analysis is shown in Figure 5. The approach is used to control the execution of the FE model and to derive the temperature and contact pressure distribution. This includes conditions of drag brake operating, material properties, structural constraints, thermal properties and constraints, the method used to rotate the disc, the contact scheme and its solution, the cooling parameters and the sensitivity of the model output to mesh density. Initial simulation results from the vanned type brake disc are presented followed by those from the lattice typed disc brake model. Results are presented to show the effect of the friction force and applied pressure concerning vehicle speed on the system transient temperature and contact pressure distribution which is generated within the brake system model. Subsequently, the results of the transient complex eigenvalue analysis are used to discuss the instability behaviour of the vanned type and lattice-type brake.
Generation of the transient heat drag braking is caused by the operating conditions input in drag braking respected to the temperature dependency of brake components material properties, which are followed by the brake assembly constraints (thermal and structural behaviour). Using proper criteria for the control strategy and contact elements at the contact surface of finite element models is caused to accurately and efficiently capture the thermal and structural behaviours that are presented at the interface.

A. OPERATING PARAMETERS
As shown in Figure 6 in the lattice type and vanned type brake disc models loading is caused by applying a uniform pressure at the backplates of pads. The applied uniform pressure is transmitted to the outboard pad through the paw fingers and to the inboard pad through the piston. The stationary heat source on each pad is governed by the rotational velocity and friction coefficient of the disc and is imposed until completion of the drag braking event duration. In actual brake events, the coefficient of friction varies between each brake application. However, the applied coefficient of friction is considered to be realistic and constant and was 0.3 with alternative relatively high values of 0.4. Each analysis is run for 10 seconds braking duration.

B. MATERIAL PROPERTIES
The temperature dependency properties of cast iron material in Koetniyom [73] are used to fully explore the effect of temperature dependency of the disc while those used by McKinlay [74] are employed for the friction material of the brake pad. The Handbook of ASM Specialty [75] is used to derive the material properties of the steel backplate used for the paw and piston elements. By increasing in the temperature in components of the braking system during braking, the component's material properties are automatically adjusted to reflect their dependence on temperature.

C. STRUCTURAL CONSTRAINTS
Structural constraints initiate the generation of tangential friction when the disc and pad surface are in contact. The disc can rotate freely but it is constrained at the top hat axially. To restrict the pad motion during the dragging process action by the rotating brake disc fully constrains in the in-plane direction are applied to the elements of pads, paw and piston. Applied constraints with respect to the dynamics of motion allow the disc brake system to generate heat, distributes pressure at the contact interface and to create structural deformation in components concurrently throughout the braking period over different operating drag braking conditions.

D. CONTROL STRATEGIES
Two methods of rotating which are beam element connector and rigid body connector are applied to the model. As both techniques are able to rotate the disc [76], the choice of the technique to be selected depends on the computation time needed to fulfil the simulation completely. The results of the investigation are performed using typical drag brake parameters and the penalty approach is shown that the rigid body connector results allow more than 50% reduction in the run time of the central processing unit (CPU) as well as efficient computational data storage of the generated results. Longtime of computation and huge amount of data storage are decreased significantly by using the rigid body connector; therefore, derivation of the contact pressure is implemented by this method and also is used to define the instability behaviour of the braking systems.
Thermal conditions at the contact interface are depended strongly on the partition of heat transferred between the contact disc and pad surface. The assigned finite sliding algorithm at the contact interface is calculated the frictional heat by bringing the slave nodes of the pad into contact on the entire disc master surface as the vehicle speed 3 or 10 km/h is imposed at the axis of rotation, with the assigned value of coefficient of friction 0.3, 0.4 or 0.5. As time is elapsed, the heat is moved to the neighbouring elements that are come into contact relative to the sliding interaction and this has resulted in the generation of phenomena which is called ''rotating heat source'' effect. The heat is circumfused on the surface of the brake disc and gotten away from the system based on the transient frictional heat generation during sliding contact. Also, heat can be transferred between the assembly components through the process of conduction and can be removed from the exposed surface of the disc. Also, pad transfers heat through convection and radiation.
The thermal and mechanical responses during frictional sliding can be determined by use of either the penalty or Lagrange multiplier methods. The time taken for the Lagrangian method is higher as no-slip is allowed when the sliding bodies come into contact. As more degrees of freedom are added to the model, the number of iterations are needed to complete the contact solution. This may have resulted in some convergence issues. By using the penalty method for defining the convergence on the analysis of the surface contact behaviour is easier with respect to the reduction of the computational time, so the method is used in this research. The cooling strategy is based on the SAE drag brake test matrix that involves periods of cooling. During these periods, no frictional heating takes place as the applied pressure is zero. The convection heat transfer coefficients are employed to allow the brake to cool.
The temperature-dependent values of convection and radiation coefficients are calculated and are shown in Figures 7(a) and 7(b) for the vanned type and latticetype brake systems, respectively. The radiative heat transfer values of the system at low temperature are very low in comparison to convection values. This shows that radiation is less effective at low speed and low temperature which explains why many researchers have omitted the radiation effect [41].

E. MESH SENSITIVITY
The contact interface performance depends strongly on the distribution of pressure across the contact surface. The level of mesh density in the contact area may, therefore, influence the accuracy of the results as well as the computation time.
To conduct the mesh sensitivity analysis, different mesh densities are generated and are analyzed as shown in Table 6. The baseline mesh model (BMM) is created and validated through the free-free modal analysis. The BMM is subsequently compared with three other levels of mesh density across the thickness of pads and disc. The number of elements for the coarse mesh model (CMM) is reduced to 15% compared with the baseline mesh model (BMM) whilst the number of elements is increased by 15% for the fine mesh model (FMM) and up to 30% for the very fine mesh model (VFMM). The element size variation is just applied to the contact sliding body, the pad friction material and the disc rubbing surface where the contact occurs.
A similar type of drag braking simulation together with identical thermal and structural constraints are used to illustrate the effectiveness of the different FE models and conduct the thermo-mechanical contact analysis. The element sizes for the different mesh density models are illustrated in Figure 8(a) for the vanned type model and Figure 8(b) for the lattice type model.
A mesh sensitivity analysis is performed to determine the appropriate mesh density of the FE model concerning the mesh element type and size. To conduct the mesh sensitivity analysis, different mesh densities are generated and analyzed. Two types of the element which are 8-noded brick elements and 20-noded iso-parametric brick elements are considered. The results of simulations are given in Table 3 and the accuracy of the mesh density models with respect to the CPU simulation time are compared. Table 7 demonstrates that the 8-noded brick element with fine mesh not only provides accurate results but also saves significant computational CPU simulation time compared with the finer mesh models. Consequently, the lowest integration order which produces acceptable results is used to reduce processing time. Also, it should be mentioned that in comparison with modal analysis, the CPU simulation time is considerably higher for the instability analysis concerning the contact approaches.
The effect of change in element density on the interfacial contact pressure distribution for the vanned type and latticetype models are shown in Figure 9 and Figure 10, respectively. A line plot along the central arc of the pad mean rubbing radius for the four different mesh density levels is given for the isothermal contact conditions and a stationary rotor (t = 0 seconds). It can be seen that the BMM and FMM models, as well as the VFMM model, show the negligible variation of contact distribution as the pressure is exerted on the pad surface. Since only differences in minor ranges are observed among the mesh densities, It is concluded that the baseline mesh (BMM) is sufficiently refined for capturing the transient behaviour of contact and subsequent instability analysis.
Accumulated temperatures at specific nodes in the vanned type model using the four levels of mesh density are given in Table 8. The temperature is predicted at the end of the drag braking warm-up period for the disc surface temperature   (node 2045) and pad (node 7569). It is shown that the nodal temperature is predicted either has growth by no more than 2 • C for the finer meshes (FMM and VFMM) or has reduction by only 1 • C for the coarser mesh (CMM).    The predicted temperatures from the coupled thermomechanical analysis of lattice-type model at the end of the first braking event are employed for four different levels of mesh density and are given in Table 9. Similar to the vanned type model, the mesh density at the baseline level (BMM) is adequate. Therefore is retained throughout the analysis since the different levels of mesh density are not strongly influenced the results.

VII. STABILITY ANALYSIS RESULTS WITH ASSOCIATED MODE SHAPES AND FREQUENCIES
The instability characteristics of braking systems are extracted at a specified time of braking. The unstable modes are defined respected to the adopting diametral modes of discs which are affected by the increment of frequency. Tables 10 and 11 are presented the CEA results of instability characteristics of vanned and lattice-type brake models.  The instability frequencies are predominated by the natural frequencies of the disc which are in line with published work carried out by Nishlwakl, et al. [35] and Dunlap, etal. [43]. Therefore, the instability can be raised by the coupling of disc diametral modes with the pad either twisting or bending modes. Hence, respected to the theoretical studies, the frequency of resonance can be determined or controlled by the excitation frequency and the diametral natural frequency. Also according to brake instability characteristic description, the modal spacing of the disc can be used to define the difference between low-frequency and high-frequency instabilities [43]. The modelling outcomes as described in Table 6 show that the low-frequency of instability is occurred at the second, third and fourth diametral modes for both brake disc types. On the other hand, high-frequency of instability is occurred at the fifth, sixth, seventh and eighth diametral modes. It can be observed that at the frequency of 5202 Hz, the first twisting mode is appeared which is categorized in high-frequency instability range. Respected to the results, it can be suggested that combing of the bending mode of the pad with the diametral mode of disc cause to raise low-frequency instability and that the pad twisting mode generates high-frequency instability.
This study points out that the lattice type brake disc is represented high-frequency instability range at 1349 Hz at mode 8ND. On the other hand, the vanned type is represented instability frequency at 9196 Hz at mode 8ND. Noted that in braking system instability analysis, higher range of instability frequency often observes at a frequency range between 4 to15 kHz and at higher sequence of vibrational modes which include 5 to 10 nodal diameters. The results indicate that the range of frequencies estimated by the analytical analysis is very close to the complex eigenvalue (CEA) results and have a maximum difference of 8%. Evaluation and comparison of brake disc at the exciting status and the free-free modal status are demonstrated that the CEA is done in a very close estimation on frequencies of instability and associate disc mode shapes. Therefore, it can be propounded that the natural instability frequencies and associated mode shapes of a rotational brake disc were both in VOLUME 8, 2020  the adjacency of the natural frequencies and mode shapes of a rotor with stationary status as it has achieved by other studies [32], [43], [44], [77].

VIII. SENSITIVITY ANALYSIS OF CONTACT PRESSURE CONCERNING TO LATTICE PARAMETERS
The lattice models of present interest (Figure 1) consist of arrays of struts in the I shape configuration. The geometry is characterized by three independent parameters which are: (i) the truss angle (θ) that is measured relative to the vertical plane of nodal points, (ii) the truss mass fraction factor η and (iii) the relative density ρ. Four models for lattice geometry in our analysis are considered and the geometrical descriptions of lattice models are shown in Figure 1.
The overarching objective is employed to analyze the effect of different I-type lattice parameters and structures on contact pressure area which enables stable design with high compressive strength, thermal and stability properties. The objective is achieved through a combination of lattice models and finite element analysis of lattice geometries and is characterized principally by the angle of the strut members θ, relative density and mass fraction factor. Parenthetically, there are multiple combinations of θ and ρ that yield a specified η. Thus, the lattice geometry and the lattice response cannot be couched solely in terms of ρ. Furthermore, using the finite element methodology is enabled to found out that the lattice geometry can be designed over a wide range of lattice geometries (characterized by 25 . Noted that for all subsequent simulations, only the case in which the lattice is fixed to rigid plates is considered. Accordingly, the models are based on the assumption that the loads are transmitted through the struts. It should be mentioned that all the four I-type lattice models are satisfied with the requirements to use in brake disc design and geometrically detailed in Table 12.

A. EFFECT OF CONTACT PRESSURE
Non-linear contact pressure analysis is carried out to investigate the contact pressure between the pad and brake rotor as a function of lattice parameters (θ, ρ, η) and the thermal loading conditions. The contact pressure distribution enables the effect of thermo-mechanical interactions during braking to be highlighted which is an essential element of the brake instability simulation. The extent to which the interface pressure distribution is influenced by the thermo-mechanical response of the structure is discussed. Figure 11 and Figure 12 are presented the interface pressure associated with an arc of friction material at its mean rubbing radius for the isothermal case (t = 0 second) for the vanned type and lattice-type brake disc. Also, Figure 11 and Figure 12 are demonstrated the transient thermal analysis results at zero, 1, 6.0 and 10.0 seconds concerning to the piston loaded (inboard) and paw loaded (outboard).
The influence of mass fraction factor, relative density and pressure on the distribution of contact pressure throughout the braking process for the lattice-type structure are represented in Figure 13. The results are represented that the mass fraction factor and relative density are substantially influenced the   contact pressure distribution. On the other hand, higher mass fraction factor (lower relative density) is generated more loss of contact (indicated by very high contact pressure) at the paw loaded area concerning to the trailing edge and also at the piston area regarding the leading edge. It can be observed that the centre of pressure in the case of the piston loaded immigrates by increasing time towards the pad leading edge. Subsequently, it may be deformed by thermal loading. The pads at the trailing edge are highlighted for more surface separation, especially at the paw loaded area. The results are demonstrated that the mass fraction factor has a moderate influence on the interfacial contact pressure distribution. When the mass fraction factor is reduced (ρ increased), the capability to uniformly distribute the loads from the lining to the friction interface is decreased and leading to a less uniform interfacial contact pressure distribution.
Accordingly, a truss mass fraction factor η has introduced, which leads to an expression that separates components of the core topology into those that directly contribute to the stiffness from those that are apply to achieve robust performance. However, at lower mass fraction factor (higher relative densities), the struts become increasingly thick and provide not enough opportunity to conduct the heat from the hot face sheet into the strut structures. This results in a less temperature gradient through the thickness of the core. If a coolant flows through the core, heat transfer at the metal coolant interface occurs less efficient and subsequently, less heat can transport away from the struts and the average temperature of the coolant cannot enough raising and increasing that of the lattice metal core of brake disc.
In order to compare the interface pressure distribution of the four lattice-type models with respect to the different coefficient of friction, the interface pressure distributions are associated with an arc of friction material at its mean rubbing radius and are plotted in Figure 14 and Figure 15. The isothermal case (t = 0 second) and the transient thermal analysis at 1, 4.0, 6.0 and 10.0 seconds for the piston loaded (inboard) and paw loaded (outboard) pads at two values of interface friction µ = 0.3 and µ = 0.4 are considered.
The results are highlighted that the friction coefficient substantially influences the contact pressure distribution and the heat transformation has a prominent effect. The lattice structure with higher relative density and lower mass fraction factor presents more contact loss at paw loaded area regarded to the trailing edge and at the piston area regarded to the leading edge concerning to process of comparing to the other models of lattices. Surface separation of the pads at paw loaded area respected to the trailing edge is more and it is because of the digging-in phenomenon [78]. So the piston side contact area cannot increase more meanwhile at the paw side surface separation continues concerning to the trailing edge. Therefore the lattice type with lower relative density and higher mass fraction factor parameters shows better results.
On the other hand, the variation of applied pressure is shown to change the contact pressure distribution over the brake application period as the surface of the disc becomes hotter. The higher temperature generates larger separation at the leading edge of the pad on the piston loading side and the trailing edge on the paw loading side for both vanned and lattice-type models. The interfacial contact pressure distributions over the contact surface for the lattice models at applied pressures P = 1.0 MPa, P = 2.0 MPa and P = 3.0 MPa are shown in Figure 16 for both the piston side and the paw loading side.
Applied pressure equal to 1.0 MPa is shown a more uniform pressure distribution compared to either 2.0 MPa or 3.0 MPa applied pressures. In other words, lower temperatures on the disc surface have resulted in more uniform contact pressure distributions. The higher values of applied pressure, are increased the surface temperature of the disc by creating more friction at the lining contact area.

IX. EFFECT OF DIFFERENT LATTICE CORE MASS FRACTION FACTORS AND RELATIVE DENSITIES ON LATTICE BRAKE DISC INSTABILITY
Definition of the relations between the relative density, mass fraction factor and contact pressure interfaces of the present lattice designs, are enabled us to investigate the instability propensity comprehensively concerning the main lattice parameters. Unstable frequencies are caused by coupling a diametral mode of lattice brake disc with the twisting or TABLE 13. Standard deviation of instability for different lattice core mass fraction factors and relative densities on lattice brake disc (P = 2.0MPa, µ = 0.3, v = 10km/h). bending mode of the pad. The instability frequencies are in command of the natural frequencies of disc in accordance with the published studies of Nishlwakl et al. [35] and Dunlap et al. [43]. In order to define the instability behaviour of the brake system and extracting the data, particular times are selected at the certain braking period and unstable modes are indicated through the adapting the diametral motion modes concerning to increment in frequency. The contact pressure distribution changes over the brake application period as the surface of the disc becomes hotter and this is well agreed by the results of investigation of the effect of pressure on the dynamic brake stability by Lee et al. [38] and Ripin and Bin [42]. They have shown that the increase of brake pressure leads up to increase in the braking system instability.
Liles [44] deduced that audible brake instability has a high probability of creating modes with significant instability and connected the instability propensity to the standard deviation of the instability measurements. Since instability propensity can be used to estimate by the standard deviation (SD) of the instability measurement within a frequency range [38]- [44]. The main assumption is that the larger propensity of the disc brake system to instability is estimated by the larger standard deviation of the instability measurements. Basically, the occurrence of instability has a statistical characteristic and there is a wide variation on the complex eigenvalues, also do not always bring to being unstable in reality.
In this study, the same method like as Liles [44] is used to define the instability propensity and it can be illustrated as a single number for the instability propensity of the different I-type lattice models under brake operating conditions, rather than a set of real and imaginary components. Therefore, the squeal propensity can be quantified by the equation of standard deviation as follows: where {g 1 , g 2 , . . . , g N } is the population sample of instability measurement, g is the average value of the whole population and N 1 is the total population of instability measurements. It should be mentioned that in the present study g = 0, so the conjugate pairs of the complex eigenvalues are considered. Noticed that eigenvalue pairs may be consist of positive and negative real parts equally. Therefore, the sensitivity of the brake instability to mass fraction factor and relative density of lattice structure is investigated concerning to the average standard deviations of the contact pressure force. The results are presented in Table 13. Sensitivity studies illustrate that the standard deviation of instability measurements and that of the contact pressure loads are both related to the mass fraction factor and relative density, and the effect is considerable. This indicates that decreasing the mass fraction factor or on the other hand increasing the relative density leads to less stability in the system that is more prone to instability. The eigenvalue plots of the instability measurement and the instability distribution measurements versus natural frequency for the different models are shown in Figure 17.
It can be observed that both the quantity and magnitude of instability measurements of unstable modes for the different mass fraction factor models, are increased slightly while the mass fraction factor is decreased. Also, an increment has occurred only within a relatively high-frequency range between 8 and 14 kHz and the system is becoming very noisy at η = 0.84, ρ = 0.45. Concerning the stiffness of the lattice brake disc (Eq. (5)), it can be seen that the instability is substantially influenced by brake disc stiffness and consequently the mass fraction factor and the contact pressure are applied to the system. Furthermore, highly unstable frequencies are predicted at lower mass fraction factor and higher relative density. Hence, the system is tended to show low instability frequencies at higher mass fraction factor and lower relative density (at the η = 0.87, ρ = 0.15 instability propensity is the lowest).
However, at lower relative densities, the struts are becoming increasingly slender and are provided ample opportunity to conduct the heat from the hot face sheet into the strut structures [20], [22], [31], [79]- [81]. This has resulted in a temperature gradient through the thickness of the core. The coolant flow is conveyed through the core, the heat transferring at the metal coolant interface is occurred and the heat is transported away from the struts and raised the average temperature of the coolant flow. Subsequently, temperature reduction at the lattice metal core of the brake disc has happened. Paying attention to Figure 16 it can be found that lower relative density is led to lower surface temperature and instability propensity.

X. CONCLUSIONS
In this study, the lattice parameters sensitivity analysis is used to investigate the instability analysis of the lattice brake disc. Analytical models and finite element calculations of the compressive response of I-type lattices are submitted for a broad parameter space. Concerning the lattice structure capability to increase the thermal properties, the system performance in term of rotor surface temperature, peak contact pressure and instability properties is shown to be sensitive to the shape of lattice, relative density and mass fraction factor (incorporated into the rotor design).
The higher relative density of lattice structure concerning the lattice strut angle is led to a higher temperature at the disc and pad surface. On the other hand, a higher mass fraction factor is led to a lower temperature at the disc and pad surface.
Concerning lattice design, higher values of mass fraction factor, are stiffened the disc structure and increased the stability of the area over which heat can be gone away from the disc. The lattice structure in the brake disc has decreased the surface temperature of the disc and consequently has reduced the peak contact pressure that experienced on both the piston and paw loading sides. The maximum contact pressure is observed in the model with less mass fraction factor and this is located at the leading edge of the piston side and the trailing edge of the paw side. More uniform pressure distributions are observed on both paw and piston sides when the mass fraction factor is increased. The results are leaded to define the optimal lattice designs which the inclination angle is of intermediate value θ = 35 • and the relative density of ρ = 0.31 and mass fraction factor of η = 0.85. CHIN AN TAN received the B.Sc. degree in mechanical engineering from UC Berkeley, the M.S. degree in aeronautics from the California Institute of Technology, Caltech, and the Ph.D. degree in mechanical engineering from UC Berkeley. He was the Associate Dean for the Academic Affairs in the college. He is currently a Professor in mechanical engineering at Wayne State University. His research skills and interests are in the fields of structural dynamics, vibration, and control; systems analysis; and mechatronics. He has many years of experience working on research problems with applications to automotive engineering (disc brake dynamics and advanced manufacturing). He is a Fellow of the American Society of Mechanical Engineers. YUANLONG WANG received the B.Sc. degree in material science and technology, and the M.Sc. and Ph.D. degrees in vehicle engineering from the Nanjing University of Science and Technology, Nanjing, China. He joined the College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, in 2017. His research interests include simulation and control of vehicle dynamic performance and dynamics of negative Poisson's ratio structure and dielectric elastomers. VOLUME 8, 2020