Analysis of Brake Disc for Regenerative Braking System

The assignment requires students to analyze a regenerative braking system, which is a critical component of hybrid electric vehicles. The student must calculate the average and maximum power required for the system, as well as the heat flux into the brake disc. Additionally, they must determine the temperature rise of the brake disc during a single stop and a fade stop. The assignment also requires students to research and discuss relevant mathematical models and simulations used in the analysis of regenerative braking systems.

BRAKING SYSTEM WITH MATLAB
Introduction
The goal of this exercise has been to create a matlab program to represent a reasonably accurate model of a
racing car travelling around an arbitrary circuit. This allows experimentation to take place with a variety of factors that are
likely to affect the performance.
A number of assumptions have been made thus far that are likely to affect the outcome, to a greater or lesser
degree:
i) The coefficient of friction of the tyre remains constant at all times
- in practice, this tails off as the load placed on the tyre increases, and hence brings into play the design of the suspension
in terms of the effects of weight transfer, and the possible benefits of adding downforce.
ii) The sample points on the circuit are close enough together to calculate the car’s performance statically at each
point. Acceleration between sample points is taken to be constant.
iii) The driver is capable of getting the maximum performance from the tyres during braking and cornering (or is able
to achieve a fixed percentage of the maximum possible forces).
iv) The circuit is completely flat – no hills or dips are present.
– in practice the presence of hills will affect the ability of the car to brake and accelerate, as
well as the grip available from the tyres for cornering.
Side Slip, Deceleration * Angular Velocity Matlab Code
close all
clear
m = 300; %mass of the vehicle
g = 9.81; %gravitation force
r = 0.32;
I = 1;
v0 = 38.88889; %mps of 140kmph velocity
w0 = 100;
desiredSlip = 0.2;
%normal force
FN = m*g;
%temporary Mb
Mb = 1500;
v(1) = v0;
w(1) = w0;
%controller
Kp = 1000;
P0 = -50;
%counter for vectors
vecPos = 1;
%step size
h = 0.01;
for t = 1:h:10

%calculate slip ratio
%slip ratio = the spin velocity divided by its actual world velocity.
%A slip ratio of -1 means full braking lock;
%a ratio of 0 means the tire is spinning at the exact same rate as the road is disappearing
below it.
%A slip ratio of 1 means it's spinning.
slipRatio(vecPos) = (v(vecPos) - w(vecPos)*r) / v(vecPos);
%adjust slip if it exceeds its limits
if(slipRatio(vecPos) > 1)
slipRatio(vecPos) = 1;
elseif(slipRatio(vecPos) < 0)
slipRatio(vecPos) = 0;
end
%calculate the error
slipError(vecPos) = abs(slipRatio(vecPos)-desiredSlip);
%adjust Mb
if(slipRatio(vecPos) > desiredSlip)
Mb = Mb - Kp*slipError(vecPos) + P0;
elseif(slipRatio(vecPos) < desiredSlip)
Mb = Mb + Kp*slipError(vecPos) + P0;
end
%calculate friction using Pacejka's Magic Formula
u = sin(1.9*atan(10*slipRatio(vecPos) - 0.97*(10*slipRatio(vecPos) -
atan(10*slipRatio(vecPos)))));
%u = 1.28*(1-exp(-23.99*slipRatio(vecPos)) - 0.52*slipRatio(vecPos));
%update wheel angular velocity
dw = (r*u*FN - Mb) / I;
w(vecPos+1) = w(vecPos) + h*dw;
%update car velocty
dv = -u*FN/m;
v(vecPos+1) = v(vecPos ) + dv*h;
%if car velocity is zero or lower, stop simulation
if(v(vecPos+1) < 0)
v(vecPos+1) = 0;
break
end
%count array
vecPos = vecPos + 1;
end
plot(v)

title('Deceleration')
xlabel('t') % x-axis label
ylabel('v') % y-axis label
figure
plot(w)
title('Angular Velocity')
xlabel('t') % x-axis label
ylabel('w') % y-axis label
figure
plot(slipRatio)
title('Side Slip')
xlabel('t') % x-axis label
ylabel('Slip Ratio') % y-axis label
%ends here
Braking Operations
clear, clc;
%% Conversion Factors
ft2Meters = 0.3048; %1ft = 0.3048m
Miles2Meters = 1609.344;
mph2ms = 1609.344/3600;
%% vehicle Constants
vehicleMass = 300; %kg
pusher_accel_rate_g = 2.4; %g - from linearly interpolated graph.
pusher_accelrate = pusher_accel_rate_g * 9.81; %m/s^2
track_distance = 5000*ft2Meters; %m (Assume 5000ft track)
pusher_distance = 1600*ft2Meters; %m (Assume 1600ft pusher accel phase)
setbraking_distance = 400*ft2Meters; %m
%setbraking_distance = ft2Meters * input('Set braking distance (feet) = '); %m
coast_distance = track_distance - pusher_distance - setbraking_distance; %m
%--Velocity
Vmax1 = 140 * mph2ms; %m/s - Max velocity after pusher releases (constant accel)
Vkmax = 140 * mph2ms; %m/s - Max velocity before hitting kantrowitz limit, full carbon shell
%% Time Calcs
%--On pusher (estimation!!)
timepusherend1 = Vmax1 / pusher_accelrate;
%--Coasting
coastingtime1 = coast_distance / Vmax1; %m
%--Braking
brakingtime1 = setbraking_distance / Vmax1; %sec
%--Total Operation Time
TotalTime1 = timepusherend1 + coastingtime1 + brakingtime1; %sec

%% Graph
%Correcting Time Intervals
coastingtime1_int = coastingtime1 + timepusherend1; %sec
brakingtime1_int = brakingtime1 + coastingtime1_int; %sec
% Test Figure 1 - Compute speed and position 1000 samples linearly spaced
t = linspace(0,TotalTime1,1000); Interval = t(2)-t(1);
% Zero acceleration
Acceleration = zeros(size(t));
% Less than Time to Maximum Velocity +ve Accn
Acceleration(t <= timepusherend1) = pusher_accelrate;
% greater than time after the measured mile -ve Accn
Acceleration(t > coastingtime1_int + brakingtime1_int) = -pusher_accelrate;
% Speed = Integral of Acceleration
Speed = cumsum(Acceleration)*Interval;
% Position = Integral of Speed
Position = cumsum(Speed)*Interval;
% Plotting commands - Test Figure 1
figure('name','Test Figure 1','numbertitle','off')
title('Test Figure 1 - Acceleration/Speed profile'); grid;
subplot(311)
plot(t,Acceleration,'LineWidth',2)
xlabel('Time (sec)'); ylabel('Acceleration (m/s^2)'); axis('normal')
title('UWHL Test Figure 1 - Accel/Speed profile'); grid;
subplot(312)
plot(t,Speed,'LineWidth',2)
xlabel('Time (sec)'); ylabel('Speed (m/s)'); grid; axis('normal')
subplot(313)
plot(t, Position,'LineWidth',2);
xlabel('Time (sec)'); ylabel('Position (meters)'); grid; axis('normal')
-Total Braking Distance and Maximum Deceleration
tr=1.2;%reaction time
ta=0.1;%application time
tb=0.18;%deceleration rise
vo=140;%max. velocity
dt=120; %total braking distance
dtot=0;
dmax=10.0;
%dmax(max. deceleration) when total braking distance=120 m
while(dtot<120)
d1=vo*(tr+ta);
d2=vo*tb-(dmax/6)*(tb^2);
d3=(vo^2)/(2*dmax)-vo*tb/2+dmax*(tb^2)/8;
dtot=d1+d2+d3; %total braking distance
dmax=dmax-0.0001;
end
% dmax and corresponding total braking distance
dmax

dtot
-Actuating Force, Braking Torque, Normal Pressure, Equivalent Radius, Force Location
h_max = 1.4; %[m] Max. Height of the Vehicle(ASSUMPTION)
mu = 0.85; %[] Coefficient of Friction Between Tires and Road(ASSUMPTION)
g = 9.81; %[m/s^2] Gravitational Coefficient
f_r = 0.055; %[] Rolling Resistance Coefficient(ASSUMPTION)
i_t_1 = 3.818; %[] First Gear Ratio
i_d = 4.412; %[] Differential Ratio
NRD=14*24.5; %[mm] Nominal Rim Diameter
NSW=175; %[mm] Nominal Section Width
PHI=0.60;
%[] Aspect Ratio
K=0.96; %[] Dimensionless Constant,0.96 for Radial Automobile Tires
m_unl = 532; %[kg] Unladen Mass
gama = 1.03+0.0016*(i_t_1*i_d)^2; %[] Rotating Mass Factor(ASSUMPTION)
m_eq_unl = m_unl*(1+gama); %[] Equivalent Mass
m_pass = 80*2; %[kg] Passenger Masw(95%Percentile(ASSUMPTION)
v_tank = 32; %[l] Fuel Tank Volume
ro_RON95 = 0.7431; %[kg/l] Density of RON95 [THUMMADETSAK]
m_fuel_max = v_tank*ro_RON95; %[kg] Fuel Mass (Full Tank)(ASSUMPTION)
m_lug = 30; %[kg] Mass of Additional Luggage(ASSUMPTION)
m_add = m_pass+m_fuel_max+m_lug; %[kg] Tot.Mass of Add. Factors(ASSUMPTION)
m_tot=m_eq_unl+m_add; %[kg] Total Mass
i=0.6701; %[]
Brake Force Distribution Factor
W=m_tot*g; %[N]
Weight of the Vehicle-Laden
dmax=-4; %[m/s^2] Maximum Deceleration
d=-(dmax/g); %[]
Dimensionless Deceleration
F_f=i*W*d %[N]
Front Braking Force
THETA=60*pi/180; %[rad] Caliper Angle Difference
R_w=K*((NRD/2)+PHI*NSW) %[mm] Tire Rolling Radius
T_bf=(F_f*R_w/1000)/2 %[Nm] Front Braking Torque
T=T_bf/2 %[Nm] Braking Torque on One Side of the Wheel
MU=0.40; %[]Coefficient of Friction for Rigid Molded Asbestos Pads
(ASSUMPTION) for 0.31-0.49
for R_r=122.5:2.5:135 %[mm] Rotor Radius for 14 in Rim Diameter
R_o=R_r-20 %[mm] Outer Radius of Pad! 20mm(ASSUMPTION)
R_i=R_o-35 %[mm] Inner Radius of Pad! 40mm(ASSUMPTION)
R_eff=(R_o+R_i)/2 %[mm] Effective Radius
F_d=(2*T)/(R_eff/1000) %[N]
F=T/(MU*(R_eff/1000)) %[N]
Actuating Force
P_a=F/(THETA*(R_o-R_i)*R_i) %[MPa] Largest Normal Pressure
P_hyd=9000*10^3 %[Pa]
Hydraulic Pressure
Pis_Num=2 %[]
Piston Number
Tot_Area=(F/P_hyd)*10^6 %[mm^2] Total Piston Area
Pis_Area=Tot_Area/2 %[mm^2] Allowable Piston Area for Ro-Ri=35 mm
R_pis=sqrt(Pis_Area/pi) %[mm]
Piston Radius
end

-Thermal Considerations, Temperature Rise
%% Thermal Analysis
k=1.15 %[] Correction Factor for Rotating Masses (Assumption)(1.05-1.15
for passenger cars in high gears)
V1=140 %[m/s]
Initial Velocity
m=m_tot %[kg]
Laden Mass of Vehicle
Eb=(k*m*V1^2)/2 %[J]
Energy Absorbed by Brake
ro_disc=7800 %[kg/m^3] Density of Disc Mat.(Steel(Assumption(7750-8050)
R_disc=130 %[mm]
Radius of Brake Disc (Assumption)
t=30 %[mm]
Thickness of Brake Disc (Assumption)(Minimum value
of the thickness is 28.1mm for ventilation brakes)
m_disc=ro_disc*pi*R_disc^2*t*10^-9 %[kg]
Mass of the Brake Disc
Cp=500 %[J/kg*C] Specific Heat Capacity
delta_T=Eb/(m_disc*Cp) %[C]
Temperature Rise
Tamb=25 %[C]
Temperature of Ambient
%% Tmax-Tamb=300
t1=60^2/24;
hr=27.5;
%[C] First Assumption for Temperature Rise
%[s] Brake was used 24 times per hour (Assumption)
%[W/m^2*C] Radiation Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
hc=7.5; %[W/m^2*C] Convective Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
MAS=12; %[m/s]
Mean Air Speed !!!Assumption
fv=6; %[]
Multiplying Factor(Shigley Figure 16.24b)
hcr=hr+fv*hc; %[W/m^2*C] Overall Heat Transfer Coefficient
A=4*pi*R_disc^2*10^-6; %[m^2]
Lateral Surface Area
W=m_disc; %[kg]
Mass of the Brake Disc
beta=(hcr*A)/(W*Cp); %[1/s]
Tmax=Tamb+(delta_T/exp(-beta*t1)); %[C]
Maximum Temperature
Temp_Rise=Tmax-Tamb; %[C]
Calculated Temperature Rise % First Assumption is not valid.
%% Tmax-Tamb=143.5295 %[C] Second Assumption for Temperature Rise
t1=60^2/24; %[s] Brake was used 24 times per hour (Assumption)
hr=12; %[W/m^2*C] Radiation Component of Heat Transfer Coefficient

hc=7.2; %[W/m^2*C] Convective Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
MAS=12; %[m/s] Mean Air Speed (Assumption)(Shigley Example 16.5)
fv=6; %[]
Multiplying Factor(Shigley Figure 16.24b)
hcr=hr+fv*hc; %[W/m^2*C] Overall Heat Transfer Coefficient
A=4*pi*R_disc^2*10^-6; %[m^2] Lateral Surface Area
W=m_disc; %[kg]
Mass of the Brake Disc
beta=(hcr*A)/(W*Cp); %[1/s]
Tmax=Tamb+(delta_T/exp(-beta*t1)); %[C]
Maximum Temperature
Temp_Rise=Tmax-Tamb;
%[C]
Calculated Temperature Rise
%Second Assumption is not valid
%% Tmax-Tamb=131.3444 %[C]
Third Assumption for Temperature Rise
t1=60^2/24; %[s]
Brake was used 24 times per hour(Assumption)
hr=11.1; %[W/m^2*C]
Radiation Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
hc=7.1; %[W/m^2*C]
Convective Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
MAS=12; %[m/s] Mean Air Speed (Assumption)(Shigley Example 16.5)
fv=6; %[]
Multiplying Factor(Shigley Figure 16.24b)
hcr=hr+fv*hc; %[W/m^2*C] Overall Heat Transfer Coefficient
A=4*pi*R_disc^2*10^-6; %[m^2]
Lateral Surface Area
W=m_disc; %[kg]
Mass of the Brake Disc
beta=(hcr*A)/(W*Cp); %[1/s]
Tmax=Tamb+(delta_T/exp(-beta*t1)); %[C]
Maximum Temperature
Temp_Rise=Tmax-Tamb; %[C]
Calculated Temperature Rise %Third Assumption is not valid
%% Tmax-Tamb=130.3379
t1=60^2/24
hr=11;
%[C] Fourth Assumption for Temperature Rise
%[s] Brake was used 24 times per hour(Assumption)
%[W/m^2*C] Radiation Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
hc=7.1; %[W/m^2*C] Convective Component of Heat Transfer
Coefficient(Shigley Figure 16.24a)
MAS=12; %[m/s] Mean Air Speed !!!Assumption(Shigley Example 16.5)
fv=6;%[]
Multiplying Factor(Shigley Figure 16.24b)
hcr=hr+fv*hc %[W/m^2*C] Overall Heat Transfer Coefficient
A=4*pi*R_disc^2*10^-6 %[m^2]
Lateral Surface Area
W=m_disc %[kg]
Mass of the Brake Disc

beta=(hcr*A)/(W*Cp) %[1/s]
Tmax=Tamb+(delta_T/exp(-beta*t1)) %[C]
Maximum Temperature
Temp_Rise=Tmax-Tamb
%[C] Calculated Temperature Rise
%Fourth Assumption is valid because the assumption and the calculated temperature rise is close enough to each other.
%Actual temperature rise=130.2711
Braking Operations
After accelerating down a straight on the racetrack, sooner or later, the car is going to come to a corner. If the
corner is tight enough, the driver will need to slow the car down in order to get around the corner. The speed of the car will
need to be reduced to the maximum corner speed by the time that the car arrives at the corner.
The initial calculation of maximum corner speed and forward acceleration can be done by progressing forwards
through each sector in turn, feeding the exit velocity from one sector into the entry velocity of the next. To calculate the
best possible braking performance, we first calculate the acceleration from rest and maximum corner speeds, then apply
the braking forces on a second pass, this time going backwards through the sectors.
For the first pass, it will be clear that we have to slow down for corners in some way or other. The simplest way to
do this seems to be as follows. When considering a new sector, the entry speed is compared against the maximum corner
speed for the sector. If the entry speed is higher, then the car would have needed to brake to this point. We reset the
corner entry to a reasonable value – the maximum corner speed.
Maximum corner speed
It is useful for us to be able to calculate the fastest speed the car could travel though a particular sector of a known radius.
We can then determine the braking effort required for the corner, and detect when the driver should balance the car on the
throttle at the maximum corner speed.
We assume that each sector is of a constant radius, and at the maximum ‘limiting’ speed, the speed will remain constant
through the corner. The maximum speed will be determined by the ability of the tyres to turn the car through the corner
and simultaneously push it forward through the air against drag to maintain a constant speed.
The maximum force that can be transmitted through the tyre without slippage is given by:
F = μ R
Where μ is the coefficient of friction of the tyre, and R is the ‘normal force’ acting on the tyre.
At rest, the normal force will be equal to the weight of the car (Mass x gravity). At speed the normal force will be affected
by the aerodynamics of the car – if it produces downforce, this will lead to the normal force being larger, and if it generates
lift at speed (a bad thing), then the normal force will be reduced.
The forces that must be transmitted through the tyre as it travels around the corner are:
Description
When the car is cornering at the maximum speed, the total force Ft will be equal to the maximum force that can be
transmitted through the tyres:
Ft = √Fc 2 + Fd 2 = μ R
We use a coefficient of friction for a racing tyre of μ = 1.4.

Braking Calculations
Mass of Vehicle – 300 kg (approx.)
Coefficient of friction of dry road- 0.75
Outer Diameter of brake disc- 1100mm
Inner Diameter of brake disc(d)- 700mm
Maximum speed of Vehicle(V)- 140km/hr
Pedal Ratio-4:1
Coefficient of friction between brake pads and disc- 0.4
Maximum stopping distance allowed -20m
Diameter of Master Cylinder- .025 m
Area of Brake Pad = 0.0018 m^2
Volume of Brake Disc= 94123.57 mm^3
Density of Brake Disc Material= 7850 Kg/m^3
Specific Heat Capacity of Brake Disc Material(c)=475 J/Kg/K
Thermal Conductivity of Brake Disc Material(k)= 44.5 W/(m K)
Ambient Temprature (Tamb.)= 22 (C)
Calculations :
Deceleration of Vehicle – Coeff. Of .fric. * g = 0.75 * 9.81 = 7.35 m/ s^2
Stopping Distance – ( V^2) / 2*coeff. Of fric. *g = (13.89)^2 / 2*0.75 * 9.81 = 13.125 m
Kinetic Energy of Vehicle – ½ * 300 * 13.89 *13.89= 28935.2.6 J
Stopping Time – Stopping Distance / Deceleration = 13.125 / 7.35 = 1.89 seconds
Braking Force – coeff. Of fric. * mass of vehicle * g = 0.75 * 300 * 9.81 = 2205 N
Effective disc Radius =1/3* (D^3 –d^3)/(D^2-d^2)= 74.83 mm
Braking Torque – (Braking Force ) * (Radius of Brake Disc ) = 1102.5 * .7483= 82.497 N m
Clamping Force- 2205 / 0.4 = 5512.5 N
Required Brake Pressure = (Clamping Force/Brake Pad area )=5512.5 / 0.0018= 3062.5 K Pa
Maximum Brake Pressure Applicable in Brake Disc – 2300 K Pa
Since allowable Pressure is greater than the Pressure applicable in current case therefore the brake disc chosen for
braking is suitable.
Average Power (P)= Kinetic Energy (E)/Stopping Time (T) =14467.6/1.89=7654.815 W
Maximum power at the onset of braking = 2* Average Power
Heat Flux into one side of the Disc= 4P/3.14*(D^2-d^2)= 4*7654.815 /3.14* (.19^2- .1^2)
=37.3615 W/cm^2
Maximum Temprature=(.527 *Q *t^1/2) /(Density of brake disc Material *c* k)+ T(amb.)
(Single Stop) = (.527*373615.2*1.38)/(7850*44.5*475)^1/2 +22(C)
=43.1 (C)
Fade Stop Temprature Rise = P*t / (Density of Brake disc material *c *V)
= 7654.815*1.89 /(7850*475*.0000941)
= 41.233 (C)

References
Ahn, J.K., Jung, K.H., Kim, D.H., Jin, H.B., Kim, H.S. and Hwang, S.H., 2009. Analysis of a regenerative braking system
for hybrid electric vehicles using an electro-mechanical brake. International Journal of Automotive Technology, 10(2),
pp.229-234.
Nian, X., Peng, F. and Zhang, H., 2014. Regenerative braking system of electric vehicle driven by brushless DC
motor. IEEE Transactions on Industrial Electronics, 61(10), pp.5798-5808.
Xiu-qin, Z., Bo, Y., Chao, Y. and Guan-neng, X., 2012. Research on ABS of multi-axle truck based on ADAMS/Car and
Matlab/Simulink. Procedia Engineering, 37, pp.120-124.
Pugi, L., Malvezzi, M., Allotta, B., Banchi, L. and Presciani, P., 2004. A parametric library for the simulation of a Union
Internationale des Chemins de Fer (UIC) pneumatic braking system. Proceedings of the institution of mechanical
engineers, part F: journal of rail and rapid transit, 218(2), pp.117-132.
Peng, D., Zhang, Y., Yin, C.L. and Zhang, J.W., 2008. Combined control of a regenerative braking and antilock braking
system for hybrid electric vehicles. International Journal of Automotive Technology, 9(6), pp.749-757.
Haber, A., Kolassa, C., Manhart, P., Nazari, P.M.S., Rumpe, B. and Schaefer, I., 2013, January. First-class variability
modeling in Matlab/Simulink. In Proceedings of the Seventh International Workshop on Variability Modelling of Software-
intensive Systems (p. 4). ACM.
Kim, D., Hwang, S. and Kim, H., 2008. Vehicle stability enhancement of four-wheel-drive hybrid electric vehicle using rear
motor control. IEEE Transactions on Vehicular Technology, 57(2), pp.727-735.
Li, J., Yu, F., Zhang, J.W., Feng, J.Z. and Zhao, H.P., 2002. The rapid development of a vehicle electronic control system
and its application to an antilock braking system based on hardware-in-the-loop simulation. Proceedings of the Institution
of Mechanical Engineers, Part D: Journal of Automobile Engineering, 216(2), pp.95-105.