Skip to content

My Site Has Moved to ConsultKeithYoung

Over the past month or so I have been migrating my site to my own hosted site.

Consult Keith Young

Here is a preview of new content on the site that isn’t here. Sorry I’m not posting here any more, but it’s time to move on to bigger and better things. I hope you enjoy the new site even more than this old one!



Calculating Boundary Layer Thickness to Choose Pitot Tube Height on F1 Car

First, here are some Prereq articles:

Bernoullis Equation for a Pitot Tube Example on a Formula 1 Race Car

Pitot Tube and Pressure Explanation With F1 and Pikes Peak Example

Introduction to Boundary Layers, Viscous Flow, and Velocity Profiles

The answer to the question in the Pikes Peak article will be below the examples here.

Lets begin.

We want to see what happens if we place the Pitot Tube too low, and we want to find the Boundary Layer Thickness to know how high the Pitot Tube must be in order to be outside the Boundary Layer, and therefore be measuring the free stream velocity.

If the Pitot Tube is too low, then it will be inside the Boundary Layer and therefore will not be measuring the free stream velocity.

Example 1 – Finding a Pitot Tube height that will ensure it’s measuring the vehicle speed and not a lower speed due to the boundary layer

Lets say we want to know the airspeed of an F1 car, and we want to put the Pitot Tube some distance up the nose. How high up does the Pitot Tube need to be mounted in order to be measuring the free stream velocity outside the boundary layer?

  • Pitot Tube is 1m back from the nose
  • Air Density is 1.225 kg/m^3
  • Air Viscosity is 1.8*10^-5 kg/ms
  • Speed range from 10 m/s to 100 m/s (360km/h or 224 mph)


  1. Infinately wide Flat Plate
  2. Incompressible
  3. Instant transition from laminar to turbulent flow
  4. Time Averaged Mean Flow
  5. Steady State
  6. Turbulent Transition at Re=500,000

Lets first introduce the equations we’ll be needing.

Re_x = \frac{\rho * V * x}{\mu}

\delta_{lam} = \frac{5*x}{\sqrt{Re_x}}

\delta_{tur} = \frac{.3747*x}{(Re_x)^{0.2}}

x_{tr} = \frac{Re_{tr}}{\frac{\rho * V}{\mu}}

To read the rest of the article please go to Calculating Boundary Layer Thickness to Choose Pitot Tube Height on F1 Car using Theory and Pitot Tube Arrays. I am in the process of migrating my website there, which is why I’m putting partial articles where my viewers can still find them.

Validating My Theory Based on Couette-Poiseuille Flow For Race Car Underbodies

In an older post about Couette Poiseuille Flow Using Aerodynamic Navier Stokes Equations I derived a theoretical equation to simulate the behavior of airflow underneath a highly idealized race car underfloor. In this writeup I intend to validate that theory by checking the Boundary Conditions of the equation representing the velocity profile of the airflow, and checking the mass flow rate using the Continuity Equation.

Here is the equation I derived last time:

u(y)=\frac{1}{2} [\frac{1}{\mu} \frac{dP}{dx}] y^{2} + \frac{V}{2h} y + \frac{\mu V - \frac{dP}{dx} h^{2}}{2\mu}

What the equation represents is the velocity of the air flow in the x direction as a function of the y direction. Let me explain that briefly, I can make a more complete post on velocity profiles, boundary layers, and shear later. The air directly touching the ground is stationary since it is in contact with the ground. The air just above the ground has “friction” with the air below it. The same principle is in effect for the bottom of the race car. The air touching the bottom of the car is stationary with respect to the car, and “friction” effects the nearby air molecules. I say “friction” because in reality the air molecules are bouncing off each other. This creates a velocity profile.

So in saying all that jibber jabber above, it simply means the air at different heights between the ground and underfloor of the race car has different speeds, and we are taking that in to account with the equation. So if this article is over your head for now, bookmark it, and come back. Like I said, I’ll make a more detailed post on this later. Promise.

So the first method I will use to validate the theory, or mostly to make sure it isn’t dead wrong, is to double check the Boundary Conditions. I said in a quick sentence in my last post that I had done it, but I didn’t describe how the method works.

We know the air speed on the ground, and on the underbody. We can take the derived equation and make sure the speed is in agreement.

Lets assume we have a race car and atmosphere with the following specs:

  • Ride height =  .1m
  • Speed          = 40m/s
  • Density        = 1.225kg/m^3
  • Viscosity      = 1.8*10^-5 kg/(s*m)
  • dP/dx=0

Lets plug and chug.


u(y)=\frac{1}{2} [\frac{1}{\mu} \frac{dP}{dx}] y^{2} + \frac{V}{2h} y + \frac{\mu V - \frac{dP}{dx} h^{2}}{2\mu}

u(y)=\frac{1}{2} [\frac{1}{1.8*10^3} *0]*[-.05]^{2} + \frac{40}{2*.05}*[-.05] + \frac{1.8*10^3 40 - 0}{2*1.8*10^3}=0

u(y)=\frac{1}{2} [\frac{1}{1.8*10^3} *0]*[.05]^{2} + \frac{40}{2*.05}*[.05] + \frac{1.8*10^3 40 - 0}{2*1.8*10^3}=40

So at the very least, we solved correctly for the Boundary Conditions.

Now, lets check the mass flow rate using the Continuity Equation.

\int \int \rho \vec{V} \cdot \vec{n} dA = 0


Since we are treating this as a 2 dimensional problem we can replace the first integral with a unit depth since there is no velocity change in the z direction due to the assumptions made in the previous article. Also, since the velocity only has an x component we don’t have any extra mess from the dot product.

\int_{-h}^{h}\rho[\frac{1}{2} [\frac{1}{\mu} \frac{dP}{dx}] y^{2} + \frac{V}{2h} y + \frac{\mu V - \frac{dP}{dx} h^{2}}{2\mu}]dy=\frac{h\rho[3\mu V-2h^2*\frac{dp}{dx}}{3\mu}] 102

\frac{h\rho[3\mu V-2h^2*\frac{dp}{dx}]}{3\mu}

Since I am slowly transitioning to my new site, please read the rest here…

Validating My Theory Based on Couette-Poiseuille Flow For Race Car Underbodies | Consult Keith Young

Pitot Tube 2

First of all, I am slowly migrating to my new site about F1 Race Car Technology.

In my previous Pitot Tube Example Using Bernoulli’s Equation I made some simplifying assumptions and did not explain things as clearly as I could have.

Before  you can really understand how Pitot Tubes work, whether it’s for aircraft or a race car, you need to understand the different “types” of pressure. You may have seen words like Dynamic Pressure, Static, Total, Absolute, Stagnation, Gage,  and others.

First, I will describe Static Pressure. For a more official description you can read a book, but the way I think of static pressure is the pressure I feel in stationary air. At the bottom of a pool that pressure goes up, up in a mountain that pressure goes down. What happens though when you are in a race car traveling at speed? The static pressure in this case would be the pressure of the air measured by an observer traveling with the air, not against it.

Static Pressure is created when the random motions of air molecules interact with a body. As the particles bounce off the surface, they transfer momentum to the surface and subjects it to a force.


In the picture above, the random motion paths are shown with the red lines. A few collisions are observed and they are imparting small forces on the body shown with black arrows.

Dynamic Pressure is used in finding lift and drag. It takes in to account the density and velocity of the air, while basically ignoring the little random motions of the air that make up the static pressure.


The image above shows that Dynamic Pressure doesn’t care about the random motions of the particles, only the momentum transfer due to the density of the air and the freestream velocity.


The image above shows how an object with more mass will impart a larger force than a smaller object with equal velocity.


Additionally, the faster of the two objects with equal mass will impart a larger force on the body.

q = \frac{1}{2} \rho V^{2}

Based on the equation above, for Dynamic Pressure, it is clear that changes in freestream velocity have a greater effect than does the density of the free stream.

Stagnation Pressure and Total Pressure are redundant terms. Stagnation pressure is the pressure measured if the moving airflow is brought to a stop Isentropically. Intuitively it makes sense. If air flow is hitting the front of an F1 car the pressure is higher since that flow had to slow down. Basically stagnation pressure adds up the static pressure from random molecular motion, and the dynamic pressure made up of the actual freestream velocity of the air.

P_0 = P + \frac{1}{2} \rho V^{2} = P + q

That’s all great. How does this apply to Pitot Tubes? A Pitot Tube is used to find the velocity of a fluid, in the case of a race car that fluid is air. To do this, a Pitot Tube compares the Static Pressure to the Stagnation Pressure.


Once both the Static Pressure and the Stagnation Pressures are known, the velocity of the race car relative to the air can be determined. This can be done using Bernoulli’s Equation.

\frac{P_1}{\rho} + \frac{1}{2} V_1^{2} + g h_1 = \frac{P_2}{\rho} + \frac{1}{2} V_2^{2} + g h_2

Since the Formula 1 car speed is much less than the speed of sound, we can make some assumptions such as constant air density. We can also assume that the difference in height of the fluid is negligible. Therefore we can rewrite the equation.

\frac{P_1}{\rho} + \frac{1}{2} V_1^{2} = \frac{P_2}{\rho} + \frac{1}{2} V_2^{2}

This is better, but still a little confusing. It can be made much simpler when we multiply both sides by the density.

P_1 + \frac{1}{2} \rho V_1^{2} = P_2 + \frac{1}{2} \rho V_2^{2}

Now let’s designate the variables. Let’s have the first values relate to the Stagnation portion and the second set of values relate to the Static values. We know the velocity at the Pitot Tube inlet is zero, hence the Stagnation Pressure. The Stagnation Pressure is measured, and therefore known. Next, we have also measured the Static Pressure. Lastly, we should know the air density. The only thing we don’t know is the Velocity of the freestream as it passes by the hole, or port. Since the Velocity is the only unknown, we can solve.

P_{stag} = P_{stat} + \frac{1}{2} \rho V_2^{2}

With a little bit of algebra, we arrive at the equation for the velocity.

V = \sqrt{\frac{[2 [P_{stag} - P_{stat}]]}{\rho}}

F1 Car Pitot Tube Example, Pikes Peak International Hill Climb:

But let’s make it a little more difficult than that. An F1 car runs up Pikes Peak. Toward the top of the mountain, what Velocity is required to make the Stagnation/Static Pressure Ratio equal to 1.0082?

According to Wolfram Alpha:

Height of Pikes Peak: 4,300m

Air Density at 4,300m: .79kg/m^3

Air Pressure at 4,300m: 59kPa

P_1 + \frac{1}{2} \rho V_1^{2} = P_2 + \frac{1}{2} \rho V_2^{2}

P_{stag} = P_{stat} + \frac{1}{2} \rho V_2^{2}

\frac{P_{stag}}{P_{stat}} = 1 + \frac{\frac{1}{2} \rho V_2^{2}}{P_{stat}}

V_2 = \sqrt{2 \frac{P_{stat}}{\rho} [\frac{P_{stag}}{P_{stat}}-1]}

V_2 = \sqrt{2 \frac{59000}{.79} [1.0082-1]} = 35 \frac{m}{s}

P_{Static} = 59{kPa}

P_{Dynamic} = \frac{1}{2} \rho V^{2} = 0.484{kPa}

P_{Stagnation} = P_{Static} + P_{Dynamic} = P + \frac{1}{2} \rho V^{2} = 59.484{kPa}

Hope this helped you understand these things a little better. This article tried to explain the differences between Static Pressure, Dynamic Pressure, and Stagnation Pressure aka Total Pressure. It also showed how to use Bernoulli’s Equation to find the velocity of the freestream by using a Pitot Tube. Please comment and rate. I am fairly new to writing and need to know if the general reader thinks this material and example was too complicated or too simple.

Think you get it? What Stagnation/Static Pressure ratio is required when the racing car is at ground level at the same speed? You can see if you got it correct in my next article.

Racing Brakes Numerical Analysis 2 – Initial Mathematical Model


    This is the second writing in a series of papers on the numerical analysis of brake rotors used in race cars. The purpose of this writing is to describe the mathematical model for this simplest initial form of the code before it is updated. The model was based on the conservation of energy. The equations will be described, as well as the algorithms and the method for error analysis.


    The goal of this paper was to introduce the mathematical model of the oversimplified brake heat model. The conservation of energy method was used. The solution to a chosen node was detailed, and the solutions to other nodes with different boundary condition types were included in the appendices. There were in total 11 different types of boundary conditions used in this current math model.

    It was found that improvements for the math model include replacing loops with matrix operations to speed iterations, and increasing the order of accuracy from one to two.



    The goal of this writing was to introduce the mathematical model used in the previous writing, which was the simplest form used mainly to troubleshoot the code. The purpose of this code and the equations in it, at the time, were not meant to represent brake rotor heating in a true sense, but rather to provide a troubleshooting platform for the overall structure of the code. Despite this, it was decided to release the math model due to its usefulness outside of braking as well as introduce the current error analysis method.


    Numerical methods are common in today’s engineering. As problems become more and more complicated, purely analytical solutions become more difficult to achieve. Given the simultaneous improvement in computing power and reduced costs, numerical analysis is quickly becoming more important.

    Numerical analysis is a huge topic. It can be used to simulate air flow around a car, heat exchange through the radiators, deflections in vehicle structure, and in this case the heating of brake rotors. Many methods exist, such as Newton’s Method, RK4, use of the Taylor Series, and others. The types of mathematical problems that can be solved include sets of equations, nonlinear equations, interpolation and curve fitting, numerical integration, and all the way up to higher order partial differential equations.

    Numerical methods can also be intuitive. In fact in 8th grade I independently invented several methods, such as numerical integration and the secant method. Other people intuitively do this to. For example, a driver may try to estimate distance driven. They could add up how much time they spent at 60mph, 30mph, and stopped. With the ability to perform very simple arithmetic a value can be estimated. The basic goal of numerical methods in computing is to make more slices of time, more velocity measurements, carry more decimal places, and arrive at a more accurate answer. This is an oversimplified explanation, but may help the reader understand the nature of numerical computing.


Mathematical model

    First, the basic equation of heat transfer will be introduced. Basically what it says is that at any point the sum of heat entering minus the sum of heat leaving equals the change of internal energy at that point in a given amount of time. In this case, the model was assumed as steady state, which by definition means that rate of change of temperature with respect to time is zero [1].

\frac{\partial}{\partial x} (K \frac{\partial T}{\partial x})+ \frac{\partial}{\partial y} (K \frac{\partial T}{\partial y})+ \dot{q}= \rho c_p \frac{\partial T}{\partial t} = 0

\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\dot{q}}{K} = \frac{1}{\alpha} \frac{\partial T}{\partial t} = 0

    The method used was the conservation of energy.

\dot{E_{in}} = \dot{E_{out}}

    In order to realize this goal at each node, heat from every direction had to be accounted for. One assumption was that no heat traveled in the Z direction. The reason for this was that given the length chosen of 1m, the heat transfer out the ends of this length would be negligible compared to the heat flow out the cross sectional area chosen.

    The example node that will be analyzed in this case is the only node which observes the 9th boundary condition according to Figure 1. This particular node has heat entering over only half of its exposed surface. It has two instances of conduction through half the usual cross sectional area, as well as convection and radiation. This makes it a good representative of the many conditions that can be faced by each node in the 11 different boundary types used.

Figure 1 – The cross section of the model and its 11 boundary condition types

Figure 2 – Detail of the 9th type of boundary condition

Figure 3 – Discretization method showing rows and columns, X’s and Y’s

    Using Figure 2 above all the heats encountered can be identified, quantified, and then summed. First, the short hand will be explained. Any X or Y is an indicator of the column or row respectively in the matrix in which the temperature is being sought. A 1 or 2 immediately after the letter indicates the node is +1 or -1 in that direction. For example, Y2X indicates row Y-1 and column X. The reason the Y and X are swapped from their usual position is that Matlab first wants the row then the column input (Fig. 3).

q_{(y,x-1) \to (y,x)} = K (\Delta y * 1m) \frac{T_{(y,x-1)} - T_{(y,x)}}{\Delta x}

    But \Delta y = \Delta x = \Delta s = ds .

q_{(y,x-1) \to (y,x)} = K [T_{(y,x-1)} - T_{(y,x)}]

q_{(y+1,x) \to (y,x)} = \frac{1}{2} K [T_{(y+1,x)} - T_{(y,x)}]

q_{(y-1,x) \to (y,x)} = \frac{1}{2} K [T_{(y-1,x)} - T_{(y,x)}]

q_{(q) \to (y,x)} = \frac{1}{2} \frac{ds}{padheight} q''

q_{(\infty) \to (y,x)} = h ds [T_{\infty} - T_{(y,x)}]

q_{(rad) \to (y,x)} = ds \epsilon \sigma T_{(y,x)}^4

K [T_{(y,x-1)} - T_{(y,x)}] + \frac{1}{2} K [T_{(y+1,x)} - T_{(y,x)}] +\frac{1}{2} K [T_{(y-1,x)} - T_{(y,x)}] + \frac{1}{2} \frac{ds}{padheight} q'' + h ds [T_{\infty} - T_{(y,x)}] + ds \epsilon \sigma T_{(y,x)}^4 = 0

    Solving for the temperature at the node after substituting a different name to the same temperature that is used to calculate radiation (“thing” was used) gets

T_{(y,x)} = \frac{K padh [2 T_{(y,x-1)} + T_{(y+1,x)} + T_{(y-1,x)}] + 2 h ds padh T_{\infty} + ds [q'' - 2 padh \epsilon \sigma T_{(y,x)}^4]}{2 padh (h ds + 2 K)}

    See the appendices for the solution to other boundary condition types.

    In order to find the error a rather simple method was used. It was placed after the computational loop in order to now slow the solution of the temperature field. In order to find error, every exposed area of surface was analyzed for heat escape, and all of these values were summed. They were then compared to the heat input from the brake pads, which was a known in this case.



Original Writer: Keith A Young

Date: 03/24/2012

Latest Update: 03/26/2012


Email: k dot young at consultkeithyoung dot com

This code may only be used under the following conditions:

Non profit use – This code may be used or modified by anyone for non profit

use only if credit is given to the original writer for any code updates,

modifications, or results. Any release of code originating from this must

give its users the same terms as listed in this note, and most include this


For profit – Permission must first be granted. For permission, contact keith dot

a dot young at comcast dot net.

Please do not remove this note. If you prefer, it is easiest to move it to

the bottom of the code.

All other rights are held by the original writer, Keith A Young.

All units are in Metric %}


qu = 0;

qu = qu + (t(y,x)-tinf)/(1/(hcon*.5*ds))+.5*ds*sig*t(y,x)^4;

for x = 2:sizex

qu = qu + (t(y,x)-tinf)/(1/(hcon*ds))+ds*(sig*t(y,x)^4);


for y=2:sizey

qu = qu + (t(y,x)-tinf)/(1/(hcon*ds))+ds*(sig*t(y,x)^4);


for x = 2:sizex-1

qu = qu + (t(y,x)-tinf)/(1/(hcon*ds))+ds*(sig*t(y,x)^4);



qu = qu + (t(y,x)-tinf)/(1/(hcon*.5*ds))+.5*ds*sig*t(y,x)^4;

qdotin = qdotin;

qu = qu;

error = (qdotin – qu)/qdotin


Results and discussion

    The mathematical model was sufficient for the simplified treatment given. However the use of so many loops didn’t utilize the strength of Matlab, which is performing matrix operations. Also, it focused on one of the weaknesses of Matlab, being the speed at which it can perform loops.

    The method used achieved first order accuracy, where second order of accuracy would be more ideal.

Conclusion and recommendations

    It was concluded that speed may be increased by utilizing matrix operations rather than loops. Also, the order of accuracy isn’t as good as it could be, being only being first order.

    Recommendations for the future include replacing as many loop operations as possible with matrix operations. This could help increase the number of operations per second and therefore may reduce solution time to a certain residual. Another improvement would be increasing the order of accuracy from one to two.


    Thanks for all the input from the forums at

    Thread –


[1] Incropera et all, Fundamentals of Heat and Mass Transfer 6th Ed., Wiley

[2] C. Gerald, P Wheatley, Applied Numerical Analysis 7th Ed., Pearson Addison Wesley

[3] D. Pitts, L. Sissom, Heat Transfer 2nd Ed., Schaum’s Outlines, McGraw Hill

[4] R.W. Hamming, Numerical Methods for Scientists and Engineers 2nd Ed., Dover Publications

[5] F. Talati, S. Jalalifar, Investigation of heat transfer phenomena in a ventilated disk brake rotor with straight radial rounded vanes, Journal of Applied Sciences 8 (20) (2008) 3583-3592

[6] J. Zhang, Numerical Accuracy,,


    Solutions for each boundary condition type


Boundary condition types









Racing Brakes Numerical Analysis 1 – Introduction of the problem


    This is the first writing in what will be a series of papers on the numerical analysis of brake rotors used in race cars. The purpose of this writing is to describe the goal of this series, to introduce the problem, its background, and the method that will be used to solve it. The most basic initial model will be introduced, along with a list of simplifying assumptions that were used, and the initial results it produced. Being the first writing in a series, the topic is covered more broadly, and many aspects of the topic are introduced at once. This resulted in a longer writing than what is anticipated for future writings.

    A detail of the mathematical model will be released at a later date.


    The goal of this paper was to introduce the topic of the series and release preliminary results. The limits of the current code were mentioned, and the effects were discussed. Iterations, error, maximum temperature, and temperature distribution were observed with changes in relaxation, step size, and residual.

    The maximum max temperature was 1,166K while the minimum max temperature was 1,158.7, a deviation of only .6%. It was observed that adding a relaxation factor significantly improved performance. For a residual of E-5 and a step size of .00025m the version without a relaxation factor took 188,642 iterations while with a relaxation factor only 9.905 iterations were necessary. The initial temperature guess also had a large effect. By guessing 1000K across the rotor before the first iteration, 22,294 iterations were necessary, while using the average temperature of each row as the guess for that row before the first iteration required only 16,264 iterations.

    Error was somewhat erratic without relaxation, but with relaxation it was found that error was directly related to step size, and the size of the grid. It was found that the solution was first order accurate, and that as residual was decreased error exactly reached a value expected of a first order accurate system. The exact error was reached at a residual of E-6 for every step size.

    The study was highly limited due to the initial assumption that the brake rotor was a slab. Though it provided a benefit in simplicity for troubleshooting the initial version of the code, it undoubtedly compromised the quantative results. Other assumptions were limiting factors, such as constant thermal conductivity with temperature.



    Besides introducing the topic, the goal is to obtain a base numerical analysis in order to have a means of later comparison. Values such as maximum and minimum temperatures are not to be treated as indicative of those seen in real life this early in the development of the engineering code. Rather, they will be used as a benchmark to observe how much a given change alters the results. Until the system has been converted to a radial, rather than the current model which must be understood by the reader to be an oversimplified model, any quantitative results of this initial analysis are only for later comparison.

    The effects of changes in grid size and residual will be observed, as will the effects of a relaxation factor. Of particular interest will be the effects on the number of iterations, the time of the solution, the changes in maximum temperature, and the amount of error when compared to the conservation of energy. The temperature distribution of the model will also be analyzed.


    Brakes are used in automobiles to convert kinetic energy in to internal or thermal energy. In race cars, the purpose is to slow the car as quickly as possible to a speed that will enable the race vehicle to negotiate the looming corner. Due to the need to eliminate kinetic energy as quickly as possible, the levels of heat created can be unwieldy. It quickly becomes apparent that the braking system must be made to properly absorb all this heat, and improvements are sought constantly.

Figure 1 – LMP Brakes

    There are several ways to improve a braking system’s ability to absorb heat, some of which will be mentioned before the method investigated by this series is described. The two most obvious methods to improve brake performance from a heat standpoint are to create materials more capable of coping with the high temperatures encountered in racing, and improving the cooling of the brakes. Some other methods involve segregating heat from the braking fluid, as well as raising the boiling point of the brake fluid.

    This series will focus on the heat distribution through the brake rotor, and aim to reduce that temperature (particularly the maximum temperature) through geometry and airflow. Certainly, cooling can be increased through increased heat radiation, and other properties such as improved thermal conductivity, however these aspects are beyond the current scope of this series of writings.

    Papers written for this series will follow a common format. Generally, when a new paper is released in the series its purpose will be to eliminate a previous simplifying assumption, fix an error, or improve the speed and accuracy of the code. There are many assumptions in this initial model which severely limit the results, and despite the slow elimination of these assumptions they will be prevalent in the writings. Mistakes will also be common. Some posts will contain no direct changes, but instead may detail the effects things like of grid size, and speed of convergence.

    The purpose of this series is for the writer to start from a very basic and broad undergraduate engineering background, and slowly improve beyond the typical undergraduate simplified methods in a way that shows the reader the full process with as much transparency as possible. In no other conceivable way is this writing unique. The problems posed have already been solved by individuals with more experience in the field than I ever expect to achieve. The mathematical methods used already exist. In fact the only methods used that aren’t already well known will be fudge factoring and seat of the pants thinking. Even some of those methods though will likely already exist.

    The reader is cautioned to keep all these points in consideration when deciding how to apply the content of this series.


Physical model

    In order to simplify initial development and troubleshooting of the engineering code the physical model currently wasn’t a disk, but a giant block or slab. The cross sectional area of the block was the same as the cross section of the brake rotor. This cross section was then split in half at the midplane between the two planar faces of the rotor. This was done due to the assumed symmetry of the system, which is not the worst assumption and will remain with the code for the foreseeable future.

    This boundary of symmetry was treated as adiabatic. Such a boundary condition is common practice, since it reduces the size of the model, makes for simpler equations, and does not increase error. The boundary conditions on the other three faces allow for convection and radiation. Radiation was assumed to contribute enough to the results to be included, and it was not difficult to add. Finally, there exists the region on the outer planar face where the brake pad creates friction with the rotor and generates heat. It was assumed that no heat transmits through the brake pad, and that the brake pad does not block escaping heat.

Figure 2 – Physical Model

    The incoming heat was estimated from observations of the hairpin at Sebring International Raceway during the 12 Hours of Sebring 2012. It was estimated that LMP cars can decelerate at approximately 3 G’s, begin braking at 80m/s, and stop braking at 22.7m/s. The time taken, considering constant deceleration, was 1.95 seconds. It was assumed that 100% of the kinetic energy of the vehicle was converted to internal energy in the brake rotors. This amount of energy was then divided by ten, since it was estimated that the duty cycle of such braking was approximately 10%. Heat per rotor was then found by dividing the value by 4, which made the assumption that each brake rotor absorbed equal heat. In reality, the front rotors would encounter more heat than the rear rotors in general. It was decided that for this initial test, the incoming heat to the representative block would be 16,975 Watts. It was assumed that the heat generation would be constant across the contact area.

    The brake rotor was discretized with equal step sizes in both directions (see Fig. 3) in the cross sectional area of the brake rotor, with a 1m depth for simplicity. The vertical steps symbolized infinitesimal steps in radius, and horizontal steps symbolized infinitesimal steps in the axial direction. This grid of nodes was represented in Matlab as a matrix, with i and j representing the location of the node, and the value in the ith column and jth row represented the temperature at that node.

Figure 3 – Discretization

    The math model was based on the conservation of energy. Since the model was treated as steady state, the sum of energy entering and exiting any region of the brake rotor must equal zero.

Results and discussion

    The effectiveness of the relaxation factor was profound. The longest duration analysis that was run with and without a relaxation factor was a step size of .00025m (this is why most figures are for data collected with a step size of .00025), which resulted in a grid size of 40,000 nodes. For a residual of E-5 (which is the maximum residual for non relaxation factor runs) the solution took 450s and 188,642 iterations to converge (see Fig. 4). On the other hand, with an optimized relaxation factor the solution took 35s (almost 13 times less time) and 9,905 iterations (19 times fewer iterations). 419 iterations per second were achieved without the relaxation factor, whereas only 283 iterations a second were achieved with the relaxation factor. This makes sense, since including the relaxation factor increases the complexity of the equations, and therefore more calculations per iteration. However, due to the vast improvement in iteration count to convergence, using the relaxation factor clearly improves the speed of convergence.



Figure 4 – Iterations as a function of Residual, taken at .00025m step size


Figure 5 – Iterations as a function of step size, taken at E-5 Residual

    Interestingly there was an improvement in error when using a relaxation factor for moderate residuals. However as the residual was reduced further, the error increased whereas under the same conditions the error decreased without the relaxation factor. It was found that the solution was first order accurate; meaning that halving the step size had the effect of halving the error. This was a desirable outcome since it validated the numerical method used was functional and none of the long equations used were accidently mistyped.


Figure 6 – Error as a function of Residual, taken at step size of .00025m


Figure 7 – Error as a function of step size, taken at maximum Residual run (E-5 to E-7)

    The maximum temperature behaved in general as was expected. With large step sizes, the maximum temperature was high, and with smaller step sizes the maximum temperature was reduced. The absolute maximum temperature predicted was 1,166K and the minimum predicted by the code was 1,158.7. This made for a deviation of only .6% which is pretty good considering that there were two methods of solving (relaxation and without) as well as a step size range of .001 to .0001, and a residual range from E-3 to E-10. Figure 8 shows the maximum temperature as related by the step size, and each point was taken at maximum residual for that run, which was not constant. Figure 9 illustrates the change in maximum temperature prediction with respect to the residual.


Figure 8 – Maximum temperature as a function of step size, taken at minumum residual


Figure 9 – Maximum temperature as a function of residual, taken at .00025m step size.

The temperature distribution predicted by the code was as expected. The centerline shows the adiabatic boundary condition, where the rate of change in temperature with respect to movement in the horizontal was zero. The maximum temperature was found near the center of the pad. Outside of the direct influence of the heat flux from the brake pads, it can be seen that convection on the planar faces of the rotor is providing a cooling effect. This can be seen most notably in the lower corners of Figure 10.

Figure 10 – Temperature distribution of the slab cut in cross section

    The initial temperature guess had a large effect on the number of iterations, as expected. Three methods were used to provide the initial temperature distribution guess. First (Fig. 11) used a simple 1000K guess for the whole rotor, which resulted in 22,294 iterations. Second, the average temperature from a previous run was applied across the entire rotor (Fig. 12), which took 21,092 iterations to converge. Last, the average temperature at each vertical location (representing changes in radius) was distributed across each row. This method resulted in only 16,264 iterations (Fig. 13). It was found that good prediction of temperature distribution can reduce iterations substantially. In this case, the iterations were reduced by 27%.

Figure 11 – Number of iterations resulting from 1000K initial guess

Figure 12 – Number of iterations resulting from Tave initial guess

Figure 13 – Number of iterations from the vertical gradient of average temperatures for each row    

Conclusion and recommendations

    It was concluded that the code itself is functional and capable of producing the results of the values of interest. Methods to reduce error and run time were observed, and will be implemented in to future code updates. Though the values predicted by the code are likely not indicative of real life values, they will still provide a benchmark for the near future when observing changes produced by differences in duty cycle, or by making material properties such as thermal conductivity vary with temperature as they do in real life.

    Limits of this writing are numerous. First and foremost is the fact that the brake rotor for the time being is treated as a block or slab rather than a cylindrical body. Still another is the solution time for runs without a relaxation factor prohibited gathering data at step sizes smaller than .00025m as well as residuals beyond E-5. This made comparisons between relaxation and non relaxation less than ideal.

    Future treatments of this subject should revise the geometry in order to properly represent the cylindrical nature of brake rotors. Time could be spent running the long simulations necessary to gather all data for non relaxation solutions, however since the model has so many other shortcomings this should be low on the priority list. While the model is this simple, making thermal conductivity and emissivity variable with temperature could be easy changes to implement in order to observe the effects on the slab.


Figure 14 – Relaxation data

Figure 15 – No relaxation data


    Thanks for all the input from the forums at



[1] Incropera et all, Fundamentals of Heat and Mass Transfer 6th Ed., Wiley

[2] C. Gerald, P Wheatley, Applied Numerical Analysis 7th Ed., Pearson Addison Wesley

[3] D. Pitts, L. Sissom, Heat Transfer 2nd Ed., Schaum’s Outlines, McGraw Hill

[4] R.W. Hamming, Numerical Methods for Scientists and Engineers 2nd Ed., Dover Publications

[5] F. Talati, S. Jalalifar, Investigation of heat transfer phenomena in a ventilated disk brake rotor with straight radial rounded vanes, Journal of Applied Sciences 8 (20) (2008) 3583-3592

[6] J. Zhang, Numerical Accuracy,,

Preliminary – Numerical Analysis of Racing Brakes

I have decided to release a preliminary post detailing some of the goals and anticipated main sources of my Numerical Analysis series. I am still working on a format for the series of posts so they can be clearer. The topic itself is complicated enough that I prefer to not add to any confusion through poor formatting and inconsistency between posts.

First, I want to make something very clear. I am an undergraduate student getting a B.S. in Mechanical Engineering with a B.S. Minor in Aerospace Engineering. I am not an expert on numerical analysis or brakes. My posts are not intended to represent the cutting edge of computer aided engineering. The intended overall goal of my posts is first and foremost to learn numerical analysis and how to write about it. My second intended goal is to document my path, complete with all the mistakes I will certainly make, in order to help the reader see what goes in to this.

Every post I make will be of another step in an incomplete project. Every time I update the Matlab code I wrote there will be simplifying assumptions that reduce the accuracy of the simulation, sometimes including wrong assumptions and mistakes on top of that. I absolutely welcome and request your criticism; however I hope you keep these things in mind.

Goals for the software (some of which may never be completed):

  • Model simple solid brake rotors operating in steady state
  • Model realistic (for example vented) brake rotors in steady state
  • Eliminate as many assumptions as possible
  • Model simple solid brake rotors operating in transient state
  • Model realistic (for example vented) brake rotors in transient state
  • Optimize brake rotor cooling in steady state
  • Optimize brake rotor cooling in transient state
  • Account for different material compositions
  • User friendly interface

My main sources may include:

Fundamentals of Heat and Mass Transfer 6th (sixth) Edition

Schaum’s Outline of Heat Transfer, 2nd Edition)

Numerical Methods for Scientists and Engineers

Applied Numerical Analysis (7th Edition)

Introduction to Linear Algebra, Fourth Edition


Once again, I welcome all constructive comments, recommendations, and criticism. You can always email me but I prefer to receive your input through comments so everyone can see any questions posed and mistakes pointed out. The idea is for this to be as transparent to everyone as possible.