%This program was developed to calculate the ignition temperature during
%heating of n powder layers on a wire filament coated with Mg powder.
%Written by Trent Ward
%Last revision April 27, 2005
%Simplified version for thesis attachment
clear all; % clear all variables in memory
close all; %close all open windows/figures
format compact; %set display format
format short g; %set display format
tic; %start timer
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%% Define Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Calculation control parameters %%%%%%%%%%%%%%%%%%%%
dtemp = 0.01; % K, Temperature change per iteration step
dt_ceil = 1e-4; % s, maximum allowed time per iteration step
dt_floor = 1e-9; % s, minimum allowed time per iteration step
tol = 2*dtemp; % K, Temp tolerance in which melting is allowed to occur
%%%%%%%%%%%%%%%%%%%%%%%%%%Known constants%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global g
g=9.81; %gravitational constant, m/s^2
R=8.31447E-3; %Universal gas constant, kJ/mol K
sigma=5.6704E-8; %Stefan-Boltzmann constant, W/m^2 K^4
%%%%%%%%%%%%%%%%%%%%%%%%%%Particle Properties%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global n D Dwire
D=9.71E-6; %particle diameter, m
filament_diameter=492E-6; %filament diameter
layer_thickness=(603E-6-filament_diameter)/2; %thickness of the coating calculated from the measured coating diameter of 603 micron, m
n=round(layer_thickness/D); %number of layers
sa=pi()*(D^2); %surface area of particle, m^2
saa=pi()*(D^2-12/4*(D/50)^2); %surface area used to estimate the actual area exposed to radiation per particle, m^2
V=pi()*D^3/6; %volume of particle, m^3
thermal_diffusivity=2.296E-7; %thermal diffusivity experimental, m^2/s
Mg_density=1740; %Mg_density of powder, kg/m^3
mass = V*Mg_density; %mass of single particle, kg
bulk_density=1259.7; %measured bulk density from laser flash diffusivity sample, kg/m^3
density_ratio=bulk_density/(.74*Mg_density); %ratio between the experimental density and the theoretical density for the HCP packing
contacts=3*density_ratio; %number of contacts between layers
Mg_heat_capacity=1024; %heat capacity, J/kg K
Tmelting=923; %melting point of Mg, K
Hf=8.7/24.305*1000^2; %Mg heat of fusion,kJ/mol--->J/kg
contact_resistance=(bulk_density*(6)^.5)/(thermal_diffusivity*(Mg_density*.74)^2*Mg_heat_capacity*D); %calculation of the individual contact resistance between particles, K/W
Z1=10E9; %Pre-exponent solid,kg/m^2 s
E1=215; %activation energy solid, kJ/mol
dH1=24.7E6; %specific heat of oxidation reaction solid,kJ/mol--->J/kg
Z2=Z1; %Pre-exponent liquid,kg/m^2 s
E2=E1; %activation energy liquid, kJ/mol
dH2=25.35E6; %specific heat of oxidation reaction liquid,kJ/mol--->J/kg
cur_HR(1:n) = 1; %specify vector for current heating rate
filament_contacts=1; %number of contacts between first layer and the filament
%%%%%%%%%%%%%%%%%%%%%%%%%%Radiation Properties%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
particle_emissivity=.75; %emissivity of particle
S=.5; %view factor between particles
%%%%%%%%%%%%%%%%%%%%%%%%%%Filament Properties%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filament_length=.0467; %length of filament, m
filament_cs_area=0.25*pi()*filament_diameter^2; %Cross section area of filament, m^2
filament_thermal_conductivity=12; %conductivity of filament, W/m K
filament_density=8400; %density of filament, kg/m^3
filament_resistivity=112E-8; %resistivity of filament, ohm m
filament_emissivity=.75; %emissivity of filament
m=75; %numer of nodes on filament
mm=round(m/2); %approximate center of filament used for filament temperature
dx=filament_length/m; %nodal spacing for filament mesh, m
coating_length=.0095; %powder coating length, m
powder_nodes=round(coating_length/dx); %number of nodes coated with powder
pyrometer_location=round(mm+(.0095/dx)/2+.002/dx); %location of the pyrometer based on experimental locations of 2 mm from edge of coating
%%%%%%%%%%%%%%%%%%%%%%%Layer Diameter Calculation%%%%%%%%%%%%%%%%%%%%%%%%%%
global Dlayer
Nlayer=1:n; %vector of layer numbers
Dlayer=filament_diameter+D.*(1+(8/3).^.5.*(Nlayer-1)); %diameter of each layer, m
%%%%%%%%%%%%%%%%%%%%%%%%%%% Air Properties%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global Tsurrounding;
global kair_p vair_p Prair_p mu_T;
Tsurrounding=298; %surrounding temperature, K
ignitiontemp = []; % declare output array so we can later append to it
temp=[20,50,100,200,300,400,500,600,700,800,900,1000,1500,2000]+273; %temperature for given properties below
kair=[.02514,.02735,.03095,.03779,.04418,.05015,.05572,.06093,.06581,.07037,.07465,.07868,.09599,.11113]; %thermal conductivity of air, W/m K
[kair_p,S_,mu_T]=polyfit(temp,kair,3); %polynomial fit to the tabulated data for the air thermal conductivity
vair=[1.516,1.798,2.306,3.455,4.765,6.219,7.806,9.515,11.33,13.26,15.29,17.41,29.22,42.7]*1E-5; %kinematic viscosity of air, m/s^2
[vair_p,S_,mu_T]=polyfit(temp,vair,3); %polynomial fit to the tabulated data for the air kinematic viscosity
Prair=[.7309,.7228,.7111,.6974,.6935,.6948,.6986,.7037,.7092,.7149,.7206,.7260,.7478,.7539]; %Prandtl number of air
[Prair_p,S_,mu_T]=polyfit(temp,Prair,7); %polynomial fit to the tabulated data for the air Prandtl number
%%%%%%%%%%%%%%%%%%%%%%%%%%%Filament Nusselt Calculation%%%%%%%%%%%%%%%%%%%%
Tnuss=298:1:2000; %temperature vector for the nusselt number calculations below
Tfilmnuss=(Tnuss+Tsurrounding)/2; %temperature vector for the film temperature
Ra_filament=g.*(1./Tfilmnuss).*(Tnuss-Tsurrounding).*filament_diameter^3.*interp1(temp,Prair,Tfilmnuss)./(interp1(temp,vair,Tfilmnuss)).^2; %Rayleigh number for the natural convection around the filament
Nu_filament=(0.6+(.387.*(Ra_filament).^(1/6))./(1+(.559./interp1(temp,Prair,Tfilmnuss)).^(9/16)).^(8/27)).^2; %Nusselt number for filament convection
[Nu_filament_p,S_,mu_Tn]=polyfit(Tnuss,Nu_filament,23); %polynomial fit to the filament Nusselt data
%%%%%%%%%%%%%%%%%%%%%%%Particle Nusselt Calculation%%%%%%%%%%%%%%%%%%%%%%%
Ra_particle=g.*(1./Tfilmnuss).*(Tnuss-Tsurrounding).*(Dlayer(n)+D).^3.*interp1(temp,Prair,Tfilmnuss)./(interp1(temp,vair,Tfilmnuss)).^2; %Rayleigh number for the natural convection around the coating
Nu_particle=(0.6+(.387.*(Ra_particle).^(1/6))./(1+(.559./interp1(temp,Prair,Tfilmnuss)).^(9/16)).^(8/27)).^2; %Nusselt number for coating convection
[Nu_particle_p,S_,mu_Tn]=polyfit(Tnuss,Nu_particle,23); %polynomial fit to the coating Nusselt data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Start of temperature calculation routine %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for HR=1:5 %heating rate index, 1 through 5 is for the five heating rates matched to the experimental, this can be changed to run only one heating rate or all at the same time
%the voltage and resistance can be varied to calculate
%the ignition temperature for any heating rate
switch HR
case 1 %lowest heating rate parameters to match with the experimental
voltage=12.3;
external_resistance=2.08;
case 2 %low heating rate parameters to match with the experimental
voltage=12.3;
external_resistance=1.27;
case 3 %middle heating rate parameters to match with the experimental
voltage=12.3;
external_resistance=.635;
case 4 %high heating rate parameters to match with the experimental
voltage=12.3;
external_resistance=.18;
case 5 %highest heating rate parameters to match with the experimental
voltage=24;
external_resistance=.187;
end
iteration=400000; %maximum number of particle iterations
dt = dt_ceil; %set initial time step to the ceiling value to start the computation
dtfil = 100 * dt; %set the initial filament time step to 100 times the particle time step
indexfil = round(iteration/100); %specifiy the size of the filament array
timefil(1:3)=0; %declare filament time vector
time(1:iteration)=0; %declare particle time vector
f=1; % time counter IC
g=1; %filament counter IC
T(1:iteration,n)=0; %Define particle temperature matrix
T(1,:)=Tsurrounding; %assign initial conditions at time=0
Tf(1:indexfil,1:m)=0; %define filament temperature matrix
Tf = Tf + Tsurrounding; %set temperature of filament to the surrounding temperature
current = zeros(indexfil,m-1); %declare an array for the filament electric current
Qcondwire = zeros(indexfil,m-1); %declare an array for the filament heat conduction
Qconvwire = zeros(indexfil,m-1); %declare an array for the filament heat convection
Qradwire = zeros(indexfil,m-1); %declare an array for the filament radiation
Gwire = zeros(indexfil,m-1); %declare an array for the filament heat generation
Qcond(1:n) = 0; %declare an vector for the particle conduction
Qconv(1:n) = 0; %declare an vector for the particle convection
Qrad(1:n) = 0; %declare an vector for the particle radiation
Qchem(1:n) = 0; %declare an vector for the particle chemical heat
Qphys(1:n) = 0; %declare vector to sum up "physical" heat flow terms
Qmelt(1:n) = 0; %declare an vector for the heat accumulated during melting
Qexcess(1:n) = 0; % this is how far we've overshot the heat req'd to melt a particle
cumulative_Qchem(1:n) = 0; % cumulative term: integrated heat production due to reaction
consumed(1:n) = 0; % cumulative term: fraction of particle that has been oxidized
sustained(1:n) = 0; % current ratio of "chemical" to "physical" heat flow terms
step_1 = false; % flag used to write data to 'ignitiontemp' array when Qc/Qp == 1
step_10 = false; % flag used to write data to 'ignitiontemp' array when Qc/Qp == 10
temp_OK = false; % flag used to indicate that a temperature jump is imminent
stop = false; % flag used to indicate that ignition has been detected
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%Main Solving Block%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~stop; % time loop; stop will be true if ignition is detected.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Filament calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate filament temperature when particle /time/
% approaches that of most recent filament T calculation
if timefil(g)-time(f) < 2 * dt || g==1 %statement to enter filament equations when particle time is near the last filament time
dtfil=100*dt; % aim to update filament every 100 particle iterations
Tf(g+1,1:m) = Tsurrounding; % add one row to Tf array, initialize with ambient temperature
for z=2:m-1 %nodal loop for filament
if z=round(mm+powder_nodes/2-1) %calculate for all nodes not in contact with powder
current(g,z)=voltage/(external_resistance+resistance_coefficient(Tf(g,z))*filament_resistivity*filament_length/filament_cs_area); %calculation of current through circuit, this is tempeature dependent
Qcondwire(g,z)=(filament_thermal_conductivity*filament_cs_area/dx)*(Tf(g,z-1)+Tf(g,z+1)-2*Tf(g,z)); %conduction heat in filament
T_ = (Tsurrounding+Tf(g,z))/2;
Qconvwire(g,z)=polyval(Nu_filament_p,T_,[],mu_Tn)*polyval(kair_p,T_,[],mu_T)*dx*pi()*(Tsurrounding-Tf(g,z)); %convection heat lost from filament
Qradwire(g,z)=filament_emissivity*sigma*dx*pi()*filament_diameter*(Tsurrounding^4-Tf(g,z)^4); %radiation to surroundings
Gwire(g,z)=current(g,z)^2*resistance_coefficient(Tf(g,z))*(filament_resistivity*dx/filament_cs_area); %heat generated in filament
Tf(g+1,z)=Tf(g,z)+(dtfil/(filament_density*filament_cs_area*dx*Cp_filament(Tf(g,z))))*(Qcondwire(g,z)+Qconvwire(g,z)+Qradwire(g,z)+Gwire(g,z)); %Calculate nodal temperature for filament at next time step
else %calculate filament temperature for nodes in contact with powder
current(g,z)=voltage/(external_resistance+resistance_coefficient(Tf(g,z))*filament_resistivity*filament_length/filament_cs_area); %calculation of current through circuit, this is tempeature dependent
Qcondwire(g,z)=(filament_thermal_conductivity*filament_cs_area/dx)*(Tf(g,z-1)+Tf(g,z+1)-2*Tf(g,z))+pi()*Dlayer(1)*dx/D^2*(T(f,1)-Tf(g,z))/contact_resistance*filament_contacts; %conduction heat in filament
Qradwire(g,z)=filament_emissivity*sigma*dx*pi()*filament_diameter*(T(f,1)^4-Tf(g,z)^4); %radiation to particles
Gwire(g,z)=current(g,z)^2*resistance_coefficient(Tf(g,z))*(filament_resistivity*dx/filament_cs_area); %heat generated in filament under the coating
Tf(g+1,z)=Tf(g,z)+(dtfil/(filament_density*filament_cs_area*dx*Cp_filament(Tf(g,z))))*(Qcondwire(g,z)+Qradwire(g,z)+Gwire(g,z)); %Calculate nodal temperature for filament nodes coated with powder at next time step
end % if
end % for z loop
timefil(g+1)=timefil(g)+dtfil; % update time for filament calc
g=g+1; % increase counter
end % if
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Particle temperature calculations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:n % loop over particle layers in coating, n layers
% Actual calculation of temperature and heat flow follows
% Go ahead if T either hasn't been calculated or if current
% calculation would result in a T jump and T needs to be recalculated
while ~temp_OK;
% [f, j, dt]
% Here we calculate the temperature of particle /j/ at
% iteration /f+1/ based on the temperature of its
% neighbors at iteration /f/
%
% We won't increase /f/ unless we are satisfied with
% the calculation
if j==1 %calculation for particle contacting the filament
if f==1 %for the first iteration, the filament temperature is set to ambient, avoids negative indexing in linear interpolation
Tf_=298;
else
Tf_=(time(f)-timefil(g))/(timefil(g-1)-timefil(g))*(Tf(g-1,mm)-Tf(g,mm))+Tf(g,mm); % linear interpolation between the calculated filament temperature points
end
Qcond(j) = 1/contact_resistance * (filament_contacts*(Tf_-T(f,j))+contacts*(T(f,j+1)-T(f,j))); %conduction heat
Qrad(j) = particle_emissivity*sigma*saa*S*(Tf_^4+T(f,j+1)^4-2*T(f,j)^4); %radiation heat
if T(f,j) <= Tmelting
Qchem(j) = sa*dH1*Z1*exp(-E1/(R*T(f,j))); %heat from chemical reaction before melting
elseif T(f,j) > Tmelting
Qchem(j) = sa*dH2*Z2*exp(-E2/(R*T(f,j))); %heat from chemical reaction after melting
end
if T(f,j)Tmelting %temperature calculation for next time step before and after melting
T(f+1,j)=T(f,j)+(dt/(Cp_part(T(f,j))*mass))*(Qcond(j)+Qrad(j)+Qchem(j)); %Temperature at next time step
elseif T(f,j)>=Tmelting-tol %start of melting given some range to know that melting is about to start
T(f+1,j)=Tmelting; %constant temperature condition during melting
Qmelt(j)=Qmelt(j) + dt * (Qcond(j)+Qrad(j)+Qchem(j)); %cumulative heat taken to melt particle for melting period
if Qmelt(j) > Hf*mass
Qexcess(j)=Qmelt(j) - Hf*mass; %calculate heat accumulated after the particle has completely melted
T(f+1,j)=T(f,j)+(dt/(Cp_part(T(f,j))*mass))*(Qexcess(j)); %temperature at next time step considering melting
end
end % if
elseif (j>1) & (j Tmelting
Qchem(j) = sa*dH2*Z2*exp(-E2/(R*T(f,j))); %heat from chemical reaction after melting
end
if T(f,j)Tmelting %temperature calculation for next time step before and after melting
T(f+1,j)=T(f,j)+(dt/(Cp_part(T(f,j))*mass))*(Qcond(j)+Qrad(j)+Qchem(j)); %Temperature at next time step
elseif T(f,j)>=Tmelting-tol %start of melting given some range to know that melting is about to start
T(f+1,j)=Tmelting; %constant temperature condition during melting
Qmelt(j)=Qmelt(j) + dt * (Qcond(j)+Qrad(j)+Qchem(j)); %cumulative heat taken to melt particle for melting period
if Qmelt(j) > Hf*mass
Qexcess(j)=Qmelt(j)-Hf*mass; %calculate heat accumulated after the particle has completely melted
T(f+1,j)=T(f,j)+(dt/(Cp_part(T(f,j))*mass))*(Qexcess(j)); %temperature at next time step considering melting
end
end
elseif j==n %calculation for outer particle in contact with boundary layer
Qcond(j) = contacts*1/contact_resistance*(Dlayer(j-1)/Dlayer(j))*(T(f,j-1)-T(f,j)); %conduction heat
Qrad(j) = particle_emissivity*sigma*saa*S*(T(f,j-1)^4+Tsurrounding^4-2*T(f,j)^4); %radiation heat
Qconv(j) = (Mg_density*.74)/bulk_density*D^2/(Dlayer(n))*polyval(Nu_particle_p,(T(f,j)+Tsurrounding)/2,[],mu_Tn)*polyval(kair_p,(T(f,j)+Tsurrounding)/2,[],mu_T)*(Tsurrounding-T(f,j)); %convection heat lost from the last particle layer to the surrounding
if T(f,j) <= Tmelting
Qchem(j) = sa*dH1*Z1*exp(-E1/(R*T(f,j))); %heat from chemical reaction before melting
elseif T(f,j) > Tmelting
Qchem(j) = sa*dH2*Z2*exp(-E2/(R*T(f,j))); %heat from chemical reaction after melting
end
if T(f,j)Tmelting %temperature calculation for next time step before and after melting
T(f+1,j)=T(f,j)+(dt/(Cp_part(T(f,j))*mass))*(Qcond(j)+Qrad(j)+Qchem(j)+Qconv(j)); %Temperature at next time step
elseif T(f,j)>=Tmelting-tol %start of melting given some range to know that melting is about to start
T(f+1,j)=Tmelting; %constant temperature condition during melting
Qmelt(j)=Qmelt(j) + dt * (Qcond(j)+Qrad(j)+Qchem(j) + Qconv(j)); %cumulative heat taken to melt particle for melting period
if Qmelt(j) > Hf*mass
Qexcess(j)=Qmelt(j)-Hf*mass; %calculate heat accumulated after the particle has completely melted
T(f+1,j)=T(f,j)+(dt/(Cp_part(T(f,j))*mass))*(Qexcess(j)); %temperature at next time step considering melting
end
end % if
end %
% now we're ready to move on if the newly calculated
% temperature is not more than 2* the targeted
% temperature interval - or - if we've already hit the
% shortest allowed time step.
if T(f+1,j)-T(f,j) < 2 * dtemp || dt/2 <= dt_floor
temp_OK = true; % this really needs a repeat..until construct
continue
end % if
dt = dt / 2; % repeat temperature calculation with decreased time step
end % while ~temp_OK
temp_OK = false; % reset the flag so that we calculate T in next iteration
Qphys(j) = (Qcond(j) + Qrad(j) + Qconv(j)); % sum up "physical" heat flow terms for convenience
cumulative_Qchem(j) = cumulative_Qchem(j) + Qchem(j)*dt; %sum up the chemical heat
consumed(j) = cumulative_Qchem(j) / (dH1 * mass); %calculate the fraction of the particle that has been oxidized
if Qphys(j) ~= 0; sustained(j) = Qchem(j) / Qphys(j); else; sustained(j) = 0; end; %write the ratio of chemical heat to physical heat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Ignition Criteria %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (1363)<=T(f,j) %satisfied if the particle temperature is above the boiling temperature of Mg
disp('The ignition temperature is:')
disp([interp1(timefil(1:g),Tf(1:g,pyrometer_location),time(f)),f,time(f)]) %display ignition temperature
for HR_fit_time=1:2:1000 %find time range for a 20K temperature interval
if (Tf(g,pyrometer_location)-Tf(g-HR_fit_time,pyrometer_location))>=20
break
end % if
end % for
real_HR = (Tf(g,pyrometer_location)-Tf(g-HR_fit_time,pyrometer_location))/(timefil(g)-timefil(g-HR_fit_time)); %calculate the heating rate near ignition
ignitiontemp=[ignitiontemp; HR, interp1(timefil(1:g),Tf(1:g,pyrometer_location),time(f)), time(f), j, real_HR, T(f,j), consumed(j), sustained(j),interp1(timefil(1:g),Tf(1:g,mm),time(f))]; %store data for this heating rate in an array
stop = true; %end loop for this heating rate
break %break from for loop
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Calculate next time step based on temperature history %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if f < 10
cur_HR(j) = 1; % low heating rate to force dt_floor
else
cur_HR(j) = (T(f,j)-T(f-4,j))/(time(f)-time(f-4)); % dT/dt determined from last 5 points [K/s]
end % if
end % for j=1:n
if stop %break from while loop
break
end
time(f+1)=time(f)+dt; %index time counter
% calculate next time step from max. current heating rate
dt = dtemp / (0.0001 + max(cur_HR(1:n))); % "+.0001" == cheap insurance against div/zero error
if dt > dt_ceil; dt = dt_ceil; end; % upper limit
if dt < dt_floor; dt = dt_floor; end; % lower limit
if (rem(f,5000)==0) %condition to output data every 5000 time steps
disp([dt,f,T(f,:)]); %display current parameters
end
f=f+1; %increment particle calculation to next time step
end %end while particle loop
end %for loop HR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Write and display data files %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dateandtime=fix(clock); %determine current time to specify run output
filename=['run_' sprintf('%04d%02d%02d%02d%02d', dateandtime(1:5)) '.csv']; %define the filename for the output
diary(filename);
disp(sprintf('%g ,pre-exponent',Z1));
disp(sprintf('%g ,activation energy',E1));
disp(sprintf('%g ,thermal thermal_diffusivity',thermal_diffusivity));
disp(sprintf('%d ,number of particle layers', n));
disp(sprintf('%g ,particle diameter', D));
disp(sprintf('%g ,measured Mg_density', bulk_density));
disp(sprintf('%g ,contacts between layers', contacts));
disp(sprintf('%g ,powder covering length', coating_length));
disp('setting,filament temperature,time,igniting layer,heating rate,particle temperature,fraction of particle consumed,Qchem/Qphys');
diary off;
dlmwrite(filename,[ignitiontemp],'-append'); %write the file
disp(ignitiontemp); %display results for all heating rates
toc %stop timer and output elapsed time