Filling in values of a formula into another formula.
4 views (last 30 days)
Show older comments
Sue Ellen de Nijs
on 3 Mar 2021
Edited: Alan Stevens
on 11 Mar 2021
Hello,
I'm working on a heat transfer in Matlab on a concrete floor during a fire. During the fire the gas temperature in the compartment, the density of the concrete and the specific heat of the concrete will change every second. To be able to calculate the surface temperature over time, the initial surface temperature is calculated at t=1. Than the formula can be filled in again with this value and so on.
But, when I'm filling the initial surface temperature into the formula for dOc, the gas temperature at t=1 is used again, while it must continue on t=2. I also want to create a loop, because the calculation will consists out of 5400 timesteps after the initial time step.
Does anyone know how I can do this? I also placed the Matlab script below.
Thank you!
%Exposed surface
L=0.24; %m, floor thickness
n=10; %number of simulations
T0=20; %deg c initial temperature floor
T2s=20; %deg c unexposed surface temperature
k_initial=1.4875; %thermal conductivity
p0_initial=2300; %initial density concrete
Cp_initial=900; %initial specific heat concrete
alfa=25; %W/m2degC convection transfer coeffiecient
dx=L/n; %, node thickness
%Terms
er=0.5; %emissivity of the boundary surface and combustoion gases
dt=1; %s fixed time step
t_final=10; %s, simulation time
x=dx/2:dx:L-dx/2;
T=ones(n,1)*T0;
dTdt=zeros(n,1);
t=dt:t_final;
%thermal diffusivity
a=k_initial/(p0_initial*Cp_initial);
%Gas temperature
Tt_celcius= 20+345*log10(8*t/60+1); %degC gas temperature in compartment
Tt_kelvin=Tt_celcius+273;
%Change in temperature in de the section over time interval dt
dOc_initial=alfa/(p0_initial*Cp_initial)*(Tt_kelvin-293)*dt; %(s)
%Surface temperature after dt
Ts_initial=dOc_initial+293;
%density
if Ts_initial <=338;
p0=20;
elseif Ts_initial <= 673;
p0=20*(0.98-0.03*((Ts_initial-200)/200));
elseif Ts_initial <= 1473;
p0=20*(0.95-0.07*((Ts_initial-400)/800));
end
%specific heat
if Ts_initial <=373;
Cp=900;
elseif Ts_initial <=573;
Cp=900+(Ts_initial-100);
elseif Ts_initial<=673;
Cp=1000+((Ts_initial-200)/2);
elseif Ts_initial<=1473;
Cp=1100;
end
dOc=alfa/(p0*Cp)*(Tt_kelvin-Ts_initial)*dt;
Ts=dOc+293;
2 Comments
darova
on 4 Mar 2021
Looks like you are trying to use difference scheme (partial diff equations)
Do you have original formulas?
Accepted Answer
Alan Stevens
on 11 Mar 2021
Edited: Alan Stevens
on 11 Mar 2021
Your equation depends only on t, not x. The following is a modification of your code where I've removed the parameters that are unused, but it shows how to incorporate temperature dependence into your density and specific heat.
(You don't seem to have included the mass term).
n=10; %number of timesteps
T0=20; %deg c initial temperature floor
k_initial=1.4875; %thermal conductivity
p0_initial=2300; %initial density concrete
Cp_initial=900; %initial specific heat concrete
alfa=25; %W/m2degC convection transfer coeffiecient
%Terms
dt=1; %s fixed time step
t_final=10; %s, simulation time
t=dt:t_final;
%thermal diffusivity
a=k_initial/(p0_initial*Cp_initial);
% Make Gas temps functions of t
%Gas temperature
Tt_celcius=@(t) 20+345*log10(8*t/60+1); %degC gas temperature in compartment
Tt_kelvin=@(t) Tt_celcius(t)+273;
%Change in temperature in de the section over time interval dt
dOc_initial=alfa/(p0_initial*Cp_initial)*(Tt_kelvin(t(1))-293)*dt; %(s)
%Surface temperature after dt
Ts_initial=dOc_initial+293;
% Initialize Ts
Ts = zeros(n,1);
Ts(1) = Ts_initial;
% Loop through timesteps
for i = 1:n-1
p0 = density(Ts(i));
Cp = spheat(Ts(i));
dOc=alfa/(p0*Cp)*(Tt_kelvin(t(i))-Ts(i))*dt;
Ts(i+1)=Ts(i) + dOc; % Update temperature
end
plot(t,Ts),grid
xlabel('time'),ylabel('Ts')
% Make density and specific heat functions of Ts
%density
function p0 = density(Ts)
if Ts <=338
p0=20;
elseif Ts <= 673
p0=20*(0.98-0.03*((Ts-200)/200));
elseif Ts <= 1473
p0=20*(0.95-0.07*((Ts-400)/800));
end
end
%specific heat
function Cp = spheat(Ts)
if Ts <=373
Cp=900;
elseif Ts <=573
Cp=900+(Ts-100);
elseif Tsl<=673
Cp=1000+((Ts-200)/2);
elseif Ts<=1473
Cp=1100;
end
end
0 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!