In :# Configure Jupyter so figures appear in the notebook %matplotlib inline # Configure Jupyter to display the assigned value after an assignment %config InteractiveShell.ast_node_interactivity='last_expr_or_assign' # import functions from the modsim.py module from modsim import *
The goal of this case study is to use simulation to guide the design of a car intended to accelerate from standing to 100 kph as quickly as possible. The world record for this event, using a car that meets the competition requirements, is 1.513 seconds.
We'll start with a simple model that takes into account the characteristics of the motor and vehicle:
The motor is an Emrax 228 high voltage axial flux synchronous permanent magnet motor; according to the data sheet, its maximum torque is 240 Nm, at 0 rpm. But maximum torque decreases with motor speed; at 5000 rpm, maximum torque is 216 Nm.
The motor is connected to the drive axle with a chain drive with speed ratio 13:60 or 1:4.6; that is, the axle rotates once for each 4.6 rotations of the motor.
The radius of the tires is 0.26 meters.
The weight of the vehicle, including driver, is 300 kg.
To start, we will assume no slipping between the tires and the road surface, no air resistance, and no rolling resistance. Then we will relax these assumptions one at a time.
First we'll add drag, assuming that the frontal area of the vehicle is 0.6 square meters, with coefficient of drag 0.6.
Next we'll add rolling resistance, assuming a coefficient of 0.2.
Finally we'll compute the peak acceleration to see if the "no slip" assumption is credible.
We'll use this model to estimate the potential benefit of possible design improvements, including decreasing drag and rolling resistance, or increasing the speed ratio.
I'll start by loading the units we need.
In :radian = UNITS.radian m = UNITS.meter s = UNITS.second minute = UNITS.minute hour = UNITS.hour km = UNITS.kilometer kg = UNITS.kilogram N = UNITS.newton rpm = UNITS.rpm
And store the parameters in a
In :params = Params(r_wheel=0.26 * m, speed_ratio=13/60, C_rr=0.2, C_d=0.5, area=0.6*m**2, rho=1.2*kg/m**3, mass=300*kg)
make_system creates the initial state,
init, and constructs an
interp1d object that represents torque as a function of motor speed.
In :def make_system(params): """Make a system object. params: Params object returns: System object """ init = State(x=0*m, v=0*m/s) rpms = [0, 2000, 5000] torques = [240, 240, 216] interpolate_torque = interpolate(Series(torques, rpms)) return System(params, init=init, interpolate_torque=interpolate_torque, t_end=3*s)
In :system = make_system(params)
The relationship between torque and motor speed is taken from the Emrax 228 data sheet. The following functions reproduce the red dotted line that represents peak torque, which can only be sustained for a few seconds before the motor overheats.
In :def compute_torque(omega, system): """Maximum peak torque as a function of motor speed. omega: motor speed in radian/s system: System object returns: torque in Nm """ factor = (1 * radian / s).to(rpm) x = magnitude(omega * factor) return system.interpolate_torque(x) * N * m
In :compute_torque(0*radian/s, system)
In :omega = (5000 * rpm).to(radian/s) compute_torque(omega, system)
Plot the whole curve.
In :xs = linspace(0, 525, 21) * radian / s taus = [compute_torque(x, system) for x in xs] plot(xs, taus) decorate(xlabel='Motor speed (rpm)', ylabel='Available torque (N m)')
In :def slope_func(state, t, system): """Computes the derivatives of the state variables. state: State object t: time system: System object returns: sequence of derivatives """ x, v = state r_wheel, speed_ratio = system.r_wheel, system.speed_ratio mass = system.mass # use velocity, v, to compute angular velocity of the wheel omega2 = v / r_wheel # use the speed ratio to compute motor speed omega1 = omega2 / speed_ratio # look up motor speed to get maximum torque at the motor tau1 = compute_torque(omega1, system) # compute the corresponding torque at the axle tau2 = tau1 / speed_ratio # compute the force of the wheel on the ground F = tau2 / r_wheel # compute acceleration a = F/mass return v, a
slope_func at linear velocity 10 m/s.
In :test_state = State(x=0*m, v=10*m/s)
In :slope_func(test_state, 0*s, system)
Now we can run the simulation.
In :results, details = run_ode_solver(system, slope_func) details
And look at the results.
After 3 seconds, the vehicle could be at 40 meters per second, in theory, which is 144 kph.
In :v_final = get_last_value(results.v)
In :def plot_position(results): plot(results.x, label='x') decorate(xlabel='Time (s)', ylabel='Position (m)') plot_position(results)
In :def plot_velocity(results): plot(results.v, label='v') decorate(xlabel='Time (s)', ylabel='Velocity (m/s)') plot_velocity(results)
In :def event_func(state, t, system): """Stops when we get to 100 km/hour. state: State object t: time system: System object returns: difference from 100 km/hour """ x, v = state # convert to km/hour factor = (1 * m/s).to(km/hour) v = magnitude(v * factor) return v - 100
In :results, details = run_ode_solver(system, slope_func, events=event_func) details
Here's what the results look like.
In :subplot(2, 1, 1) plot_position(results) subplot(2, 1, 2) plot_velocity(results) savefig('figs/chap11-fig02.pdf')
According to this model, we should be able to make this run in just over 2 seconds.
In :t_final = get_last_label(results) * s
At the end of the run, the car has gone about 28 meters.
In :state = results.last_row()
If we send the final state back to the slope function, we can see that the final acceleration is about 13 $m/s^2$, which is about 1.3 times the acceleration of gravity.
In :v, a = slope_func(state, 0, system) v.to(km/hour)
In :g = 9.8 * m/s**2 (a / g).to(UNITS.dimensionless)
It's not easy for a vehicle to accelerate faster than
g, because that implies a coefficient of friction between the wheels and the road surface that's greater than 1. But racing tires on dry asphalt can do that; the OEM team at Olin has tested their tires and found a peak coefficient near 1.5.
So it's possible that our no slip assumption is valid, but only under ideal conditions, where weight is distributed equally on four tires, and all tires are driving.
Exercise: How much time do we lose because maximum torque decreases as motor speed increases? Run the model again with no drop off in torque and see how much time it saves.
In this section we'll see how much effect drag has on the results.
Here's a function to compute drag force, as we saw in Chapter 21.
In :def drag_force(v, system): """Computes drag force in the opposite direction of `v`. v: velocity system: System object returns: drag force """ rho, C_d, area = system.rho, system.C_d, system.area f_drag = -np.sign(v) * rho * v**2 * C_d * area / 2 return f_drag
We can test it with a velocity of 20 m/s.
In :drag_force(20 * m/s, system)
Here's the resulting acceleration of the vehicle due to drag.
In :drag_force(20 * m/s, system) / system.mass
We can see that the effect of drag is not huge, compared to the acceleration we computed in the previous section, but it is not negligible.
Here's a modified slope function that takes drag into account.
In :def slope_func2(state, t, system): """Computes the derivatives of the state variables. state: State object t: time system: System object returns: sequence of derivatives """ x, v = state r_wheel, speed_ratio = system.r_wheel, system.speed_ratio mass = system.mass omega2 = v / r_wheel * radian omega1 = omega2 / speed_ratio tau1 = compute_torque(omega1, system) tau2 = tau1 / speed_ratio F = tau2 / r_wheel a_motor = F / mass a_drag = drag_force(v, system) / mass a = a_motor + a_drag return v, a
And here's the next run.
In :results2, details = run_ode_solver(system, slope_func2, events=event_func) details
The time to reach 100 kph is a bit higher.
In :t_final2 = get_last_label(results2) * s
But the total effect of drag is only about 2/100 seconds.
In :t_final2 - t_final
That's not huge, which suggests we might not be able to save much time by decreasing the frontal area, or coefficient of drag, of the car.
Next we'll consider rolling resistance, which the force that resists the motion of the car as it rolls on tires. The cofficient of rolling resistance,
C_rr, is the ratio of rolling resistance to the normal force between the car and the ground (in that way it is similar to a coefficient of friction).
The following function computes rolling resistance.
In :system.set(unit_rr = 1 * N / kg)
In :def rolling_resistance(system): """Computes force due to rolling resistance. system: System object returns: force """ return -system.C_rr * system.mass * system.unit_rr
The acceleration due to rolling resistance is 0.2 (it is not a coincidence that it equals
In :rolling_resistance(system) / system.mass
Here's a modified slope function that includes drag and rolling resistance.
In :def slope_func3(state, t, system): """Computes the derivatives of the state variables. state: State object t: time system: System object returns: sequence of derivatives """ x, v = state r_wheel, speed_ratio = system.r_wheel, system.speed_ratio mass = system.mass omega2 = v / r_wheel * radian omega1 = omega2 / speed_ratio tau1 = compute_torque(omega1, system) tau2 = tau1 / speed_ratio F = tau2 / r_wheel a_motor = F / mass a_drag = drag_force(v, system) / mass a_roll = rolling_resistance(system) / mass a = a_motor + a_drag + a_roll return v, a
And here's the run.
In :results3, details = run_ode_solver(system, slope_func3, events=event_func) details
The final time is a little higher, but the total cost of rolling resistance is only 3/100 seconds.
In :t_final3 = get_last_label(results3) * s
In :t_final3 - t_final2
So, again, there is probably not much to be gained by decreasing rolling resistance.
In fact, it is hard to decrease rolling resistance without also decreasing traction, so that might not help at all.
The gear ratio 13:60 is intended to maximize the acceleration of the car without causing the tires to slip. In this section, we'll consider other gear ratios and estimate their effects on acceleration and time to reach 100 kph.
Here's a function that takes a speed ratio as a parameter and returns time to reach 100 kph.
In :def time_to_speed(speed_ratio, params): """Computes times to reach 100 kph. speed_ratio: ratio of wheel speed to motor speed params: Params object returns: time to reach 100 kph, in seconds """ params = Params(params, speed_ratio=speed_ratio) system = make_system(params) system.set(unit_rr = 1 * N / kg) results, details = run_ode_solver(system, slope_func3, events=event_func) t_final = get_last_label(results) a_initial = slope_func(system.init, 0, system) return t_final
We can test it with the default ratio:
In :time_to_speed(13/60, params)
Now we can try it with different numbers of teeth on the motor gear (assuming that the axle gear has 60 teeth):
In :for teeth in linrange(8, 18): print(teeth, time_to_speed(teeth/60, params))
Wow! The speed ratio has a big effect on the results. At first glance, it looks like we could break the world record (1.513 seconds) just by decreasing the number of teeth.
But before we try it, let's see what effect that has on peak acceleration.
In :def initial_acceleration(speed_ratio, params): """Maximum acceleration as a function of speed ratio. speed_ratio: ratio of wheel speed to motor speed params: Params object returns: peak acceleration, in m/s^2 """ params = Params(params, speed_ratio=speed_ratio) system = make_system(params) a_initial = slope_func(system.init, 0, system) * m/s**2 return a_initial
Here are the results:
In :for teeth in linrange(8, 18): print(teeth, initial_acceleration(teeth/60, params))
As we decrease the speed ratio, the peak acceleration increases. With 8 teeth on the motor gear, we could break the world record, but only if we can accelerate at 2.3 times the acceleration of gravity, which is impossible without very sticky ties and a vehicle that generates a lot of downforce.
In :23.07 / 9.8
These results suggest that the most promising way to improve the performance of the car (for this event) would be to improve traction.