%% Steady-State 2D Light Bulb Model

% Authors; Chad Clogston, Ryan Cox, Nathan Wilde

clear;

clc;

% Code Description:

% This script solves the steady-state 2D heat conduction problem for a

% small cross-section of a 40 W T10 tubular incandescent bulb using a

% 15-node network representing the tungsten filament, the inert bulb gas

% (argon/nitrogen, modeled with air-like properties), and the glass

% envelope. Energy balances at each node include conduction between

% neighboring nodes, internal heat generation in the filament, a fixed

% base temperature, and combined convection and linearized radiation at

% the outer glass surface; these balances are assembled into a linear

% system A*T = C whose solution gives the steady-state nodal

% temperatures and corresponding temperature distribution plots.

%% Geometry/grid

dx = 0.002; % m

dy = dx; % m

b = 0.01; % m

% Material properties

k_fil = 30; % W/m-K tungsten

k_gas = 0.026; % W/m-K gas/air region

k_glass = 1.4; % W/m-K soda-lime glass

% Heat generation in filament region

P_bulb = 40; % W, total power

N_fil = 3; % number of filament nodes

V_node = b dx dy; % m^3 volume per node

V_gen = N_fil * V_node; % total generating volume

q_dot = P_bulb / V_gen; % W/m^3 applied to nodes 1–3

%% 4. External convection & radiation

T_inf = 295; % K, ambient temperature

h_conv = 10; % W/m^2-K, convection coefficient

epsilon_glass = 0.90; % glass emissivity

sigma = 5.67e-8; % W/m^2-K^4

T_ref = T_inf; % linearization temperature for radiation

h_rad = epsilon_glass sigma (T_ref + T_ref) * (T_ref^2 + T_ref^2);

h_total = h_conv + h_rad; % convection + radiation

T_base = T_inf;

%% 5. Matrix A and Vector C

A = zeros(15,15);

C = zeros(15,1);

%% 6. Nodal energy-balance equations

% Node 1

A(1,1) = -(3*k_fil + k_gas);

A(1,2) = k_fil;

A(1,7) = k_gas;

C(1) = -2*k_fil*T_base - q_dot*dx^2;

% Node 2

A(2,1) = k_fil;

A(2,2) = -(2*k_fil + k_gas);

A(2,3) = k_fil;

A(2,8) = k_gas;

C(2) = -q_dot*dx^2;

% Node 3

A(3,2) = k_fil;

A(3,3) = -(k_fil + 2*k_gas);

A(3,4) = k_gas;

A(3,9) = k_gas;

C(3) = -q_dot*dx^2;

% Node 4

A(4,3) = k_fil;

A(4,4) = -(k_gas + k_glass + k_fil);

A(4,5) = k_glass;

A(4,10) = k_gas;

C(4) = 0;

% Node 5

A(5,4) = k_gas;

A(5,5) = -(k_gas + 2*k_glass);

A(5,6) = k_glass;

A(5,11) = k_glass;

C(5) = 0;

% Node 6

A(6,5) = k_glass;

A(6,6) = -(k_glass + 2*h_total*dx);

C(6) = -2*h_total*dx*T_inf;

% Node 7

A(7,1) = k_fil;

A(7,7) = -(k_fil + 3*k_gas + k_glass);

A(7,8) = k_gas;

A(7,12) = k_glass;

C(7) = -2*k_gas*T_base;

% Node 8 (interior gas / glass)

A(8,2) = k_fil;

A(8,7) = k_gas;

A(8,8) = -(2*k_gas + k_fil + k_glass);

A(8,9) = k_gas;

A(8,13) = k_glass;

C(8) = 0;

% Node 9

A(9,3) = k_fil;

A(9,8) = k_gas;

A(9,9) = -(2*k_gas + k_fil + k_glass);

A(9,10) = k_gas;

A(9,14) = k_glass;

C(9) = 0;

% Node 10

A(10,4) = k_gas;

A(10,9) = k_gas;

A(10,10) = -(k_glass + 3*k_gas);

A(10,11) = k_gas;

A(10,15) = k_glass;

C(10) = 0;

% Node 11

A(11,5) = k_glass;

A(11,10) = k_gas;

A(11,11) = -(k_gas + k_glass + 2*h_total*dy);

C(11) = -2*h_total*dy*T_inf;

% Node 12

A(12,7) = k_gas;

A(12,12) = -(3*k_gas + k_glass + h_total*dx);

A(12,13) = k_glass;

C(12) = -h_total*dy*T_inf - 2*k_gas*T_base;

% Node 13 (

A(13,8) = k_gas;

A(13,12) = k_glass;

A(13,13) = -(k_gas + 2*k_glass + h_total*dx);

A(13,14) = k_glass;

C(13) = -h_total*dy*T_inf;

% Node 14

A(14,9) = k_gas;

A(14,13) = k_glass;

A(14,14) = -(k_gas + 2*k_glass + h_total*dx);

A(14,15) = k_glass;

C(14) = -h_total*dy*T_inf;

% Node 15

A(15,10) = k_gas;

A(15,14) = k_glass;

A(15,15) = -(k_glass + k_gas + 2*h_total*dy);

C(15) = -2*h_total*dy*T_inf;

%% 7. Solve linear system

T = A \ C;

fprintf('Nodal temperatures for Light Bulb (K):\n');

for n = 1:15

fprintf('T%-2.0f = %8.2f K\n', n, T(n));

end

%% 8. Map to nodal layout and plot contour

Tplot = zeros(3,6);

Tplot(3,1:4) = [T(12) T(13) T(14) T(15)];

Tplot(2,1:5) = [T(7) T(8) T(9) T(10) T(11)];

Tplot(1,1:6) = [T(1) T(2) T(3) T(4) T(5) T(6)];

[x,y] = meshgrid(1:6,1:3);

[xq,yq] = meshgrid(linspace(1,6,100), linspace(1,3,100));

Tinterp = griddata(x, y, Tplot, xq, yq, 'linear');

figure;

contourf(xq,yq,Tinterp,30,'LineColor','none');

hold on;

contour(xq,yq,Tinterp,15,'k');

hold off;

colorbar; colormap(jet);

title('Steady-State Temperature Distribution (K)');

xlabel('Node Position (x)');

ylabel('Node Position (y)');

set(gca,'YDir','normal');

axis equal; grid on;