Part A

We have a steady state non-dimensional heat equation given by \[\frac{d^2 \theta}{d\xi^2} - \beta \theta = -\beta \theta_f \label{eq:ss}\] Where \(\theta_f\) does not depend on \(\xi\). The exact solution is a function of the parameters \(\beta\) and \(\theta_f\). The solution of this equation takes the form: \[\theta(\xi) = c_1 e^{\sqrt{\beta} \xi} + c_1 e^{-\sqrt{\beta} \xi} - \frac{\theta_f}{\beta}\] The exact analytic solution can be found for a particular pair of boundary conditions \(\theta(\xi_1)\) and \(\theta(\xi_2)\) by solving the system given by: \[\begin{bmatrix} e^{\sqrt{\beta} \xi_1} & e^{-\sqrt{\beta} \xi_1} \\ e^{\sqrt{\beta} \xi_2} & e^{-\sqrt{\beta} \xi_2} \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} \theta(\xi_1) \\ \theta(\xi_2) \end{bmatrix}\] For the coefficients \(c_1\), \(c_2\).

We can obtain a numerical approximation of [eq:ss] at a grid point \(i\) by replacing the derivative term with a finite difference approximation over a grid with spacing \(h\): \[\frac{d^2 \theta_i}{d\xi^2} = \frac{ \theta_{i-1} - 2 \theta_i + \theta_{i+1} }{h^2} + \mathcal{O}(h^2) \label{eq:finitediff2}\] Substituting [eq:finitediff2] into [eq:ss] yields: \[\begin{aligned} \frac{\theta_{i-1}}{h^2} + \frac{-2 \theta_i}{h^2} + \frac{\theta_{i-1}}{h^2} + \mathcal{O}(h^2) - \beta{\theta_i} &= -\beta \theta_{i,f}\\ \left(\frac{-1}{h^2}\right) \theta_{i-1} + \left(\frac{2}{h^2} + \beta \right) \theta_{i} + \left(\frac{-1}{h^2}\right) \theta_{i+1} + \mathcal{O}(h^2) &= \beta \theta_{i,f} \label{eq:Amatrixeqn} \end{aligned}\] Using this equation, we can obtain a numerical approximation for each interior grid point \(i\). For a grid with step size \(h=\frac{1}{6}\), we have \(n=6\) grid points, \(n=4\) of which are internal grid points. We will use a zero-indexed notation for the grid points. In this case, \(\theta\) is known at the endpoints, \(i=0\) and \(i=5\). We can construct the matrix system \([A] \{\theta\} = \{b\}\) as follows. \[\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{h^2} & \frac{2}{h^2} + \beta & -\frac{1}{h^2} & 0 & 0 & 0 \\ 0 & -\frac{1}{h^2} & \frac{2}{h^2} + \beta & -\frac{1}{h^2} & 0 & 0 \\ 0 & 0 & -\frac{1}{h^2}& \frac{2}{h^2} + \beta & -\frac{1}{h^2} & 0 \\ 0 & 0 & 0 & -\frac{1}{h^2} & \frac{2}{h^2} + \beta & -\frac{1}{h^2} \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \theta_3 \\ \theta_4 \\ \theta_5 \end{bmatrix} = \begin{bmatrix} \theta(\xi_0) \\ \theta_f \\ \theta_f \\ \theta_f \\ \theta_f \\ \theta(\xi_5) \end{bmatrix} \label{eq:ss_Axb}\]

This tridiagonal system can be solved quickly using the Thomas algorithm (TDMA) to obtain the full vector of temperatures. The solution for this system for several parameters is shown in 1

We compare the finite difference approximation to the exact solution for various grid sizes in 2 and verify that the error is \(\mathcal{O}(h^2)\).

Part 5

Consider a control volume that encloses the fin. We have two areas through which heat can pass: the fin surface, exposed to the fluid, and the fin ends, which are attached to some solid substrate. The heat flux through the fin surface occurs via convection, and the heat flux through the fin ends occurs via conduction. If we assume that no energy is generated inside of the fin, the total heat flux is given by: \[Q = Q_{\text{cond}} + Q_{\text{conv}} = 0\] Then, we can say that \[Q_{\text{cond}} = -Q_{\text{conv}}\] That is, depending on the fluid temperature, heat travels from the substrate, through the fin, and into the fluid, consistent with a cooling application. To evaluate the transfer by convection, we have that: \[Q_{\text{conv}} = \int_{0}^{L} h (T-T_f) dx\] Where \(h\) is the coefficient of convection (not grid spacing!). We can non-dimensionalize to obtain the conduction term: \[\beta \int_{0}^{1}(\theta - \theta_f) d\xi \label{eq:convterm}\] The convection term is modeled by Newton’s law of cooling: \[k A (\left. \frac{dT}{dx} \right|_{0} - \left. \frac{dT}{dx} \right|_{L}) = 0\] which can be non-dimensionalized: \[\left. \frac{d\theta}{d\xi} \right|_{0} - \left. \frac{d\theta}{d\xi} \right|_{1} \label{eq:condterm}\] If the energy balance is to hold, we have then that convection out of the control volume ([eq:convterm]) equals conduction into the control volume ([eq:condterm]): \[0 = \left. \frac{d\theta}{d\xi} \right|_{0} - \left. \frac{d\theta}{d\xi} \right|_{1} + \beta \int_{0}^{1}(\theta - \theta_f) d\xi\] The convection term can be obtained via integration of the temperature profile. We can use the analytic solution to ([eq:ss]) to obtain the temperature profile. Then, we can numerically integrate.

Integration to solve convection was carried out using the trapezoid rule for \(h=\frac{1}{8}\), \(h=\frac{1}{16}\), and \(h=\frac{1}{32}\). Additionally, a \(\mathcal{O}(h^4)\) and \(\mathcal{O}(h^6)\) romberg method was used.

Meanwhile, differentiation to solve conduction was kept steady at a 2nd-order method for each boundary point.

Then, a series of problems with varying \(\beta\) values were solved. Because the convection term depends on \(\beta\), the accuracy of an integration methods will depend not only on \(h\), but also on \(\beta\). The results are shown in 3.

Part B

Now, we consider a problem where the free stream temperature is a function of its station along the fin: \(\theta_f = f(\xi)\). We consider a gaussian "plume" of high temperature in the free stream, given by the function: \[f(\xi) = A \exp{\left( \frac{-(\xi - \xi_0)^2}{2 \sigma^2} \right)}\] The plume has magnitude \(A\), is centered at \(\xi_0\), and has a width corresponding to \(\frac{1}{\sigma ^2}\).

Although finding the analytic solution of the \(\theta_f=f(\xi)\) problem may be difficult or impossible, the numerical solution is easy: the procedure is the same regardless of the fluid temperature profile. We need only substitute the \(b\) vector.

Part C

Now we change the conditions so that the fluid temperature is also a function of time. \[\theta_f(\xi, \tau) = A \sin{\omega \tau} \exp{\left(\frac{-(\xi-\xi_0)^2}{2 \sigma ^2}\right)} \label{eq:thetaf_dynamic}\] This represents a fluid with a plume centered at \(\xi_0\), the temperature of which oscillates between \(A\) and \(-A\) in a sinusoidal fashion with a period \(\omega\).

We can rewrite our governing equation. The \(\theta_{fluid \hspace{0.1cm} i,j}\) term, rather than being a steady-state function only of \(\xi\), is now a function both of time and space, governed by [eq:thetaf_dynamic].

Explicit Solution to Time-Dependent Fin Temperature

\[\frac{\partial \theta}{\partial \tau} = \frac{\partial^2 \theta}{\partial \xi^2} + \beta \theta_f - \beta \theta \label{eq:governing_time}\] We can give the RHS of this equation a similar treatment as in 1: \[\begin{aligned} \frac{\partial \theta}{\partial \tau} &= \frac{\theta_{i-1} -2 \theta_i + \theta_{i-1}}{h^2}+ \beta \left(\theta_{f,i} - \theta_i \right) + \mathcal{O}(h^2) \\ \frac{\partial \theta}{\partial \tau} &= \frac{\theta_{i-1}}{h^2} + \frac{-2 \theta_i}{h^2} + \frac{\theta_{i-1}}{h^2}+ \beta \left(\theta_{f,i} - \theta_i \right) + \mathcal{O}(h^2) \\ \frac{\partial \theta}{\partial \tau} &= \left(\frac{1}{h^2}\right) \theta_{i-1} + \left(\frac{-2}{h^2} - \beta \right) \theta_i + \left(\frac{1}{h^2}\right) \theta_{i+1} + \left(\beta\right) \theta_{f,i,j} + \mathcal{O}(h^2) \label{eq:explicit_time} \end{aligned}\] This produces a matrix equation of the form: \[\frac{\partial \theta}{\partial \tau} = -[A] \theta + b\] Upon expanding the terms of this equation, we can see that they are identical to those of 1. In fact, when \(\frac{\partial \theta}{\partial \tau} = 0\), this equation is identical to [eq:ss_Axb]. When there is no time derivative in the sytem, we are solving a steady-state equation.

With this formulation, we have the form \(\frac{\partial \theta}{\partial \tau} = f(\theta, \tau)\), which we can use any explicit method to solve. We have three different boundary conditions: in order to increment time, we need the first state to be a fully fixed (steady state) condition, which can be ascertained by solving [eq:ss_Axb] for \(\left. \theta_f \right|_{\tau=0}\). The second and third boundary conditions are the left and right fin ends: \(i=0\) and \(i=n-1\), corresponding to the fixed temperature of the sink at \(\xi=0\) and \(\xi=1\). The solution is physically accurate for the step size \(t=0.49 h^2\). We can see in 5 that the behavior of the fin follows our physical intuition: the higher frequency fluid temperature oscillation does not allow enough time for the fin to absorb the heat, and hence the fin temperature is steadier than the lower-frequency temperature oscillation.

The same explicit method cannot be used to accurately solve the problem for a step size of \(t>0.5 h^2\), because of stability issues in the explicit method. This can prove to be a significant obstacle in both runtime and memory, especially where the dimensionality of the problem increases, for example when a 2-D or 3-D temperature profile must be calculated.

To obtain the solution shown in 5, 166 steps were needed for the short period \(\omega=400\) case. However, 16415 timesteps were needed to compute the long-period \(\omega=4\) case, because the timestep size had to be kept small. The required small timestep makes the solution quite accurate, but the runtime can be a serious obstacle.

A runtime analysis of the explicit method can be found in 6. These severe stability downsides are addressed by the Crank-Nicolson method in the following section.

"Semi-Implicit" Crank-Nicolson Solution to Time Dependent Fin Temperature

Rather than solving the governing equation explicitly, we can use the implicit Crank-Nicolson method, which alleviates the stability problems we saw in 3.1. We begin with the matrix form of our governing equation: \[\frac{\partial \theta}{\partial \tau} = -[A] \theta + b\] We can obtain an approximation for the equation at the next timestep by formulating the quantity \(\theta_{j+1}\) with the mid-point rule, which is the combination of the forward and backward Euler methods yielding \(\mathcal{O}(h^2)\) error. Because the backward Euler method has an implicit dependence on the solution at \(j+1\), solving this equation means solving a linear system \(\theta_{j+1} = B \theta_{j} + c\) at each timestep. \[\begin{aligned} \frac{\theta_{j+1} - \theta_{j}}{t} &= \frac{1}{2} \left[ \left(\frac{\partial \theta}{\partial \tau}\right)_{j} + \left(\frac{\partial \theta}{\partial \tau}\right)_{j+1}\right] \\ \frac{\theta_{j+1} - \theta_{j}}{t} &= \frac{1}{2} \left[ \left(-A \theta + b\right)_{j} + \left(-A \theta + b\right)_{j+1} \right] \\ \theta_{j+1} &= \theta_{j} + \frac{t}{2} \left[ \left(-A \theta + b\right)_{j} + \left(-A \theta + b\right)_{j+1} \right] \\ \theta_{j+1} &= \theta_{j} - \frac{t}{2} A \left( \theta_{j+1} + \theta_{j} \right) + \frac{1}{2} b \left( \theta_{j+1} + \theta_{j} \right) \\ \left( I + \frac{t}{2} A \right) \theta_{j+1} &= \left( I - \frac{t}{2} A \right) \theta_{j} + \frac{t}{2} \left( b_{j+1} + b_{j} \right) \label{eq:cn} \end{aligned}\] Since \(A\) is tridiagonal, the matrix \(\left( I + \frac{t}{2} A \right)\) is also tridiagonal, and hence can be solved quickly. Although the method is iterative, because the value \(\theta_{j+1}\) implicitly considered, it does not suffer from the same stability issues as the explicit method. The downside is the computational cost at each timestep (inverting the matrix, rather than a simple matrix multiplication). However, because the timestep can be made much larger than explicit RK2 method, Crank-Nicolson can give us solutions to otherwise infeasible problems. The solutions in 7 took only 156 timesteps to compute, both for the \(\omega=4\) case and the \(\omega=400\) case – significantly better than the explicit case above.