# Development of numerical model for analysis of biomass based pottery furnace

Jha Rahul Binod

Third Year Undergraduate, Department of Mechanical Engineering, National Institute of Technology Sikkim, Ravangla 737139, Sikkim

Prof. M. R. Ravi

Department of Mechanical Engineering, Indian Institute of Technology Delhi, Hauz Khas, New Delhi 110016

## Abstract

Furnaces are in use since a long time for different applications like bell metal craft, pottery making, brick making, bangle manufacturing, etc. Generally, in today’s world, industries use forced draft furnaces because of the availability of resources and electricity. These forced draft furnaces have good efficiency since the air flow rate and thus the power rating of the furnace can be controlled during operation. But most of the small artisans in rural areas use simple low investment furnaces which do not require electricity and are based on natural draft. These traditional furnaces do not have good efficiency and have very high emissions, observing which scientific community started working on development of models to gain insights into the energy distribution of these systems so as to improve their efficiency. The performance of a natural draft furnace is dependent on its configuration. Irrespective of the applications, the governing equations involved in the analysis of any furnace would be conservation of mass, momentum and energy. Thus, students at IIT-Delhi developed a generic model, which can be used for simulating furnaces of different configurations to determine different parameters like energy utilization, temperature distribution, etc. As a part of Summer Research Fellowship project, I along with another summer fellow, Raj Joshi have developed a numerical model for simulating a biomass based pottery furnace installed at village Khurrampur in Haryana. For this we first separated the whole furnace in different control volume, then we used the discretised form of conservation equation of mass, momentum and energy for every control volume. For numerical computation of equations we have made a matlab code which gives the temperature and pressure of different control volume with time taking mass flow rate as an input.

Keywords: natural draft, firing chamber, pottery chamber, MATLAB code, simulation

## Abbreviations

 Abbreviations P Pressure (Pa) ρ Density ( $kg/m^3$) T Temperature (K) Q Heat flux $(W/m^2)$ v Velocity (m/s) g Acceleration due to gravity $(m/s^2)$ A Cross sectional area $(m^2)$ ​ $\epsilon$ Emissivity of air/gas ​ $\sigma$ Stephan-Boltzmann constant $(W/m^2/T^4)$​ R Gas constant (J/kg/K) k' Effective conductivity of middle part of wall (W/m/K) H Heat transfer coefficient $(W/m^2/K)$ Dp Effective Particle Diameter (m) μ Coefficient of dynamic viscosity(kg/m/s) D Diameter (m) dx Length(m) of each control volume in Pottery Chamber CV Calorific Value (J/kg) Cv Specific heat capacity at constant volume(J/kg/K) Cp Specific heat capacity at constant Pressure (J/kg/K) V Volume $(m^3)$ Φ Porosity m Mass floe rate (kg/s) k Thermal conductivity (W/m/K) t time (s) K Co-efficient of Pressure loss Subscripts in Inlet out Outlet g / gas Flue gas ref Reference F Fuel air Air amb / a Ambient i No of control volume Superscripts t At time t t+Δt At time t+Δt

## Background/Rationale

Pottery Industries are in use since ancient times. Small potters in developing nations like India still use traditional pottery kilns to make pots and different ceramic components. These pottery kilns uses both biomass or LPG for firing. In villages small potters do not use LPG due to its inavailability or higher cost. Also pottery furnances can be both natural draft or forced draft. Since the natural draft furnance systems do not require any electrical supply they are widely used in rural area for numerous application. Comprehensive understanding of functioning of forced draft furnace is available to the scientific community but such widespread insights about natural draft furnance system is not available. Owing to the widespread use and poor energy efficiency of natural draft furnace systems, scientific community started working on gaining insights into the energy distribution and performance of these systems; engineers started making designs to improve the fuel efficiency after analyzing the results of model’s simulations. Considerable amount of work have been done to gain insights about the natural draft furnance system in IIT Delhi. Biomass based pottery furnace have been installed by IIT Delhi at a village khurrampur in Haryana. As a part of the initiative taken by IIT Delhi, this project is about to develop a numerical model for the furnance.

Schematic of front view of pottery furnace

Based on the previous literary work done by different teams on natural draft furnace system, this model has been developed which, considering the pottery furnace as a combination of different control volume gives the temperature and pressure of these control volumes as output taking mass flow rate of air flow and fuel feed rate as an input.

## Objectives of the Research

The project intends to create a model that can estimate the energy distribution (i.e. temperature and pressure) in different control volume of the biomass based pottery furnace. The model estimates the temperature at different points taking firing rate and air fuel rate as an input.

To estimate the effect of providing no. of holes for air flow on the performance of the pottery furnance.

To validate the model by comparing the estimates predicted by the model with the experimental observations available due to the research held in IIT Delhi.

## LITERATURE REVIEW

Considerable efforts have been made by students at IIT Delhi to gain insight into the function of natural draft furnance. Oza et al. [1] developed a generic mathematical model of natural draft furnace systems with respect to fuel, application and geometry which estimates energy distribution and parameters such as fuel consumption as function of time by taking system’s configuration as input. The model does not take firing rate as an input and estimates the mass flow rate in the system by using the pressure difference available across the system. The model has provision for components like Combustion chamber, stack, pipe, blower and damper. Rishabh et al. [2] modified the model in order to bear close resemblance to the field practices and validated the project with a diverse set of experimental data and also created a user friendly interface.

In order to develop the numerical model of biomass based pottery furnance which is also a natural draft type system, we have used some of their model which is common in our project and modied them according to our requirement so that it become suitable for the analysis of pottery furnace. The special features of the furnace concerned here like the rat trap bond construction of wall and a separate existance of pottery chamber needed to be modelled separately which has been looked into the project.

## METHODOLOGY

For the analysis of pottery furnace, we have divided the furnace into N control volumes, which includes inlet holes, grate, firing chamber and pottery chamber. Then the conservation equation of mass, momentum, and energy has been applied to each and every control volume. Taking mass flow rate of fuel (constant), mass flow rate of air for one hole (which changes for different time steps) as an inlet the optimum value of mass flow rate through other holes has been found by applying boundary condition of pressure as a constraint.

Then the temperature distribution at different control volume of the pottery furnace is found solving the set of equations for different control volumes in MATLAB for a certain time step. For this we first wrote the conservation equation of mass, energy and momentum for each control volume, which results in a total of 8N+12 equations. We wrote all these equations in explicit form and then solved it simultaneously using non-linear multi variable Newton Raphson Method [6]. The iteration is then carried out repeatedly for all time steps which is equal to the total time during which the experiment of firing the pottery furnace was completed divided by the magnitude of each time step. The value of pressure and temperature obtained for different control volume in each time step is recorded and is used as an input for another time step. The value of temperature and pressure for the last time step is then taken as the result.

The result obtained through the numerical computation is then verified with the experimental data obtained at the time of firing the pottery furnance. If the obtained result differ from the obtained experimental data then the simulation is again continued varying the mass flow rate of air, untill the desired output is obtained.

The density has been considered dependent only on temperature and not on pressure for the ease of computations.

Gas/Air constant =

Density of Gas/Air =

## Modelling of Components

Inlet

At the bottom of the pottery chamber, an inlet is provided for charging the fuel and to allow passage of air. This inlet cross-section is separated into two parts due to the presence of fuel bed, thus allowing air two ways to pass through. Also three holes are provided at some distance below the grate to provide sufficient air necessary for complete combustion of fuel.

Various losses associated with the flow of air through the inlet are, as in [3]

Entrance loss (Contraction) =   ...(1)

Expansion and Contraction losses at inlet

Exit loss (Expansion) =     ...(2)

Pressure loss due to bend =       ...(3)

Fuel Bed

The fuel bed separates the inlet of the pottery chamber into two parts, therefore the air rushes into the pottery chamber from two different paths. The air which enters through the hole below the fuel bed has to pass through the fuel bed which acts like a porous media, due to which the pressure of the flow decreases.

Ergun's equation as by Oza et al. [1] is used to estimate the pressure drop. Thus

...(4) Where Dp is the effective fuel particle diameter. It is estimated using the equation while, porosity of the fuel bed is taken as input.

...(5)

Firing Chamber

It has been assumed that as soon as air enters the firing chamber, combustion of fuel happens inside the chamber and flue gases are formed instantly without any delay. Inside the firing chamber, the exhaust gas gains heat due to the combustion of fuel and loses heat to the walls of firing chamber via convection and radiation. Radiation is treated as a surface phenomenon for the walls of the firing chamber whereas it is treated as a volumetric phenomenon for the exhaust gas in the firing chamber. The first law of thermodynamics equation used is:

...(6)

Grate

Grate is provided in the pottery furnace above the inlet and the firing chamber on which the pottery to be baked is kept. The exhaust gases produced in the combustion chamber passes through the grate to the pottery chamber. The pressure of the exhaust gases changes while passing through the grate due to sudden contraction and expansion, this loss in pressure due to the presence of grate is found using given expression.

Grate in Pottery Furnance at Khurrampur

Pressure loss due to grate contraction =   ...(7)

if

else =

Pressure loss due to grate expansion =     ...(8)

Pottery Chamber

The compartment above the grate where potteries are kept is termed as pottery chamber. The exhaust gases produced in the combustion chamber comes in contact with the pottery surface in the pottery chamber. We have divide the entire pottery chamber in different control volumes. Each control volume is considered as a porous medium composed of spheres of the size of pottery. The porosity of medium is given by the formula;

...(9)

The height of each control volume is given by dx = $\frac{H-{h}_{grate}}{N}$  ...(10)

The first law equation for pottery in pottery chamber is;

...(11)

The energy equation for gas present in control volume is;

...(12)

The pressure drop in each control volume is evaluated using ergun's equation

as used by Oza et al.[1]

​ ...(13)

Wall

The wall of the entire pottery chamber is made in the form of a rat trap bond construction to provide better thermal insulation.

Rat trap bond construction [7]

The heat transfer equation through the wall is made considering both the conductive resistance of bricks and the radiative resistance of the the space between the two layers of bricks. To obtain the governing equation for heat transfer via conduction in the wall, it has been assumed that the difference in cross sectional surface area of innermost and outermost walls arising due to wall thickness is negligible. Heat transfer across the walls is assumed to happen in 1 dimension only. It has been assumed that the properties of the walls (conductivity, diffusivity etc.) remain same everywhere.

Electrical anology of thermal resistance of rat trap bond construction

The governing equation of heat transfer through different section of wall are as follows:

$\frac{K\;A_{wall\;}(T_{w1}\;-T_{w2})}t\;=\;(H_r+H_c)\;A_{wall}\;(T_g-T_{w1})$  ...(14)

$\frac{}{}$K Awall(Tw1-Tw2)t=K' Awall(Tw2-Tw3)t ...(15)

$\frac{{K}^{\text{'}}{A}_{wall}\left({T}_{w2}-{T}_{w3}\right)}{t}=\frac{K{A}_{wall}\left({T}_{w3}-{T}_{w4}\right)}{t}$  ...(16)

...(17)

Pressure Evaluation

The pressure evaluation for every control volume has been done using the following equation

..(18)

Determination of heat transfer co-efficient

H =   ...(19)

## RESULTS AND DISCUSSION

From reviewing the literature, appropriate heat transfer and fluid flow equations for the required model have been framed and discretised using finite difference method - forward in time and central in space scheme. All the equations have been written in implicit form and presented as a program in MATLAB and the jacobian computed, is solved by multi variable Newton-Raphson Method. The program is being analysed for results. The work is being continued by Raj Joshi, another summer research fellow of the academies for his remaining course of fellowship.

## Conclusion

Review of heat transfer processes was done and applied for the analysis of furnace. Conduction, Convection and Radiation were studied and accounted for it.

Compressible fluid flow was studied and accounted for change in pressure,velocity and density.

Numerical methods are used to solve all the above equations thereby developing a numerical model for the updraft natural draft furnace in MATLAB.

## REFERENCES

[1] Oza, H.I., Sai Anurag, M., Gaba, V., 2015, “Development of Generic Model to Analyze Natural Draft Furnace Systems,” B.Tech project thesis, Indian Institute of Technology, Delhi.

[2] Baid, R., Ravi Teja, K., 2017, "Development of a Generic Simulation Model to Analyze Natural Draft Furnace Systems", B.Tech project thesis, Indian Institute of Technology, Delhi.

[3] White, F.M. , 2012, Fluid Mechanics, McGraw-Hill, New Delhi.

[4] Incropera, F.P., Dewitt, D.P., Bergman, T.L., Lavine, A.S., 2007, Fundamentals of Heat and Mass Transfer, Wiley India, New Delhi.

[5] Bejan, A., 2006, Convection Heat Transfer, Wiley India, New Delhi.

[6] Chapra, S., Canale, R., 2010, Numerical Methods For Engineers, McGraw-Hill, New Delhi.

[7] https://architecturelive.in/of-bricks-and-bonds-the-rat-trap-bond/

## APPENDICES

Appendix 1: Input Data
 Parameter Values Atmosphere Ambient Temperature (K) 307.45 air constant 359.5 Hext 5 Cp for air (J/kg/K) 1005.7 Fuel Rgas 292 Calorific Value (J/kg) 45700000 Reference Temperature (K) 298 eps 0.1 Cp for flue gas (J/kg/K) 1160 Viscosity (kg/m/s) a 3E-08 mu=aT+b b 0.000007 gas constant 400 Thermal Conductivity (W/m/K) a 0.00007 b 0.0031 Pr gas 0.7 Fuel bed Porosity 0.3 Fuel bed thickness (m) 0.1 Fuel particle effective diameter (m) 0.007 Grate No of holes 24 Area of hole 0.016875 Geometry H 2.573 h channel 0.225 h inlet 0.675 D 1.61 width inlet 0.45 h grate 1.143 N 5 Time dt (s) 4 time_combustion (min) 92 Pottery Cp pot 800 rho pot 259 phi pot 0.6621 Wall n air hole 0 area above hole 0.84375 m1 0.01 k 0.75 k eff 0.4421 Asf total 20.285273 Dp pot 0.29295
Appendix 2: Fuel Feed Rate
 Variation of fuel flow rate with time time(min) mass flow rate(kg/s) 35 0.004915714 60 0.006866667 70 0.016716667 78 0.020416667 82 0.042 85 0.056111111 89 0.042083333 92 0.055555556 96 0.0415 99 0.056555556 103 0.042208333 106 0.056722222 113 0.023809524 116 0.055555556 119 0.055555556 128 0.014814815

## Appendix 3

The code prepared until now and under scrutiny for results is as follows:

global all

time_firing=all(29);

dt=all(28);

N=all(27);

D=all(24);

h_grate=all(26);

n=all(19);

a=all(20);

T_amb=all(1);

air_constant=all(2);

gas_constant=all(12);

Rgas=all(5);

Cp_gas=all(9);

Cv_gas=Cp_gas-Rgas;

Cp_air=all(4);

Cp_pot=all(30);

rho_pot=all(31);

CV=all(6);

g=9.81;

H=all(21);

dx=(H-h_grate)/N;

eps=all(8);

sig=5.67*10^(-8);

phi=all(32);

Hext=all(3);

k=all(36);

k_dash=all(37);

Pr_gas=all(15);

Asf=all(38)/N;

width_inlet=all(25);

h_inlet=all(23);

h_grate=all(26);

h_channel=all(22);

n_holes=all(33);

A1_ent=width_inlet*h_channel;

A2_ent=width_inlet*(h_inlet-h_channel);

A_holes=all(34);

phi_fuel=all(16);

Dp=all(39);

Dp_fuel=all(18);

h=h_grate-h_channel;

t=all(40);

ke_grate=(1-n*a/(pi/4*D^2))^2;

if n*a/(pi/4*D^2)<0.76^2

kc_grate=0.42*(1-n*a/(pi/4*D^2));

else

kc_grate=(1-n*a/(pi/4*D^2))^2;

end

F=zeros(8*N+12,1);

X=zeros(8*N+12,1);

X_old=zeros(8*N+12,1);

Z=zeros(8*N+12);

for i=1:5+6*N

%X(i)=T_amb;

X_old(i)=T_amb;

end

% for i=5+6*N+1:length(X)

% X(i)=0.1;

% end

p = xlsread('main.xlsx','fuel mass flow rate');

time = p(:,1);

Massf = p(:,2);

rr=1;

for i=1:length(time)

if rr*dt/60 <= time(i)

Mf=Massf(i);

break;

end

end

for j=1:10

Hr=eps*sig*(X(2)+X(1))*(X(2)^2+X(1)^2);

k_gas=all(13)*X(1)+all(14);

Hc=8.235*k_gas/D;

F(1)=X(5+6*N+2)*Cp_gas*X(1)/2+(Hr+Hc)*pi*D*h*X(1)-gas_constant*X_old(1)*pi*D^2*h*Cv_gas/4/dt/X(1)-(Hr+Hc)*pi*D*h*X(2)+X(5+6*N+2)*Cp_gas*X(6)/2-X(5+6*N+2)*Cp_gas*T_amb+gas_constant*Cv_gas*pi*D^2*h/4/dt-Mf*CV;

Z(1,1)=X(5+6*N+2)*Cp_gas/2+(Hr+Hc)*pi*D*h-gas_constant*X_old(1)*pi*D^2*h*Cv_gas/4/dt/X(1)^2;

Z(1,2)=-(Hr+Hc)*pi*D*h;

Z(1,6)=X(5+6*N+2)*Cp_gas/2;

Z(1,5+6*N+2)=Cp_gas*X(1)/2+Cp_gas*X(6)/2-Cp_gas*T_amb;

Hr=eps*sig*(X(7)+X(6))*(X(7)^2+X(6)^2);

k_gas=all(13)*X(6)+all(14);

Hc=8.235*k_gas/D;

Hsf=6*k_gas/D;

F(6)=-X(5+6*N+2)*Cp_gas*X(1)/2+X(5+6*N+3)*Cp_gas*X(6)/2+(Hr+Hc)*pi*D*dx*X(6)+Hsf*Asf*X(6)-gas_constant*X_old(6)*phi*pi*D^2*dx*Cv_gas/4/dt/X(6)-X(5+6*N+2)*Cp_gas*X(6)/2-Hsf*Asf*X(7)-(Hr+Hc)*pi*D*dx*X(8)+X(5+6*N+3)*Cp_gas*X(12)/2+X(5+6*N+2)*Cp_gas*T_amb-X(5+6*N+3)*Cp_gas*T_amb+gas_constant*Cv_gas*phi*pi*D^2*dx/4/dt;

Z(6,1)=-X(5+6*N+2)*Cp_gas/2;

Z(6,6)=X(5+6*N+3)*Cp_gas/2+(Hr+Hc)*pi*D*dx+Hsf*Asf-gas_constant*X_old(6)*phi*pi*D^2*dx*Cv_gas/4/dt/X(6)^2-X(5+6*N+2)*Cp_gas/2;

Z(6,7)=-Hsf*Asf;

Z(6,8)=-(Hr+Hc)*pi*D*dx;

Z(6,12)=X(5+6*N+3)*Cp_gas/2;

Z(6,5+6*N+2)=-Cp_gas*X(1)/2-Cp_gas*X(6)/2+Cp_gas*T_amb;

Z(6,5+6*N+3)=Cp_gas*X(6)/2+Cp_gas*X(12)/2-Cp_gas*T_amb;

for i=3:N

Hr=eps*sig*(X(5+6*(i-2)+2)+X(5+6*(i-2)+1))*(X(5+6*(i-2)+2)^2+X(5+6*(i-2)+1)^2);

k_gas=all(13)*X(5+6*(i-2)+1)+all(14);

Hc=8.235*k_gas/D;

Hsf=6*k_gas/D;

F(5+6*(i-2)+1)=-X(5+6*N+i)*Cp_gas*X(5+6*(i-3)+1)/2+X(5+6*N+i+1)*Cp_gas*X(5+6*(i-2)+1)/2+(Hr+Hc)*pi*D*dx*X(5+6*(i-2)+1)+Hsf*Asf*X(5+6*(i-2)+1)-gas_constant*X_old(5+6*(i-2)+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5+6*(i-2)+1)-X(5+6*N+i)*Cp_gas*X(5+6*(i-2)+1)/2-Hsf*Asf*X(5+6*(i-2)+2)-(Hr+Hc)*pi*D*dx*X(5+6*(i-2)+3)+X(5+6*N+i+1)*Cp_gas*X(5+6*(i-1)+1)/2+X(5+6*N+i)*Cp_gas*T_amb-X(5+6*N+i+1)*Cp_gas*T_amb+gas_constant*Cv_gas*phi*pi*D^2*dx/4/dt;

Z(5+6*(i-2)+1,5+6*(i-3)+1)=-X(5+6*N+i)*Cp_gas/2;

Z(5+6*(i-2)+1,5+6*(i-2)+1)=X(5+6*N+i+1)*Cp_gas/2+(Hr+Hc)*pi*D*dx+Hsf*Asf-gas_constant*X_old(5+6*(i-2)+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5+6*(i-2)+1)^2-X(5+6*N+i)*Cp_gas/2;

Z(5+6*(i-2)+1,5+6*(i-2)+2)=-Hsf*Asf;

Z(5+6*(i-2)+1,5+6*(i-2)+3)=-(Hr+Hc)*pi*D*dx;

Z(5+6*(i-2)+1,5+6*(i-1)+1)=X(5+6*N+i+1)*Cp_gas/2;

Z(5+6*(i-2)+1,5+6*N+i)=-Cp_gas*X(5+6*(i-3)+1)/2-Cp_gas*X(5+6*(i-2)+1)/2+Cp_gas*T_amb;

Z(5+6*(i-2)+1,5+6*N+i+1)=Cp_gas*X(5+6*(i-2)+1)/2+Cp_gas*X(5+6*(i-1)+1)/2-Cp_gas*T_amb;

end

Hr=eps*sig*(X(5+6*(N-1)+2)+X(5+6*(N-1)+1))*(X(5+6*(N-1)+2)^2+X(5+6*(N-1)+1)^2);

k_gas=all(13)*X(5+6*(N-1)+1)+all(14);

Hc=8.235*k_gas/D;

Hsf=6*k_gas/D;

F(5+6*(N-1)+1)=-X(5+6*N+N+1)*Cp_gas*X(5+6*(N-2)+1)/2+X(5+6*N+N+2)*Cp_gas*X(5+6*(N-1)+1)-X(5+6*N+N+1)*Cp_gas*X(5+6*(N-1)+1)/2+(Hr+Hc)*pi*D*dx*X(5+6*(N-1)+1)+Hsf*Asf*X(5+6*(N-1)+1)-gas_constant*X_old(5+6*(N-1)+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5+6*(N-1)+1)-Hsf*Asf*X(5+6*(N-1)+2)-(Hr+Hc)*pi*D*dx*X(5+6*(N-1)+3)+X(5+6*N+N+1)*Cp_gas*T_amb-X(5+6*N+N+2)*Cp_gas*T_amb+gas_constant*Cv_gas*phi*pi*D^2*dx/4/dt;

Z(5+6*(N-1)+1,5+6*(N-2)+1)=-X(5+6*N+N+1)*Cp_gas/2;

Z(5+6*(N-1)+1,5+6*(N-1)+1)=X(5+6*N+N+2)*Cp_gas-X(5+6*N+N+1)*Cp_gas/2+(Hr+Hc)*pi*D*dx+Hsf*Asf-gas_constant*X_old(5+6*(N-1)+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5+6*(N-1)+1)^2;

Z(5+6*(i-2)+1,5+6*(N-1)+2)=-Hsf*Asf;

Z(5+6*(i-2)+1,5+6*(N-1)+3)=-(Hr+Hc)*pi*D*dx;

Z(5+6*(i-2)+1,5+6*N+N+1)=-Cp_gas*X(5+6*(N-2)+1)/2-Cp_gas*X(5+6*(N-1)+1)/2+Cp_gas*T_amb;

Z(5+6*(i-2)+1,5+6*N+N+2)=Cp_gas*X(5+6*(N-1)+1)-Cp_gas*T_amb;

for i=2:N+1

k_gas=all(13)*X(5+6*(i-2)+1)+all(14);

Hsf=6*k_gas/D;

F(5+6*(i-2)+2)=-Hsf*Asf*X(5+6*(i-2)+1)+Hsf*Asf*X(5+6*(i-2)+2)+rho_pot*Cp_pot*(1-phi)*pi*D^2*dx/4/dt*X(5+6*(i-2)+2)-rho_pot*Cp_pot*(1-phi)*pi*D^2*dx/4/dt*X_old(5+6*(i-2)+2);

Z(5+6*(i-2)+2,5+6*(i-2)+1)=-Hsf*Asf;

Z(5+6*(i-2)+2,5+6*(i-2)+2)=Hsf*Asf+rho_pot*Cp_pot*(1-phi)*pi*D^2*dx/4/dt;

end

Hr=eps*sig*(X(2)+X(1))*(X(2)^2+X(1)^2);

k_gas=all(13)*X(1)+all(14);

Hc=8.235*k_gas/D;

F(2)=-(Hr+Hc)*X(1)+(Hr+Hc)*X(2)+k/t*X(2)-k/t*X(3);

Z(2,1)=-Hr-Hc;

Z(2,2)=Hr+Hc+k/t;

Z(2,3)=-k/t;

F(3)=k*X(2)-(k+k_dash)*X(3)+k_dash*X(4);

Z(3,2)=k;

Z(3,3)=-k-k_dash;

Z(3,4)=k_dash;

F(4)=k_dash*X(3)-(k+k_dash)*X(4)+k*X(5);

Z(4,3)=k_dash;

Z(4,4)=-k-k_dash;

Z(4,5)=k;

F(5)=k/t*X(4)-k/t*X(5)-Hext*X(5)+Hext*T_amb;

Z(5,4)=k/t;

Z(5,5)=-k/t-Hext;

for i=2:N+1

Hr=eps*sig*(X(5+6*(i-2)+2)+X(5+6*(i-2)+1))*(X(5+6*(i-2)+2)^2+X(5+6*(i-2)+1)^2);

k_gas=all(13)*X(5+6*(i-2)+1)+all(14);

Hc=8.235*k_gas/D;

F(5+6*(i-2)+3)=-(Hr+Hc)*X(5+6*(i-2)+1)+(Hr+Hc)*X(5+6*(i-2)+3)+k/t*X(5+6*(i-2)+3)-k/t*X(5+6*(i-2)+4);

Z(5+6*(i-2)+3,5+6*(i-2)+1)=-Hr-Hc;

Z(5+6*(i-2)+3,5+6*(i-2)+3)=Hr+Hc+k/t;

Z(5+6*(i-2)+3,5+6*(i-2)+4)=-k/t;

F(5+6*(i-2)+4)=k*X(5+6*(i-2)+3)-(k+k_dash)*X(5+6*(i-2)+4)+k_dash*X(5+6*(i-2)+5);

Z(5+6*(i-2)+4,5+6*(i-2)+3)=k;

Z(5+6*(i-2)+4,5+6*(i-2)+4)=-k-k_dash;

Z(5+6*(i-2)+4,5+6*(i-2)+5)=k_dash;

F(5+6*(i-2)+5)=k_dash*X(5+6*(i-2)+4)-(k+k_dash)*X(5+6*(i-2)+5)+k*X(5+6*(i-2)+6);

Z(5+6*(i-2)+5,5+6*(i-2)+4)=k_dash;

Z(5+6*(i-2)+5,5+6*(i-2)+5)=-k-k_dash;

Z(5+6*(i-2)+5,5+6*(i-2)+6)=k;

F(5+6*(i-2)+6)=k/t*X(5+6*(i-2)+5)-k/t*X(5+6*(i-2)+6)-Hext*X(5+6*(i-2)+6)+Hext*T_amb;

Z(5+6*(i-2)+6,5+6*(i-2)+5)=k/t;

Z(5+6*(i-2)+6,5+6*(i-2)+6)=-k/t-Hext;

end

F(5+6*N+1)=X(5+6*N+1)-X(5+6*N+N+3)-X(5+6*N+N+4)-X(5+6*N+N+5)-Mf;

Z(5+6*N+1,5+6*N+1)=1;

Z(5+6*N+1,5+6*N+N+3)=-1;

Z(5+6*N+1,5+6*N+N+4)=-1;

Z(5+6*N+1,5+6*N+N+5)=-1;

F(5+6*N+2)=pi*D^2*h*gas_constant/4/dt/X(1)-X(5+6*N+1)+X(5+6*N+2)-pi*D^2*gas_constant*h/4/dt/X_old(1);

Z(5+6*N+2,1)=pi*D^2*h*gas_constant/4/dt/X(1)^2;

Z(5+6*N+2,5+6*N+1)=-1;

Z(5+6*N+2,5+6*N+2)=1;

for i=3:N+1

F(5+6*N+i)=phi*pi*D^2*dx*gas_constant/4/dt/X(5+6*(i-3)+1)+X(5+6*N+i)-X(5+6*N+i-1)-phi*pi*D^2*dx*gas_constant/4/dt/X_old(5+6*(i-3)+1);

Z(5+6*N+i,5+6*(i-3)+1)=phi*pi*D^2*dx*gas_constant/4/dt/X(5+6*(i-3)+1)^2;

Z(5+6*N+i,5+6*N+i-1)=-1;

Z(5+6*N+i,5+6*N+i)=1;

end

F(5+6*N+N+2)=phi*pi*D^2*dx*gas_constant/4/dt/X(5+6*(N-1)+1)-X(5+6*N+N+1)+X(5+6*N+N+2)-phi*pi*D^2*dx*gas_constant/4/dt/X_old(5+6*(N-1)+1);

Z(5+6*N+N+2,5+6*(N-1)+1)=phi*pi*D^2*dx*gas_constant/4/dt/X(5+6*(N-1)+1)^2;

Z(5+6*N+N+2,5+6*N+N+1)=-1;

Z(5+6*N+N+2,5+6*N+N+2)=1;

mu_gas=all(10)*X(1)+all(11);

F(5+6*N+N+3)=X(5+6*N+N+5+1)-X(5+6*N+N+5+2)-g*h*gas_constant/X(1)-0.5*((X(1)+X(6))/2/gas_constant)*(X(5+6*N+2)/(pi/4*D^2))^2+0.5*T_amb/air_constant*(X(5+6*N+N+3)/A1_ent)^2+0.5*T_amb/air_constant*(X(5+6*N+N+4)/A2_ent)^2+0.5*n_holes*T_amb/air_constant*(X(5+6*N+N+5)/A_holes)^2-0.5*0.5*air_constant*(X(5+6*N+N+3)*T_amb/A1_ent/air_constant)^2/T_amb-(1-A1_ent/(D*h))^2*0.5*(X(5+6*N+N+3)*T_amb/A1_ent/air_constant)^2*air_constant/T_amb-0.5*(ke_grate+kc_grate)*(X(5+6*N+N+3)/(n*a))^2/air_constant*X(1)-0.25*0.5*(X(5+6*N+N+3)/(n*a))^2/air_constant*X(1)-150*mu_gas*(1-phi_fuel)^2*(X(5+6*N+N+3)*X(1)*4/gas_constant/pi/D^2)^2*h_inlet/(phi_fuel^3*Dp_fuel^2)-1.75*(1-phi_fuel)/gas_constant*X(1)*(X(5+6*N+N+3)*4/pi/D^2)^2*h_inlet/(phi_fuel^3*Dp_fuel)-0.5*(ke_grate+kc_grate)*(X(5+6*N+2)/(n*a))^2/gas_constant*(X(1)+X(6))/2;

Z(5+6*N+N+3,1)=-g*h*gas_constant/X(1)^2-0.5/2/gas_constant*(X(5+6*N+2)/(pi/4*D^2))^2-0.5*(ke_grate+kc_grate)*(X(5+6*N+N+3)/(n*a))^2/air_constant-0.25*0.5*(X(5+6*N+N+3)/(n*a))^2/air_constant-150*mu_gas*(1-phi_fuel)^2*X(5+6*N+N+3)*4/gas_constant/pi/D^2*h_inlet/(phi_fuel^3*Dp_fuel^2)*2*X(5+6*N+N+3)-1.75*(1-phi_fuel)/gas_constant*(4/pi/D^2)^2*h_inlet/(phi_fuel^3*Dp_fuel)*2*X(5+6*N+N+3);

Z(5+6*N+N+3,6)=-0.5/2/gas_constant*(X(5+6*N+2)/(pi/4*D^2))^2-0.5*(ke_grate+kc_grate)*(X(5+6*N+2)/(n*a))^2/gas_constant/2;

Z(5+6*N+N+3,5+6*N+2)=-0.5*((X(1)+X(6))/2/gas_constant)*(1/(pi/4*D^2))^2*2*X(5+6*N+2)-0.5*(ke_grate+kc_grate)*(1/(n*a))^2/gas_constant*(X(1)+X(6))/2*2*X(5+6*N+2);

Z(5+6*N+N+3,5+6*N+N+3)=0.5*T_amb/air_constant*(1/A1_ent)^2*2*X(5+6*N+N+3)-0.5*0.5*air_constant*(T_amb/A1_ent/air_constant)^2/T_amb*2*X(5+6*N+N+3)-(1-A1_ent/(D*h))^2*0.5*(T_amb/A1_ent/air_constant)^2*air_constant/T_amb*2*X(5+6*N+N+3)-0.5*(ke_grate+kc_grate)*(1/(n*a))^2/air_constant*X(1)*2*X(5+6*N+N+3)-0.25*0.5*(1/(n*a))^2/air_constant*X(1)*2*X(5+6*N+N+3)-150*mu_gas*(1-phi_fuel)^2*(X(1)*4/gas_constant/pi/D^2)^2*h_inlet/(phi^3*Dp_fuel^2)-1.75*(1-phi)/gas_constant*X(1)*(4/pi/D^2)^2*h_inlet/(phi^3*Dp)*2*X(5+6*N+N+3);

Z(5+6*N+N+3,5+6*N+N+4)=0.5*T_amb/air_constant*(1/A2_ent)^2*2*X(5+6*N+N+4);

Z(5+6*N+N+3,5+6*N+N+5)=0.5*n_holes*T_amb/air_constant*(1/A_holes)^2*2*X(5+6*N+N+5);

Z(5+6*N+N+3,5+6*N+N+5+1)=1;

Z(5+6*N+N+3,5+6*N+N+5+2)=-1;

F(5+6*N+N+4)=X(5+6*N+N+4)^2*(0.5*0.5*air_constant*(T_amb/A2_ent/air_constant)^2/T_amb+(1-A2_ent/(D*h))^2*0.5*(T_amb/A2_ent/air_constant)^2*air_constant/T_amb+0.25*0.5*(T_amb/A2_ent/air_constant*A2_ent/(D*h))^2*air_constant/X(1))-X(5+6*N+N+3)^2*(0.5*0.5*air_constant*(T_amb/A1_ent/air_constant)^2/T_amb+(1-A1_ent/(D*h))^2*0.5*(T_amb/A1_ent/air_constant)^2*air_constant/T_amb+0.5*(ke_grate+kc_grate)*(X(1)/air_constant/(n*a))^2*air_constant/X(1)+0.25*0.5*(X(1)/air_constant/(n*a))^2*air_constant/X(1));

Z(5+6*N+N+4,1)=X(5+6*N+N+4)^2*0.25*0.5*(T_amb/A2_ent/air_constant*A2_ent/(D*h))^2*air_constant/X(1)^2-X(5+6*N+N+3)^2*(0.5*(ke_grate+kc_grate)*(1/air_constant/(n*a))^2*air_constant+0.25*0.5*(1/air_constant/(n*a))^2*air_constant);

Z(5+6*N+N+4,5+6*N+N+3)=-2*X(5+6*N+N+3)*(0.5*0.5*air_constant*(T_amb/A1_ent/air_constant)^2/T_amb+(1-A1_ent/(D*h))^2*0.5*(T_amb/A1_ent/air_constant)^2*air_constant/T_amb+0.5*(ke_grate+kc_grate)*(X(1)/air_constant/(n*a))^2*air_constant/X(1)+0.25*0.5*(X(1)/air_constant/(n*a))^2*air_constant/X(1));

Z(5+6*N+N+4,5+6*N+N+4)=2*X(5+6*N+N+4)*(0.5*0.5*air_constant*(T_amb/A2_ent/air_constant)^2/T_amb+(1-A2_ent/(D*h))^2*0.5*(T_amb/A2_ent/air_constant)^2*air_constant/T_amb+0.25*0.5*(T_amb/A2_ent/air_constant*A2_ent/(D*h))^2*air_constant/X(1));

mu_gas=all(10)*X(1)+all(11);

F(5+6*N+N+5)=X(5+6*N+N+5)^2*(0.5*0.5*air_constant*(T_amb/A_holes/air_constant)^2/T_amb+(1-A_holes/(D*h))^2*0.5*(T_amb/A_holes/air_constant)^2*air_constant/T_amb+0.25*0.5*(T_amb/A_holes/air_constant*A_holes/(D*h))^2*air_constant/X(1))-n_holes^2*(X(5+6*N+N+3)*150*mu_gas*(1-phi_fuel)^2*X(1)*4/gas_constant/pi/D^2*(h_grate-h_channel)/(phi_fuel^3*Dp^2)+X(5+6*N+N+3)^2*(0.5*0.5*air_constant*(T_amb/A1_ent/air_constant)^2/T_amb+(1-A1_ent/(D*h))^2*0.5*(T_amb/A1_ent/air_constant)^2*air_constant/T_amb+0.5*(ke_grate+kc_grate)*(X(1)/air_constant/(n*a))^2*air_constant/X(1)+0.25*0.5*(X(1)/air_constant/(n*a))^2*air_constant/X(1)+1.75*(1-phi_fuel)*gas_constant/X(1)*(X(1)*4/gas_constant/pi/D^2)^2*(h_grate-h_channel)/(phi_fuel^3*Dp)));

Z(5+6*N+N+5,1)=X(5+6*N+N+5)^2*0.25*0.5*(T_amb/A_holes/air_constant*A_holes/(D*h))^2*air_constant/X(1)^2-n_holes^2*(X(5+6*N+N+3)*150*mu_gas*(1-phi_fuel)^2*4/gas_constant/pi/D^2*(h_grate-h_channel)/(phi_fuel^3*Dp^2)+X(5+6*N+N+3)^2*(0.5*(ke_grate+kc_grate)*(1/air_constant/(n*a))^2*air_constant+0.25*0.5*(1/air_constant/(n*a))^2*air_constant+1.75*(1-phi_fuel)*gas_constant*(4/gas_constant/pi/D^2)^2*(h_grate-h_channel)/(phi_fuel^3*Dp)));

Z(5+6*N+N+5,5+6*N+N+3)=n_holes^2*(150*mu_gas*(1-phi_fuel)^2*X(1)*4/gas_constant/pi/D^2*(h_grate-h_channel)/(phi_fuel^3*Dp^2)+2*X(5+6*N+N+3)*(0.5*0.5*air_constant*(T_amb/A1_ent/air_constant)^2/T_amb+(1-A1_ent/(D*h))^2*0.5*(T_amb/A1_ent/air_constant)^2*air_constant/T_amb+0.5*(ke_grate+kc_grate)*(X(1)/air_constant/(n*a))^2*air_constant/X(1)+0.25*0.5*(X(1)/air_constant/(n*a))^2*air_constant/X(1)+1.75*(1-phi_fuel)*gas_constant/X(1)*(X(1)*4/gas_constant/pi/D^2)^2*(h_grate-h_channel)/(phi_fuel^3*Dp)));

Z(5+6*N+N+5,5+6*N+N+5)=2*X(5+6*N+N+5)*(0.5*0.5*air_constant*(T_amb/A_holes/air_constant)^2/T_amb+(1-A_holes/(D*h))^2*0.5*(T_amb/A_holes/air_constant)^2*air_constant/T_amb+0.25*0.5*(T_amb/A_holes/air_constant*A_holes/(D*h))^2*air_constant/X(1));

F(5+6*N+N+5+1)=X(5+6*N+N+5+1)-air_constant*g*H;

Z(5+6*N+N+5+1,5+6*N+N+5+1)=1;

mu_gas=all(10)*X(6)+all(11);

F(5+6*N+N+5+2)=X(5+6*N+N+5+2)-X(5+6*N+N+5+3)+0.5*(X(1)+X(6))/2/gas_constant*(X(5+6*N+2)/(phi*pi/4*D^2))^2-0.5*(X(6)+X(12))/2/gas_constant*(X(5+6*N+3)/(phi*pi/4*D^2))^2-g*dx*gas_constant/X(6)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(6)*(X(5+6*N+2)+X(5+6*N+3))/(2*gas_constant*phi*pi/4*D^2))^2-1.75*(1-phi)/phi*3/Dp*dx*X(6)/gas_constant*((X(5+6*N+2)+X(5+6*N+3))/(2*phi*pi/4*D^2))^2;

Z(5+6*N+N+5+2,1)=0.5/2/gas_constant*(X(5+6*N+2)/(phi*pi/4*D^2))^2;

Z(5+6*N+N+5+2,6)=0.5/2/gas_constant*(X(5+6*N+2)/(phi*pi/4*D^2))^2-0.5/2/gas_constant*(X(5+6*N+3)/(phi*pi/4*D^2))^2-g*dx*gas_constant/X(6)^2-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*((X(5+6*N+2)+X(5+6*N+3))/(2*gas_constant*phi*pi/4*D^2))^2*2*X(6)-1.75*(1-phi)/phi*3/Dp*dx/gas_constant*((X(5+6*N+2)+X(5+6*N+3))/(2*phi*pi/4*D^2))^2;

Z(5+6*N+N+5+2,12)=-0.5/2/gas_constant*(X(5+6*N+3)/(phi*pi/4*D^2))^2;

Z(5+6*N+N+5+2,5+6*N+2)=0.5*(X(1)+X(6))/2/gas_constant*(1/(phi*pi/4*D^2))^2*2*X(5+6*N+2)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(6)/(2*gas_constant*phi*pi/4*D^2))^2*2*(X(5+6*N+2)+X(5+6*N+3));

Z(5+6*N+N+5+2,5+6*N+3)=-0.5*(X(12)+X(6))/2/gas_constant*(1/(phi*pi/4*D^2))^2*2*X(5+6*N+3)-1.75*(1-phi)/phi*3/Dp*dx*X(6)/gas_constant*(1/(2*phi*pi/4*D^2))^2*2*(X(5+6*N+2)+X(5+6*N+3));

Z(5+6*N+N+5+2,5+6*N+N+5+2)=1;

Z(5+6*N+N+5+2,5+6*N+N+5+3)=-1;

for i=3:N+1

mu_gas=all(10)*X(5+6*(i-2)+1)+all(11);

F(5+6*N+N+5+i)=X(5+6*N+N+5+i)-X(5+6*N+N+5+i+1)+0.5*(X(5+6*(i-3)+1)+X(5+6*(i-2)+1))/2/gas_constant*(X(5+6*N+i)/(phi*pi/4*D^2))^2-0.5*(X(5+6*(i-2)+1)+X(5+6*(i-1)+1))/2/gas_constant*(X(5+6*N+i+1)/(phi*pi/4*D^2))^2-g*dx*gas_constant/X(5+6*(i-2)+1)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(5+6*(i-2)+1)*(X(5+6*N+i)+X(5+6*N+i+1))/(2*gas_constant*phi*pi/4*D^2))^2-1.75*(1-phi)/phi*3/Dp*dx*X(5+6*(i-2)+1)/gas_constant*((X(5+6*N+i)+X(5+6*N+i+1))/(2*phi*pi/4*D^2))^2;

Z(5+6*N+N+5+i,5+6*(i-3)+1)=0.5/2/gas_constant*(X(5+6*N+i)/(phi*pi/4*D^2))^2;

Z(5+6*N+N+5+i,5+6*(i-2)+1)=0.5/2/gas_constant*(X(5+6*N+i)/(phi*pi/4*D^2))^2-0.5/2/gas_constant*(X(5+6*N+i+1)/(phi*pi/4*D^2))^2-g*dx*gas_constant/X(5+6*(i-2)+1)^2-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*((X(5+6*N+i)+X(5+6*N+i+1))/(2*gas_constant*phi*pi/4*D^2))^2*2*X(5+6*(i-2)+1)-1.75*(1-phi)/phi*3/Dp*dx/gas_constant*((X(5+6*N+i)+X(5+6*N+i+1))/(2*phi*pi/4*D^2))^2;

Z(5+6*N+N+5+i,5+6*(i-1)+1)=-0.5/2/gas_constant*(X(5+6*N+i+1)/(phi*pi/4*D^2))^2;

Z(5+6*N+N+5+i,5+6*N+i)=0.5*(X(5+6*(i-3)+1)+X(5+6*(i-2)+1))/2/gas_constant*(1/(phi*pi/4*D^2))^2*2*X(5+6*N+i)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(5+6*(i-2)+1)/(2*gas_constant*phi*pi/4*D^2))^2*2*(X(5+6*N+i)+X(5+6*N+i+1));

Z(5+6*N+N+5+i,5+6*N+i+1)=-0.5*(X(5+6*(i-1)+1)+X(5+6*(i-2)+1))/2/gas_constant*(1/(phi*pi/4*D^2))^2*2*X(5+6*N+i+1)-1.75*(1-phi)/phi*3/Dp*dx*X(5+6*(i-2)+1)/gas_constant*(1/(2*phi*pi/4*D^2))^2*2*(X(5+6*N+i)+X(5+6*N+i+1));

Z(5+6*N+N+5+i,5+6*N+N+5+i)=1;

Z(5+6*N+N+5+i,5+6*N+N+5+i+1)=-1;

end

F(5+6*N+N+5+N+2)=X(5+6*N+N+5+N+2);

Z(5+6*N+N+5+N+2,5+6*N+N+5+N+2)=1;

X=X-Z\F;

X

end

More
Written, reviewed, revised, proofed and published with