%% Light Bulb Transient – Explicit Finite Difference

% Authors; Chad Clogston, Ryan Cox, Nathan Wilde

% This script computes the transient temperature response of a 2D

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

% explicit finite-difference method on a 15-node network representing

% the tungsten filament, inert bulb gas (argon/nitrogen, modeled with

% air-like properties), and the glass envelope. Starting from a uniform

% initial temperature, each timestep updates the nodal temperatures

% based on the transient energy balance with conduction, filament heat

% generation, a fixed base temperature, and combined convection and

% linearized radiation at the outer glass surface; the stored results

% are used to generate time histories, a final contour plot, and an

% animation of the warm-up process.

clear;

clc;

close all;

%% 1. Material properties

% Filament (tungsten)

k_fil = 30; % W/m-K

rho_fil = 19300; % kg/m^3

cp_fil = 134; % J/kg-K

% Bulb gas (argon/nitrogen, modeled as air-like)

k_gas = 0.026; % W/m-K

rho_gas = 1.2; % kg/m^3

cp_gas = 1000; % J/kg-K

% Glass envelope (soda-lime glass)

k_glass = 1.4; % W/m-K

rho_glass= 2500; % kg/m^3

cp_glass = 840; % J/kg-K

%% 2. Geometry and boundary conditions

dx = 0.002; % m, grid spacing in x

dy = dx; % m, grid spacing in y

b = 0.01; % m, depth into page

% Representative bulb (40 W T10 tubular)

P_bulb = 40; % W, power applied

D_outer = 0.032; % m, 1.25 in diameter

L_bulb = 0.13; % m, 5 in length

% Thermal boundary conditions

T_inf = 295; % K, ambient temperature

T_base = 295; % K, fixed base temperature

h_conv = 10; % W/m^2-K, convection at outer glass

epsilon_glass = 0.90;

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

T_ref = T_inf;

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

h_t = h_conv + h_rad; % combined convection + radiation coefficients

%% 3. Volumetric heat generation in filament (nodes 1–3)

N_fil = 3; % 3 filament nodes

V_node = b dx dy; % node control volume

V_gen = N_fil * V_node;

q_dot = P_bulb / V_gen; % W/m^3

%% 4. Time step based on stability

alpha_fil = k_fil /(rho_fil*cp_fil);

alpha_gas = k_gas /(rho_gas*cp_gas);

alpha_glass = k_glass/(rho_glass*cp_glass);

alpha_max = max([alpha_fil alpha_gas alpha_glass]);

dt_max = 0.25 * dx^2 / alpha_max;

dt = 0.5*1e-5; % conservative choice (half of max)

t_final = 10; % seconds

nt = round(t_final/dt);

fprintf('alpha_max = %.3e m^2/s\n', alpha_max);

fprintf('dt_max = %.3e s\n', dt_max);

fprintf('dt = %.3e s, nt = %d\n', dt, nt);

%% 5. Initial condition

T0 = T_inf; % K

T = T0 * ones(15,1); % 15-node temperature vector

T_history = zeros(15,nt+1);

T_history(:,1) = T;

%% 6. Explicit time stepping loop

for p = 1:nt

Told = T;

% Node 1

T(1) = Told(1) + dt/(rho_fil*cp_fil*dx^2) * ( ...

k_fil*(Told(2)-Told(1)) + ...

k_gas*(Told(7)-Told(1)) + ...

2*k_fil*(T_base - Told(1)) + ...

q_dot*dx^2 );

% Node 2

T(2) = Told(2) + dt/(rho_fil*cp_fil*dx^2) * ( ...

k_fil*(Told(1)-Told(2)) + ...

k_fil*(Told(3)-Told(2)) + ...

k_gas*(Told(8)-Told(2)) + ...

q_dot*dx^2 );

% Node 3

T(3) = Told(3) + dt/(rho_fil*cp_fil*dx^2) * ( ...

k_fil*(Told(2)-Told(3)) + ...

k_gas*(Told(4)-Told(3)) + ...

2*k_gas*(Told(9)-Told(3)) + ...

q_dot*dx^2 );

% Node 4

T(4) = Told(4) + dt/(rho_gas*cp_gas*dx^2) * ( ...

k_fil*(Told(3)-Told(4)) + ...

k_glass*(Told(5)-Told(4)) + ...

k_gas*(Told(10)-Told(4)) );

% Node 5

T(5) = Told(5) + dt/(rho_glass*cp_glass*dx^2) * ( ...

k_gas*(Told(4)-Told(5)) + ...

k_glass*(Told(6)-Told(5)) + ...

2*k_glass*(Told(11)-Told(5)) );

% Node 6

T(6) = Told(6) + dt/(rho_glass*cp_glass*dx^2) * ( ...

k_glass*(Told(5)-Told(6)) - ...

2*h_t*dx*(Told(6)-T_inf) );

% Node 7

T(7) = Told(7) + dt/(rho_gas*cp_gas*dx^2) * ( ...

k_fil*(Told(1)-Told(7)) + ...

k_gas*(Told(8)-Told(7)) + ...

k_glass*(Told(12)-Told(7)) + ...

2*k_gas*(T_base - Told(7)) );

% Node 8

T(8) = Told(8) + dt/(rho_gas*cp_gas*dx^2) * ( ...

k_fil*(Told(2)-Told(8)) + ...

k_gas*(Told(7)-Told(8)) + ...

k_gas*(Told(9)-Told(8)) + ...

k_glass*(Told(13)-Told(8)) );

% Node 9

T(9) = Told(9) + dt/(rho_gas*cp_gas*dx^2) * ( ...

k_fil*(Told(3)-Told(9)) + ...

k_gas*(Told(8)-Told(9)) + ...

k_gas*(Told(10)-Told(9)) + ...

k_glass*(Told(14)-Told(9)) );

% Node 10

T(10) = Told(10) + dt/(rho_gas*cp_gas*dx^2) * ( ...

k_gas*(Told(4)-Told(10)) + ...

k_gas*(Told(9)-Told(10)) + ...

k_glass*(Told(11)-Told(10)) + ...

k_glass*(Told(15)-Told(10)) );

% Node 11

T(11) = Told(11) + dt/(rho_glass*cp_glass*dx^2) * ( ...

2*k_glass*(Told(5)-Told(11)) + ...

k_gas*(Told(10)-Told(11)) - ...

2*h_t*dy*(Told(11)-T_inf) );

% Node 12

T(12) = Told(12) + dt/(rho_glass*cp_glass*dx^2) * ( ...

k_gas*(Told(7)-Told(12)) + ...

k_glass*(Told(13)-Told(12)) - ...

h_t*dx*(Told(12)-T_inf) + ...

2*k_glass*(T_base - Told(12)) );

% Node 13

T(13) = Told(13) + dt/(rho_glass*cp_glass*dx^2) * ( ...

k_gas*(Told(8)-Told(13)) + ...

k_glass*(Told(12)-Told(13)) + ...

k_glass*(Told(14)-Told(13)) - ...

h_t*dx*(Told(13)-T_inf) );

% Node 14

T(14) = Told(14) + dt/(rho_glass*cp_glass*dx^2) * ( ...

k_gas*(Told(9)-Told(14)) + ...

k_glass*(Told(13)-Told(14)) + ...

k_glass*(Told(15)-Told(14)) - ...

h_t*dx*(Told(14)-T_inf) );

% Node 15

T(15) = Told(15) + dt/(rho_glass*cp_glass*dx^2) * ( ...

k_gas*(Told(10)-Told(15)) + ...

k_glass*(Told(14)-Told(15)) - ...

2*h_t*dy*(Told(15)-T_inf) );

T_history(:,p+1) = T;

end

%% 7. Final nodal temperatures

disp('Final transient nodal temperatures:');

for n = 1:15

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

end

%% 8. Transient line plot

figure(1);

time = (0:nt)*dt;

plot(time,T_history(1,:), ...

time,T_history(3,:), ...

time,T_history(6,:), ...

time,T_history(13,:), ...

time,T_history(15,:), ...

'LineWidth',1.5);

grid on;

xlabel('Time (s)');

ylabel('Temperature (K)');

title('Transient Temperature Response');

legend('Node 1','Node 3','Node 6','Node 13','Node 15','Location','best');

%% 9. Final temperature contour

figure(2);

Tfinal = T_history(:,end);

Tplot = zeros(3,6);

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

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

Tplot(1,1:6) = [Tfinal(1) Tfinal(2) Tfinal(3) Tfinal(4) Tfinal(5) Tfinal(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');

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

colorbar; colormap(jet);

title('Final Transient Temperature Distribution');

xlabel('Node Position (x)');

ylabel('Node Position (y)');

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

axis equal; grid on;

hold on; contour(xq,yq,Tinterp,15,'k'); hold off;

%% 10. Transient movie

videoName = 'LightBulb_Transient_Explicit_Final.mp4';

v = VideoWriter(videoName,'MPEG-4');

v.FrameRate = 30;

open(v);

figure(3);

frameCount = 0;

t_movie = 6;

p_max = round(t_movie/dt);

p_max = min(p_max, nt+1);

skip = max(1, round(nt/200));

for p = 1:skip:p_max;

Tnow = T_history(:,p);

Tplot = zeros(3,6);

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

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

Tplot(1,1:6) = [Tnow(1) Tnow(2) Tnow(3) Tnow(4) Tnow(5) Tnow(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');

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

colorbar; colormap(jet);

title(sprintf('Transient Temperature Distribution, t = %.4f s',(p-1)*dt));

xlabel('Node Position (x)');

ylabel('Node Position (y)');

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

axis equal; grid on;

hold on; contour(xq,yq,Tinterp,15,'k'); hold off;

drawnow;

frame = getframe(gcf);

writeVideo(v,frame);

frameCount = frameCount + 1;

end

close(v);

fullpath = fullfile(pwd, videoName);

fprintf('Movie saved as: %s\n', fullpath);

fprintf('Movie saved as: %s\n', videoName);