Topic: Newton’s Law of Cooling and Euler’s Method to IVP ODE
Subject: Physics and Numerical Methods
Tool: Scilab
This article presents a solution as to how long it will take for an object to cool down to a specified temperature and thereby to ambient temperature. The object is being taken out of the oven at \(230^\circ C\) and put in place where the ambient temperature is \(30^\circ C\). Neglecting the heat transfer due to conduction where the object is in contact with, the time to cool it down to a particular temperature due to natural convection and radiation can be estimated using the ordinary differential equation (ODE) representing the Newton’s Law of Cooling shown in Equation (1).
$$\frac{dT}{dt}\;=\;-k \left(T(t) - T_{env} \right) \tag{1}$$
Where:
The analytical solution to the ODE
representing Newton’s law of cooling can be found by integrating both sides of
the equation as shown in Equation (2).
$$\int_{T_{0}}^{T_{t}} dT\;=\;-k \int_{0}^{t} \left(T(t)-T_{env} \right) dt \tag{2}$$
For this article, a numerical solution
using the Euler method of approximating the solution of a first-order differential
equation is presented with the use of Scilab. Euler methods require an ODE to
be in the form of Equation (3).
$$\frac{dT}{dt}\;=\;f\left(t,T\right) \tag{3}$$
With the given initial condition,
$$\left(t_{0},\;T_{0}\right) \tag{4}$$
The approximate solution using Euler's method is shown in Equation (5).
$$\frac{\Delta T}{\Delta t}\; \approx \; \frac{T_{n+1}-T_{n}}{h}\; \approx \; f\left(t,T \right) \tag{5}$$
Or by rearranging and solving for \(T_{n+1}\),
$$T_{n+1}\;=\;T_{n}+h \;f\left(t,T\right) \tag{6}$$
The \( f(t,T)\) requires finding the cooling
constant, \(k\), for a specific material. If this property seems to be unavailable or difficult to find, experimental data can be used as the basis for
calculating the \(k\) using Equation (7).
$$k\;=\;- \frac{1}{t} \ln \left[ \frac{T(t) - T_{env}}{T_{0} - T_{env}} \right] \tag{7}$$
Figure (1) is the illustration of the
condition to be given a solution using Newton’s law of cooling by applying
Euler's method of solving ODE.
Figure 1. An object at \(230^\circ C\) exposed in an environment at \(30^\circ C\).
To provide a simulation of the given
condition in Figure (1), experimental data were gathered to calculate the cooling
constant, \(k\), from the given temperature and time. The cooling constant, \(k =
0.08 \;\text{per min}\), was calculated from the given parameters and Scilab script.
// Copyright (C) 2024 - Gani Comia
// Experimental Cooling Constant
clear;clc;
// (1) defining cooling constant function
function k=f(t, Tt, T0, Tenv)
// where:
// t - min, time consideration
// Tt - deg C, temperature at time t
// T0 - deg C, initial temperature
// Tenv - deg C, environment temperature
k = -(1/t).*log((Tt - Tenv)./(T0 - Tenv));
endfunction
// (2) experimental data
t = 30; // min, time between T0 and Tt
T0 = 230; // deg C, initial temp
Tt = 50; // deg C, temp at time t
Tenv = 30; // deg C, ambient temp
// (3) cooling constant calculation
k = f(t,Tt,T0,Tenv)
mprintf("Cooling Constant, k = %4.2f per min\n", k)
The first block of the script defines the
physical parameter under consideration.
clear;clc;
// (1) physical parameters
T0 = 230; // deg C, initial temp of the object
Tenv = 30; // deg C, ambient temp
k = 0.08; // per min, experimental cooling constant
dt = 0.05; // min, time step
totalTime = 80; // min, total simulation timeThe time domain in the simulation is
divided into smaller step sizes for accuracy of solution. Initial time and
temperature are given on the next block of script.
// (2) initial conditions and derived parameters
time = 0:dt:totalTime;
nSteps = length(time);
T = zeros(nSteps, 1);
T(1) = T0;Equation (6), which is the numerical
solution to the ODE using the Euler method, is transformed in the Scilab script. It uses the \(for-loop\) function to obtain the solution over the desired range.
// (3) numerical solution to ODE using Euler Method
for i = 1:nSteps-1
dTdt = -k * (T(i) - Tenv); // deg C/min, rate of temp change
T(i+1) = T(i) + dTdt * dt; // deg C, temp at the next time step
endThe Scilab \(plot()\) function is used to
visualize the solution from the Euler method of approximation. Finding the
temperature for a given time can be estimated using interpolation by the use of the \(interp1()\) function.
// (4) graph for visualization of solution
clf;
plot(time, T, "b","linewidth",2.5);
xlabel("Time, t [ min ]","fontsize",2);
ylabel("Temperature, T [ °C ]","fontsize",2);
title("Newton''s Law of Cooling Simulation","fontsize",3);
note = ["T0="+string(T0)+"°C, Tenv="+string(Tenv)+"°C, k="+string(k)]
legend(note,with_box = %F);
xgrid(3,1);
// (5) interpolating temp, Tx, for the particular time, tx
format(6)
tx = 17.5; // min, sample time
Tx = interp1(time, T, tx); // deg C, interpolated temp
disp(tx, Tx)
// (6) plotting lines for tx and Tx
xptVer = [tx tx];
yptVer = [0 Tx];
plot(xptVer, yptVer, "r--","linewidth",1.5)
xptHor = [0 tx];
yptHor = [Tx Tx];
plot(xptHor, yptHor, "r--","linewidth",1.5)
plot(tx,Tx,"ro");
note1 = ["t (min) ="+string(tx)+"; T (°C) ="+string(Tx)];
xstring(tx,Tx,note1);
xstring(45,200,["https://gani-mech-toolbox.blogspot.com"])
eqn1=["$\LARGE\frac{dT}{dt}\;=\;{-k}\left(T-T_\mathbf{{env}}\right)$"];
xstring(45,166,eqn1);
a = gca();
a.children.children.mark_background = 9;
This Scilab toolbox is a powerful tool
for visualizing solution data and calculating specific conditions.
Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.
