Abstract

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.

Summary

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.

Introduction

Objectives

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.

Background

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.

Model

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

note.

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 %}

y=1;x=1;

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);

end

for y=2:sizey

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

end

for x = 2:sizex-1

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

end

x=1;

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.

Acknowledgements

Thanks for all the input from the forums at F1Technical.net.

References

[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, tamu.edu, https://ceprofs.civil.tamu.edu/jzhang/cven-302/chap04b.ppt

Appendices

Solutions for each boundary condition type

Boundary condition types

BC1

BC2

BC3

BC7

BC8

BC9

BC10

BC11