Thermal Element Formulation and Stiffness

TL;DR

You'll learn how to model heat transfer in materials by breaking them into smaller pieces called thermal elements. We'll focus on creating equations that describe heat flow within these elements and how they connect to form a system. This approach lets you calculate temperature distributions and heat fluxes in complex structures.

1. The Mental Model

Imagine you want to know how heat spreads through a metal bar. Instead of solving one giant, complex equation for the whole bar, you cut it into many small, connected slices. Each slice is a "thermal element," and we figure out how heat moves in and out of each slice, then combine them to understand the whole bar.

2. The Core Material

When we talk about thermal element formulation, we're essentially developing the mathematical equations that describe heat transfer within a single, small piece of a material. This piece is our "element." The "stiffness" part then refers to how these elements resist (or conduct) heat flow.

Discretization and Nodal Temperatures

Photo by Matheus Bertelli on Pexels

The first step is called discretization. This means we take a continuous object and divide it into a finite number of discrete elements connected at specific points called nodes. For thermal problems, the primary unknown we're trying to find at these nodes is temperature.

Consider a 1D rod. We could discretize it into several 1D line elements. Each element has two nodes, one at each end.

graph TD
    Rod["Continuous Rod"] --> Discretization["Discretization (into elements)"]
    Discretization --> Element1["Element 1"]
    Discretization --> Element2["Element 2"]
    Discretization --> ElementN["... Element N"]
    Element1 --> NodeA["Node A (T_A)"]
    Element1 --> NodeB["Node B (T_B)"]
    Element2 --> NodeB
    Element2 --> NodeC["Node C (T_C)"]
    NodeA -- "Heat Flux" --> NodeB
    NodeB -- "Heat Flux" --> NodeC

Here, T_A, T_B, T_C are the temperatures at the nodes.

Heat Transfer Mechanisms

Photo by Nic Wood on Pexels

Thermal elements typically account for:
* Conduction: Heat transfer through direct contact. This is usually the dominant mechanism we model within the element's stiffness matrix.
* Convection: Heat transfer between a solid surface and a moving fluid (like air or water). This often appears as a boundary condition or can be integrated into the element formulation for specific elements.
* Radiation: Heat transfer via electromagnetic waves. This is often non-linear and handled separately or as a boundary condition.

For basic thermal element stiffness, we usually focus on conduction.

Element Stiffness Matrix (1D Conduction)

Photo by The Six on Pexels

Let's look at a simple 1D line element with two nodes, i and j, at positions x_i and x_j. The length of the element is L = x_j - x_i.
Assume:
* Constant cross-sectional area A.
* Constant thermal conductivity k.
* No internal heat generation.
* Steady-state conditions (temperatures don't change with time).

The heat flow q through the element is given by Fourier's Law of Conduction:
q = -kA (dT/dx)

If we make a linear approximation for temperature variation within the element (T(x) = N_i T_i + N_j T_j, where N_i and N_j are shape functions), we can derive the element stiffness matrix [k_e].

For a 1D element, the relationship between nodal heat fluxes (q_i, q_j) and nodal temperatures (T_i, T_j) is:

{q_e} = [k_e] {T_e}

Where q_e is the vector of nodal heat fluxes, T_e is the vector of nodal temperatures, and [k_e] is the element thermal stiffness matrix.

For a 1D two-node element, [k_e] is:

[k_e] = (kA/L) * [[ 1, -1 ], [ -1, 1 ]]

Let's break that down:
* kA/L: This term represents the thermal conductance of the element. A higher k or A means more heat can flow, so it's "stiffer" thermally. A longer L means more resistance, so it's "less stiff."
* The [[ 1, -1 ], [ -1, 1 ]] part ensures that if T_i = T_j, there's no net heat flow. If T_i > T_j, heat flows from i to j (positive q_i, negative q_j).

Global Stiffness Matrix Assembly

Photo by Markus Spiske on Pexels

Once you have the stiffness matrix for each individual element, the next step is to combine them into a global stiffness matrix [K] for the entire system. This process is called assembly.

The assembly process is essentially summing the contributions of each element's stiffness matrix into the correct positions in the larger global matrix. For example, if element 1 connects nodes 1 and 2, and element 2 connects nodes 2 and 3, the stiffness contribution at node 2 from both elements gets added together.

The global system equation looks like this:

{Q} = [K] {T}

Where:
* {Q} is the global vector of applied heat loads or fluxes at each node.
* [K] is the global thermal stiffness matrix.
* {T} is the global vector of unknown nodal temperatures.

Boundary Conditions

To solve for the unknown temperatures, you need to apply boundary conditions. These define how the system interacts with its surroundings. Common types include:
* Prescribed temperature (Dirichlet boundary condition): You know the temperature at a specific node (e.g., T_1 = 100°C).
* Prescribed heat flux (Neumann boundary condition): You know the heat entering or leaving a node (e.g., Q_2 = 50 W). A perfectly insulated surface means Q = 0.
* Convection boundary condition: Heat transfer to or from an ambient fluid (e.g., h A (T_surface - T_ambient)). This term can be added to the nodal heat flux vector and often modifies the stiffness matrix itself.

Solving the System

Once the global system of equations is assembled and boundary conditions are applied, you solve [K] {T} = {Q} for the unknown temperatures T. This usually involves matrix inversion or other numerical methods.

3. Worked Example

Let's consider a simple steel rod, 0.2 m long, with a constant cross-sectional area A = 0.001 m^2.
Its thermal conductivity k = 50 W/(m·K).

We'll discretize it into two elements of equal length.
* Element 1: Nodes 1 and 2.
* Element 2: Nodes 2 and 3.

Each element length L = 0.1 m.

1. Element Stiffness Matrices:

For Element 1 (connecting nodes 1 and 2):
k_e1 = (kA/L) = (50 * 0.001 / 0.1) = 0.5 W/K
[k_e1] = 0.5 * [[ 1, -1 ], [ -1, 1 ]] = [[ 0.5, -0.5 ], [ -0.5, 0.5 ]] (for nodes 1, 2)

For Element 2 (connecting nodes 2 and 3):
k_e2 = (kA/L) = (50 * 0.001 / 0.1) = 0.5 W/K
[k_e2] = 0.5 * [[ 1, -1 ], [ -1, 1 ]] = [[ 0.5, -0.5 ], [ -0.5, 0.5 ]] (for nodes 2, 3)

2. Global Stiffness Matrix Assembly:

We need a 3x3 global stiffness matrix for nodes 1, 2, and 3.
The general form is:
[[ K_11, K_12, K_13 ], [ K_21, K_22, K_23 ], [ K_31, K_32, K_33 ]]

• Element 1 contributes to K_11, K_12, K_21, K_22.
• Element 2 contributes to K_22, K_23, K_32, K_33.

[K] = [[ 0.5, -0.5, 0 ],
[ -0.5, (0.5+0.5), -0.5 ],
[ 0, -0.5, 0.5 ]]

[K] = [[ 0.5, -0.5, 0 ],
[ -0.5, 1.0, -0.5 ],
[ 0, -0.5, 0.5 ]]

3. Apply Boundary Conditions:

Let's say:
* Node 1 is held at T_1 = 100 °C.
* Node 3 is exposed to a convection boundary where h = 10 W/(m^2·K) and T_ambient = 20 °C. (This is simplified; usually convection is added as a flux or to the matrix for edge elements). For simplicity here, let's assume Q_3 due to convection. A more advanced treatment would add a term hA to K_33 and hA * T_ambient to Q_3.
* Let's simplify and say Q_3 = hA(T_ambient - T_3)
A more typical and simpler boundary condition for Q_3 for this example would be Q_3 = 0 (insulated) or a fixed flux.

Let's use a simpler boundary condition for the example:
* T_1 = 100 °C (fixed temperature).
* Q_3 = 0 (node 3 is perfectly insulated).
* Q_2 = 0 (no external heat input at node 2).

Our system equation:
[[ 0.5, -0.5, 0 ],
[ -0.5, 1.0, -0.5 ],
[ 0, -0.5, 0.5 ]] * [[ T1 ], [ T_2 ], [ T_3 ]] = [[ Q_1 ], [ Q_2 ], [ Q_3 ]]`

We know T_1 = 100. We know Q_2 = 0 and Q_3 = 0. Q_1 is unknown – it's the reaction heat flux required to maintain T_1.

Let's modify the system to incorporate T_1 = 100.
From the first equation: 0.5*T_1 - 0.5*T_2 = Q_1
From the second equation: -0.5*T_1 + 1.0*T_2 - 0.5*T_3 = Q_2
From the third equation: -0.5*T_2 + 0.5*T_3 = Q_3

Substitute T_1 = 100, Q_2 = 0, Q_3 = 0:

0.5*(100) - 0.5*T_2 = Q_1 (1)
-0.5*(100) + 1.0*T_2 - 0.5*T_3 = 0 (2)
-0.5*T_2 + 0.5*T_3 = 0 (3)

From (3): 0.5*T_3 = 0.5*T_2, so T_3 = T_2.
Substitute T_3 = T_2 into (2):
-50 + 1.0*T_2 - 0.5*T_2 = 0
-50 + 0.5*T_2 = 0
0.5*T_2 = 50
T_2 = 100 °C

Since T_3 = T_2, then T_3 = 100 °C.

So, the temperature distribution is T_1 = 100 °C, T_2 = 100 °C, T_3 = 100 °C.
This makes sense: if one end is held at 100°C and the other is insulated with no internal heat generation, and no heat is lost from the sides, the entire rod will eventually reach a uniform temperature of 100°C.

We can also find Q_1 using equation (1):
0.5*(100) - 0.5*(100) = Q_1
50 - 50 = Q_1
Q_1 = 0 W
This also makes sense: once the rod is at a uniform temperature, no more heat needs to flow in to maintain T_1.

4. Key Takeaways

• You break down complex shapes into simple elements for easier analysis (discretization).
• Each element has a stiffness matrix that defines how it resists/conducts heat flow based on material properties and geometry.
• The element stiffness matrices are combined into a global stiffness matrix, representing the entire system.
• Boundary conditions (known temperatures or heat fluxes) are essential to solve for unknown temperatures.
• The solution gives you the temperature at each node, from which you can find heat fluxes within elements.

Common Mistakes to Avoid:
- Incorrect assembly: Make sure element contributions are added to the correct global matrix positions, especially at shared nodes.
- Wrong units: Always be careful with units (e.g., thermal conductivity W/(m·K), length m, area m^2).
- Missing boundary conditions: Without sufficient boundary conditions, your system of equations won't have a unique solution.
- Ignoring convection/radiation: For many real-world problems, these heat transfer modes are significant and must be included, often as boundary conditions or modified stiffness terms.