Filling in values of a formula into another formula.

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;

[darova]

Looks like you are trying to use difference scheme (partial diff equations)

Do you have original formulas?

[Sue Ellen de Nijs]

I'm sorry for the late response. The goal is to fill in the formula below for timesteps of 1 second, in the Matlab file called %change in temperature in the section over time interval dt. For calculating the following second, the properties of the concrete due to the temperature of the previous timestep should be used. The properties of the concrete like the density and the specific heat can be calculated with the formulas above. Those are the original formulas from the Eurocode. The problem is that I do not know how I can combine the formulas in such a way that the right numbers are used.

Answer 1

[Alan Stevens]

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