Rico

5.2  2D PLPM Thermal Analysis of the Experiment

Introducing the Experiment

The experimental apparatus by P. A. Hilton (Hilton, 2006; Hilton, 2011) is designed to study radial heat transfer. In figure 1, a photo of the H112B unit (Hilton, 2011) is shown. In figure 2, a cross section of the apparatus is shown. In the center we see a heater that supplies heat to the surface of an insulated annular cylinder made of brass, which we will call the "disk." Along the outer edge of the disk, we see water cooling. Along a ray from the center, 6 thermocouples are embedded. Their radial positions from the center of the disk are as follows:

  • 7 mm
  • 10 mm
  • 20 mm
  • 30 mm
  • 40 mm
  • 50 mm
Figure 1: Photo of the experimental apparatus (Hilton, 2011)
Figure 2: Diagram of the experimental apparatus (Hilton, 2011)

The heater is an electrical resistance that is supplied with variable AC power from the service unit H112 (Hilton, 2006), controlled by a potentiometer (knob labeled "VOLTAGE CONTROLLER") on the front panel, as shown in figure 3. The readings of the current \(i_{\text{RMS}}\) and voltage \(v_{\text{RMS}}\) on the front panel are in RMS, so the average power dissipated in the resistive heater is \(P = i_{\text{RMS}} v_{\text{RMS}}\). All of the dissipated power (units W) is heat transferred through the rod; that is, \(\dot{q} = P\). The maximum voltage is 240 V, although the automatic power shutdown that prevents overheating will begin to interfere with the experiment at higher power levels.

Figure 3: Photo of the service unit (Hilton, 2006)

Although we won't use this feature, the H112B thermocouples can be connected to the sockets and corresponding temperature readings can be displayed on the front panel of the service unit.

Instead of capturing the thermocouple readings with the service unit, we will use two NI 9211 data loggers, shown in figure 4, connected via USB to a PC running LabVIEW.

Figure 4: Photo of the NI data loggers.

Analyzing the Experiment

We will use a 2D partial lumped-parameter model (PLPM) to analyze the experiment. 2D conduction can be modeled analytically using the heat equation, which is a partial differential equation (PDE) that describes the distribution of heat in a given region over time. However, the heat equation is difficult to solve analytically for complex geometries and boundary conditions, and PDE solutions are beyond the scope of most undergraduate courses. Using PLPM, we can derive a set of linear ordinary differential equations (ODEs) that describe the temperature distribution in the rod as a function of time. These can be solved analytically, but we will instead use an easy-to-use numerical ODE solver available in Python's SciPy library.

The Finite Elements

Due to the symmetry of the disk, we can expect identical temperature distributions along any ray from the center of the disk to its edge. Therefore, we can model the disk in a manner similar to that of the rod in the previous two chapters. We divide the disk into \(m\) even sectors, each with interior angle \(\theta = 2\pi/m\). We segment a sector of the disk into \(n\) finite elements, each with identical radial length \(r_{k+1} - r_k\) and height \(\ell\), as shown in figure 5, which shows a single sector. Note that the finite elements increase in volume as we move away from the center of the disk, which means that the thermal capacitance increases as well. The volume \(V_k\) of an element \(k\) is given by \[ V_k = \pi (r_k^2 - r_{k-1}^2) \ell / m. \] The volume of an element affects its thermal capacitance, as we will consider in a moment. Similarly, the contact surface area at radius \(r_k\) of adjacent elements in a sector is given by \[ A_k = r_k \theta \ell. \] This interface area affects the thermal resistance between adjacent elements in a sector, which we will also consider in a moment.

Figure 5: Representative finite elements of the disk in a single sector. The center of the disk corresponds to the left boundary, and the outer edge corresponds to the right boundary.

We will model each finite element as a lumped capacitance \(C_i\) connected by thermal resistors \(R_i\) and, in most cases, \(R_{i-1}\), to neighboring elements. We assume that all heat transfer occurs through the radial direction, which is to say that the temperature is assumed to be constant along cylindrical surfaces. As we have already noted, this is valid due to the symmetry of the problem. The leftmost element has the heat transfer rate source \(\dot{q}_S\) entering the system. The rightmost element has the temperature source \(T_S\) (sink, as it were), which will model for us the free-stream temperature of the cooling water. The rightmost conduction resistance \(R_n\) has only half the radial length of the other resistors. The resistance \(R_{n+1}\) models the convection between the rod and the cooling water.

The Linear Graph Model

We will use a linear graph modeling approach. Consider the linear graph shown in figure 6.

Figure 6: Linear graph model with highlighted normal tree.

Variable Definitions

We define variables as follows:

  • Primary Variables: \(T_{C_1}\), \(T_{C_2}\), ..., \(T_{C_n}\), \(T_{R_{n+1}}\), \(\dot{q}_{R_1}\), \(\dot{q}_{R_2}\), ..., \(\dot{q}_{R_n}\)
  • Secondary Variables: \(\dot{q}_{C_1}\), \(\dot{q}_{C_2}\), ..., \(\dot{q}_{C_n}\), \(\dot{q}_{R_{n+1}}\), \(T_{R_1}\), \(T_{R_2}\), ..., \(T_{R_n}\)
  • State Variables: \(T_{C_1}\), \(T_{C_2}\), ..., \(T_{C_n}\)
  • Vectors: \[\bm{x} = \begin{bmatrix} T_{C_1} \\ T_{C_2} \\ \vdots \\ T_{C_n} \end{bmatrix}, \qquad \bm{u} = \begin{bmatrix} \dot{q}_S \\ T_S \end{bmatrix}, \qquad \bm{y} = \bm{x}.\]

Elemental, Continuity, and Compatibility Equations

Now we are ready to write the equations for the system. The elemental equations are as follows:

\[\begin{align} C_1 &: \dot{T}_{C_1} = \frac{1}{C_1} \dot{q}_{C_1} \\ C_2 &: \dot{T}_{C_2} = \frac{1}{C_2} \dot{q}_{C_2} \\ \vdots &\phantom{:} \vdots \\ C_n &: \dot{T}_{C_n} = \frac{1}{C_n} \dot{q}_{C_n} \\ R_1 &: \dot{q}_{R_1} = \frac{1}{R_1} T_{R_1} \\ R_2 &: \dot{q}_{R_2} = \frac{1}{R_2} T_{R_2} \\ \vdots &\phantom{:} \vdots \\ R_n &: \dot{q}_{R_n} = \frac{1}{R_n} T_{R_n} \\ R_{n+1} &: T_{R_{n+1}} = \dot{q}_{R_{n+1}} R_{n+1} \end{align}\]

A continuity equation is required for each passive branch in the normal tree. They are as follows:

\[\begin{align} C_1 &: \dot{q}_{C_1} = \dot{q}_S - \dot{q}_{R_1} \\ C_2 &: \dot{q}_{C_2} = \dot{q}_{R_1} - \dot{q}_{R_2} \\ \vdots &\phantom{:} \vdots \\ C_n &: \dot{q}_{C_n} = \dot{q}_{R_{n-1}} - \dot{q}_{R_n} \\ R_{n+1} &: \dot{q}_{R_{n+1}} = \dot{q}_{R_n} \end{align}\]

A compatibility equation is required for each passive link (element not in the normal tree). They are as follows:

\[\begin{align} R_1 &: T_{R_1} = T_{C_1} - T_{C_2} \\ R_2 &: T_{R_2} = T_{C_2} - T_{C_3} \\ \vdots &\phantom{:} \vdots \\ R_n &: T_{R_n} = T_{C_n} - T_S - T_{R_{n+1}} \end{align}\]

This is a set of \(4 n + 2\) linear algebraic and differential equations and \(4 n + 2\) unknown variables.

Reducing to the State-Space Model

We can systematically reduce the system of equations to a set of \(n\) first-order differential equations in the state variables of \(\bm{x}\) by eliminating non-state and non-input variables, as follows.

We begin by substituting the continuity and compatibility equations into the elemental equations, eliminating a secondary variable with each substitution. This will eliminate half of the unknown variables and equations. The substitutions yield the following set of equations:

\[\begin{align} C_1 &: \dot{T}_{C_1} = \frac{1}{C_1} (\dot{q}_S - \dot{q}_{R_1}) \\ C_2 &: \dot{T}_{C_2} = \frac{1}{C_2} (\dot{q}_{R_1} - \dot{q}_{R_2}) \\ \vdots &\phantom{:} \vdots \\ C_n &: \dot{T}_{C_n} = \frac{1}{C_n} (\dot{q}_{R_{n-1}} - \dot{q}_{R_n}) \\ R_1 &: \dot{q}_{R_1} = \frac{1}{R_1} (T_{C_1} - T_{C_2}) \\ R_2 &: \dot{q}_{R_2} = \frac{1}{R_2} (T_{C_2} - T_{C_3}) \\ \vdots &\phantom{:} \vdots \\ R_n &: \dot{q}_{R_n} = \frac{1}{R_n} (T_{C_n} - T_S - T_{R_{n+1}}) \\ R_{n+1} &: T_{R_{n+1}} = \dot{q}_{R_n} R_{n+1} \end{align}\]

The equations for \(C_1\) to \(C_n\) can become the state equations by eliminating \(q_{R_i}\) via substitution of the \(R_1\) to \(R_n\) equations. Most of these substitutions are straightforward, but the substitution for \(\dot{q}_{R_n}\) requires two steps. Substituting the \(R_{n+1}\) equation into the equation for \(R_n\) yields can eliminate \(T_{R_{n+1}}\) as follows:

\[ \dot{q}_{R_n} = \frac{1}{R_n} (T_{C_n} - T_S - \dot{q}_{R_n} R_{n+1}). \]

Solving for \(\dot{q}_{R_n}\) yields:

\[ \dot{q}_{R_n} = \frac{1}{R_n + R_{n+1}} (T_{C_n} - T_S). \]

The substitutions are now straightforward and yield the following set of state equations:

\[\begin{align} \dot{T}_{C_1} &= \frac{1}{C_1} \left(\dot{q}_S - \frac{1}{R_1} (T_{C_1} - T_{C_2}) \right) \\ \dot{T}_{C_2} &= \frac{1}{C_2} \left( \frac{1}{R_1} (T_{C_1} - T_{C_2}) - \frac{1}{R_2} (T_{C_2} - T_{C_3}) \right) \\ \vdots &\phantom{:} \vdots \\ \dot{T}_{C_n} &= \frac{1}{C_n} \left( \frac{1}{R_{n-1}} (T_{C_{n-1}} - T_{C_n}) - \frac{1}{R_n + R_{n+1}} (T_{C_n} - T_S) \right) \end{align}\]

Combining terms yields

\[\begin{align} \dot{T}_{C_1} &= \frac{1}{C_1} \left(\dot{q}_S - \frac{1}{R_1} (T_{C_1} - T_{C_2}) \right) \\ \dot{T}_{C_2} &= \frac{1}{C_2} \left( \frac{1}{R_1} T_{C_1} - \left(\frac{1}{R_1} + \frac{1}{R_2}\right) T_{C_2} - \frac{1}{R_2} T_{C_3} \right) \\ \vdots &\phantom{:} \vdots \\ \dot{T}_{C_n} &= \frac{1}{C_n} \left( \frac{1}{R_{n-1}} T_{C_{n-1}} - \left( \frac{1}{R_{n-1}} + \frac{1}{R_n + R_{n+1}} \right) T_{C_n} + \frac{1}{R_n + R_{n+1}} T_S \right) \end{align}\]

The State-Space Model in Standard Form

The state equation has standard form \(\dot{\bm{x}} = A \bm{x} + B \bm{u}\), where \(A\) and \(B\) are matrices containing system parameters. We can transcribe our state equations into this standard form as follows:

\[\begin{align} \begin{bmatrix} \dot{T}_{C_1} \\ \dot{T}_{C_2} \\ \\ \vdots \\ \\ \dot{T_{C_n}} \end{bmatrix} &= \begin{bmatrix} \frac{-1}{R_1 C_1} & \frac{1}{R_1 C_1} & 0 & 0 & \cdots & 0 & 0 & 0 \\ \frac{1}{R_1 C_2} & \frac{-1}{C_2}\left(\frac{1}{R_1} + \frac{1}{R_2}\right) & \frac{1}{R_2 C_2} & 0 & \cdots & 0 & 0 & 0 \\ & & & & & & & \\ \vdots & & \ddots & \ddots & \ddots & & & \vdots \\ & & & & & & & \\ 0 & 0 & 0 & 0 & \cdots & 0 & \frac{1}{R_{n-1} C_n} & \frac{-1}{C_n} \left( \frac{1}{R_{n-1}} + \frac{1}{R_n + R_{n+1}} \right) \end{bmatrix} \begin{bmatrix} T_{C_1} \\ T_{C_2} \\ \\ \vdots \\ \\ T_{C_n} \end{bmatrix} + \\ &+ \begin{bmatrix} \frac{1}{C_1} & 0 \\ 0 & 0 \\ & \\ \vdots & \vdots \\ & \\ 0 & \frac{1}{C_n(R_n + R_{n+1})} \end{bmatrix} \begin{bmatrix} \dot{q}_S \\ T_S \end{bmatrix}. \end{align}\]

The Output Equation in Standard Form

The output equation has standard form \(\bm{y} = C \bm{x} + D \bm{u}\), where \(C\) and \(D\) are matrices containing system parameters. We have said that the output variables are simply \(\bm{y} = \bm{x}\), where \(\bm{x}\) is the state vector, so the output matrix \(C = I_{n \times n}\) is the identity matrix, and the direct transmission matrix \(D = 0_{n \times 2}\) is a zero matrix.

Determining the Parameters

The parameters used in the state-space model are the thermal resistances \(R_i\) and capacitances \(C_i\) of the system. We will try to determine these parameters from the material properties and geometry of the system.

Thermal Capacitance

The lumped thermal capacitance of each element is given by \(C_i = m_i c_i\), where \(m_i\) is the mass of each element and \(c_i\) is the specific heat capacity of the material. For the experiment we will run, the parameters are as follows:

  • Mass of each element: \(m_k = \rho_k V_k = \rho_i \pi (r_k^2 - r_{k-1}^2) \ell / m\), where \(\rho_k\) is the density of the material, \(\rho_k = 8.6 \cdot 10^3 \text{ kg}/\text{m}^3\) (brass), and \(V_k\) is the volume of the element, which we derived above.
  • Specific heat capacity of the material: \(c_k = 380 \text{ J}/(\text{kg}\cdot\text{K})\) for brass.

Conductive Thermal Resistance

The conductive heat transfer between two adjacent elements in the system has resistance \(R_k\) for \(k = 1, 2, \ldots, n\). Let the radial length of each element be \(\delta r = (r_n - r_0)/n\), where \(r_0\) and \(r_n\) are the inner and outer radii of the disk, respectively. Then the length between element centers is \(\delta r\). Let the angle of a sector be \(\theta = 2\pi/m\). For each element \(k\), we allow a potentially different thermal conductivity \(\kappa_k\). Finally, we approximate the cross-sectional area between elements \(k\) and \(k+1\) as \(A_k = r_k \theta \ell\), as described above. Therefore the thermal resistance \(R_k\) is given by \[R_k = \frac{\delta r}{\kappa_k A_k} = \frac{\delta r}{\kappa_k r_k \theta \ell},\] which has units \(\text{K}/\text{W}\). It makes sense that the thermal resistance increases with \(\delta r\) and decreases with the thermal conductivity and cross-sectional area. For the experiment we will run, the parameters are as follows:

  • Disk outer radius: \(r_n = 55\) mm
  • Disk inner radius: \(r_0 = 7\) mm
  • Thickness of the disk: \(\ell = 3.2\) mm
  • Thermal conductivity (brass): \(\kappa_i = 121 \text{ W}/(\text{m}\cdot\text{K})\)

Convective Thermal Resistance

While most surfaces are insulated, the surface of the disk that is exposed to the flow of water is intentionally subject to convective heat transfer. The average water temperature is given by \(T_S\). The convective thermal resistance \(R_{n+1}\) is given by \[R_{n+1} = \frac{1}{h a / m},\] where \(h\) is the convective heat transfer coefficient and \(a\) is the surface area of the disk exposed to water. For the experiment we will run, the parameters are as follows:

  • Convective heat transfer coefficient: \(h = 20{,}000 \text{ W}/(\text{m}^2\cdot\text{K})\), based on the flowrate of water and the geometry of the setup
  • Surface area of the disk exposed to water: \(a = 2\pi r_n \ell\)
  • Temperature of the water: Record in the lab, but typically around \(20^\circ\text{C}\).

Simulating the Experiment

This script simulates transient heat conduction in a brass disk using finite element analysis. The disk is heated from the center and loses heat to the surrounding environment through convection along the outer surface. Due to symmetry, we model this as a 1D heat conduction problem with a heat source inner boundary condition and a convective outer boundary condition. Begin by importing necessary libraries as follows:

import numpy as np
import matplotlib.pyplot as plt
import control

System Parameters

Define the physical properties of the brass disk and the experimental conditions. These parameters are based on typical values for brass and the experimental setup.

rn = 55e-3  # (m) Outer radius of the disk
r0 = 7e-3  # (m) Inner radius of the disk
l = 3.2e-3  # (m) Thickness of the disk
kappak = 121  # (W/(m*K)) Thermal conductivity of the disk (brass)
ck = 380  # (J/(kg*K)) Specific heat capacity of the disk (brass)
rhok = 8.6e3  # (kg/m^3) Density of the disk (brass)
h = 20_000  # (W/(m^2*K)) Convective heat transfer coefficient
a = 2 * np.pi * rn * l  # (m^2) Surface area of the disk exposed to convective heat transfer
v = 120  # (V) RMS heater voltage
i = 0.246  # (A) RMS heater current
qS_total = v * i  # (W) Power generated by the heater
Tinf = (66 - 32) * 5 / 9   # (C) Water temperature

Finite Element Discretization

Divide the disk into finite elements and calculate the thermal properties for each element. The finite element method allows us to solve the partial differential equation numerically.

m = 50  # Number of sectors
theta = 2 * np.pi / m  # Angle of each sector
n = 50  # Number of finite elements per sector
dr = (rn - r0) / n  # Finite element length
print(f"Finite element radial length δr: {dr * 1e3:0.3g} mm")
r_array = np.linspace(r0, rn, n + 1)  # Radii of interfaces between finite elements
C_array = np.zeros(n)
R_array = np.zeros(n + 1)
for k in range(0, n):  # Loop over finite elements
    Vk = np.pi * (r_array[k+1] ** 2 - r_array[k] ** 2) * l / m
    mk = rhok * Vk  # Mass of the finite element
    C_array[k] = mk * ck
    Ak = r_array[k+1] * theta * l
    R_array[k] = dr / (kappak * Ak)
    if k == n - 1:
        R_array[k] = R_array[k] / 2  # Half width distance to the edge
R_array[n] = 1 / (h * a / m)
Finite element radial length δr: 0.96 mm

State-Space Model Construction

Create the state-space matrices A, B, C, and D for the thermal system. This represents the system of ordinary differential equations after spatial discretization.

A = np.zeros((n, n))
for k in range(0, n):
    if k == 0:
        A[0, 0] = -1 / (R_array[k] * C_array[k])
        A[0, 1] = 1 / (R_array[k] * C_array[k])
    elif k == n - 1:
        A[n - 1, n - 2] = 1 / (R_array[k-1] * C_array[k])
        A[n - 1, n - 1] = -1 / C_array[k] * (1/R_array[k-1] + 1/(R_array[k] + R_array[k+1]))
    else:
        A[k, k - 1] = 1 / (R_array[k-1] * C_array[k])
        A[k, k] = -1 / C_array[k] * (1/R_array[k-1] + 1/R_array[k])
        A[k, k + 1] = 1 / (R_array[k] * C_array[k])

Input and Output Matrix Definition

Define how external inputs (heating power and ambient temperature) affect the system and how we observe the system outputs (temperature at each element).

B = np.zeros((n, 2))
B[0, 0] = 1 / (C_array[0])
B[-1, 1] = 1 / (C_array[n-1] * (R_array[n-1] + R_array[n]))

C = np.eye(n)  # Identity matrix of size n
D = np.zeros((n, 2))

Create Control System Model

Assemble the state-space model using the control systems library. This allows us to easily simulate the system response to various inputs.

sys = control.ss(A, B, C, D)

Simulation Function Definition

Create a flexible simulation function that can handle different heating and cooling durations. This allows us to match experimental conditions or explore different heating scenarios.

def simulate(heat_time, cool_time, qS, nt=200):
    """
    Simulate the system response with specified heat and cool times.

    Parameters:
    - heat_time: Duration of heating phase (seconds)
    - cool_time: Duration of cooling phase (seconds)
    - nt: Number of time points for simulation

    Returns:
    - t: Time array (seconds)
    - y_sim: Temperature response array (C), shape (n_elements, n_time)
    """
    total_time = heat_time + cool_time
    t = np.linspace(0, total_time, nt)

    # Create input array
    u = np.zeros((2, nt))
    heat_end_index = int(np.floor((heat_time / total_time) * nt))
    u[0, :heat_end_index] = qS  # Heat input for specified duration
    u[1, :] = Tinf  # Ambient temperature

    # Initial conditions (all elements start at ambient temperature)
    T0 = np.ones(n) * Tinf

    # Run simulation
    y_sim = control.forced_response(sys, t, u, x0=T0).outputs

    return t, y_sim

Temperature Interpolation Function

Create an interpolation function to estimate temperature at any position and time. This is useful for comparing simulation results with experimental sensor data that may not be located at finite element centers.

def simulation_interpolate(time, position, t_sim, y_sim):
    """
    Interpolate the simulation results at given time and position arrays.
    Uses sequential 1D interpolation with linear extrapolation for positions
    outside the finite element grid.

    Parameters:
    - time: Time point(s) for interpolation (scalar or array)
    - position: Position(s) along disk for interpolation (scalar or array)
    - t_sim: Time array from simulation
    - y_sim: Temperature response array from simulation

    Returns:
    - Interpolated temperature values
    """
    from scipy.interpolate import interp1d

    # Define the grid for interpolation
    # Position grid: center of each finite element
    r_grid = np.linspace(r0 + dr/2, rn - dr/2, n)

    # Handle scalar vs array inputs
    time_is_scalar = np.isscalar(time)
    position_is_scalar = np.isscalar(position)

    # Convert to arrays for consistent processing
    time_array = np.atleast_1d(time)
    position_array = np.atleast_1d(position)

    # Sequential 1D interpolation
    result = np.zeros((len(time_array), len(position_array)))

    # Interpolate for each time point
    for i, t in enumerate(time_array):
        # Step 1: Interpolate in time to get temperature profile at desired time
        if t <= t_sim[0]:
            # Use first time point
            temp_profile = y_sim[:, 0]
        elif t >= t_sim[-1]:
            # Use last time point
            temp_profile = y_sim[:, -1]
        else:
            # Interpolate in time
            time_interp = interp1d(t_sim, y_sim, axis=1, kind='linear',
                                 bounds_error=False, fill_value='extrapolate')
            temp_profile = time_interp(t)

        # Step 2: Interpolate in space with extrapolation
        space_interp = interp1d(r_grid, temp_profile, kind='linear',
                               bounds_error=False, fill_value='extrapolate')
        result[i, :] = space_interp(position_array)

    # Return in appropriate format based on input types
    if time_is_scalar and position_is_scalar:
        return float(result[0, 0])
    elif time_is_scalar:
        return result[0, :]
    elif position_is_scalar:
        return result[:, 0]
    else:
        return result

Sensor Position Definition

Define the physical positions of temperature sensors along the disk. These positions match the experimental setup and will be used for predictions.

sensor_positions = 1e-3 * np.array([7, 10, 20, 30, 40, 50])
n_sensors = len(sensor_positions)

Example Simulation

Run a demonstration simulation with typical experimental parameters. This shows the temperature evolution during heating and cooling phases.

# Example simulation parameters
heat_time_ = 4 * 60  # (s) heating time
cool_time_ = 2 * 60  # (s) cooling time

# Run simulation
t_, y_ = simulate(heat_time_, cool_time_, qS_total / m)

Temperature Predictions at Sensor Locations

Use interpolation to predict temperatures at the actual sensor positions. This demonstrates how the simulation would compare with experimental measurements.

# Predict temperatures at sensor positions
T_sensor_predicted = np.zeros((len(t_), n_sensors))
for i, t in enumerate(t_):
    T_sensor_predicted[i, :] = simulation_interpolate(t, sensor_positions, t_, y_)

Visualization of Sensor Predictions

Plot the predicted temperature evolution at sensor locations. This shows what we would expect to measure experimentally.

fig, ax = plt.subplots(figsize=(10, 6))

# Define colors for each sensor
colors = plt.cm.tab10(np.linspace(0, 1, n_sensors))

# Plot temperature predictions at sensor positions
for i in range(n_sensors):
    position_mm = sensor_positions[i] * 1000  # Convert to mm
    ax.plot(t_ / 60, T_sensor_predicted[:, i], color=colors[i],
            label=f'Sensor {i+1} ({position_mm:.0f} mm)', linewidth=2)

ax.axvline(heat_time_ / 60, color='red', linestyle='--', alpha=0.7, label='Heating ends')
ax.set_xlabel('Time (min)')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Predicted Temperature Evolution at Sensor Locations')
ax.grid(True, alpha=0.3)
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
Predicted temperature evolution at sensor locations during heating and cooling phases

Temperature Profile Visualization

Show how the temperature profile along the disk changes at different times. This provides insight into the spatial temperature distribution.

fig, ax = plt.subplots(figsize=(10, 6))

# Select specific time points to show profiles
time_points = [0, heat_time_/4, heat_time_/2, heat_time_,
               heat_time_ + cool_time_/2, heat_time_ + cool_time_]
time_labels = ['Start', 'Heat 25%', 'Heat 50%', 'Heat End', 'Cool 50%', 'Cool End']

# Define element positions for plotting
element_positions = np.linspace(r0 + dr/2, rn - dr/2, n)

for i, (t_point, label) in enumerate(zip(time_points, time_labels)):
    t_idr = np.argmin(np.abs(t_ - t_point))
    ax.plot(element_positions * 1000, y_[:, t_idr], 'o-', label=label)

ax.set_xlabel('Position (mm)')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Temperature Profiles at Different Times')
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
Temperature profiles along the disk at different times during heating and cooling

Sensor Layout Visualization

Show the relationship between finite elements and sensor positions. This helps understand how interpolation maps simulation results to measurement locations.

fig, ax = plt.subplots(figsize=(12, 6))

# Show finite element positions
ax.scatter(element_positions * 1000, np.zeros_like(element_positions),
           marker='|', s=100, c='lightblue', alpha=0.6, label=f'{n} Finite Elements')

# Show sensor positions with different colors
for i, pos in enumerate(sensor_positions):
    ax.scatter(pos * 1000, 0, marker='o', s=20, c=[colors[i]],
               label=f'Sensor {i+1}', alpha=0.9, edgecolors='black', linewidth=1)

ax.set_xlabel('Position (mm)')
ax.set_title('Finite Element Grid and Sensor Positions')
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_ylim(-0.5, 0.5)
ax.set_ylabel('')
ax.set_yticks([])
plt.tight_layout()
plt.show()
Finite element grid and sensor positions showing spatial discretization and measurement locations

2D Temperature Distribution Heatmaps

Create a 2x3 grid of heatmaps showing the temperature distribution at key moments during heating and cooling phases. The heatmaps reconstruct the 2D temperature field from the 1D simulation using cylindrical symmetry.

def create_2d_temperature_field(time_point, n_angular=100):
    """
    Create a 2D temperature field at a specific time using cylindrical symmetry.

    Parameters:
    - time_point: Time at which to create the temperature field
    - n_angular: Number of angular divisions for 2D field (not used, kept for compatibility)

    Returns:
    - x, y: 2D coordinate grids
    - T_2d: 2D temperature field
    """
    # Create 2D coordinate system
    x = np.linspace(-rn, rn, 201)
    y = np.linspace(-rn, rn, 201)
    X, Y = np.meshgrid(x, y)
    R = np.sqrt(X**2 + Y**2)

    # Create 2D temperature field using simulation_interpolate
    T_2d = np.zeros_like(R)
    mask_valid = (R >= r0) & (R <= rn)  # Between inner and outer radius
    T_2d[mask_valid] = simulation_interpolate(time_point, R[mask_valid], t_, y_)
    T_2d[~mask_valid] = np.nan  # Outside valid region (inner hole and outer boundary)

    return X, Y, T_2d
fig, axes = plt.subplots(2, 3, figsize=(10, 6))

# Define key time points for visualization
heating_times = [heat_time_/4, heat_time_/2, heat_time_*0.995]  # 25%, 50%, 100% of heating
cooling_times = [heat_time_ + cool_time_/4, heat_time_ + cool_time_/2, heat_time_ + cool_time_*0.995]  # 25%, 50%, 100% of cooling

all_times = heating_times + cooling_times
time_labels = ['25% Heating', '50% Heating', '100% Heating',
               '25% Cooling', '50% Cooling', '100% Cooling']

# Define element positions for temperature field creation
element_positions = np.linspace(r0 + dr/2, rn - dr/2, n)

# Find global temperature range for consistent colorbar across all simulation data
vmin = np.min(y_)
vmax = np.max(y_)

# Create normalization for consistent colorbar
from matplotlib.colors import Normalize
norm = Normalize(vmin=vmin, vmax=vmax)

# Create heatmaps
im_list = []
for i, (t_point, label, ax) in enumerate(zip(all_times, time_labels, axes.flat)):
    # Create 2D temperature field
    X, Y, T_2d = create_2d_temperature_field(t_point)

    # Create heatmap with shared normalization
    im = ax.contourf(X*1000, Y*1000, T_2d, levels=50, cmap='plasma', norm=norm)
    im_list.append(im)

    # Add contour lines
    contours = ax.contour(X*1000, Y*1000, T_2d, levels=10, colors='white', alpha=0.3, linewidths=0.5)

    # Add inner and outer boundaries
    circle_inner = plt.Circle((0, 0), r0*1000, fill=False, color='white', linewidth=2)
    circle_outer = plt.Circle((0, 0), rn*1000, fill=False, color='white', linewidth=2)
    ax.add_patch(circle_inner)
    ax.add_patch(circle_outer)

    # Formatting
    ax.set_aspect('equal')
    ax.set_xlim(-rn*1000*1.1, rn*1000*1.1)
    ax.set_ylim(-rn*1000*1.1, rn*1000*1.1)
    ax.set_xlabel('x (mm)')
    ax.set_ylabel('y (mm)')
    ax.set_title(f'{label}\n(t = {t_point/60:.1f} min)')
    ax.grid(True, alpha=0.3)

plt.tight_layout()

# Add single shared colorbar for all subplots
cbar = fig.colorbar(im_list[0], ax=axes, shrink=0.8, aspect=30, pad=0.02)
cbar.set_label('Temperature (°C)', fontsize=12)
plt.show()
2D temperature distribution heatmaps showing key moments during heating (upper row) and cooling (lower row) phases. The heatmaps reconstruct the temperature field from the 1D simulation using cylindrical symmetry.

3D Surface Temperature Distribution Snapshots

Create 3D surface plots showing multiple snapshots of temperature distributions. Left subplot shows heating phase, right subplot shows cooling phase with semitransparent meshgrid plots.

def create_3d_temperature_wireframe(time_point, n_radial=20, n_angular=30):
    """
    Create 3D temperature wireframe data at a specific time using polar coordinates.

    Parameters:
    - time_point: Time at which to create the temperature wireframe
    - n_radial: Number of radial grid lines
    - n_angular: Number of angular grid lines

    Returns:
    - X, Y: 2D coordinate grids in polar pattern
    - Z: 2D temperature field (with NaN for invalid regions)
    """
    # Create polar coordinate system
    r = np.linspace(r0, rn, n_radial)
    theta = np.linspace(0, 2*np.pi, n_angular)
    R, Theta = np.meshgrid(r, theta)

    # Convert to cartesian coordinates
    X = R * np.cos(Theta)
    Y = R * np.sin(Theta)

    # Create temperature field using simulation_interpolate
    # Flatten R for interpolation, then reshape back to original shape
    R_flat = R.flatten()
    Z_flat = simulation_interpolate(time_point, R_flat, t_, y_)
    Z = Z_flat.reshape(R.shape)

    return X, Y, Z
fig = plt.figure(figsize=(10, 6))

# Define time snapshots for 3D visualization
heating_snapshots = [heat_time_/4, heat_time_/2, heat_time_*0.995]
cooling_snapshots = [heat_time_ + cool_time_/4, heat_time_ + cool_time_/2, heat_time_ + cool_time_*0.995]

# Colors for different snapshots
heating_colors = ['red', 'orange', 'darkred']
cooling_colors = ['blue', 'cyan', 'darkblue']
alphas = [0.4, 0.6, 0.8]  # Increasing opacity for time progression

# Left subplot: Heating phase
ax1 = fig.add_subplot(121, projection='3d')

for i, (t_snap, color, alpha) in enumerate(zip(heating_snapshots, heating_colors, alphas)):
    X, Y, Z = create_3d_temperature_wireframe(t_snap)

    # Create wireframe plot
    wire = ax1.plot_wireframe(X*1000, Y*1000, Z, color=color, alpha=alpha,
                             linewidth=1.5, antialiased=True)

ax1.set_xlabel('x (mm)')
ax1.set_ylabel('y (mm)')
ax1.set_zlabel('Temperature (°C)')
ax1.set_title('Heating Phase\nTemperature Evolution', fontsize=12, fontweight='bold')
ax1.view_init(elev=15, azim=45)

# Right subplot: Cooling phase
ax2 = fig.add_subplot(122, projection='3d')

for i, (t_snap, color, alpha) in enumerate(zip(cooling_snapshots, cooling_colors, alphas)):
    X, Y, Z = create_3d_temperature_wireframe(t_snap)

    # Create wireframe plot
    wire = ax2.plot_wireframe(X*1000, Y*1000, Z, color=color, alpha=alpha,
                             linewidth=1.5, antialiased=True)

ax2.set_xlabel('x (mm)')
ax2.set_ylabel('y (mm)')
ax2.set_zlabel('Temperature (°C)')
ax2.set_title('Cooling Phase\nTemperature Evolution', fontsize=12, fontweight='bold')
ax2.view_init(elev=15, azim=45)

# Set consistent z-axis limits for both subplots
z_min, z_max = vmin, vmax
ax1.set_zlim(z_min, z_max)
ax2.set_zlim(z_min, z_max)

# Add time legends
heating_labels = [f'Heating {100 * t / heat_time_:.1f}%' for t in heating_snapshots]
cooling_labels = [f'Cooling {100 * (t - heat_time_) / cool_time_:.1f}%' for t in cooling_snapshots]

# Create custom legend elements
from matplotlib.patches import Rectangle
heating_legend = [Rectangle((0,0),1,1, facecolor=color, alpha=alpha, label=label)
                 for color, alpha, label in zip(heating_colors, alphas, heating_labels)]
cooling_legend = [Rectangle((0,0),1,1, facecolor=color, alpha=alpha, label=label)
                 for color, alpha, label in zip(cooling_colors, alphas, cooling_labels)]

ax1.legend(handles=heating_legend, loc='upper left', bbox_to_anchor=(0, 1))
ax2.legend(handles=cooling_legend, loc='upper left', bbox_to_anchor=(0, 1))

plt.tight_layout()
plt.show()
3D wireframe plots showing temperature distribution snapshots during heating (left) and cooling (right) phases. Polar grid pattern with constant radius and angle lines shows the cylindrical symmetry. Multiple semitransparent wireframes represent different time points, with increasing opacity indicating time progression.

Bibliography

  1. Hilton, P. A.. "Optional Radial Heat Conduction Unit H112B: Experimental Operating and Maintenance Manual".
  2. Hilton, P. A.. "Heat Transfer Service Unit H112: Experimental Operating and Maintenance Manual".