OPEN ACCESS
This paper describes several simulations of the seawater cooling phenomenon using a theoretical approach based on analytical method, numerical method of Euler’s and RungeKutta of fourthorder (RK4). These research objectives, i.e. to get the energy balance equations, to use the solution equations, and to do the simulation processes. The methods used, i.e. (i) do the completion of mathematical equations to get the constants for energy balance based on the ordinary differential equations (ODEs), (ii) do the made of solving equations for simulation, and (iii) does the simulation processes assisted by a spreadsheet application and result in analysis. The results are in the form of (a) the constants for energy balance, i.e. b1 is 0.9134∙10^{3} sec^{1} and b2 is 0.31∙10^{6} sec^{1}, (b) produced three solving equations for the simulation, and (c) obtained curves of the temperature changes. The results of the simulation processes based on a spreadsheet application have obtained the results, that the time span 1800 seconds (30 minutes) can cause changes in the temperature of the fluid from 25℃ to (i) 4.823770℃ with a rate of change of 0.092%, if used of the analytical method; (ii) 4.819088℃ with a rate of change of 0.092%, if used the numerical method of Euler’s; and (iii) 5.600404℃ with a rate of change of 0.083%, if used the numerical method of RK4’s. The conclusion in this paper, that all the curves of changes in the temperature of the fluid are the nonlinear curves. Even though the final value of temperature is the highest, but the RK4 more thoroughly. Suggestions for future work, that the simulation of its phenomena must be begun with making the mathematical models with the analytical and/or numerical method and implemented them into the computer application.
analytical method, numerical methods of Euler’s and RK4’s, ordinary differential equations, simulation on the seawater cooling phenomena, single rectangular platefin
The closedform or analytic solutions to mathematical problems and applications known as actual or exact solutions [14]. All mathematical problems cannot be solved by exactly, because the completion processes are often quite complicated and take a lot of time, making its inefficient [5]. In this case, most of the forms are lowdegree polynomials which are solved by the exact solution [6], in the other words, analytic solutions are not easy to find [7, 8]. The limitations of analytical solutions contained in these mathematical problems solved by numeric, although no exact value generated, at least the expectation values has been approached [713]. The nature of numerical solutions is only in the form of approximation solutions, but errors from approximation solutions can always calculated [2, 7, 1417]. A condition where there are several conditions and methods used for obtaining a good settlement result, so if the problem has been a problem that has a high complexity, so even numerical solutions cannot with a good solution, then the simulation method must be used [2, 7, 1620].
The cooling process of solid or fluid or solid to fluid is a physical operation and there is a temperature difference and the heat has removed [2023]. The state of the cooling process is a form of the phenomenon that has explained by the implementation of Newton's Law [11, 19, 2025]. Newton’s Law described that the rate of change of the temperature of an object is proportional to the difference between its temperature and the ambient temperature of its surroundings [2123]. The statement about an instantaneous rate of change of the temperature in Newton's Law will then be a function equation that tracks the complete record of the temperature over time [18, 24, 25]. It will see that when it translates this verbal statement into a differential equation. In the other sentence, the object’s temperature can be modeled to a differential equation. It is arriving at a differential equation as the solution to the cooling process phenomena [2125].
The phenomena of a seawater cooling process are also with a temperature difference and a removing the heat. The temperature difference that is completely recorded over time and the heat that has removed from the process is a physical operation [11, 2125] that has calculated using the ordinary differential equations (ODEs) [18, 20, 2628]. The general form of ODEs is an equation with variables $x, y$, and derivatives of $y$ to $x$ that has written with the Eq. (1) [2628] while the form of linear equations with ordern has written as the Eq. (2) [20, 2628].
$F\left(x, y, \frac{d}{d x} y, \frac{d^{2}}{d x^{2}} y, \cdots, \frac{d^{n}}{d x^{n}} y\right)=0$ (1)
$a_{0}(x) \frac{d^{n}}{d x^{n}} y+a_{1}(x) \frac{d^{n1}}{d x^{n1}} y+\cdots+a_{n1}(x) \frac{d}{d x} y$ $+a_{n}(x) y=f(x)$ (2)
where, $a_{0}, a_{1}, \cdots, a_{n1}, a_{n}$, and $f$ is a free variable function of $x$ and $a_{0} \neq 0$.
The form of ODEs such as Eq. (2) has solved by methods of analytical (exact) and/or numerical (approximation) [2, 7, 11, 12, 14, 16, 18, 20, 29] used the numerical methods of Euler’s [2, 7, 12, 29], or its development, and/or RungeKutta fourthorder method (RK4) [2, 7, 11, 12, 16, 30]. Analytic solution based on the Eq. (2) [2628] for the seawater cooling process [24, 25], if:
#i) with $n=1$, has known as linear differential equations of order1, as shown in the Eq. (3), i.e.
$\frac{d}{d x} y+P(x) y=Q(x)$ (3)
#ii) with $Q(x)=0$, hence $\frac{d}{d x} y+P(x) y$ the order of homogeneous linear differential equations of order1 has used a general solution as shown in the Eq. (4) [2628], i.e.
$y=K \cdot e^{\int P(x) d x}$ (4)
#iii) with $Q(x) \neq 0$, hence $\frac{d}{d x} y+P(x) y=Q(x)$ is known as a nonhomogeneous differential equation with a general solution has shown in Eq. (5) [2628], i.e.
$y=K \cdot e^{\int P(x) d x}+e^{\int P(x) d x} \cdot \int e^{\int P(x) d x} \cdot Q(x)$$\cdot d x$ (5)
where, $K$ as an integration constant according to boundary conditions.
The numerical solution was been done through the simplest linear differential equations solving, i.e. solving with Euler's numerical based on the description of the Taylor series [2, 7, 11, 12, 14, 16, 18, 20, 2629]. The explanation of the Taylor series of a function in mathematics is a sum of terms in an infinite condition that is expressed in terms of the derivatives of a function at a single point. The general form of Taylor series as shown in Eq. (6).
$\begin{aligned} f\left(x_{i+1}\right)=f\left(x_{i}\right) &+f^{\prime}\left(x_{i}\right) \frac{\Delta x}{1 !}+f^{\prime \prime}\left(x_{i}\right) \frac{\Delta x^{2}}{2 !} \\ &+f^{\prime \prime \prime}\left(x_{i}\right) \frac{\Delta x^{3}}{3 !}+\cdots \\ &++f^{n}\left(x_{i}\right) \frac{\Delta x^{(n)}}{n !}+R_{n} \end{aligned}$ (6)
Based on Eq. (6), it can be explained that if a function $f(x)$ is known at point $x_{i}$ and all derivatives of $f$ with respect to $x$ are known at that point, then the Taylor series can be expressed the value of $f$ at point of $x_{i+1}$ which is located at a distance $\Delta x$ from point of $x_{i}$. Another form of the Taylor series is Eq. (7).
$f(x+\Delta x)=f(x)+\frac{d}{d x} y \cdot \Delta x+\frac{d^{2}}{d x^{2}} y \cdot \Delta x^{2}+\cdots$ (7)
where: $\frac{d}{d x} y=\frac{f(x+\Delta x)+f(x)}{\Delta x}$.
Furthermore, Eq. (7) can be expressed by Eq. (8).
$f(x+\Delta x)=\sum_{n=0}^{\infty} \frac{y^{(n)}\left(x_{0}\right)}{n !} \Delta x^{n}$ (8)
Or written as Eq. (9).
$f(x)=\sum_{n=0}^{\infty} \frac{y^{(n)}\left(x_{0}\right)}{n !}\left(xx_{0}\right)^{n}$ (9)
The results of ordinary differential equations models for observing the phenomena of temperature changes on a single rectangular platefin [20] refers to the book by Bird et al. [21]. The results are about the changes of temperature value on the copper bar and the fluid [20]. A shape of a single rectangular platefin with a thick much smaller than length (B << L) [21, 20] as shown in Figure 1.
Figure 1. A shape of a single rectangular platefin with a thick much smaller than length (B << L)
Based on Figure 1 can be explained, that the mathematical equation on the copper bar has shown in Eq. (10) [20].
$T_{c u}=\frac{\cosh \left[\sqrt{\left.\frac{h_{c u} \quad\cdot L^{2}}{k_{c u} \quad\cdot B} \cdot\left(1\frac{Z}{L}\right)\right]}\right.}{\cosh \sqrt{\frac{h_{c u} \quad\cdot L^{2}}{k_{c u} \quad\cdot B}}} \cdot\left(T_{w}T_{f}\right)+T_{f}$ (10)
The temperature value on the fluid has influenced by the time has shaped by the exponential curve with an initial temperature value of 20 degrees that based on the Eq. (10) as shown in Eq. (11) [20].
$\frac{d}{d t} T_{f}=\frac{h_{c u} \cdot A_{c u}}{m_{f} \cdot C_{p_{f}}} \cdot\left(T_{c u}T_{f}\right)+\frac{U_{w} \cdot A_{w}}{m_{f} \cdot C_{p_{f}}}$ $\cdot\left(T_{f}T_{w}\right)$ (11)
Based on several descriptions about the current status that have been published, then the determined three research objectives. The first, to get the constants of energy balance through completion of mathematical equations based on ODEs. The second, to made the three solving equations for simulation in the form of analytic, numerical of Euler’s, and RK4 solution. The third, to do the simulation processes using the spreadsheet application and result in analysis after simulation has been done.
2.1 Materials of the simulation
The equation of energy balance in the form of the single rectangular platefin (fin beams simple single rod) has done through the application of the calculation of cooling efficiency $[21,20]$. Solving with other numerical methods, i.e. the solution using the RK4 method which is a development of the basic Euler's method. The basic Euler's method and the Euler's method which has improved and modified, both of which are onestep methods, because they use initial information at a starting point, for example, $\left(t_{0}, y_{0}\right)$ to get estimates on $y\left(t_{1}\right)[2,7,11,12,16,29,30]$. Furthermore, the estimated value of $y\left(t_{1}\right)$ will use for estimation of $y\left(t_{2}\right)$. And so on, where at each stage the value of $t_{1}$ has advanced by $h_{1}$ and the value of $y$ will change by $\partial y_{i}[2,7,11,12,16,29,30]$.
Use mathematic, the process can be formulated by Eq. (12) [29, 30].
$y_{i+1}=y_{i}+e s t .\left[\partial y_{i}\right]$ (12)
Based on Eq. (12), the value of est. $\left[\partial y_{i}\right]$ has determined. The basic Euler’s method for determining value of est. $\left[\partial y_{i}\right]$ as shown in the Eq. (13) [2, 7, 11, 12, 16, 29], i.e.
$\partial y_{i}=\int_{t_{1}}^{t_{i+1}} y_{i}^{\prime} d t \approx y_{i}^{\prime}\left(t_{i+1}t_{i}\right)$ (13)
Eq. (12) has substituted into Eq. (13), so Eq. (14) has obtained [2, 7, 11, 12, 16, 29].
$y_{i+1}=y_{i}+\left(t_{i+1}t_{i}\right) y_{i}^{\prime}=y_{i}+\Delta t \cdot y_{i}^{\prime}$ (14)
The basic Euler's method has improved and modified through the application for estimation of $\partial y_{i}$ uses a value of $y_{i}{ }^{\prime}$ is smaller than true value and vice versa, if the value of $\partial y_{i}$ is estimated with a value of $y_{i+1}^{\prime}$, then the value of $\partial y_{i}$ has estimated will be greater. Euler's method improved through $y_{i}{ }^{\prime}$ and $y_{i+1}{ }^{\prime}$ for the determination the value of $\partial y_{i}$ as shown in Eq. (15) $[2,7,11,12,16,29]$, i.e.
$\partial y_{i}=\frac{h}{2}\left(y_{i}^{\prime}+y_{i+1}^{\prime}\right)$ (15)
Eq. (15) described as shown in Eq. (16).
$\partial y_{i}=\frac{h}{2}\left[f\left(t_{i}, y_{i}\right)+f\left(t_{i+1}, y_{i+1}\right)\right]$ (16)
Based on Eq. (16), the value of $\partial y_{i}$ cannot calculated, because the value of $y_{i+1}$ unknown. The revised and modified application of the basic Euler method requires the value of $y_{i+1}$ to calculated first using the basic Euler’s method according to Eq. (14), so that Eq. (17) obtained [2, 7, 11, 12, 16, 29].
$\begin{aligned} y_{i+1}=y_{i}+\frac{h}{2}\{& f\left(t_{\mathrm{i}}, y_{i}\right) \\ &\left.+f\left[t_{i+1}, y_{i+1}+h \cdot f\left(t_{\mathrm{i}}, y_{i}\right)\right]\right\} \end{aligned}$ (17)
Based on Eq. (17) shown, that the improved and modified Euler’s method is more accurate than the basic Euler’s method.
The RungeKutta method uses the average weight approach to determine the value of $\partial y_{i}$ with the aim that the acquisition of $\partial y_{i}$ value with an accuracy approach that can be achieved through the Taylor series [2, 7, 11, 12, 16, 30]. The value of $\partial y_{i}$ calculated by Eq. (18).
$\partial y_{i}=\left[a_{1}\left(k_{1}\right)_{i}+a_{2}\left(k_{2}\right)_{i}+\cdots+a_{r}\left(k_{r}\right)_{i}\right] \cdot h$ (18)
where:
$a_{1}, a_{2}, \cdots, a_{r}$ are constants
$k$ is a function of $f(t, y)$ given in the form as Eq. (19).
$\left(k_{1}\right)_{i}=f\left(t_{i}, y_{i}\right)\left(k_{2}\right)_{i}=f\left(t_{i}+p_{1} h, y_{i}+\right.$$\left.q_{1,1} k_{1} h\right)$
$\vdots$
$\left(k_{r}\right)_{i}=f\left(\begin{array}{c}t_{i}+p_{r1} h, y_{i}+q_{r1,1} k_{1} h+q_{r1,2} k_{2} h \\ +\cdots+q_{r1, r12} k_{r1} h\end{array}\right)$ (19)
The $r$ denotes the level of the RungeKutta method. For $r=1$, the $\mathrm{RK}$ method called the RK1 method. The RK2 method obtained from $r=2$, as well as the RK4 method obtained from $r=4$.
Using the RK2 method, it is necessary to determine the values $a_{1}, a_{2}, p_{1}$, and $q_{1,1}$ according to Eqns. (18) and (19), then the value of $y_{i+1}$ given by Eq. (20) $[2,7,11,12,16,29]$.
$y_{i+1}=y_{i}+a_{1} k_{1} h+a_{2} k_{2} h$ (20)
Eq. (20) described as shown in Eq. (21).
$\begin{aligned} y_{i+1}=y_{i}+a_{1} f &\left(t_{i}, y_{i}\right) h \\ &+a_{2} f\left(t_{i}+p_{i} h, y_{i}+q_{11} k_{1} h\right) h \end{aligned}$ (21)
The third term on the right side of Eq. (21), which is $f\left(t_{i}+p_{i} h, y_{i}+q_{11} k_{1} h\right)$ can expressed in the Taylor series so that it obtained like Eq. (22) [2, 7, 11, 12, 16, 30].
$\begin{aligned} f\left(t_{i}+p_{i} h, y_{i}+q_{11} k_{1} h\right) & \\ &=f\left(t_{i}, y_{i}\right)+p_{i} h \frac{\partial}{\partial t} f \\ &+q_{1,1} k_{1} h \frac{\partial}{\partial t} f+0\left(h^{2}\right) \end{aligned}$ (22)
Eq. (22) substituted into Eq. (21) with stages, i.e.
$\begin{aligned} y_{i+1}=y_{i}+a_{1} f(&\left.t_{i}, y_{i}\right) \\ &+a_{2}\left[f\left(t_{i}, y_{i}\right)+p_{i} h \frac{\partial}{\partial t} f\right.\\ &\left.+q_{1,1} k_{1} h \frac{\partial}{\partial t} f+0\left(h^{2}\right)\right] h \end{aligned}$
$\begin{aligned} y_{i+1}=y_{i}+a_{1} f(&\left.t_{i}, y_{i}\right) h+a_{2} f\left(t_{i}, y_{i}\right) h \\ &+a_{2}\left(p_{i} h \frac{\partial}{\partial t} f\right) h \\ &+a_{2}\left(q_{1,1} k_{1} h \frac{\partial}{\partial t} f\right) h \\ &+a_{2}\left(0\left(h^{2}\right)\right) h \end{aligned}$
Eq. (23) is obtained.
Referring to the equation $\frac{\partial}{\partial t} f=\frac{\partial}{\partial t} f(t, y)$ and $\frac{\partial}{\partial y} f=$ $\frac{\partial}{\partial y} f(t, y)$ compared to Eq. (23), the value of $y_{i+1}$ from the Taylor series method obtained as Eq. (24) $[2,7,11,12,16,30]$.
$\begin{aligned} y_{i+1}=y_{i}+\left[a_{1} f\right.&\left.\left(t_{i}, y_{i}\right)+a_{2} f\left(t_{i}, y_{i}\right)\right] h \\ &+\left[a_{2} p_{1} \frac{\partial}{\partial t} f\right.\\ &\left.+a_{2} q_{1,1} f\left(t_{i}, y_{i}\right) \frac{\partial}{\partial t} f\right] h^{2} \\ &+a_{2} 0\left(h^{3}\right) \end{aligned}$ (23)
$y_{i+1}=y_{i}+\left.h \frac{\partial}{\partial t} y\right_{i}+\left.\frac{1}{2} h^{2} h \frac{\partial^{2}}{\partial t^{2}} y\right_{i}+0\left(h^{3}\right)$ (24)
where:
$\left.\frac{\partial}{\partial t} y\right_{i}=f\left(t_{i}, y_{i}\right)$ (25)
$\left.\frac{\partial^{2}}{\partial t^{2}} y\right_{i}=f^{2}\left(t_{i}, y_{i}\right)=\frac{\partial}{\partial t} f+\frac{\partial}{\partial t} f+\frac{\partial}{\partial t} y$ (26)
The substitution of Eqns. (25) and (26) to Eq. (24) obtained as Eq. (27).
$y_{i+1}=y_{i}+f\left(t_{i}, y_{i}\right) h+\left(\frac{\partial}{\partial t} f+\frac{\partial}{\partial t} f+\frac{\partial}{\partial t} y\right) \frac{h^{2}}{2}$ $+0\left(h^{3}\right)$ (27)
Comparison between Eq. (23) and Eq. (27), the values are obtained as Eq. (28).
$a_{1}+a_{2}=1 ; a_{2} \cdot p_{1}=0.5 ;$ and $a_{2} \cdot q_{1,1}=0.5$ (28)
For four unknown variables from the three equations, one variable must be assumed $[2,7,11,12,16,30] .$ For $a_{2}$ it assumed to be $0.5$, then other constants obtained $a_{1}=0.5$ and $p_{1}=q_{1,1}=1$. This method known as the Heun's method (RK2). The constants that have obtained, substituted for Eq. (27), then the equation for the solution using the RK4 method obtained as Eq. (29) $[2,7,11,12,16,30]$.
$y_{i+1}=y_{i}+\frac{h}{2} f\left(t_{i}, y_{i}\right)+\frac{h}{2} f\left[t_{i}+h, y_{i}+h f\left(t_{i}, y_{i}\right)\right]$ (29)
2.2 Methods of the simulation
In this simulation, we assume a cooling chamber and a rectangular finshaped copper rod connected to the copper plate covering the cooling chamber. Cross section of the cooling container for simulation as shown in Figure 2.
Figure 2. Cross section of the cooling container for simulation
The assumption values have explained in the following descriptions.
#a) A rectangular box of seawater with the dimensions:
#a1) length, $l_{\text {space }}$ = 15 cm = 0.15 m;
#a2) width, $w_{\text {Space }}$ = 15 cm = 0.15 m;
#a3) height, $h_{\text {Space }}$ = 20 cm = 0.2 m; and
#a4) cooling container volume, vol $_{\text {space }}=(0.15 \cdot 0.15 \cdot 0.2)$ m^{3} = 0.0045 m^{3}.
#b) The cover of the box is maintained at a temperature of 0 ^{0}C and in the center of the cover is given a rectangular copper rod, $5 \mathrm{~cm} \cdot 5 \mathrm{~cm}$ and a length of 10 cm.
#c) The initial temperature of the seawater is at 25℃.
#d) A heat insulator is made of cork with a thickness of 5 cm.
#e) Additional information:
#e1) Density of copper, $\rho_{c u}$ = 8.9 x 10^{3} [kg/m^{3}];
#e2) Density of cork (for wall), $\rho_{w}$ = 160 [kg/m^{3}];
#e3) Density of seawater, $\rho_{f}$ = 1.00 x 10^{3} [kg/m^{3}];
#e4) Heat capacity of copper, $C_{\rho, c u}$ = 0.3831 [kJ/kg.K];
#e5) Heat capacity of cork (for walls), $C_{\rho, w}$ = 0.35 [kJ/kg.K];
#e6) Heat capacity of water (fluid), $C_{\rho, f}$ = 4.2 [kJ/kg.K];
#e7) Heat loss coefficient (assumed), $U_{w}$ = 4.8 [watts/m^{2}. ℃] = 0.01758 [watts/m^{2}K];
#e8) Heat coefficient of water (fluid), $h_{f}$ = 1.0 [watts/m^{2}. ℃];
#e9) Initial Condition (IC): $t=0 ; T_{c u}=0^{0} \mathrm{C}$; and (fluid temperature) $T_{f}$ = 25℃;
#e10) Boundary Condition (BC): $t \geq 0$; (air temperature), $T_{a}$ = constant = 28℃; and (wall temperature), $T_{w}$ = constant = 25℃.
An assumption is also complemented, that no additional energy loss in the study.
The steps in the simulation are: (a) making assumptions about the length, width, and height of the cooling room made of cork; (b) provision of preliminary conditions and limits; (c) implementing the stages of solving energy balance and obtaining mathematical equations for simulation (using analytical and numerical methods of Euler’s and RK4); and (d) implementation of simulations assisted by a spreadsheet application of the equations obtained, both using the analytical and numerical methods of Euler’s and RK4.
This chapter explains the objectives of the simulation that includes the completion of mathematical equations to get the constants for energy balance, solving equations for simulation, and simulation processes assisted by a spreadsheet application and result in analysis.
3.1 The constants for energy balance
The solution to the energy balance based on Eq. (11) obtained Eq. (30).
$\frac{d}{d t} T_{f}=b 1 \cdot\left(T_{c u}T_{f}\right)+b 2 \cdot\left(T_{f}T_{w}\right)$ (30)
where: $b 1=\frac{h_{c u} \cdot A_{c u}}{m_{f} \cdot C_{p_{f}}}$ and $b 2=\frac{U_{w} \cdot A_{w}}{m_{f} \cdot C_{p_{f}}}$.
Further solution to Eq. (30) obtained Eq. (31).
$\frac{d}{d t} T_{f}=b 1 \cdot T_{c u}+b 1 \cdot T_{f}+b 2 \cdot T_{f}b 2 \cdot T_{w}$ (31)
Eq. (31) simplified to Eq. (32).
$\frac{d}{d t} T_{f}=(b 1+b 2) \cdot T_{f}b 1 \cdot T_{c u}b 2 \cdot T_{w}$ (32)
Another form of Eq. (32) is Eq. (33).
$\frac{d}{d t} T_{f}(b 1+b 2) \cdot T_{f}=b 1 \cdot T_{c u}b 2 \cdot T_{w}$ (33)
The form of Eq. (33) identical to the form of Eq. (3), so that:
$\# \frac{d}{d t} T_{f}=\frac{d}{d x} y$,
$\#(b 1+b 2) \cdot T_{f}=P(x) y$, and
$\#b 1 \cdot T_{c u}b 2 \cdot T_{w}=Q(x)$.
Based on the value of $Q(x) \neq 0$, the analytical solution to Eq. (32) must be in the form of Eq. (5), so that a general solution obtained such as Eq. (34).
$y=K \cdot e^{\int(b 1+b 2) d t}+e^{\int(b 1+b 2) d t} \cdot \int e^{\int(b 1+b 2) d t}$
$\cdot\left(b 1 \cdot T_{c u}b 2 \cdot T_{w}\right) \cdot d t$ (34)
Further solution to Eq. (34) is obtained Eq. (35).
$\begin{aligned} T_{f}(t)=K & \cdot e^{(b 1+b 2) \cdot t}+e^{(b 1+b 2) \cdot t} \\ & \cdot \frac{\left(b 1 \cdot T_{c u}b 2 \cdot T_{w}\right)}{(b 1b 2)} \\ & \cdot e^{(b 1+b 2) \cdot t} \end{aligned}$ (35)
Simplifying Eq. (35) obtained Eq. (36).
$T_{f}(t)=K \cdot e^{(b 1+b 2) \cdot t}\frac{\left(b 1 \cdot T_{c u}b 2 \cdot T_{w}\right)}{(b 1+b 2)}$ (36)
Using Eq. (32) and $\frac{d}{d t} T_{f}=\frac{T_{f(i+1)}T_{f}}{\Delta t}$, then it obtained Eq. (37).
$\frac{T_{f(i+1)}T_{f}}{\Delta t}=(b 1+b 2) \cdot T_{f}b 1 \cdot T_{c u}b 2 \cdot T_{w}$ (37)
Further simplification of Eq. (37) obtained Eq. (38).
$T_{f(i+1)}=T_{f}+\left[(b 1+b 2) \cdot T_{f}b 1 \cdot T_{c u}b 2 \cdot T_{w}\right]$$\cdot \Delta t$ (38)
The solution using the numerical method of RK4 based on Euler's numerical method, so that Eqns. (38) and (29) used to obtain solving equation, such as Eq. (39) through some stages and detailed derivation process in the appendix.
$T_{f(i+1)}=T_{f}+\frac{4}{6}\left\{\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]^{2}\right.$ $\cdot \Delta t\}$ (39)
In achieving many stages, three basic forms of equations for the solution have generated, i.e. #i) the analytical method uses Eq. (36), #ii) the Euler's numerical uses Eq. (38), and #iii) the RK4 numerical uses Eq. (39). Based on the three equations, i.e. Eqns. (36), (38), and (39), it is necessary to determine the values of values $b_{1}$ and $b_{2}$. The calculation to get the values of $b_{1}$ and $b_{2}$ have done through:
(*) the fluid used is seawater, then:
$m_{f}=\rho_{f} \cdot v o l_{f}=\rho_{f}\left(v o l_{s p a c e}v o l_{c u}\right) \quad\left[\left(\mathrm{kg} / \quad \cdot \quad \mathrm{m}_{m}^{3}\right]=\right.$ $4.35625 \mathrm{~kg} ;$
$C_{p, f}=1[\mathrm{~kJ} / \mathrm{kg} . \mathrm{K}] ;$ and
$m_{f} \cdot C_{p, f}=4,35625[\mathrm{~kJ} / \mathrm{K}]$;
(**) fin simple block form is made of copper, then the value according to Adams and Rogers in 1973 [18] as like Eq. (40).
$h_{c u}=C_{p, c u} \cdot k_{c u} \cdot \theta^{\prime}(0) \cdot \frac{4}{3} \cdot L^{3} / 4 / L$ (40)
with:
value of $\theta^{\prime}(0)$ 0.5046;
value of $C_{p, c u}=$ 0.3831= 0.3831 [kJ/kg.K];
value of $k_{c u}=$ 386 [watts/m.K];
length of copper rod, L = 0.1 m;
so that the value of $h_{c u}$ contained in Eq. (39) is:
$h_{c u}=(0.3831) \cdot 386 \cdot(0.5046) \cdot \frac{4}{3} \cdot 0.1^{3 / 4} / 0.1$
$h_{c u}=(0.3831) \cdot 386 \cdot(0.5046) \cdot 1.333 \cdot 1.1778$
$h_{c u}=176.851\left[\mathrm{watt} / \mathrm{m}_{\omega}^{2} \cdot \mathrm{K}\right] .$
Besides, it is also obtained:
$A_{c u}=(0.05 \cdot 0.05)+4(0.05 \cdot 0.1)=0.0225 \mathrm{~m}^{2}$
$b 1=\frac{h_{c u} \cdot A_{c u}}{m_{f} \cdot C_{p_{f}}}$ = ($176.851 \cdot 0.0225$)/4.35625 = 3.9792 [watts/K]/4.35625 [kJ/K] = 0.9134 [watts/kJ] = $0.9134 \cdot 10^{3}$sec^{1}.
$U_{w}=0.01758[$ watts $/ \mathrm{m} 2 \mathrm{~K}]$;
$A_{w}=(0.15 \cdot 0.15)+4(0.15 \cdot 0.2)+[(0.15 \cdot 0.15)$
$(0.05 \cdot 0.050)]=0.7625 \mathrm{~m}^{2}$;
$b 2=\frac{U_{w} \cdot A_{w}}{m_{f} \cdot C_{p_{f}}}=(0.017 \cdot 0.7625) / 4.35625=0.01365[\mathrm{watts} / \mathrm{K}] /$$4.35625[\mathrm{~kJ} / \mathrm{K}]=0.31 \cdot 10^{6} \mathrm{~s}^{1}$.
3.2 The solving equations for simulation
Based on the values of $b 1=0.9134 \cdot 10^{3} \mathrm{sec}^{1}$ and $b 2=$ $0.31 \cdot 10^{6} \mathrm{sec}^{1}$, a mathematical equation for the rate of change in fluid temperature with time obtained using the analytical and numerical methods of Euler's and RK4. The use of analytical methods for obtaining mathematical equations based on Eq. (36) and substituting the values of $b 1=$ $0.9134 \cdot 10^{3} \mathrm{sec}^{1}$ and $b 2=0.31 \cdot 10^{6} \mathrm{~s}^{1}$ with stages:
$T_{f}(t)$$=K \cdot e^{\left(0.9134 \cdot 10^{3}+0.31 \cdot 10^{6}\quad\right) \cdot t}$ $\frac{\left(0.9134 \cdot 10^{3} \cdot T_{c u}0.31 \cdot 10^{6}\quad \cdot \cdot T_{w}\right)}{\left(0.9134 \cdot 10^{3}+0.31 \cdot 10^{6}\quad\right)}$
$T_{f}(t)$ $=K \cdot e^{0.91309 \cdot 10^{3} \quad\cdot t}$ $\frac{\left(0.9134 \cdot 10^{3} \cdot T_{c u}0.31 \cdot 10^{6} \quad\cdot \cdot T_{w}\right)}{0.91309 \cdot 10^{3}}$
So that the Eq. (36) turns into Eq. (41).
$T_{f}(t)=K \cdot e^{0.91309 \cdot 10^{3} \quad\cdot t}+1.000339507 \cdot T_{c u}$ $0.000339507 \cdot T_{w}$ (41)
For the solution to Eq. (41) used IC and BC that have determined. Substituting IC and BC into Eq. (41) the value of K obtained is:
$25=K \cdot e^{0.91309 \cdot 10^{3} \cdot 0}+1.000339507 \cdot 0$
$0.000339507 \cdot 0$
$25=K \cdot e^{0}+00>>>25=K \cdot 1 \gg>>25=K$ $K=25$.
Substituting the value of K into Eq. (41), the equation for the cooling rate obtained using an analytical method such as Eq. (42).
$T_{f}(t)=25 \cdot e^{0.91309 \cdot 10^{3} \cdot t}+1.000339507 \cdot T_{c u}$ $0.000339507 \cdot T_{w}$ (42)
The final equation form for the analytical solution given by Eq. (42), that change in temperature value (in seawater fluid) as a function of the time influenced by the initial fluid temperature, the copper rod temperature, and the wall temperature.
The use of the numerical methods of Euler's for obtaining mathematical equations based on Eq. (38), substituting the values of $b 1=0.9134 \cdot 10^{3} \mathrm{sec}^{1}$ and $b 2=0.31 \cdot 10^{6} \mathrm{sec}^{1}$, and used IC and BC that have determined, so that the Eq. (38) turns into Eq. (43), i.e. from,
$T_{f(i+1)}=T_{f}+\left(0.91309 \cdot 10^{3} \cdot T_{f}+0.9134\right.$ $\left.\cdot 10^{3} \cdot T_{c u}0.31 \cdot 10^{6} \cdot T_{w}\right)$ $\cdot \Delta t$
become:
$T_{f(i+1)}$ $=T_{f}\left(\begin{array}{c}0.0009134 \cdot T_{f}+0.0009134 \cdot T_{c u} \\ 0.00000031 \cdot T_{w}\end{array}\right) \Delta t$ (43)
The final equation form for the numerical methods of Euler’s solution given by Eq. (43), that change in temperature value (in seawater fluid) as a function influenced by the previous temperature of fluid, minus the sum or difference in temperature of the fluid, copper rod, and wall, and multiplied by the change in time.
The use of the numerical methods of RK4's for obtaining mathematical equations based on Eq. (39), substituting the values of $b 1=0.9134 \cdot 10^{3} \mathrm{sec}^{1}$ and $b 2=0.31 \cdot 10^{6} \mathrm{sec}^{1}$, and used $\mathrm{IC}$ and $\mathrm{BC}$ that have determined, so that the Eq. (39) turns into Eq. (44), i.e. from,
$T_{f(i+1)}=T_{f}+\frac{4}{6}\left\{\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]^{2}\right.$$\cdot \Delta t\}$
become:
$T_{f(i+1)}$$=T_{f}+\frac{4}{6}$$\cdot\left(\begin{array}{c}0.9134 \cdot 10^{3} \cdot T_{c u}0.91309 \cdot 10^{3} \cdot T_{f} \\ 0.31 \cdot 10^{6} \cdot T_{w}\end{array}\right)^{2} \cdot \Delta t$ (44)
The final equation form for the numerical methods of RK4’s solution given by Eq. (44), that change in temperature value (in seawater fluid) as a function influenced by the previous temperature of fluid, plus foursixths of the square of the difference between the temperatures of the copper rod, fluid, and wall, and multiplied by the change in time.
3.3 The simulation processes assisted by spreadsheet application and result in analysis
Using Eqns. (42), (43), and (44) and assisted by a spreadsheet application, the curve of temperature change in the fluid obtained. Changes in temperature on the fluid as shown in Figure 3.
Figure 3. Changes in temperature on the fluid
Based on Figure 3 can explained, that the results of the simulation processes based on a spreadsheet application have obtained the results, that the time span 1800 seconds (30 minutes) can cause changes in the temperature of the fluid from 25℃ to (i) 4.823770℃ with a rate of change of 0.092%, if used of the analytical method; (ii) 4.819088℃ with a rate of change of 0.092%, if used the numerical method of Euler’s; and (iii) 5.600404℃ with a rate of change of 0.083%, if used the numerical method of RK4’s. All the curves of changes in the temperature of the fluid are the nonlinear curves. Even though the final value of temperature is the highest, but the RK4 more thoroughly.
Based on results and discussions can concluded according to the research objectives. The completion of mathematical equations to get the constants for energy balance, i.e. the values of constant $b 1$ is $0.9134 \cdot 10^{3} \mathrm{sec}^{1}$ and constant $b 2$ is $0.31 \cdot 10^{6} \mathrm{sec}^{1}$. Constants b1 and b2 are important factors for making the solving equations for simulation related to the analytical method, the numerical method of Euler’s, and the numerical method of RK4’s. All nonlinear curves that obtained based on an exponential mathematical equation. Predicting the final temperature value for a time span of 1800 seconds, using the RK4 numerical method is the best because the error value is the smallest, even though the temperature value is the highest. Suggestions for future work, for all the simulation of its phenomena and/or the others that use the ODE’s begin with making the mathematical models and implemented them into the computer application with the analytical method or numerical methods of Euler’s and/or RK4 has chosen.
Submitting of thanks to Prof. Dr. Kamaruddin Abdullah who had given the science and knowledge of the transport phenomena subjects in the lecturing on several years ago.
ODEs 
ordinary differential equations 
RK 
RungeKutta 
$x$ 
independent variable; $x$axis 
$y$ 
dependent variable; $y$ axis 
$F, P, Q, f$ 
functions 
K 
an integration constant according to boundary conditions 
d 
differentiating with one variable 
$\partial$ 
differentiating with two or more variables 
$e$ 
exponential number = 2.71828… 
∫ 
integrating 
$f\left(x_{i}\right)$ 
function at point of $x_{i}$ 
$f\left(x_{i+1}\right)$ 
function at point of $x_{i+1}$ 
$f^{\prime}, f^{\prime \prime}, \cdots f^{(n)}$ 
the first, second, ..., n^{th} derivative of the function 
$\Delta x$ 
distance between $x_{i}$ and $x_{i+1}$ 
$R_{n}$ 
truncation error 
$\sum$ 
summing 
$\Delta$ 
the segment of a 
n 
constant 
$!$ 
factorial operator 
B 
fin thickness 
L 
fin lengthiness 
W 
fin widths 
z 
z axis 
q 
heat flux at the surface 
h 
heat transfer coefficient 
T 
absolute temperature 
K 
thermal conductivity, W.m^{1}. K^{1} 
$C_{p}$ 
specific heat, J. kg^{1}. K^{1} 
$a_{1}, a_{2}, \cdots, a_{r}$ 
contants 
$k_{1}, k_{2}, \cdots, k_{r}$ 
contants 
$t_{0}$ 
initial time for information at a starting point 
$y_{0}$ 
initial dependent variable for information at a starting point 
$t_{1}$ 
estimate time for information 
$y_{1}$ 
estimate dependent variable for information 
$y_{t+1}$ 
nextorder in time for estimate dependent variable 
$t$ 
time 
U 
overall heat transfer coefficient 
A 
area, m^{2} 
IC 
initial conditions 
BC 
boundary conditions 
$v o l_{f}$ 
volume of fluid 
vol $_{\text {space }}$ 
volume of cooling space 
$v o l_{c u}$ 
volume of copper bar 
Greek symbols 

$\Theta$ 
dimensionless temperature 
$\zeta$ 
dimensionless distance 
N 
dimensionless heat transfer coefficient 
Superscripts 

$1,2, n, n1$ 
sequence to 
$(n)$ 
nth derivative 
Subscripts 

$1,2, n$ 
sequence to 
met. 
etal 
cu 
copper 
$a$ 
air 
$f$ 
fluid 
w 
wall 
z 
in zaxis 
The detailed of derivation process to get Eq. (39).
$T_{f(i+1)}=T_{f}+$$\frac{1}{6}\left(\begin{array}{c}{\left\{\begin{array}{c}{\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]} \\ +4\left[\frac{b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}}{2}\right] \\ +\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]\end{array}\right\}} \\ {\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right] \Delta t}\end{array}\right)$
$T_{f(i+1)}=T_{f}+$$\frac{1}{6}\left(\begin{array}{c}\left\{\begin{array}{c}2\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]+ \\ 4\left[\frac{b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}}{2}\right]\end{array}\right\} \\ {\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right] \Delta t}\end{array}\right)$
$T_{f(i+1)}=T_{f}+$$\frac{1}{6}\left(\begin{array}{c}\left\{\begin{array}{c}2\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]+ \\ 2\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]\end{array}\right\} \\ {\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right] \Delta t}\end{array}\right)$
$\begin{aligned} T_{f(i+1)}=T_{f}+\frac{1}{6} &\left\{4\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]\right.\\ &\left.\cdot\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right] \Delta t\right\} \end{aligned}$
$\begin{aligned} T_{f(i+1)}=T_{f}+\frac{4}{6} &\left\{\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right]\right.\\ &\left.\cdot\left[b_{1} T_{c u}\left(b_{1}b_{2}\right) T_{f}b_{2} T_{w}\right] \Delta t\right\} \end{aligned}$
[1] Borwein, J.M., Crandall, R.E. (2013). Closed forms: What they are and why we care. Notices of the American Mathematical Society, 60(1): 5065. https://doi.org/10.1090/noti936
[2] Süli, E. (2010). Numerical solution of ordinary differential equations. Mathematical Institute, University of Oxford. https://people.maths.ox.ac.uk/suli/nsodes.pdf.
[3] Mex, L., CruzVillar, C.A., Peñuñuri, F. (2015). Closedform solutions to differential equations via differential evolution. Discrete Dynamics in Nature and Society, pp. 111. https://dx.doi.org/10.1155/2015/910316
[4] Meyer, J.F. (1982). Closedform solutions of performability. IEEE Transactions on Computers, 31(7): 648657. https://dx.doi.org/10.1109/TC.1982.1676062
[5] Co, T.B. (2013). Analytical solutions of ordinary differential equations. Methods of Applied Mathematics for Engineers and Scientists, Cambridge: Cambridge University Press, pp. 235272. https://doi.org/10.1017/CBO9781139021821.010
[6] McNamee, J.M., Pan, V.Y. (2013). Lowdegree polynomials. In Studies in Computational Mathematics, 16: 527556. https://doi.org/10.1016/b9780444527301.000060
[7] Atkinson, K.E., Han, W., Stewart, D.E. (2009). Numerical Solution of Ordinary Differential Equations. John Wiley & Sons.. http://dx.doi.org/10.1002/9781118164495
[8] Christodoulou, D.M., Katatbeh, Q.D. (2017). A powerful diagnostic tool of analytic solutions of ordinary secondorder linear homogeneous differential equations. Advances in Difference Equations, 2017(1): 19. https://doi.org/10.1186/s1366201712725
[9] Ardourel, V., Jebeile, J. (2017). On the presumed superiority of analytical solutions over numerical methods. European Journal for Philosophy of Science, 7(2): 201220. https://doi.org/10.1007/s1319401601522
[10] Brownlee, J. (2018). Analytical vs numerical solutions in machine learning. Machine Learning Mastery Website. https://machinelearningmastery.com/analyticalvsnumericalsolutionsinmachinelearning.
[11] Dios, A.Q., Encinas, A.H., Vaquero, J.M., del Rey, A.M., Pérez, J.J.B., Sánchez, G.R. (2015). How engineers deal with mathematics solving differential equation. Procedia Computer Science, 15: 19771982. https://doi.org/10.1016/j.procs.2015.05.462
[12] Ming, C.Y. (2017). Solution of differential equations with applications to engineering problems. Dynamical SystemsAnalytical and Computational Techniques, pp. 233264. https://dx.doi.org/10.5772/67539
[13] Liu, M., Xue, W., Makarov, S.B., Qi, J., Li, B. (2019). Comparative study of analytic solution and numerical solution of baseband symbol signal based on optimal generic function. Mathematical Problems in Engineering, 2019: 8045217. https://doi.org/10.1155/2019/8045217
[14] Shawagfeh, N., Kaya, D. (2004). Comparing numerical methods for the solutions of systems of ordinary differential equations. Applied Mathematics Letters. 17(3): 323328. https://doi.org/10.1016/S08939659(04)900705
[15] Devries, P.L., Hasbun, J.E. (2011). A First Course in Computational Physics. Second edition. Jones and Bartlett Learning, Burlington, Massachusetts.
[16] Hsu, T.R. (2018). Numerical solution methods for engineering analysis. Applied Engineering Analysis, John Wiley & Sons, Hoboken, NJ, pp. 339380.
[17] Anandrao, P., Koushal, R. (2018). An overview of numerical and analytical methods for solving ordinary differential equations. International Journal of Advance Research and Innovative Ideas in Education (IJARIIE). 4(2): 50905095. https://doi.org/10.48550/arXiv.2012.07558
[18] Adams, J.A., Rogers, D.F. (1973). Computer Aided Heat Transfer Analysis. McGrawHill, New York.
[19] Joel, A.S., ElNafaty, U.A., Makarfi, Y.I., Mohammed, J., Musa, N.M. (2017). Study of effect of fin geometry on cooling process of computer microchips through modelling and simulation. International Journal of Industrial and Manufacturing Systems Engineering, 2(5): 4856. https://doi.org/10.11648/j.ijimse.20170205.11
[20] Goeritno, A. (2021). Ordinary differential equations models for observing the phenomena of temperature changes on a single rectangular plate fin. Mathematical Modelling of Engineering Problems, 8(1): 8994. https://doi.org/10.18280/mmep.080111290
[21] Bird, R.B., Stewart, W.E, Lightfoot, E.N. (1960). Transport Phenomena. John Wiley & Sons, New York.
[22] Glasgow, L.A. (2010). Transport Phenomena: An Introduction to Advanced Topics. John Wiley & Sons, Inc., Hoboken, New Jersey, 83100.
[23] Lienhard, I.V., John, H. (2005). A Heat Transfer Textbook. Phlogiston Press. Cambridge, Massachusetts, USA.
[24] Besson, U. (2012). The history of the cooling law: When the search for simplicity can be an obstacle. Science & Education, 21: 10851110. https://doi.org/10.1007/s1119101093241
[25] Zheng, L., Zhang, X. (2017). Modeling and analysis of modern fluid problems. Academic Press. pp. 361455. https://doi.org/10.1016/B9780128117538.000086
[26] Coddington, E.A., Levinson, N. (1987). Theory of Ordinary Differential Equations (9th Reprinted). Tata McGrawHill Publishing Co. Ltd., New Delhi.
[27] Archibald, T., Fraser, C., GrattanGuinness, I. (2005). The history of differential equations, 16701950. Oberwolfach Reports, 1(4): 27292794. https://doi.org/10.4171/OWR/2004/51
[28] Logan, J.D. (2011). A First Course in Differential Equations, second edition. Springer Science+Business Media, New York, NY.
[29] Barker, C.A. (2017). Numerical methods for solving differential equations: Euler’s method. Computer Laboratory, Mathematics and Science Learning Center, San Joaquin Delta College, Stockton, CA.
[30] Islam, M.A. (2015). Accurate solutions of initial value problems for ordinary differential equations with fourth order RungeKutta method. Journal of Mathematics Research, 7(3): 4145. https://doi.org/10.5539/jmr.v7n3p41