8.00: Physical Problems for Ordinary Differential Equations

Lesson 1: General Engineering

Summary

To find the temperature of a heated ball as a function of time, a first-order ordinary differential equation must be solved.

Figure $\PageIndex{1.1}$. A heated ball is cooled to a specified temperature in a quenching bath.

Modeling

You are working for a ball bearing company - Ralph’s Bearings. For a long lasting life of some spherical bearings made in this company, after heating them to a high temperature of $1000\ \text{K}$, they need to be quenched for $30$ seconds in water that is maintained at a room temperature of $300\ \text{K}$. However, it takes time to take the ball from the furnace to the quenching bath, and its temperature falls. The temperature of the furnace is $1200\ \text{K}$, and it takes $10$ seconds to take the ball to the quenching bath. Does the ball need to be in the air for less or more than $10$ seconds to get to a temperature of $1000\ \text{K}$?

Mathematical model

To keep the mathematical simple, one can assume the spherical bearing to be a lumped mass system. What does a lumped system mean? It implies that the internal conduction in the sphere is large enough that the temperature throughout the ball is uniform. This assumption allows us to assume that the temperature is only a function of time and not of the location in the spherical ball. This dependence of temperature only on time for a lumped system means that if a differential equation governs this physical problem, it would be an ordinary differential equation. For a nonlumped system, a partial differential equation would govern. In your heat transfer course, you will learn when a system can be considered lumped or nonlumped. In simplistic terms, this distinction is based on the material, geometry, and heat exchange factors of the ball with its surroundings.

Now assuming a lumped-mass system, let us develop the mathematical model for the above problem. When the ball is taken out of the furnace at the initial temperature of $\theta_{0}$ and is cooled by radiation to its surroundings at the temperature of $\theta_{a}$, the rate at which heat is lost to radiation is

$$\text{Rate of heat lost due to radiation} = A \epsilon \sigma\left( \theta^{4} - {\theta_{a}}^{4} \right)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.1}) \nonumber$$

where

$$A = \text{surface area of the ball,}\ \text{m}^{2} \nonumber$$

$$\epsilon = \text{emittance} \nonumber$$

$$\sigma = \text{Stefan-Boltzmann constant,}\ 5.67 \times 10^{- 8}\ \frac{\text{J}}{\text{s}\cdot \text{m}^{2}\cdot \text{K}^{4}} \nonumber$$

$$\theta = \text{temperature of the ball at a given time,}\ \text{K} \nonumber$$

In the above equation for radiation heat loss, emittance $\epsilon$ is defined as the total radiation emitted divided by total radiation that would be emitted by a blackbody at the same temperature. The emittance is always between $0$ and $1$. A black body is a body that emits and absorbs at any temperature the maximum possible amount of radiation at any given wavelength.

The Stefan-Botzmann constant, $\sigma$ was discovered by two Austrian scientists – J. Stefan and L. Boltzmann. Stefan found it experimentally in 1879, and Boltzmann derived it theoretically in 1884.

The energy stored in the mass is given by

$$\text{Energy stored by mass} = mC\theta\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.2}) \nonumber$$

where

$$m = \text{mass of the ball,}\ \text{kg} \nonumber$$

$$C = \text{specific heat of the ball,}\ \text{J}/\left( \text{kg} \cdot \text{K} \right) \nonumber$$

From an energy balance,

$$\text{Rate at which heat is gained} - \text{Rate at which heat is lost} = \text{Rate at which heat is stored} \nonumber$$

gives

$$- A \epsilon \sigma\left( \theta^{4} - {\theta_{a}}^{4} \right) = mC\frac{d \theta}{dt}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.3}) \nonumber$$

Given that

$$\text{Radius of the ball,}\ r = 2.0\ \text{cm} \nonumber$$

$$\text{Density of the ball,}\ \rho = 7800\ \text{kg/m}^{3} \nonumber$$

$$\text{Specific heat,}\ C = 420\ \text{J}/(\text{kg}\cdot \text{K}) \nonumber$$

$$\text{Emittance,}\ \epsilon = 0.85 \nonumber$$

$$\text{Stefan-Boltzmann Constant,}\ \sigma = 5.67 \times 10^{-8}\ \frac{\text{J}}{\text{s}\cdot \text{m}^{2}\cdot \text{K}^{4}} \nonumber$$

$$\text{Initial temperature of the ball,}\ \theta\left( 0 \right) = 1200\ \text{K}, \nonumber$$

$$\text{Ambient temperature,}\ \theta_{a} = 300\ \text{K}, \nonumber$$

we have

$$\begin{split} \text{Surface area of the ball: } A &= 4 \pi r^{2} \\ &= 4 \pi (0.02)^{2} \\ &= 5.02654 \times 10^{-3} \ \text{m}^{2} \end{split}$$

$$\begin{split} \text{Mass of the ball: } M &= pV \\ &= \rho \left[ \frac{4}{3} \pi r^{3} \right] \\ &= 7800 \times \left( \frac{4}{3} \right) \pi (0.02)^{3} \\ &= 0.261380 \ \text{kg} \end{split}$$

Hence

$$- A \epsilon \sigma(\theta^{4} - {\theta_{a}}^{4}) = mC\frac{d \theta}{dt}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.3} \text{ - repeated}) \nonumber$$

reduces to

$$- (5.02654 \times 10^{- 3})(0.85)(5.67 \times 10^{- 8})(\theta^{4} - 300^{4}) = (0.261380)(420)\frac{d \theta}{dt} \nonumber$$

$$\frac{d\theta}{dt} = - 2.20673 \times 10^{- 12}\ (\theta^{4} - 81 \times 10^{8})\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.4}) \nonumber$$

An improved mathematical model

The heat from the ball can also be lost due to convection. The rate of heat lost due to convection is

$$\text{Rate of heat lost due to convection} = hA\left( \theta - \theta_{a} \right),\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.5}) \nonumber$$

where

$$h = \text{the convective cooling coefficient},\ \left\lbrack \text{W}/\left( \text{m}^{2} \cdot \text{K} \right) \right\rbrack. \nonumber$$

Hence the heat is lost is due to convection as well as radiation and is given by

$$\text{Rate of heat lost due to convection and radiation} = A \epsilon \sigma\left( \theta^{4} - \theta_{a}^{4} \right) + hA\left( \theta - \theta_{a} \right) \nonumber$$

$$\text{Energy stored by mass} = mC\theta\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.6}) \nonumber$$

From an energy balance,

$$\text{Rate at which heat is gained} - \text{Rate at which heat is lost} = \text{Rate at which heat is stored} \nonumber$$

gives

$$- A \epsilon \sigma\left( \theta^{4} - {\theta_{a}}^{4} \right) - hA\left( \theta - \theta_{a} \right) = mC\frac{d \theta}{dt}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.7}) \nonumber$$

Given $h = 350 \ \frac{\text{J}}{\text{s} \cdot \text{m}^{2} \cdot \text{K}}$ and considering both convection and radiation, and substituting the values of the constants given before, we have

$$- (5.02654 \times 10^{- 3})(0.85)(5.67 \times 10^{- 8})(\theta^{4} - 300^{4}) \nonumber$$

$$- (350)(5.02654 \times 10^{- 3})(\theta - 300) = (0.261380)(420) \frac{d\theta}{dt} \nonumber$$

$$\frac{d \theta}{dt} = - 2.20673 \times 10^{- 13}(\theta^{4} - 81 \times 10^{8}) - 1.60256 \times 10^{- 2}(\theta - 300)\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{1.8}) \nonumber$$

The solution to the above ordinary differential equation with the initial condition of $\theta\left( 0 \right) = 1200\ \text{K}$ would give us the temperature of the ball as a function of time. We can then find the time at which the ball temperature drops to $1000\ \text{K}$.

Questions

(1) Note that the ordinary differential equation given by Equation $(\PageIndex{1.8})$ is non-linear. Is there is an exact solution to the problem?

(2) You are asked to solve the inverse problem, that is, when the dependent variable temperature will be equal to $1000\ \text{K}$. How would you solve the inverse problem using different numerical methods such as Euler’s and Runge-Kutta’s methods?

(3) Can you find if convection or radiation can be neglected? How would you quantify the effect of neglecting one or the other?

(4) Find the following at $t = 5\text{ s}$

1. Rate of change of temperature,
2. Rate at which heat is lost due to convection,
3. Rate at which heat is lost due to radiation,
4. Rate at which heat is stored in the ball.

Lesson 2: Mechanical Engineering

Summary

A physical problem of finding how much time it would take a trunnion to cool down in a refrigerated chamber is presented. To estimate this time, the problem would be modeled as an ordinary differential equation.

Modeling

To make the fulcrum (Figure $\PageIndex{2.1}$) of a bascule bridge, a long hollow steel shaft called the trunnion is shrunk-fit into a steel hub. The resulting steel trunnion-hub assembly is then shrunk-fit into the girder of the bridge.

Figure $\PageIndex{2.1}$. Trunnion-Hub-Girder (THG) assembly.

This assembly is done by first immersing the trunnion in a cold medium such as liquid nitrogen. After the trunnion reaches the steady-state temperature, that is, the temperature of the cold medium, the trunnion outer diameter contracts. The trunnion is taken out of the medium and slid though the hole of the hub (Figure $\PageIndex{2.1})$.

When the trunnion heats up, it expands and creates an interference fit with the hub. Because of the large thermal shock that the trunnion undergoes, it is suggested that the trunnion be cooled in stages. First, put the trunnion in a refrigerated chamber, then dip it in a dry-ice/alcohol mixture and then immerse it in a bath of liquid nitrogen. However, this approach will take more time. One is mainly concerned about the cooling time in the refrigerated chamber. Assuming that the room temperature is $27^{\circ} \text{C}$ and the refrigerated chamber is set is $-33^{\circ} \text{C}$, how much time would it take for the trunnion to cool down to $-33^{\circ} \text{C}$?

Let us do a simplified problem of a solid trunnion and assume the trunnion can be treated as a lumped-mass system. What does a lumped system mean?

Let us develop the mathematical model for the above problem. When the trunnion is placed in the refrigerated chamber, the trunnion loses heat to its surroundings by convection.

Figure $\PageIndex{2.2}$. Trunnion slides through the hub after contracting

$$\text{Rate of heat lost due to convection} = h\left( \theta \right)A\left( \theta - \theta_{a} \right). \;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.1}) \nonumber$$

where

$$h\left( \theta \right) = \text{the convective cooling coefficient,}\ \text{W/m}^{2}\cdot{\text{K}}\ \text{and is a function of temperature} \nonumber$$

$$A =\text{surface area} \nonumber$$

$$\theta_{a} =\text{ambient temperature of the refrigerated chamber} \nonumber$$

The energy stored in the mass is given by

$$\text{Energy stored by mass} = mC\theta\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.2}) \nonumber$$

where

$$m = \text{mass of the trunnion,}\ \text{kg} \nonumber$$

$$C = \text{specific heat of the trunnion,}\ \text{J/(kg} \cdot \text{K)} \nonumber$$

From an energy balance,

$$\text{Rate at which heat is gained} - \text{Rate at which heat is lost} = \text{Rate at which heat is stored} \nonumber$$

gives

$$- h\left( \theta \right)A\left( \theta - \theta_{a} \right) = mC\frac{d \theta}{dt}\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.3}) \nonumber$$

Note that the convective cooling coefficient is a function of temperature. Other material parameters such as density (affecting mass), specific heat, and thermal conductivity (affecting whether the system can be considered lumped or not) of the trunnion material are functions of temperature as well, but within our temperature range of $27^{\circ} \text{C}$ to $-33^{\circ} \text{C}$ these vary by only $10\%$ from the room temperature values. So we assume these parameters to be constant.

Let us now determine the constants needed for the above ordinary differential equation.

Convection coefficient of air

The convection coefficient of air, $h\left( \theta \right)$ is not a constant function of temperature and is given by

$$h = \frac{Nu \times k}{D}\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.4}) \nonumber$$

where

$${Nu}\ \text{is the Nusselt number,} \nonumber$$

$$k\ \text{is the thermal conductivity of air,} \nonumber$$

$$D\ \text{is the characteristic diameter and is taken as the outer diameter of the trunnion.} \nonumber$$

The Nusselt number for a vertical cylinder is given by the empirical formula^[1]:

$$Nu = \left( 0.825 + \frac{0.387Ra^{\frac{1}{6}}}{\left(1+\left(\frac{0.492}{\Pr} \right)^{\frac{9}{16}} \right)^{\frac{8}{27}}}\right)^2\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.5}) \nonumber$$

where

$Pr$ is the Prandtl number given by:

$$Pr = \frac{v_{k}}{\alpha}\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.6}) \nonumber$$

$$\nu_{k}\ \text{is the kinematic viscosity of the fluid,} \nonumber$$

$$\alpha\ \text{is the thermal diffusivity,} \nonumber$$

$Ra$ is the Rayleigh number given by:

$$Ra = Gr \times Pr\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.7}) \nonumber$$

$Gr$ is the Grashoff number given by:

$$Gr = \frac{g\beta(T_{\text{wall}} - T_{\text{fluid}}) D^{3}}{{\nu_{k}}^2}\;\;\;\;\;\;\;\;\;\;\;\; (\PageIndex{2.8}) \nonumber$$

where

$$g\ \text{is the gravitational constant,} \nonumber$$

$$\beta\ \text{is the coefficient of volumetric thermal expansion,} \nonumber$$

$$T_{\text{wall}}\ \text{is the temperature of the wall,} \nonumber$$

$$T_{\text{fluid}}\ \text{is the temperature of the fluid.} \nonumber$$

To calculate the convection coefficient, we use the following values for kinematics viscosity $\nu$, thermal conductivity $k$, thermal diffusivity $\alpha$, and volumetric coefficient $\beta$ of air as a function of temperature, $\theta$. [Ref. 1]

Table $\PageIndex{2.1}$. Properties of air as a function of temperature.

$T_{lookup}$	$\nu$	$k$	$\alpha$	$\beta$
$^{\circ} \text{C}$	$\text{m}^2/\text{s}$	$\text{W}/(\text{m} \cdot \text{K})$	$\text{m}^2/\text{s}$	$1/\text{K}$
$-173.15$	$2.00 \times 10^{- 6}$	$9.34 \times 10^{- 3}$	$2.54 \times 10^{- 6}$	$1.00 \times 10^{- 2}$
$-123.15$	$4.43 \times 10^{- 6}$	$1.38 \times 10^{- 2}$	$5.84 \times 10^{- 6}$	$6.67 \times 10^{- 3}$
$-73.15$	$7.59 \times 10^{- 6}$	$1.81 \times 10^{- 2}$	$1.03 \times 10^{- 5}$	$5.00 \times 10^{- 3}$
$-23.15$	$1.14 \times 10^{- 5}$	$2.23 \times 10^{- 2}$	$1.59 \times 10^{- 5}$	$4.00 \times 10^{- 3}$
$26.85$	$1.59 \times 10^{- 5}$	$2.63 \times 10^{- 2}$	$2.25 \times 10^{- 5}$	$3.33 \times 10^{- 3}$

Other constants needed for this calculation are:

$$D = 0.25\ \text{m} \nonumber$$

$$g = 9.8\ \text{m/s}^{2} \nonumber$$

$$T_{\text{fluid}} = - 33^{\circ}\text{C} \nonumber$$

Table $\PageIndex{2.2}$. Convection coefficient of air as a function of temperature.

$T_{wall}$	$T_{\text{average}}=\frac{1}{2} \left( T_{\text{fluid}} + T_{\text{wall}} \right)$	$h$
$^{\circ} \text{C}$	$^{\circ} \text{C}$	\(\text{W{/(m^2.K)\)
$-33$	$-16.50$	$5.846 \times 10^{-2}$
$-18$	$-9.00$	$4.527$
$-8$	$-4.00$	$5.214$
$2$	$1.00$	$5.702$
$27$	$13.50$	$6.553$

The values in Table $\PageIndex{2.2}$ are interpolated from Table $\PageIndex{2.1}$ where temperatures are chosen as the average value for the wall and ambient temperature. The above data is interpolated as:

$$h\left( \theta \right) = - 3.69 \times 10^{- 6}\theta^{4} + 2.33 \times 10^{- 5}\theta^{3} + 1.35 \times 10^{- 3}\theta^{2} + 5.42 \times 10^{- 2}\theta + 5.59\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.9}) \nonumber$$

$$\text{Outer radius of trunnion,}\ a = 0.125\ \text{m} \nonumber$$

$$\text{Length of trunnion,}\ L = 1.36\ \text{m} \nonumber$$

$$\begin{split} \text{Surface area,} \ A &= (2\pi a)L + 2\pi a^{2}\\ &= 2\pi\left( .125 \right)\left( 1.36 \right) + 2\pi(0.125)^{2}\\ &= 1.166\ \text{m}^{2} \end{split}$$

$$\text{Length of the trunnion,}\ L = 1.36\ \text{m} \nonumber$$

$$\text{Density of trunnion material,}\ \rho = 7800 \text{kg/m}^3 \nonumber$$

$$\begin{split} \text{Mass,} \ m &= {\rho V}\\ &=\rho(\pi r^{2}L)\\ &=7800\pi\left( 0.125 \right)^{2}1.36\\ &= 520.72\ \text{kg} \end{split}$$

Other constants:

$$\text{Specific heat,}\ C = 420\ \text{J}/ (\text{kg} \cdot \text{K}), \nonumber$$

$$\text{Initial temperature of the trunnion,}\ T\left( 0 \right) = 27^{\circ} \text{C}, \nonumber$$

$$\text{Ambient temperature,}\ T_{a} = - 33{^\circ}\text{C} \nonumber$$

The first-order ordinary differential equation

$$- h\left( \theta \right)A\left( \theta - \theta_{a} \right) = mC\frac{d\theta}{dt} \nonumber$$

reduces to

$$- \left( - 3.69 \times 10^{- 6}\theta^{4} + 2.33 \times 10^{- 5}\theta^{3} + 1.35 \times 10^{- 3}\theta^{2} + 5.42 \times 10^{- 2}\theta + 5.588 \right) \times \left( 1.166 \right) \times \left( \theta + 33 \right) = \left( 520.72 \right) \times \left( 420 \right) \times \frac{d\theta}{dt} \nonumber$$

$$\frac{d\theta}{dt} = - 5.331 \times 10^{- 6}( - 3.69 \times 10^{- 6}\theta^{4} + 2.33 \times 10^{- 5}\theta^{3} + 1.35 \times 10^{- 3}\theta^{2} + 5.42 \times 10^{-2}\theta+5.588)(\theta+33)\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.10}) \nonumber$$

$$\theta\left( 0 \right) = 27{^\circ}\text{C} \nonumber$$

Is the assumption of the trunnion considered as a lumped system correct?

To determine whether a system is lumped, we calculate the Biot number $Bi$, which is defined as

$$Bi = \frac{hL}{k_{s}}\;\;\;\;\;\;\;\;\;\;\;\;(\PageIndex{2.11}) \nonumber$$

where

$$h = \text{average surface conductance,} \nonumber$$

$$L = \text{significant length dimension (volume of body/surface area),} \nonumber$$

$$k_{s} = \text{thermal conductivity of solid body.} \nonumber$$

If $B_{i} < 0.1$, the temperature in the body is uniform within $5\%$ error.

In our case:

$$h = 4.407\ \frac{\text{W}}{\text{m}^{2}\cdot{\text{K}}} \nonumber$$

$$L = 1.36\ \text{m} \nonumber$$

$$k_{s} = 81\ \frac{\text{W}}{\text{m}\cdot{\text{K}}}\;\;\;\;\;\;\;\;\;\;\;\;[\text{Ref.} \ 2] \nonumber$$

$$\begin{split} Bi &= \frac{4.407 \times 1.36}{81}\\ &= 0.074 < 0.1 \end{split}$$

The Biot number is less than $0.1$, and one can hence assume the trunnion to be a lumped mass system.

References

[1] F. B. Incropera, D. P. Dewitt, Introduction to Heat Transfer, 3rd Ed, John Wiley & Sons, Inc., New York, NY, 2000, Chap. 9.

[2] F. Kreith, M. S. Bohn, Principle of Heat Transfer, 4th Ed, Harper & Row, Publishers, New York, NY, 1993, Appendix 2.