axial temperature profile in a Burner in 1D

Good Morning !

I'm trying to monitor aproximately the temperature profile inside the burner.. ( i'm not considering diffusion in radium coordinates i want only to monitor the axial profial of temperature, so the problem is semplyfied in 1D )

For determining how the Temperature profile is varying i have made both mass and Energy balances ( i'm attaching a file where all this balances are derived )..... the file name is mass and energy balances

After determining the balances i've tryed to discretize the differentials terms of equation like dT/dz with a first order explicit finite difference method . for example dT/dz wil become (T(i+1)-T(i))/h

the problem is that after 2 iteration the solution seems to remain constant and gives me back the same value of temperature profile and i think that is probably because the term h/v*(k0*exp(-Ea/(T(i)))*CH4(i)*O2(i)^2) has an order of 10^-17.. wich is probably negligible but it shouldn't...

Does anyone has an idea why this problem occurs even if i have taken all the values of cp,Ea,deltaHra,k0 from the ansys database

The original task is :

simulate a burner with methane and air at the inlet with no premixing.

The burner has a diameter of 0.225 meters and a length of 1,8 meters.

Methane is fed at the nozzle wich has the following dimensions : diameter of the nozzle = 0,005 m with a length of 0,01 m. The nozzle is placed at the center of the burner. methane and air are fed at 300 K . volumetric flowrate of methane at the inlet is 6*10^-3 [m^3/s] with an excess of air respect to stochiometric condition of 50 % . keep the cp constat for the smulation. With a matlab code : Try to evaluate the temperature profile considering only the axial coordinate 1D.

Thanks in advance for helping!

%%
clear 
close all
clc
D = 0.225; % diameter of the burner 
R1 = D/2; % radium of the burner 
L = 1.800; % length of the burner [m] 
l = 0.010; % length of the nozzle
d = 0.005; % diameter of the nozzel from here i feed the methane 
h = d/4; % discretizetaion step 
N = L/h+1; % index of row = i 
T(N) = 0; % vector for the Temperature 
CH4(N) = 0; % vector for concentration of CH4 
O2(N) = 0;% vector for  concentration of O2
nCH4 = 0.2437; % [mol/s] molar flowrate of CH4 
nO2 = 3*nCH4; % [mol/s] molar flowrate of O2
T(1)= 300;  % inlet temperature for air-methane mixture 
portataV = 0.0857; %[m^3/s] volumetric flowrate of air
portata02_volumetrica = portataV*0.21; % [m^3/s ] volumetric flowrate of oxygen
portataVCH4 = 0.006; % [m^3/s] volumetric flowrate of CH4 
O2(1) = nO2/portata02_volumetrica; % initial concentration of O2
CH4(1) = nCH4/portataVCH4; % initial concentration of CH4 
portataCH4 = 0.003942; % massic flowrate of CH4 [kg/s]
portataAria = 1.225*portataV; % [kg/s] massic flowrate of air 
m_tot = portataAria+portataCH4; % total massic flowrate air+ methane [kg/s]
V_reattore =  L * pi*D^2/4; % volume of the burner  [m^3]
cp = 1030;% J/(kg*K)
%k0 = 2.119*10^11; % ko di anys 
k0 = 1.5*10^13;  % pre-exponental factor of  Arrhenius equation 
R = 8.3145; % gas constant 
Ea = 2.027 * 10^4; % J/mol di ansys 
EadivR = 20000 ; % activation energy/gas constant [K] 
v = 200; % [m/s] initial velocity of methane inside the nozzle of the burner 
deltaHr =  890.8*10^3; % [J/N*m^3*mol] 
% for i = 1:N-1
%     T(i+1) = T(i)+ h/(m_tot*cp)*deltaHr*k0*exp(-Ea/(R*T(i)))*CH4(i)^0.2*O2(i)^1.3*pi*R^2;
%     CH4(i+1) = CH4(i) - h*v/(k0*exp(-Ea/(R*T(i)))*CH4(i)^0.2*O2(i)^1.3);
%     O2(i+1) = O2(i) - h*v/(k0*exp(-Ea/(R*T(i)))*CH4(i)^0.2*O2(i)^1.3);
% end
for i = 1:N-1
    T(i+1) = T(i)+ h/(m_tot*cp)*deltaHr*k0*exp(-EadivR/(T(i)))*CH4(i)*O2(i)^2*pi*R1^2;
    CH4(i+1) = CH4(i) - h/v*(k0*exp(-EadivR/(T(i)))*CH4(i)*O2(i)^2);
    O2(i+1) = O2(i) - h/v*(k0*exp(-EadivR/(T(i)))*CH4(i)^1*O2(i)^2);
end

Answer 1

[vidyesh]

Hi Amin Bouzaiene,

I understand that you are unable to calculate the temperature after 2 iterations because the difference between them is very small.

By default, MATLAB uses 16 digits of precision. For higher precision, use the ‘vpa’ function in Symbolic Math Toolbox. ‘vpa’ provides variable precision which can be increased without limit.

For your code you can make the following changes.

T(N) = vpa(0,20); %increase precision to 20 digits 
%Inside the for loop
T(i+1) = T(i)+ vpa(h/(m_tot*cp)*deltaHr*k0*exp(-EadivR/(T(i)))*CH4(i)*O2(i)^2*pi*R1^2); 

You can find more information on ‘vpa’ function here:

https://www.mathworks.com/help/symbolic/vpa.html

Please note that this may increase the run time of your program, an alternative would be to focus your analysis on the difference between the temperatures rather than their absolute values.

Hope this helps.