1 Introduction

Disc brakes are the most vital safety part of train vehicles in reducing, stopping and maintaining a constant speed. The high temperature created in braking may cause several negative effects on the brake systems, including, but not limited to brake fade [1], premature wear [2], thermal cracks [3] and reduced service life due to thermal fatigue [4]. These failures are predominantly expected in light rail transit trains, because of their frequent stops, due to many stations per route. To counteract such negative impacts, genuine FEM (finite element method) modelling with less expensive computational modelling is essential. According to the review conducted by the authors on railway disc brake [5], four different FE thermal modellings were highlighted: axisymmetric, uncoupled thermomechanical, coupled thermomechanical and non-axisymmetric moving heat source.

There is a large volume of published studies on axisymmetric modelling, because it improves the computational efficiency by assuming the brake pads to cover 360° of the contact surface [6]. Computational problems have been easily managed by assuming two-dimensional geometry. In spite of these achievements, the issue of axisymmetric modelling has been a controversial in many aspects. Firstly, while pad–disc interface is exposed to heat source, its discontinued position (surface free from contact) is exposed to convective heat dissipation. This argument agrees well with other studies [7], who confirmed presence of spatial temperature variation. Secondly, the amount of disc geometry sector taken cannot predict the actual number of hot spots created on friction surface. In [3], for example, the FE-based thermal analysis on 30° disc showed two and half hot spots; meanwhile, experimental investigation on full disc displayed six hotspot areas. Thirdly, although many studies revealed spatially (radially and tangentially) varying temperature [7, 8] and cracks of different sizes on friction surface [9], it couldn’t determine specific position of crack initiation site, as it assumes circumferentially equal temperature and thermal stress on friction surface.

In contrast, uncoupled thermomechanical modelling has been performed by determining heat flux from interface pressure calculation [7]. And, directly coupled thermomechanical modelling appears realistic in considering spatial variation in heat input [10]. In spite of their wear consideration, both models are vulnerable to high computational cost and computer memory usage. Consequently, both are limited to short braking times. Besides, convergence problem due to massive CPU computational efforts and huge data storage is also observed here [10].

Due to the problems raised in previous methods, non-axisymmetric thermal modelling of disc brakes expressed in moving heat source has emerged, attracting attention of many investigators [7, 11]. Because, variations in heat input and boundary conditions with space and time variables are easily managed by programs written in many software codes, including APDL [11], ABAQUS [12] and COMSOL Multiphysics [13]. Despite its effectiveness in managing spatially varying heat and convergence problem, there has been very few researches conducted on non-axisymmetric moving heat source, and still question of spatial and temporal variation in heat has not been fully answered.

Gao and Lin, 2002 [7] carried out moving heat source, using coupled thermoelastic contact analysis that enabled them to account for space (r, θ) and time variables. This modelling is further extended in determining the source of the thermal fatigue under residual stress analysis in [11]. The main advantage of the model was its presentation of radial thermal gradient, although its extent to which it varies circumferentially remained unclear. Furthermore, its application was limited to short braking time(less than 4.274 s), due to convergence problem in uncoupled thermomechanical (thermoelastic) analysis.

Furthermore, non-axisymmetric moving heat source was modelled in [14] for assessing the influence of convective mode of heat dissipation. Similar modelling was also performed in [15] on functionally graded brake disk material. In both studies, one rotation time of the disc brake was shared between heating time and cooling time, according to the proportion of pad cover angle (60°) to disc angle (360°). Full disc surface was seen heated or cooled alternatively by their corresponding time share within different load steps. Although radial variation in heat input was considered in [14] (but not in [15]), the method was still axisymmetric for each load steps.

Likewise, Baron Saiz et al. [16] studied the thermomechanical performance of disc brakes with different geometries using moving heat source. Even though it was revealed successful in selecting the best shape of the brake rotor, using spatially constant heat flux at pad–disc interface and using time independent heat flux application for a single wheel rotation were areas requiring further improvement. Similarly, rotating heat source was represented by rotating heat source over stationary disc in [13]. However, the model didn’t consider the spatial variations in heat input. After four years, the same authors modelled moving heat source by rotating disc layer close to friction surface around fixed heat source [17]. Though the model contributed significant drop in computational time, it requires further improvement in accommodating further complexity of disc and pad geometry, including disc surface bolt holes and pad grooves. At the same years, [10] modelled modified heat flux with cosine functions, to consider its discontinuity at friction surface. Still, the model lacked spatial variation in heat flux at pad–disc interface, due to pad and disc geometry, and absence of convective heat dissipation on discontinuous area of friction surface. Consequently, the results obtained by the model were seen oscillated about the axisymmetric model.

Finally, moving heat source was modelled by Pan and Cai [12] to investigate temperature and stress analysis of disc brake. It was reported that the model had better engineering value compared to other axisymmetric and fully coupled modelling. Besides, the actual physical movement or rotation of heat on friction surface was seen obviously, compared to other types of moving heat source modelling. Despite that, spatially uniform thermal results illustrated at pad–disc interface required further improvement. Furthermore, pad groove and convective coefficient effects on friction surface were fully excluded from the modelling. To the best of authors knowledge, no model revealed in the literature, in accounting spatial and time variations in heat input and boundary condition simultaneously.

This study proposes a new non-axisymmetric moving heat source (here in after abbreviated to NAMHS) FE modelling that accounts spatial and time variation in heat, with feasible computational time and computer storage. In addition, its effectiveness in spatial temperature assessment is evaluated by comparing it with two traditional models (used as a reference): axisymmetric and non-axisymmetric moving heat source of Pan and Cai,2018 [12].The reason behind the selection of this study [12] is that it clearly displayed physical movement of heat source, compared to others studies [14,15,16]. However, it should be noted that only the modelling methods (physical rotation or movement of heat source along friction surface in [12], and circumferentially the same heat application in axisymmetric) are taken as a modelling. But, the disc geometries, thermophysical material properties (Table 2) or any other input parameter is consistent for all the three models. Furthermore, the same braking power calculations are used for all modellings. ANSYS parametric design language (APDL) is utilized to code commands and operators in defining spatial and timely varying heat input and boundary conditions.

To attain this objective, the paper is organized under five sections. Section two introduced disc input parameters, geometric modelling, boundary conditions and heat input calculations. And the NAMHS and traditional models are summarized in the third section. And then, results and discussions accompanied by validation of the model are elaborated in section four. Finally, conclusions are drawn, and recommendations are set in the fifth section.

2 Input Parameters, Heat Source and Boundary Condition Calculations

2.1 Disc Specifications and Its FE Modelling

Input parameters for geometric modelling and thermal calculations are taken from trailer bogie of Addis Ababa Light Rail Transit (AALRT), manufactured by CNR-CRC (Chinese Railway Corporation). It is a 70% low floor vehicle with grey cast iron disc material having maximum running speed of 65 kph on 18 km East to West journey. Each journey has 21 stations (braking point) and fifteen trips per day. As a result, disc brake is vulnerable to thermal variations as high as 315 times per a day. All disc and pad geometries are measured by calliper and used for modelling in ANSYS APDL (Table 1).

Table 1 70% Low floor light rail vehicle braking conditions and disc geometry parameters (AALRT)

Besides, thermo-physical material properties for grey cast iron disc consider non-linearity due to temperature for thermal conductivity and specific heat capacity, by Eqs. (1 ) and (2), respectively, where T is the temperature (°C) [18]. In contrast, asbestos friction material on the rubber–resin binder is used as pads material, in which physical properties are highlighted in Table 2 [13].

$$k = 50 - 0.0178T \;{\text{W/m}}\;^\circ {\text{C}}$$
(1)
$$c = 506.5 + 0.377T \;{\text{J/kg}}\;^\circ {\text{C}}$$
(2)
Table 2 Disc and pad material input parameters [13, 18], AALRT

The disc brake type under investigation is split type similar with [10], having six bolt holes 60° far from one another at radial distance of 141 mm. On its bottom side, 48 cooling fins or vanes are used to dissipate heat (Fig. 1a-d). Two reinforcement structures with 5 mm thickness are added around each hole vanes and collectively 12(Fig. 1d). Due to axial symmetry, only one-half section is considered for simulation (Fig. 1a-d). In contrast, pad geometry is made up of one groove width 5mm width running radially and contact area of 0.01200 m2 with 50° cover angle (Fig. 2c).

Fig. 1
figure 1

NAMHS heat flux and convection applications, at different braking times: a early, b middle, c late, d middle for bottom disc convection

Fig. 2
figure 2

Heat flux applications at mid-braking times (a axis symmetric, b Pan and Cai [12]), c pad cover angle geometry, and d modified disc geometry

This original disc is a little bit modified in geometry by removing holes and reinforcement structures (Fig. 2d). Because, the impacts of these holes and reinforcement structures on spatial thermal result variations, such as temperature, stress, and thermal fatigue are not insignificant. Even, they might dominate other factors raising spatial temperature variation, such as heat flux and convection on the friction surface [19]. As a result, non-symmetrical moving heat source is investigated for both original and modified disc geometry in this study.

Fig. 3
figure 3

a Input heat flux and convection heat dissipation, b increment-independent test (Δθ)

SOLID90 is used for the meshing and creating finite elements. It has a temperature degree of freedom (an independent variable required to define temperature), and applied to two shapes: hexahedron shape of 20 nodes (fins and enforcing structures) and tetrahedral shape of 10 nodes (disc without fins and enforcing structures) [20]. Although disc and fins are meshed differently, they have common nodes at their interface. It is well suited to model curved boundaries. Accuracy and computation feasibility are conducted by mesh convergence test by varying the number of elements from 11443 to 172715. Elements on the surface are more refined than those far away from it, for each test. Finally, the solid disc is discretized into 87527 elements (small pieces) with 143628nodes are selected to proceed to next simulation, with computation time of 6 hours on Intel (R) Xeon(R) silver 2.1 GHz, 32 GB RAM computer, because the variation in temperature between the finest (172715 elements) and the selected mesh (87527 elements) is less than 0.6%.

2.2 Heat Input and Boundary Condition Calculations

In actual braking process, all modes of heat transfer (conduction, convection and radiation) are involved. However, radiation was shown limited contribution in many studies [15]. Therefore, only conduction and convection are focused in this modelling. By description, the retardation power equals to the negative of the work done by the friction forces per unit time in the pad–disc interfaces for the whole brakes. By considering constant deceleration, the amount of kinetic energy predicted for one disc brake is comparable to the energy balance written below (Eq. 3) [21]:

$$c\frac{1}{2}Mv_{o}^{2} = \mathop \int \limits_{0}^{{t_{s} }} P\left( t \right){\text{d}}t = 2F_{{{\text{disc}}}} \mathop \int \limits_{0}^{{t_{s} }} v_{{{\text{disc}}}} \left( t \right){\text{d}}t$$
(3)

where

c is the percentage of weight distribution on single disc calculated as follows

$$c = \frac{{\text{mass on single disc}}}{{\text{total mass of vehicle}}}*100 = \frac{M/8}{M} = 12.5\%$$

The braking (friction) force and corresponding instantaneous heat flux radial and temporal variation on the single pad–disc interface are expressed in Eq. 4 and 5, respectively [21]. It is noted that pad–disc interface area (Ap) is used for non-axisymmetric moving heat input [12] and NAMHS. Similar empirical formula is used for axisymmetric disc thermal modelling, except that disc contact surface area (Ad) and mean effective radius (rm) are used instead (Eq. 6). Heat source spatial variation in [12] was constant at pad–disc interface. Moreover, the dependence of heat input on radial variation in disc surface is seen for NAMHS (Eq. 5). However, the tangential variation in heat input requires new modelling, because the empirical formula couldn’t evaluate it.

$$F_{{{\text{disc}}}} = \frac{{c\frac{1}{2}Mv_{o}^{2} }}{{2\frac{{R_{d} }}{{R_{w} }}\left( {v_{o} t_{s} - 0.5at_{s}^{2} } \right)}}$$
(4)
$$Q\left( {r,t} \right) = \frac{{\gamma F_{{{\text{disc}}}} v_{{{\text{disc}}}} \left( {r,t} \right)}}{{A_{p} }} = \frac{{\gamma F_{{{\text{disc}}}} \omega \left( t \right)r}}{{A_{p} }}$$
(5)
$$Q\left( t \right) = \frac{{\gamma F_{{{\text{disc}}}} \omega \left( t \right)r_{m} }}{{A_{d} }}$$
(6)

where heat partition coefficient [22], mean effective radius [13] and full disc friction surface area [23] are expressed according to Eqs. 7, 8 and 9, respectively.

$$\gamma = \frac{1}{{1 + \left( {\frac{{\rho_{d} c_{d} K_{d} }}{{\rho_{p} c_{p} K_{p} }}} \right)^{1/2} }}$$
(7)
$$r_{m} = \frac{2}{3}\frac{{\left( {R_{p}^{3} - r_{p}^{3} } \right)}}{{R_{p}^{2} - r_{p}^{2} }}$$
(8)
$$A_{d} = \pi \left( {R^{2} - r^{2} } \right)$$
(9)

Besides, disc brake is also exposed to convection boundary conditions at rubbing surface(hs) and vanes or fins(hf). Time-dependent coefficient of heat transfer convection is approximated for laminar (Re < 2.4 × 105) and turbulent (Re > 2.4 × 105) flow on friction surface, based on Eqs. 10 and 11, respectively [12]. But it should be noted that rubbing surface convection is applied only in NAMHS, as it was not applied in traditional models.

$$h_{s} = 0.7\left( {\frac{{k_{a} }}{{D_{d} }}} \right){\text{Re}}^{0.55}$$
(10)
$$h_{s} = 0.04\left( {\frac{{k_{a} }}{{D_{d} }}} \right){\text{Re}}^{0.8}$$
(11)
$${\text{Re}} = \frac{{v\left( t \right)_{{{\text{disc}}}} R_{d} }}{\nu }$$
(12)

However, for ventilated side of disc brake, timely varying convection coefficient is estimated for turbulent (Re > 104) and laminar(Re < 104) flow, according to Eqs. 13 and 14, respectively [15, 18]. Since vanes cross section varies throughout their length, averaged hydraulic diameter (dh) and cross-sectional area (Ac) are determined at mid-length of vanes. Air kinematic viscosity and thermal conductivity are taken at 1 atm and 25 °C. Moreover, it must be noted that for all convection calculations made, Prandtl number is taken a constant value of 0.7.

$$h_{f} = 0.023\left[ {1 + \left( {\frac{{d_{h} }}{l}} \right)^{0.67} } \right]{\text{Re}}^{0.8} {\text{Pr}}^{0.33} \frac{{k_{a} }}{{d_{h} }}$$
(13)
$$h_{f} = 1.86\left( {{\text{RePr}}} \right)^{1/3} \left( {\frac{{d_{h} }}{l}} \right)^{0.33} \frac{{k_{a} }}{{d_{h} }}$$
(14)

where

$${\text{Re}}_{{d_{h} }} = \frac{{d_{h} }}{\nu }V_{{{\text{average}}}}$$
(15)
$$d_{h} = \frac{{4A_{c} }}{{2 \left( {\delta_{f} + w} \right)}}$$
(16)
$$v_{{{\text{avg}}}} = \frac{{v_{{{\text{input}}}} + v_{{{\text{output}}}} }}{2}$$
(17)
$$v_{in} = 0.0158N\left( {D_{d}^{2} - d_{d}^{2} } \right)^{1/2}$$
(18)
$$v_{{{\text{out}}}} = v_{in} \frac{{A_{{{\text{in}}}} }}{{A_{{{\text{out}}}} }}$$
(19)

As a summary, calculated convective heat dissipations, heat input and train speed are displayed as a function time (Fig. 3a). Convective coefficient is the same for all modelling, except their application area at the friction surface. Furthermore, heat input is five times higher in non-axisymmetric, because pad–disc interface area and full disc surface area are used for heat flux calculations, respectively, in non-axisymmetric [12] and axisymmetric modelling. In addition, heat flux is steadily raised up to 4s, as hydraulic pressure is fully operated after 4s [24].

3 Identification of Heat Source Spatial Variation and its Modelling in Non-Axisymmetric

3.1 Moving Heat Source Leading and Trailing Edge

Disc is the rotating component, while pads are the stationary. For this type of simulation, the leading edge of a brake pad is the one that first touching the rotor as it rotates, while the trailing edge touches the same point at last [25]. For decelerating type of braking, the heat absorbed by specific point on the disc from the leading edge of the pad is expected to be higher than the heat absorbed by the same point of the disc from the trailing edge because the angular speed of the disc is higher when it enters the leading edge of the pad, and is lower when it leaves trailing edge of the pads. As a result, the leading edge of the pad could be exposed to higher temperature. This argument is also consistent with other researchers who found higher temperature at the leading edge of the pads [26]. This observation could guide us to hypothesize the presence of heat input spatial variation at pad–disc interface.

In contrast, the pads and the disc are assumed, respectively, rotating and stationary component of braking system in moving heat source modelling [13]. For clockwise rotation of the disc brake (assumed in this study), heat source rotates counterclockwise. Rotating heat source geometric size is assumed equivalent to pad cover angle [12]. In addition, heat source angular speed is assumed equivalent to disc rotational speed. For single load step, it is assumed that heat source starts to heat at specific angular displacement on the disc and continues to heat the remaining friction surface until pad cover angle is completed (one load step). This assumption is based on aforementioned disc angular velocity variation between the leading and trailing edges, for stationary pad and rotating disc. For example, for a pad cover angle distance (50° in this case), if it starts heating at 0°, it will continue to heat until it reaches the other side of the pad (50°). Thus, friction surface that receives heat earlier (0°) is the leading edge, and the one that heated later will be the trailing edge. As a result, aforementioned terminologies in rotating disc and stationary pads are interchanged.

3.2 Identifying Spatial Variation in Heat Input Between Leading and Trailing Edge of the Pad

Firstly, while heat source is rotating, its speed variation due to deceleration reduces the amount of heat input from leading to trailing edge. Thus, particular node on disc surface does not receive the same amount of heat from all contacting nodes on pad circumference. This point of view is consistent with the other investigations who reported contact pressure drop from the leading to trailing edge, for rotating disc and stationary pads [27, 28]. Thus, it is hypothesized that heat flux could vary in the same way across pad–disc interface, as it could be calculated from contact pressure.

Secondly, braking power and corresponding heat flux are directly proportional to disc radius (Eq. 5), in spite of its constant assumption in traditional modellings. Thirdly, pads are accompanied with grooves of different layouts for composite type: radial [28], tangential [24] or both [2]. Besides, it might have different geometry and the arrangement for sintered type [29]. This diversified design is reported to be the main factor which could affect spatial temperature distribution [2]. The last but not the least factor that can affect heat flux spatial variation is hydraulic pressure in early braking time. Its increase in a short period of time requires further spatial and time investigations on heat flux and temperature. Surprisingly, the simultaneous effects of these factors on spatial and temporal temperature distribution have not been comprehensively examined in moving heat source.

3.3 Non-Axisymmetric Modelling of Moving Heat Input

To consider space and time variation in heat input and boundary condition, APDL programming model is displayed (Fig. 4). The first part (the left side) is applied for calculating and storing them as array and scalar parameters for later use. Once the calculation is completed, the second part (the right side) is applied for calling and applying them on the required location at the required time. And the model is executed through APDL parameters (scalars and arrays) similar to FORTRAN written commands.

Fig. 4
figure 4

Flow diagram of heat flux and convection calculations (left) and applications (right)

Prior to calculation, train kinematic inputs (Vo, a, tb), disc geometric dimensions (Table 1) and material properties (Table 2) are defined as scalars in ANSYS database, from which heat source total rotation angle (θtot in °) can be calculated. Then circumferential heat variation is evaluated and updated based on small heat distance increment Δθ and is also saved in database as arrays. Thus, the circumferential variation in heat source is assumed to occur every Δθ, and the heat source is assumed to start rotational displacement from θo = 0°. Hence, nodes within the same Δθ are assumed to have the same angular speed (rpm), but differ in tangential speed (m/s), due to their radius. Consequently, they will have the same and varying heat flux, respectively, along circumferential and radial direction (Eq. 20). And then, θtot can be converted to total number of heat source circumferential variation (n).

Nonetheless, the exact selection of Δθ seems challenging as it can dramatically affects the circumferential heat variation (0°–360°). In previous moving heat source modelling, it was assumed to be 360° for each load step [14,15,16]. In contrary, this assumption is a little bit modified in [12], since the authors reduced circumferential thermal variation to pad cover angle (1.152 rad, corresponding to 66°). Despite that, still heat variation is expected within pad–disc interface, as it is aforementioned. To select appropriate heat distance increment, increment independent test is carried out for 25° and 10°.This is by assuming, respectively, two heat source variation and five heat source variation areas, within pad–disc interface. Besides, this test is also conducted further for 5° and 2.5° in the same way. Consequently, circumferential heat variation areas are raised further to 10 and 20, respectively. The corresponding circumferential temperature variation is plotted on mid-braking time at mean radius (Fig. 3b). It is inferred that no more increase in number of varying heat area is required beyond 2.5°, as similar result is displayed for 5° also. As a result, pad cover angle (50°) is partitioned into 20 different heat source areas, decreasing or increasing in heat and temperature from leading to trailing edge, based on braking time.

Nonetheless, radial variation in heat is performed in different manner during heat source application. For every angular increment θi, time increment (ti), friction force (Fdisc) and angular velocity ω(ti) could be calculated and stored as array parameters, by the formulae expressed in the flow diagram (Fig. 4). And then, at every ti, friction force (Eq. 4) and convection coefficients (Eqs. 1014) are also calculated and saved in ANSYS database for later calling. While the calculation is conducting, every angle θi is checked if it is equal or not equal to θtot. If θi is not equal to θtot, the process will continue for another θi. Or else, the calculation is completed for θi = θtot. As a result, heat input is expected to drop on every θi (i varied from 0 to n) until θi=n= θtot, due to deceleration.

Once the θitot is satisfied for calculations, the second part of the modelling is initiated for heat source calculation and application, and boundary condition application. This is performed according to the ti and its corresponding convection coefficient. One load step is assumed to have locally varying heat source application on pad cover angle domain [12], in addition to convection coefficient on the remaining surface nodes (hs) and fin nodes (hf). Surface nodes are selected based on heat variation distance (from θi-1 to θi, where θi - θi-1 = 2.5°) extracted from database. Nodes at location of radial pad groove (insulated) are excluded from these selections, and there are FORTRAN commands written to do so in APDL [30]. And, at bolt holes also no nodes are selected (for original disc), even when it is ordered by the commands. For specific location of θi-1 to θi, heat source is calculated from arrays extracted within database (Fdisc, ω (ti), Ap) and from inner radial location of selected nodes (Eq. 25).

Finally, the linearly varying heat source with radial location is calculated and applied from heat gradient slope [31], for every θi. And this loop will continue until pad cover angle distance is completed (θi = multiple of 50°, or i = multiple of 20) for θi < θn. Once θi assumes pad cover angle multiple (at i multiple of 20), the surface convection (hs) and fin convections (hf) are called from database and applied on their respective locations free from heat flux. For the first load step (i = 1–20, or θi = 0–50°), initial condition (25 °C) is used, while for the remaining load steps, nodal temperature at previous load step (\(t_{i} - t_{i - 20}\)) is used as initial condition. Moreover, the solution time for one load step (tls) is taken as the time required for the heat source rotation to cover a pad cover angle rotation distance. It is unique for each load step and increase from early to late braking times. Because, higher time is required to complete one pad cover angle in late than early braking time.

Following the solution of each load step, every multiple of pad cover angle θi is checked if it is equal or not equal to θn. If θi is not equal to θn, the process will continue for another load step. Otherwise, the solution is completed for θi = θn. However, this condition is very rare case to occur, and expected if and only if θn will be multiple of pad cover angle, or n will be multiple of 20.

If θn is not multiple of pad cover angle, the same procedure is followed up to the highest multiple of pad cover angle θm(where m is the highest multiple number of 20), which will be close to θn by θn−θm<pad cover angle. In another way, the highest multiple number of 20 (i = m) will be close to n by n-m<20. After θm, the heat will be applied in the same trends for θi = θm to θn. But since θn is not multiple of pad cover angle, the solution for the final load step will be executed in different loop (θi = θn satisfied), with heat source application area less than a pad cover angle distance.

For this particular braking conditions and disc dimensions, θtot = θn=27387.5°, n = 10956, m = 10940, θm = 27350°, and total load step is 548. For the first 547 load steps, heat is applied on full pad cover angle, and for the last load step, it is applied on 37.5°. The position occupied by heat flux at a given time is displaced by convection at the next time, and vice versa. As a result, friction surface is heated and cooled on small angular increments periodically, with heat flux and convective heat dissipation, respectively (Eq. 22 and Fig. 1a-c).

To observe the impact of the NAMHS on spatial heat flux variation, three braking times are randomly taken: early (0.089s), middle (8.151s) and late (15.783s) braking times (Fig. 1). The influence of radius on heat flux variation is observed at all braking times. Nonetheless, the impact of circumferential heat flux variation is prevalent with increase in braking time (Fig. 1c) because heat source displacement comparable to pad cover angle distance is covered within a very short period of time during early and mid braking time than late braking time (due to deceleration). This means circumferential speed variation among pad incremental distance will have insignificant variation in early braking and significant during late braking. As a result, the magnitude of heat source within pad cover angle decreased from leading to the trailing edge, but not seen abrupt in mid-braking time (Fig. 1b), as in late time.

$$Q\left( {r_{j} ,\theta_{i} ,t_{i} } \right) = \frac{{\gamma F_{{{\text{disc}}}} \omega \left( {\theta_{i} ,t_{i} } \right)r_{j} }}{{A_{p} }}$$
(20)

Even though pad leading edge is expected to have relatively higher heat flux, it is surprising to observe higher heat at pad trailing edge, contrary to expectations (Fig. 1a). This rather contradictory result might be due to significant increase in hydraulic pressure with time, resulting in an increase in heat flux from pad leading to the trailing edge (Fig. 3a). This dominates the reduction in heat flux from leading to trailing edge of the pad due to deceleration (reduction in speed with time). For the late braking however, heat flux decreases from leading to the trailing edge as expected, as no hydraulic pressure interferes the heat flux (Fig. 1c). In this modelling, consideration of wear heat generation is neglected as its effect is low compared to friction heat generation [11]. Besides, NAMHS model is conducted on composite pad with one radial pad groove, though its application on sintered and other composite pads (one tangential and two radial grooves) is tested and seen successful (Fig. 5).

Fig. 5
figure 5

Flexibility of NAMHS to accommodate various pad design: a sintered, b composite

Heat input model is also displayed for traditional models: axisymmetric (Fig. 2a) and Pan and Cai [12] (Fig. 2b). Heat flux is spatially constant on all friction surface, in the former one. Meanwhile, heat flux is independent of space at pad–disc interface in the latter. Heat flux is calculated at mean effective radius and varied timely in both traditional modellings. Besides, for Pan and Cai [12], friction surface area is shared for heat flux only, and the remaining surface area is insulated. Heat dissipative convection is applied on external radius (Eq. 23) and vane surface areas (Eq. 24), while fin bottom surfaces and internal radius areas are insulated for all modelling (Eq. 25 and Fig. 1d).

3.4 Friction Heating Boundary Value Problem

Non-axisymmetric boundary value problem of parabolic heat conduction equation is presented for transient analysis T (r, θ, y, t) in cylindrical coordinates of the disc brake, where, r, θ, y and t represent radial, circumferential, axial and time variables (Eq. 20) [13]. In addition, other boundary conditions, including heat conduction and convection, are also highlighted (Eq. 21-24). It should be noted that disc surface is taken as reference (0) from where other disc thickness is measured negatively along axial direction. These equations are written for one load step time (tls = ti+20-ti) and repeated for the remaining load steps (Fig. 4). And finally, Eq. 5 is a little bit modified to account spatial variation in heat (Eq. 20), where i and j are circumferential and radial variation in heat, respectively.

$$\frac{{\partial^{2} T}}{{\partial^{2} r}} + \frac{1}{r}\frac{\partial T}{{\partial r}} + \frac{1}{{r^{2} }}\frac{{\partial^{2} T}}{{\partial^{2} \theta }} + \frac{{\partial^{2} T}}{{\partial^{2} y}} = \frac{1}{{K_{d} }}\frac{\partial T}{{\partial t}} {\text{for}}\;\left\{ {\begin{array}{*{20}c} {r_{d} < r < R_{d} , \left| {\theta_{i} } \right| \le \pi ,} \\ { - \delta_{d} < y < 0} \\ {t_{i} < t < t_{i + 20} } \\ \end{array} } \right.$$
(21)
$$\left. {K_{d} \frac{\partial T}{{\partial z}}} \right|_{{y = 0^{ - } }} = \left\{ {\begin{array}{*{20}c} {Q\left( {r_{j} ,\theta_{i} ,t} \right), r_{p} < r_{j} < R_{p} ,\left| {\theta_{i} } \right| \le \theta_{o} ,t_{i} < t < t_{i + 20} } \\ {} \\ {h_{s} \left( {T_{a} - T\left( {R_{d} ,\theta_{i} ,0,t} \right)} \right),\left\{ {r_{d} < r < r_{p} ,\left| {\theta_{i} } \right| \le \pi } \right\} \vee } \\ {\left\{ { r_{p} < r < R_{p} ,\left| {\theta_{i} } \right| \ge \theta_{o} } \right\} \vee } \\ { \left\{ {R_{p} < r < R_{d} ,\left| {\theta_{i} } \right| \le \pi } \right\},t_{i} < t < t_{i + 20} } \\ \end{array} } \right.$$
(22)
$$\left. {K_{d} \frac{\partial T}{{\partial r}}} \right|_{{\left\{ {r = R_{d} } \right.}} = h_{s} \left[ {T_{a} - T\left( {R_{d} ,\theta_{i} ,y,t} \right)} \right],\left| {\theta_{i} } \right| \le \pi , - \delta_{d} < y < 0,t_{i} < t < t_{i + 20}$$
(23)
$$\left. {K_{d} \frac{\partial T}{{\partial r}}} \right|_{{\left\{ { - \left( {\delta_{f} + \delta_{d} } \right) < y < - \delta_{f} } \right.}} = h_{f} \left[ {T_{a} - T\left( {r,\theta_{i} ,y,t} \right)} \right],\left| {\theta_{i} } \right| \le \pi ,t_{i} < t < t_{i + 20}$$
(24)
$$\left. {\frac{\partial T}{{\partial r}}} \right|_{{\left\{ {r = r_{d} } \right.}} = \left. {\frac{\partial T}{{\partial r}}} \right|_{{\left\{ {y = - \left( {\delta_{f} + \delta_{d} } \right)} \right.}} = 0, r_{d} < r < R_{d} ,\left| {\theta_{i} } \right| \le \pi ,t_{i} < t < t_{i + 20} ,$$
(25)

4 Results and Discussions

4.1 Radial Transient Temperature Evaluation

For better result analysis, it would be sensible to take six radially varying points(but circumferentially 90°) on the surface: inner radius(R1 = 113 mm), inner pad radius (R2 = 122.5 mm), mean or effective radius(R3 = 146.5 mm), transitional radius(Rt = R4 = 154 mm), outer pad radius (R5 = 170 mm) and outer disc radius(R6 = 180 mm). The results of these points are illustrated for actual disc geometry on the left side of Fig. 6, while the right side of the figure displays the variation between the NAMHS and traditional models (axisymmetric and [12]). R3 is taken from Eq. 8, position at which heat flux or breaking power is calculated for both traditional models (Eq. 6). Besides, it is presented as the best position to compare thermal results, because all models will have the same radial heat flux (Eqs. 5 and 6) here.

Fig. 6
figure 6

Transient temperature at different radii: a disc inner (R1), b pad inner (R2), c mean (R3), d Transitional (R4), e pad outer (R5), f disc outer (R6)

On the contrary, Rt is the NAMHS finding from this investigation, as the NAMHS temperature starts to surpass Pan and Cai, [12] model (Fig. 6d). Despite registration of lower temperature displayed by the NAMHS from R1 to R3, this trend is not shown beyond Rt. At this point, both moving heat source models display equivalent results, except at early (0–4s) and late braking times. This might be associated with NAMHS’s hydraulic pressure increase and prevalence of surface convective coefficient on former and later braking times, respectively. Beyond Rt, the NAMHS model is shown to surpass [12], due to higher radius used in heat source calculation.

The temperature variation results at the mean effective radius (R3) demonstrated that both traditional models exaggerated surface temperature distribution as high as 0%–60% and 0–5%, respectively, for axisymmetric and moving heat source of Pan and Cai [12] (Fig. 6c). It should be noted that temperature variation calculated here is in absolute value for all radii. This much variation in axisymmetric modelling is likely attributed to huge difference in surface area shared for heat flux application. In axisymmetric, 360° surface area is exposed to heat flux, meanwhile only pad cover angle in moving heat sources. Overestimation of the results by axisymmetric is rather surprising result, as higher heat flux value is displayed for moving heat sources (Fig. 3a). The prediction of conservative fatigue life in [4] might seem due to this axisymmetric modelling. Therefore, the findings reported here suggests that the size of area on which heat applied would have a remarkable impact on thermal results, compared to any other factors affecting heat flux calculation and application.

Though the moving heat sources are observed to have similar results, still the presented value is significant (5%) in late braking time (Fig. 6c). This could be further enlightened by insulated pad groove and circumferential variation in heat flux (due to hydraulic pressure and train deceleration). Furthermore, exposure of friction surface outside the pad–disc interface to convection in NAMHS could be associated with this variation. Particularly, the rising impact of circumferential heat flux and surface convection coefficient on the late time appears responsible for gradual rise of their difference. Unfortunately, this finding is somewhat disappointing, as their share or their individual contribution on result variation is not revealed. Furthermore, this outcome is contrary to that of Pan and Cai [12] who found almost the same temperature results, for axisymmetric and moving heat source. The insulated surface area outside the pad–disc interface could explain this discrepancy, as it could raise temperature in moving heat source up to axisymmetric temperature.

The other most striking result to emerge from NAMHS is the contribution of radius in affecting temperature. This is apparently seen where pad extreme radii are used in NAMHS (Fig. 6b and e). It is displayed that the variation in NAMHS from [12] model is dropped from 10% to 0, as radius rises from R2 to Rt. The rise of heat flux with radius in NAMHS explains this observation. In contrast, from Rt to pad external radius, its variation is also shown to rise with radius, due to the same reasoning. But its variation is not substantial as in inner radius (Fig. 6e). This unexpected result might be associated with higher contribution of circumferential heat variation in reducing temperature, at outer radius other than inner radius (Figs. 1a, c and 7). Consequently, the effect of radius in rising temperature is minimized. These results suggest the NAMHS model capability in revealing the combined effect of radial and tangential variation in heat input on spatial variation in temperature.

Fig. 7
figure 7

Circumferential temperature distribution at various radius on late braking time

Particularly, the effect of radius consideration on temperature is obviously seen when equal radius is used for all models in heat flux calculation (Fig. 6c). The variation is massively reduced with gradually increased (0–5%). Hence, it could be inferred that the NAMHS is the best in revealing radial variation in temperature. Furthermore, this can highlight the convincing evidence for the argument that radius takes significant share, if not the only one, for temperature variation. In contrast, a possible explanation for 5% variation might give a clue or tentative evidence for presence of addition factors to be considered. These include circumferential heat flux variation, insulated pad groove and friction surface convection in NAMHS. Hence, this finding has drawn our attention to the significance of considering spatial variation in heat source in modelling disc brake, which couldn’t have been supported in traditional models. So, the findings reported here suggest the radial distance as the second major factor to be considered, in causing spatial temperature variation.

However, the most surprising result observed on NAMHS model is the registration of the lowest temperature for R value ranging from R2 up to Rt, in spite of higher R value used in calculating heat flux (Eq. 5). This result also offered another additional evidence on contribution of circumferential heat input variation (deceleration), friction surface convection and consideration of pad groove in thermal modelling. Collectively, these factors seen to overcome heat rise by 7.5 mm increase in radius from Rm to Rt. And hence, these findings suggested the third factors to be considered in spatial temperature variation, as deceleration, pad groove and friction surface convection.

Regardless of this investigation, some recently reviewed studies by the authors reported almost half of the disc brake thermal studies (49%) regrettably disregarded heat flux radial variation [5]. Yet, few studies reported its consideration, but their heat input calculation methods were different from the current investigation: uncoupled thermomechanical [3, 32] and coupled thermomechanical [10, 33]. It is obvious to consider spatial effect in such modellings, though convergence problem in these methods enforced us to consider different modelling.

For all presented modellings, temperature is displayed higher at middle radii (Fig. 6c and d) than pad extreme radius. Besides, disc extreme radii are observed to have the lowest temperature. This is not surprising, as these points are ‘far’ from heat source area. In accordance with the presented results, other studies also confirmed maximum temperature location at mid-surface of disc [24, 34]. On the other hand, this finding is contrary to previous studies who reported maximum temperature at outer radius of disc [23]; meanwhile, others reported extreme radius (inner and outer) as location of maximum temperature [22, 35]. These variations in conclusion seem emerge from using different pad types (cover angle, size, groove number and layout), disc types and design.

4.2 Circumferential Surface Temperature Variations

To exclude the effect of disc geometry on circumferential temperature fluctuation, the result is displayed for both modified (uniform disc thickness—Fig. 8a–c) and original (12 hole reinforcing structures around 6 holes–Fig. 8d–f) disc geometries. These results are reported for three braking times at mean radius for all modellings.

Fig. 8
figure 8

Circumferential temperature distribution for modified (ac) and original (df) disc at mean radius and different braking time: early (a,d), middle (b,e), late (c,f)

Six sinusoidally fluctuated temperatures are displayed for original disc geometry in a way that, steadily raised and then dropped, for all modellings (Fig. 8d–f, Fig. 9d–f). These results are directly related to two reinforcement structures added around each bolt hole between two consecutive fins, resulting increase in disc thickness, every 60° (Fig. 1d). This increase will lead to drop of temperature. Consequently, surface temperature variation is presented 30 °C and 42 °C by axisymmetric model on late and mid-braking times, respectively (Fig. 8e and f). Although such types of results are not concern of this study, they are useful to observe disc geometry role in complicating temperature fluctuations in moving heat sources.

Fig. 9
figure 9

Temperature plot for modified disc (NAMHS: a early, b mid, c late braking times) and original disc (d axisymmetric late, e. [12] late and f. NAMHS late braking times )

Further analysis of all braking time in moving heat sources models reveals slight rise of each peak or interface (Fig. 8), as we move further from lower to the higher angle. A possible explanation for this might be the disappearance of high temperature positions as time passed, owing to heat dissipation by either conduction [12] or both conduction and convections(NAMHS). As a result, interfaces close to pad cover angle are highlighted higher in temperature, compared to those far from it.

Since this impact of disc geometry is shown dominant, the modified geometry seems preferred to investigate the individual effect of NAMHS model at mean radius, just to investigate circumferential temperature fluctuation. Hence, two specific positions are selected on mid (28° and 263°)- and late (28° and 302°) braking times. The highest surface temperature difference is revealed by the NAMHS as 49 °C, and14 °C at mid- and late braking times, respectively. Likewise, [12] displays 42 °C and 10 °C at mid- and late braking times, respectively. Nonetheless, 0 °C and 2 °C surface temperature variations are revealed by axisymmetric model at mid- and late braking time (Fig. 8b–c). This late and early braking time variation in axisymmetric model might seem due to free mesh variation every quarter of circle and conduction along fins, respectively. Besides, the capability of all models in accommodating spatial temperature variation is also obviously revealed in Fig. 9d–f. Spatial variation is seen better displayed in NAMHS than both traditional models, and [12] model is illustrated better than axisymmetric.

During early braking time, surface temperature differences are seen clearly for both moving heat sources only, though their magnitude are not seen significant (1 °C) (Fig. 8a and d, Fig. 9a). This might be due to less temperature development in early braking time. Moreover, the maximum angular velocity attained here could lead to inconsequential surface temperature variation. In spite of these, both moving heat sources displayed gradual increase in temperature in the first half (0–180°), followed by ‘dramatic’ increase as it approached pad–disc interface area (250°–300°). Beyond this interface, they have shown ‘massive’ drop of temperature

Despite similarity in the trends of moving heat sources, their careful observation identifies their differences (Fig. 8a and d). In [12], stepwise rise in temperature is between consecutive interfaces, but not varied with time within interfaces. However, in the NAMHS, continuous rise in temperature (including, at pad–disc interface) is shown instead. One reason that could explain this observation is the circumferential increase in heat flux due to hydraulic pressure in the NAMHS model (Figs. 1a and 9a). Meanwhile, its constant assumption in traditional is related to constant temperature at the interfaces (Fig. 2a and b). One of the issues that emerged from these findings is consideration of hydraulic pressure as the fourth factor that could affect temperature evolution.

Another remarkable effect of heat flux circumferential variation in NAMHS is prevailed in terms of circumferential temperature drop, due to deceleration. At contact interface of 300°–350°, temperature variation is reported 7 °C (Figs. 8 and 9c, f); meanwhile, the traditional models couldn’t predict. Nonetheless, such trends of temperature are continued to exist at mid- braking times (Fig. 8b, e). It is attributed to relatively higher rotational speed of heat source, resulting insignificant variation.

Thus, the result presented in the NAMHS revealed encouraging outcome, as it imitated linear rise of hydraulic power and shown the best surface temperature difference, compared to traditional models. And these results uncovered two remarkable findings. Firstly, the NAMHS is successfully reported spatial variation in temperature. Secondly, disc geometry is accounted for significant amount of surface temperature difference in all models (Fig. 9d–f).

Two specific positions appear a clear gap (free from heat) in between two heat source areas in Figs. 8 and 9. Particularly, on early braking time of both figures, they are seen obviously. Two possible explanations are put forwarded: Firstly, insulation on the pad groove which might reduce temperature. Secondly, migration of surface heat from pad extreme border to neighbouring nodes, either by convection and conduction(NAMHS) or conduction only [12]. This seems jump in heat source rotation between consecutive pad cover angle, but not in actual loading. Furthermore, spatial surface temperature presented here (Fig. 9) is not observed similar with hot spots and hot bands in other studies [6], which might be due to negligence of thermoelastic contact wear between pad–disc interface, and variation in input parameters. But, compared to traditionally known models, it is believed that the NAMHS had successfully uncovered surface temperature spatial variations at any braking time.

NAMHS findings seem to be consistent with other research which found 50°C temperature difference at rubbing zones [7]. In contrast, another result reported 5% error between axisymmetric and non-axisymmetric thermal analysis [6]. This might appear due to higher initial velocity used (twice compared to this study). Because, the time required by heat source to complete one revolution is relatively small, resulting close result between axisymmetric and non-axisymmetric modelling. Thus, NAMHS model might seem recommended for low-speed trains, having frequent stop, like light rail transit. But this didn’t mean not applicable to high speed, rather, similar result would be expected with axisymmetric modelling.

Furthermore, the circumferential temperature at the interface 300°–350° is expected to decreasing from leading to trailing edge (Fig. 8c and f). But this factor is seen apparently at outer radius other than inner radius (Fig. 7). And, this didn’t mean no circumferential heat flux variation at inner radius. Rather, it could be explained by the impact of circumferential heat variation prevalence at outer radius than inner, as aforementioned in Sect. 4.1. On the contrary, interface temperature variation for original disc is depicted differently. This could be enlightened by two issues. Firstly, the proximity of the leading edge or first half of pad cover angle(301°–325°) to the bolt hole area at 300°, where low temperature is registered. Consequently, leading and trailing edge appeared to have similar temperature (Figs. 8f and 9f). Secondly, the inner radius is closer to bolt centre hole (by 28 mm) than outer radius, which is 39 mm far. Due to these issues, surprisingly the results seem to contradict the previous argument, resulting lower or similar temperature at the leading edge (Fig. 8f). Therefore, these results need to be interpreted with caution for original disc, as geometry could interfere temperature distribution. For either of the disc however, this observation of the higher temperature at the leading edge in NAMHS is consistent with previous investigation of temperature [26] and contact pressure [27, 28] (which could be the source of higher temperature) of pad–disc interface.

4.3 The Effect of Pad Groove in Circumferential Surface Temperature Variations

As we approached the interface, increasing trends of sharp drop and rise of temperature is observed throughout the braking time, due to pad groove(Figs. 7 and 8). Since temperature is not developed yet in early braking, its difference at pad groove compared to its neighbourhood nodes is seen insignificant (1 °C). In addition, at mid-braking time for mid-radius, relatively the influence of pad groove is also observed (Fig. 8b). Likewise, in late braking time, pad groove effect is illustrated influential as it approached the last peak (Figs. 7 and 8c,f). The maximum temperature variation between the groove location and its neighbouring nodes within the same pad cover angle is uncovered 44 °C and 15 °C at the mid- and late braking, respectively (Fig. 8b and c). Temperature variation due to groove is seen rising, as it approached the peak interface. This is associated with migration of high heat area to lower temperature area by conduction (to pad groove) and convection (to surrounding). Finally, temperature variation at peak pad cover angle has shown negative slope (change in temperature with in pad–disc interface to pad cover angle) of −0.122 °C/° (Fig. 8c); meanwhile, positive slope is demonstrated for early braking. Nevertheless, such results are not noticed in mid-braking times, as tangential heat variation is perceived during early and late braking times (Fig. 1b).

This effect of pad groove in spatial temperature reduction is in line with the findings of other studies [24]. Thus, further reduction in temperature is hypothesized as pad groove raised in number, and this is seen agree with other investigations [36], who found lower distortions and Mises stress for pad with groove. Hence, this result provides another support for the hypothesis that pad groove consideration could influence result variation, spatially. Consequently, the NAMHS model might help us to optimize pad design based on its numbers, shapes and location of grooves as in [2, 22, 29].

4.4 Temperature Investigation Through Disc Thickness and Validation of the Model

Further analysis by the NAMHS model shows zigzagged behaviour of temperature to have only 6.5mm depth from the surface (Fig. 10a). This cease of zigzagged behaviour is likely related to the delay of heat conduction along circumferential direction, as the disc thickness increased from the surface. This is in agreement with other studies who reported 5 mm depth [17]. Hence, the NAMHS provides the better chance to reveal spatial and temporal variation in temperature.

Fig. 10
figure 10

Temperature distribution a along disc thickness by NAMHS and b model validation

Finally, the NAMHS model is applied to experimentally and analytically studied solid disc brake [37]. All input experimental conditions or parameters (geometry, thermophysical material properties, braking times, braking powers) are extracted to APDL simulation from M. Fermér [37], and then the NAMHS model is applied. The surface temperatures at 232mm radial and random circumferential location is displayed (Fig. 10b). Though deviation from experiment and analytical study is seen during early and late braking, the modelling shows acceptable correlation. This deviation might be due to the limited number of thermocouples (only three) used in experimental investigation, while six thermocouples are required according to UIC 541-3 [24].

5 Conclusions

Non-axisymmetrical moving heat source (NAMHS) is developed for disc brake temperature analysis, and its success is evaluated with traditional models, within feasible computational time. NAMHS model is successfully implemented the major factors those affecting spatial temperature distributions, including friction surface exposed to convection and heat source, radial and circumferential heat flux variation, early braking hydraulic pressure and pad groove effects. The main findings of the model are summarized below:

  • Size and shape of area in which heat flux applied is shown remarkable effect on surface temperature variation

  • Consideration of radial distance showed surface temperature variation as high as 10% and 60% in [12] and axisymmetric, respectively.

  • Heat flux variation within pad cover angle, and partition of friction surface area to heat source and convection in NAMHS resulted maximum circumferential surface temperature variation as high as 49 °C at mean radius on mid-braking time

  • The consideration of pad groove at interface led to 44 °C difference of temperature within pad cover angle.

  • Temperature variation with in pad cover angle is shown maximum negative slope (−0.122 °C/°), positive slope and unnoticed slope at late, early and mid-braking times, respectively.

  • Reinforcement geometry around bolt hole contributed significant surface temperature variation

  • The effect of circumferential heat variation is prevalent at outer radius than inner radius

  • Hydraulic pressure in early braking time dominates the effect of deceleration on temperature within pad cover angle.

  • The model might predict close temperature with axisymmetric model for high-speed trains

In light of this remarkable findings, it would seem to provide convincing evidence for the capability of the NAMHS model in forecasting specific area of high temperature, within acceptable computation time. Furthermore, its results shed new light on extending how to implement different pad groove size, layout and shapes. Consequently, our work has led us to conclude that the model could be applied quite reliably in estimating pad optimization, thermal stress and fatigue life of disc brake. But, the generalizability of this model is subject to certain limitations, including lack of wear, CFD convection and experimental validation. Nevertheless, this study is the first step towards enhancing the observance of spatial and temporal thermal results on friction surface of disc brake, with acceptable validation. Using the developed modelling, further research into disc brake fatigue life estimation is in progress by the authors. Although it is recommended for trains with low speed with frequent stop (like light rail transit), it is also applicable for any braking types (emergency or service), disc types (solid or ventilated) and pad types(sintered or composite). Hence, we hope that our finding could influence disc brake manufacturers, researchers and maintenance personnel in disc brake damage investigation.