% Models cane diameter evolution in a glass cane drawing process
% JJ Hu, 09/21/2015

clear all;
sigma = 5.67E-8;    % Stefan-Boltzmann constant
emissivity = 0.8;   % Emissivity of glass
c_P = 2.058;    % Heat capacity in J/cm^3/K
D0_1 = 1;   % Starting diameter of segment 1 in cm
D0_2 = 0.9;   % Starting diameter of segment 2 in cm
eta0 = 4E-9;    % Pre-exponential factor for glass viscosity
Delta_E_a = 2.9;   % in eV
k_B = 8.62E-5;  % in eV/K
t = 0;      % Starting time in s
T0 = 1173;  % Starting temperature in K
ii = 1;
dt = 0.1;   % Time step in s
t_end = 10; % End time for simulation
T_1(1) = T0;
T_2(1) = T0;
D_1(1) = D0_1;
D_2(1) = D0_2;

for t = 0:dt:t_end;
    eta_1 = eta0 * exp(Delta_E_a/k_B/T_1(ii));  % Viscosity of segment 1
    eta_2 = eta0 * exp(Delta_E_a/k_B/T_2(ii));  % Viscosity of segment 2
    dL_dt_1 = 0.025 * D_2 * eta_2 / [(D_1 * eta_1) + (D_2 * eta_2)];
    dL_dt_2 = 0.025 * D_1 * eta_1 / [(D_1 * eta_1) + (D_2 * eta_2)];
    % Assume a total length change rate of 2.5% per dt time step distributed
    % between two segments
    dD_dt_1 = 1 - sqrt(1 / (1 + dL_dt_1));  % Given in terms of fractional diameter change
    dD_dt_2 = 1 - sqrt(1 / (1 + dL_dt_2));
    D_1(ii+1) = D_1(ii) * (1 - dD_dt_1);    % Diameter of segment 1 as a function of time
    D_2(ii+1) = D_2(ii) * (1 - dD_dt_2);    % Diameter of segment 2 as a function of time
    dT_dt_1 = 4 * emissivity * sigma * T_1(ii)^4 / (c_P*1E6) / (D_1(ii)/100);
    dT_dt_2 = 4 * emissivity * sigma * T_1(ii)^4 / (c_P*1E6) / (D_2(ii)/100);
    T_1(ii+1) = T_1(ii) - dT_dt_1 * dt;     % Temperature of segment 1 as a function of time
    T_2(ii+1) = T_2(ii) - dT_dt_2 * dt;     % Temperature of segment 2 as a function of time
    ii = ii + 1;
end

% Data plot
set(gca, 'box', 'on')
hold on
plot(0:dt:t_end+dt, D_1, 'b')
plot(0:dt:t_end+dt, D_2, 'r')
xlabel('Time (s)', 'Fontsize', 12)
ylabel('Cane diameter (cm)', 'Fontsize', 12)
axis([0 10 0.5 1])
figure
set(gca, 'box', 'on')
hold on
plot(0:dt:t_end+dt, T_1-273, 'b')
plot(0:dt:t_end+dt, T_2-273, 'r')
xlabel('Time (s)', 'Fontsize', 12)
ylabel('Temperature (C)', 'Fontsize', 12)
axis([0 10 700 900])