Summer Research Fellowship Programme of India's Science Academies

Development of numerical model for analysis of biomass based pottery furnace

Raj Joshi

Department of Civil Engineering, A. P. Shah Institute of Technology (University of Mumbai), Kasarvadavali, Thane 400615

Prof. M. R. Ravi

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


A natural draft furnace for firing pottery has been installed in Village Khurrampur, Tehsil Farukhnagar, District Gurgaon, Haryana. The walls of this furnace have been built in rat trap bond to reduce heat dissipated to the atmosphere and holes are provided in the wall just below the grate to place pottery, to improve the supply of air to the fuel. Since the furnace works on natural draft mechanism, electricity is not required to fire it and can be used only with the availability of biomass as fuel. The aim of the project is to evaluate the efficiency of this furnace, in terms of heat lost to atmosphere and supply of appropriate quantity of air. Previously, research has been done in Indian Institute of Technology Delhi to develop a generic model of furnaces which can be used to analyse and design furnaces of varied geometries, use and fuel. The furnace concerned here differs from the generic model in the aspects of multiple entry routes for air and placing the pottery above the grate, distinct from the combustion chamber. To accommodate these features of the above furnace, the generic numerical model is modified to be specific to this case of furnace. The furnace is divided into sections of firing chamber and pottery chamber to individually model them. The sections have been divided into small Control Volumes to ensure accuracy of the numerical model. The equations of fluid flow and heat transfer have been discretized using Finite Difference Method, forward in time and central in space implicit scheme. A program is written in MATLAB to run the numerical computations. The dimensions of furnace are being taken as input from the user in form of an MS Excel sheet and thus this model can be used for all furnaces with similar geometry and working. This work has been done along with Mr. Rahul Jha, another Summer Research Fellow of The Academies.

Keywords: natural draft furnace, rat trap bond construction, finite difference method, MATLAB


 PPressure (Pa) 
  ρ  Density ( kg/m3kg/m^3)
 TTemperature (K) 
 QHeat flux (W/m2)(W/m^2)
 vVelocity (m/s) 
 gAcceleration due to gravity (m/s2)(m/s^2)
A Cross sectional area (m2)(m^2)
 ​ ϵ\epsilonEmissivity of air/gas 
σ\sigmaStephan-Boltzmann constant (W/m2/T4)(W/m^2/T^4)​ 
 RGas constant (J/kg/K) 
 k' Effective conductivity of middle part of wall (W/m/K) 
 HHeat transfer coefficient (W/m2/K)(W/m^2/K)
 DpEffective Particle Diameter (m)  
 μCoefficient of dynamic viscosity(kg/m/s) 
 DDiameter (m) 
 dxLength(m) of each control volume in Pottery Chamber 
 CVCalorific Value (J/kg)  
 CvSpecific heat capacity at constant volume(J/kg/K)  
 CpSpecific heat capacity at constant Pressure (J/kg/K)  
 VVolume (m3)(m^3)
m Mass floe rate (kg/s) 
 kThermal conductivity (W/m/K) 
 ttime (s) 
 KCo-efficient of Pressure loss  
 out Outlet
 g / gasFlue gas 
 amb / a Ambient 
 iNo of control volume 
 tAt time t 
 t+ΔtAt time t+Δt 



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 use both biomass or LPG for firing. In villages small potters do not use LPG due to its unavailability or higher cost. Also pottery furnaces can be both natural draft or forced draft. Since the natural draft furnace systems do not require any electrical supply they are widely used in rural area for numerous applications. Comprehensive understanding of functioning of forced draft furnace is available to the scientific community but such widespread insight about natural draft furnace system is not available. Owing to the widespread use of natural draft furnace systems, scientific community in developing countries started working on gaining insights into the energy distribution and performance of these systems. If not designed properly, Natural draft furnaces tend to have poor efficiency. Hence 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 furnace system in IIT Delhi.

Biomass based pottery furnace has been installed by at village Khurrampur, tehsil Farukhnagar, district Gurgaon, Haryana. This project is about to develop a numerical model for the 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 properties of fuel and furnace and fuel feed rate as an input.

About Furnace

furnace elevation.PNG
    Elevation of Furnace

    The pottery furnace, concerned here, is in the form of a hollow cylinder, whose front view has been shown in Figure No. 1. It is surrounded by brick walls in rat trap bond construction on all sides and is open on the top. An inlet has been provided at the bottom on one side to insert and fire the fuel. A grate/channel has been provided to place the fuel. This reduces the contact between burning fuel and ground thereby saving the heat energy consumed by the ground. A grate has been provided above the firing unit to place the pottery. Few small holes have been provided in the brick wall just below the grate to enhance supply of air and ensure complete combustion.

    rat trap bond.PNG
      Rat trap bond construction

      The rat trap bond construction used here is one brick thick with plan as shown in Figure No. 2. The voids between the inner and outer layers of bricks reduce the weight of wall. Also the heat conducted by the wall to the atmosphere is reduced as the conductivity of air is less than that of bricks. Therefore the amount of heat retained by the system increases.

      Statement of the Problems

      • To create a numerical model that can estimate the energy distribution (i.e. temperature and pressure) at different intervals of the biomass based pottery furnace.
      • To estimate the effect of providing aeration holes for air flow and suggesting their numbers and size based on the performance of the pottery furnace.
      • To check for the saving in wastage of energy due to rat trap bond construction of walls.


      There was no exhaustive literature available as a generic model of natural draft furnaces; only case specific analysis had been done. So, 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, crucible, stack, pipe, blower and damper.

      Baid 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. They included dominating physical phenomena like radiation in their model. They provided different cases of crucibles so that user can choose the most resembling with the use. They provided enhanced types of fuel in the library of code thereby making the model more generic. A Graphical User Interface was made by them to run the simulations.

      In order to develop the numerical model of the above biomass based pottery furnace at Khurrampur, which is also a natural draft type system, the generic models developed by Oza et.al. and Baid et.al. needed to be modified the distinct features of this furnace which were not a part of generic models. The special features of the furnace concerned here like the rat trap bond construction of wall and a separate existence of pottery chamber needed to be modelled separately which has been looked into the project.



      The Furnace has been dividing into 2 sections: Firing Chamber and Pottery Chamber. The part below the grate to place potteries is designated as firing chamber where the fuel is stacked and fired. The part above grate is designated as pottery chamber, where hot gases pass through the stacked pottery. The pottery chamber is divided into 'N' control volumes, N being specified by the user. The flow of flue gases in pottery chamber through the large number of stacked pots is considered as flow through porous media.

      All the dimensions concerned with the furnace, the experimental rate of fuel input (in kg/s) and the properties of fuel bed and pottery chamber are to be specified by the user in the form of MS Excel sheet taken by the program as input. Example of such sheet has been attached as Appendix No. 1.

      The air coming into the firing chamber is considered from three paths. There will be a part of air entering from below the channel, passing through it and then entering the fuel bed. This has been designated as path 1. The primary air that will enter through the main entrance above the channel is designated as path 2. And the air coming from the aeration holes is considered as path 3. The ratio of the flow rates of different paths of air is estimated using the relationship between the pressure losses.


      wall sections_1.PNG
        Sections of the wall

        The walls have been divided into 3 sections as shown in Figure No. 3. Each section is of the width of the header as shown in figure No. 2. The inner and outer sections comprise only of the bricks, their thermal conductivity designated as k. The middle section consists of alternate sections of brick and air, as can be seen in Figure No. 2. The thermal conductivity of section 2 (k') is determined as a combination of conductivities of brick and air in parallel, as in Figure no. 4. The resistance for a thermal conductivity 'k' over length 'L' and cross sectional area 'A' is expressed as Lk  A\frac L{k\;A} [4]​

        effective conductivity.PNG
          schematic to calculate effective conductivity of section 2 of the wall

          Let k' is the effective conductivity and A be cross sectional area against one header,


          k=kair4+3kb4k'=\frac{k_{air}}4+\frac{3k_b}4​ ......(1)

          Heat Transfer Coefficients

          In the pottery chamber, the flow of hot gases is treated as convection through porous media. The heat transfer coefficient between flue gases and pottery is determined using the nusselt number in case of flow through porous media [5]

          HpH_p= Nu.kgasD\frac{Nu.k_{gas}}D ; Nu = 6 ......(2)

          Similarly the heat transfer coefficient between flue gas and wall is determined treating it as a case of flow through 2 large paltes [5]

          Hw=Nu.kgasDH_w=\frac{Nu.k_{gas}}D ; Nu=8.235 ......(3)

          In general Heat radiating from a substance with temperature T1 to another substance with temperature T2 is expressed as [4]

          Q=ϵ.σ.(T14T24)Q=\epsilon.\sigma.(T1^4-T2^4) ......(4)


          The thermal conductivity and dynamic viscosity of flue gas is expressed in the form of following expressions, whose coefficients have to be given by the user as input

          kg=aTg+bk_g=aT_g+b ......(5)

          μg=mTg+n\mu_g=mT_g+n ......(6)

          The density of flue gases is assumed to be negligibly affected by pressure change and thus affected only by temperature change. Thus it is expressed as

          ρg=PambTg ; Pamb=101325 Pa ......(7)

          The specific heat and universal gas constants of flue gas and air are assumed as constants, not varying with change in pressure or temperature.


          There are 5 temperature variables assigned in each control volume. The gas temperature, Grate temperature and 3 temperatures of the sections of the wall in the first control volume of the firing chamber. The entire firing chamber is considered as a single control volume since most of its region would be at same temperature because of the fuel burning in it. The grate over which pottery is placed is considered as a part of this control volume. In the remaining N control volumes of the pottery chamber, the gas temperature, pottery temperature and 3 temperatures of wall sections form the 5 variables. First Law of Thermodynamics is applied to the gas, grate, pottery and the wall sections in each control volume. Thus we have total 5(N+1) temperature variables and those many equations.

          The mass rate of flow of gas (kg/s) at inlet of each control volume is declared as a variable. Also the mass flow rate at outlet of last control volume is declared a variable. The relationship between successive inlet mass flow rates is determined using the continuity equation. Thus there are (N+2) mass flow rate variables and these many equations. Also 3 inlet air flow rates from the 3 distinct paths have been declared variables.

          As we have (N+1) control volumes, the pressure values at (N+2) nodes have been declared as variables. The pressure at the topmost node is assumed as zero, and pressure at the bottommost node as ρamb.g.H\rho_{amb}.g.H.

          Thus the total number of variables and equations is 5(N+1)+(N+2)+3+(N+2)=7N+12.


          All the equations are discretized using forward in time implicit scheme. Thereby, there is no data at any time step assumed from the previous time step barring the unsteady terms.

          All these (7N+12) equations are written in mathematically implicit form as


          These equations are solved simultaneously using multivariable iterative Newton-Raphson Method [6]

          Xi+1=Xi    Zi1.FiX_{i+1}=X_i\;-\;Z_i^{-1}.F_i ......(8)

          X  =  [x1,  x2,  x3,  ....  ,  xn]F=  [  f1,  f2,  f3,  .....  ,  fn]Z=[f1x1f1x2....f1xnf2x1......f2xn....................fnx1fnx2....fnxn]\begin{array}{l}\begin{array}{l}X\;=\;\lbrack x_1,\;x_2,\;x_3,\;....\;,\;x_n\rbrack\\F=\;\lbrack\;f_1,\;f_2,\;f_3,\;.....\;,\;f_n\rbrack\\Z=\begin{bmatrix}\frac{\partial f_1}{\partial x_1}&\frac{\partial f_1}{\partial x_2}&..&..&\frac{\partial f_1}{\partial x_n}\\\frac{\partial f_2}{\partial x_1}&..&..&..&\frac{\partial f_2}{\partial x_n}\\..&..&..&..&..\\..&..&..&..&..\\\frac{\partial f_n}{\partial x_1}&\frac{\partial f_n}{\partial x_2}&..&..&\frac{\partial f_n}{\partial x_n}\end{bmatrix}\\\end{array}\\\end{array}

          where X= vector containing the values of all variables

          F= vector containing the values of all LHS of the functions designated above

          Z= Jacobian of the functions designated above

          This process is repeated for all time steps and the temperatures, pressures and flow rates calculated successively.


          Firing Chamber

          First law of Thermodynamics is applied to the flue gases as

          ρg(1).Cvg.V.dTg(i)dt=Min(1).Cpair.(TambTref)Min(2).Cpg.(Tg(1)Tref)+Mf.CV+(Hr+Hw)Aw(Tw1(1)Tg(1))+(Hr+Hw)Agrate(Tgrate(1)Tg(1))\rho_g(1).Cv_g.V.\frac{\operatorname dT_g(i)}{\operatorname dt}=M_{in}(1).Cp_{air}.(T_{amb}-T_{ref})-M_{in}(2).Cp_g.(T_g(1)-T_{ref})+Mf.CV+(H_r+H_w)A_w(T_{w1}(1)-T_g(1))+(H_r+H_w)A_{grate}(T_{grate}(1)-T_g(1)) ......(9)

          upon discretization, LHS of above equation becomes ρg(1).Cvg.V.Tg(1)t+δtTg(1)tδt\rho_g(1).Cv_g.V.\frac{T_g(1)^{t+\delta t}-T_g(1)^t}{\delta t}

          First Law of Thermodynamics is applied to the grate as

          Mgrate.Cpgrate.dTgratedt=(Hr+Hw).(Tg(1)Tgrate)M_{grate}.Cp_{grate}.\frac{\operatorname dT_{grate}}{\operatorname dt}=(H_r+H_w).(T_g(1)-T_{grate}) ......(10)

          upon discretization, LHS of above equation becomes Mgrate.Cpgrate.Tgratet+δtTgratetδtM_{grate}.Cp_{grate}.\frac{T_{grate}^{t+\delta t}-T_{grate}^t}{\delta t}

          Pottery Chamber

          First law of Thermodynamics is applied to the flue gases as

          ρg(i).Cvg.Vg.dTg(i)dt=Min(i).Cpg.(Tin(i)Tref)Min(i+1).Cpg.(Tout(i)Tref)+Mf.CV+(Hr+Hw)Aw(Tw(i)Tg(i))+Hp.Ap.(Tp(i)Tg(i))\rho_g(i).Cv_g.V_g.\frac{\operatorname dT_g(i)}{\operatorname dt}=M_{in}(i).Cp_{g}.(T_{in}(i)-T_{ref})-M_{in}(i+1).Cp_g.(T_{out}(i)-T_{ref})+Mf.CV+(H_r+H_w)A_w(T_{w}(i)-T_g(i))+H_p.A_p.(T_p(i)-T_g(i)) ......(11)

          upon discretization, LHS of above equation becomes ​ ρg(i).Cvg.V.Tg(i)t+δtTg(i)tδt\rho_g(i).Cv_g.V.\frac{T_g(i)^{t+\delta t}-T_g(i)^t}{\delta t}

          First Law of Thermodynamics is applied to the grate as

          ρp.Vp.CppdTp(i)dt=Ap.Hp.(Tg(i)Tp(i))\rho_p.V_p.Cp_p\frac{\operatorname dT_p(i)}{\operatorname dt}=A_p.H_p.(T_g(i)-T_p(i)) ......(12)

          upon discretization, LHS of above equation becomes ​ ρp.Cpp.V.Tp(i)t+δtTp(i)tδt\rho_p.Cp_p.V.\frac{T_p(i)^{t+\delta t}-T_p(i)^t}{\delta t}


          Referring figure 3,

          ρb.Vb.Cpb.dTw1dt=(Hr+Hc).(TgTw1)Tw1Tw2t2.(1k+1k)\rho_b.V_b.Cp_b.\frac{dT_{w1}}{dt}=(H_r+H_c).(T_g-T_{w1})-\frac{T_{w1}-T_{w2}}{\frac{t}{2}.(\frac{1}{k}+\frac{1}{k'})} ......(13)

          upon discretization, LHS of above equation becomes ρb.Cpb.Vb.Tw1t+δtTw1tδt\rho_b.Cp_b.V_b.\frac{T_{w1}^{t+\delta t}-T_{w1}^t}{\delta t}

          ρb.Vb.Cpb.dTw2dt=Tw1Tw2t2.(1k+1k)Tw2Tw3t2.(1k+1k)\rho_b.V'_b.Cp_b.\frac{dT_{w2}}{dt}=\frac{T_{w1}-T_{w2}}{\frac{t}{2}.(\frac{1}{k}+\frac{1}{k'})}-\frac{T_{w2}-T_{w3}}{\frac{t}{2}.(\frac{1}{k}+\frac{1}{k'})} ......(14)

          upon discretization, LHS of above equation becomes ρb.Cpb.Vb.Tw2t+δtTw2tδt\rho_b.Cp_b.V'_b.\frac{T_{w2}^{t+\delta t}-T_{w2}^t}{\delta t}

          ρb.Vb.Cpb.dTw1dt=Tw2Tw3t2.(1k+1k)Hamb.(Tw3Tamb)\rho_b.V_b.Cp_b.\frac{dT_{w1}}{dt}=\frac{T_{w2}-T_{w3}}{\frac{t}{2}.(\frac{1}{k}+\frac{1}{k'})}-H_{amb}.(T_{w3}-T_{amb}) ......(15)

          upon discretization, LHS of above equation becomes ρb.Cpb.Vb.Tw3t+δtTw3tδt\rho_b.Cp_b.V_b.\frac{T_{w3}^{t+\delta t}-T_{w3}^t}{\delta t}

          Continuity Equations

          The inlet mass flow rate in the first control volume is expressed as,

          Min(1)=M1+M2+n.M3+MfM_{in}(1)=M_1+M_2+n.M_3+M_f ......(16)

          The inlet mass flow rate for the further control volumes are computed using the equation,

          Mint+δt(i+1)=Mint+δt(i)+Vgρgt+δtρgtδtM_{in}^{t+\delta t}(i+1)=M_{in}^{t+\delta t}(i)+V_g\frac{\rho_g^{t+\delta t}-\rho_g^t}{\delta t} ......(17)

          Pressure Equations

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

          P(i)+12ρ(i)v2(i)+ρ(i)gz(i)=P(i+1)+12ρ(i+1)v2(i+1)+ρ(i+1)gz(i+1)+Pressure  lossesP(i)+\frac12\rho(i)v^2(i)+\rho(i)gz(i)=P(i+1)+\frac12\rho(i+1)v^2(i+1)+\rho(i+1)gz(i+1)+Pressure\; losses


          The following losses are considered for the 3 air flow paths

          • Path 1 : Contraction, Expansion, Grate (contraction & expansion), Bend, Loss through fuel bed in upward direction
          • Path 2 : Contraction, Expansion, Loss through Fuel bed in horizontal direction, Bend. The 2 paths are assumed to meet at this juncture, followed by loss due to buoyancy, being intercepted by path 3 loss due to buoyancy in firing chamber = ρg(1).g.(HgrateLength  of  fuel  bed)\rho_g(1).g.(H_{grate}-Length\; of\; fuel \;bed)
          • Path 3 : Contraction, Expansion, Bend

          The various pressure losses are calculated as follows [3]

            Expansion and Contraction losses at inlet

            Sudden contraction at inlet

            δP=0.5ρ.Vg22\delta P=0.5\frac{\rho.V_g^2}{2} ......(19)

            Sudden expansion at inlet

            δP=(1AinAout)2  Vin2  ρ2\delta P =(1-\frac{A_{in}}{A_{out}})^2\;V_{in}^2\;\frac\rho2​ ......(20)


            δP=0.25ρ.Vg22\delta P=0.25\frac{\rho.V_g^2}{2} ......(21)

              bend correlation as in [3]

              Porous medium loss

              Pressure loss in porous medium, in fuel bed as well as pottery chamber has been estimated using Ergun's equation [1]

              δP=  150  μgas  (1ϕ)2  v2  Lϕ3  Dp2+1.75  ρ  (1ϕ)  v2  Lϕ3  Dp\delta P=\;\frac{150\;\mu_{gas}\;(1-\phi)^2\;v^2\;L}{\phi^3\;D_p^2}+\frac{1.75\;\rho\;(1-\phi)\;v^2\;L}{\phi^3\;D_p} ......(22)

              Dp = 6*Fuel particle volume Fuel particle surface area for fuel bed ......(23)

              For pottery chamber, Effective Particle diameter calculated considering average of spherical diameters of pots in experimental data.

              Loss at grates​

              δP=Kρ.Vg22\delta P=K\frac{\rho.V_g^2}{2} ; K=Ke+Kc ......(24)

              Kc=0.42 (1-AinAout)     if AinAout  <  0.762  \frac{A_{in}}{A_{out}}\;<\;0.76^2\;

                else = (1AinAout)2    (1-\frac{A_{in}}{A_{out}})^2\;\; ​ ......(25)

              Ke=(1AinAout)2  K_e=(1-\frac{A_{in}}{A_{out}})^2\;    ......(26)


              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. 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, implicit scheme.

              First, it was approached to solve the linear equations of temperature first, based on a guess mass rate, then solving the pressure and continuity equations based on the computed temperatures. A non-linear equation of pressures of each control volume interlinked, expressed in mathematically implicit form as a function of input mass flow rates was written and solved numerically by single-variable Newton-Raphson method to obtain correct input mass flow rate. It was observed that the mass flow rate and thus the temperatures and pressures obtained depended significantly on the initial guess. Thus this approach was needed to be modified.

              Now, all the equations have been written in mathematically 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. It is being observed that the multivariable Newton-Raphson Method is failing to converge with the given initial guess or the sequence of equations used. Further trials need to be conducted to choose an optimum initial guess and correct sequence of equations so that the jacobian is diagonally dominant.


              [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.


              We thank The Science Academies for giving us an opportunity to pursue this project.

              We thank our guide Prof. M. R. Ravi for guiding us and giving us a sense of direction at moments when we got stuck while trying to solve the problem.

              We thank Prof. Sangeeta Kohli for organising our site visit and helping us during the need of the hour.


              Appendix 1 : Input in form of MS Excel file


                  Appendix 2 : MATLAB Code

                  1. global all
                  2. all=xlsread('main.xlsx','properties','F3:F46');
                  3. time_firing=all(29);
                  4. dt=all(28);
                  5. N=all(27);
                  6. D=all(24);
                  7. h_grate=all(26);
                  8. n=all(19);
                  9. a=all(20);
                  10. T_amb=all(1);
                  11. air_constant=all(2);
                  12. gas_constant=all(12);
                  13. Rgas=all(5);
                  14. Cp_gas=all(9);
                  15. Cv_gas=Cp_gas-Rgas;
                  16. Cp_air=all(4);
                  17. Cp_pot=all(30);
                  18. rho_pot=all(31);
                  19. CV=all(6);
                  20. g=9.81;
                  21. H=all(21);
                  22. dx=(H-h_grate)/N;
                  23. eps=all(8);
                  24. sig=5.67*10^(-8);
                  25. phi=all(32);
                  26. Hext=all(3);
                  27. k=all(36);
                  28. k_dash=all(37);
                  29. Pr_gas=all(15);
                  30. Asf=all(38)/N;
                  31. width_inlet=all(25);
                  32. h_inlet=all(23);
                  33. h_grate=all(26);
                  34. h_channel=all(22);
                  35. n_holes=all(33);
                  36. A1_ent=width_inlet*h_channel;
                  37. A2_ent=width_inlet*(h_inlet-h_channel);
                  38. A_holes=all(34);
                  39. phi_fuel=all(16);
                  40. Dp=all(39);
                  41. Dp_fuel=all(18);
                  42. h=h_grate-h_channel;
                  43. t=all(40);
                  44. Cp_grate=all(41);
                  45. m_grate=all(42);
                  46. rho_brick=all(43);
                  47. Cp_brick=all(44);
                  48. ke_grate=(1-n*a/(pi/4*D^2))^2;
                  49. if n*a/(pi/4*D^2)<0.76^2
                  50. kc_grate=0.42*(1-n*a/(pi/4*D^2));
                  51. else
                  52. kc_grate=(1-n*a/(pi/4*D^2))^2;
                  53. end
                  54. F=zeros(5*(N+1)+N+5+N+2,1);
                  55. X=zeros(5*(N+1)+N+5+N+2,1);
                  56. X_old=zeros(5*(N+1)+N+5+N+2,1);
                  57. Z=zeros(5*(N+1)+N+5+N+2);
                  58. for i=1:5*(N+1)
                  59. X(i)=400;
                  60. X_old(i)=T_amb;
                  61. end
                  62. for i=5*(N+1)+1:length(X)
                  63. X(i)=1;
                  64. end
                  65. p = xlsread('main.xlsx','fuel mass flow rate');
                  66. time = p(:,1);
                  67. Massf = p(:,2);
                  68. rr=1;
                  69. for i=1:length(time)
                  70. if rr*dt/60 <= time(i)
                  71. Mf=Massf(i);
                  72. break;
                  73. end
                  74. end
                  75. for j=1:10
                  76. Hr=eps*sig*(X(3)+X(1))*(X(3)^2+X(1)^2);
                  77. k_gas=all(13)*X(1)+all(14);
                  78. Hc=8.235*k_gas/D;
                  79. F(1)=X(5*(N+1)+2)*Cp_gas*X(1)+(Hr+Hc)*pi*D*h*X(1)+(Hr+Hc)*(pi/4*D^2-n*a)*X(1)-gas_constant*X_old(1)*pi*D^2*h*Cv_gas/4/dt/X(1)-(Hr+Hc)*(pi/4*D^2-n*a)*X(2)-(Hr+Hc)*pi*D*h*X(3)-X(5*(N+1)+2)*Cp_gas*T_amb+gas_constant*Cv_gas*pi*D^2*h/4/dt-Mf*CV;
                  80. Z(1,1)=X(5*(N+1)+2)*Cp_gas+(Hr+Hc)*pi*D*h+(Hr+Hc)*(pi/4*D^2-n*a)-gas_constant*X_old(1)*pi*D^2*h*Cv_gas/4/dt/X(1)^2;
                  81. Z(1,2)=-(Hr+Hc)*(pi/4*D^2-n*a);
                  82. Z(1,3)=-(Hr+Hc)*pi*D*h;
                  83. Z(1,5*(N+1)+2)=Cp_gas*X(1)-Cp_gas*T_amb;
                  84. Hr=eps*sig*(X(6)+X(8))*(X(8)^2+X(6)^2);
                  85. k_gas=all(13)*X(6)+all(14);
                  86. Hc=8.235*k_gas/D;
                  87. Hsf=6*k_gas/D;
                  88. F(6)=-X(5*(N+1)+2)*Cp_gas*X(1)+X(5*(N+1)+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)-Hsf*Asf*X(7)-(Hr+Hc)*pi*D*dx*X(8)+X(5*(N+1)+3)*Cp_gas*X(11)/2+X(5*(N+1)+2)*Cp_gas*T_amb-X(5*(N+1)+3)*Cp_gas*T_amb+gas_constant*Cv_gas*phi*pi*D^2*dx/4/dt;
                  89. Z(6,1)=-X(5*(N+1)+2)*Cp_gas;
                  90. Z(6,6)=X(5*(N+1)+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;
                  91. Z(6,7)=-Hsf*Asf;
                  92. Z(6,8)=-(Hr+Hc)*pi*D*dx;
                  93. Z(6,11)=X(5*(N+1)+3)*Cp_gas/2;
                  94. Z(6,5*(N+1)+2)=-Cp_gas*X(1)/2+Cp_gas*T_amb;
                  95. Z(6,5*(N+1)+3)=Cp_gas*X(6)/2+Cp_gas*X(11)/2-Cp_gas*T_amb;
                  96. for i=3:N
                  97. Hr=eps*sig*(X(5*(i-1)+3)+X(5*(i-1)+1))*(X(5*(i-1)+3)^2+X(5*(i-1)+1)^2);
                  98. k_gas=all(13)*X(5*(i-1)+1)+all(14);
                  99. Hc=8.235*k_gas/D;
                  100. Hsf=6*k_gas/D;
                  101. F(5*(i-1)+1)=-X(5*(N+1)+i)*Cp_gas*X(5*(i-2)+1)/2+X(5*(N+1)+i+1)*Cp_gas*X(5*(i-1)+1)/2+(Hr+Hc)*pi*D*dx*X(5*(i-1)+1)+Hsf*Asf*X(5*(i-1)+1)-gas_constant*X_old(5*(i-1)+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5*(i-1)+1)-X(5*(N+1)+i)*Cp_gas*X(5*(i-1)+1)/2-Hsf*Asf*X(5*(i-1)+2)-(Hr+Hc)*pi*D*dx*X(5*(i-1)+3)+X(5*(N+1)+i+1)*Cp_gas*X(5*i+1)/2+X(5*(N+1)+i)*Cp_gas*T_amb-X(5*(N+1)+i+1)*Cp_gas*T_amb+gas_constant*Cv_gas*phi*pi*D^2*dx/4/dt;
                  102. Z(5*(i-1)+1,5*(i-2)+1)=-X(5*(N+1)+i)*Cp_gas/2;
                  103. Z(5*(i-1)+1,5*(i-1)+1)=X(5*(N+1)+i+1)*Cp_gas/2+(Hr+Hc)*pi*D*dx+Hsf*Asf-gas_constant*X_old(5*(i-1)+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5*(i-1)+1)^2-X(5*(i-1)+i)*Cp_gas/2;
                  104. Z(5*(i-1)+1,5*(i-1)+2)=-Hsf*Asf;
                  105. Z(5*(i-1)+1,5*(i-1)+3)=-(Hr+Hc)*pi*D*dx;
                  106. Z(5*(i-1)+1,5*i+1)=X(5*(N+1)+i+1)*Cp_gas/2;
                  107. Z(5*(i-1)+1,5*(N+1)+i)=-Cp_gas*X(5*(i-2)+1)/2-Cp_gas*X(5*(i-1)+1)/2+Cp_gas*T_amb;
                  108. Z(5*(i-1)+1,5*(N+1)+i+1)=Cp_gas*X(5*(i-1)+1)/2+Cp_gas*X(5*i+1)/2-Cp_gas*T_amb;
                  109. end
                  110. Hr=eps*sig*(X(5*N+3)+X(5*N+1))*(X(5*N+3)^2+X(5*N+1)^2);
                  111. k_gas=all(13)*X(5*N+1)+all(14);
                  112. Hc=8.235*k_gas/D;
                  113. Hsf=6*k_gas/D;
                  114. F(5*N+1)=-X(5*(N+1)+N+1)*Cp_gas*X(5*(N-1)+1)/2+X(5*(N+1)+N+2)*Cp_gas*X(5*N+1)-X(5*(N+1)+N+1)*Cp_gas*X(5*N+1)/2+(Hr+Hc)*pi*D*dx*X(5*N+1)+Hsf*Asf*X(5*N+1)-gas_constant*X_old(5*N+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5*N+1)-Hsf*Asf*X(5*N+2)-(Hr+Hc)*pi*D*dx*X(5*N+3)+X(5*(N+1)+N+1)*Cp_gas*T_amb-X(5*(N+1)+N+2)*Cp_gas*T_amb+gas_constant*Cv_gas*phi*pi*D^2*dx/4/dt;
                  115. Z(5*N+1,5*(N-1)+1)=-X(5*(N+1)+N+1)*Cp_gas/2;
                  116. Z(5*N+1,5*N+1)=X(5*(N+1)+N+2)*Cp_gas-X(5*(N+1)+N+1)*Cp_gas/2+(Hr+Hc)*pi*D*dx+Hsf*Asf-gas_constant*X_old(5*N+1)*phi*pi*D^2*dx*Cv_gas/4/dt/X(5*N+1)^2;
                  117. Z(5*N+1,5*N+2)=-Hsf*Asf;
                  118. Z(5*N+1,5*N+3)=-(Hr+Hc)*pi*D*dx;
                  119. Z(5*N+1,5*(N+1)+N+1)=-Cp_gas*X(5*(N-1)+1)/2-Cp_gas*X(5*N+1)/2+Cp_gas*T_amb;
                  120. Z(5*N+1,5*(N+1)+N+2)=Cp_gas*X(5*N+1)-Cp_gas*T_amb;
                  121. Hr=eps*sig*(X(1)+X(2))*(X(1)^2+X(2)^2);
                  122. k_gas=all(13)*X(1)+all(14);
                  123. Hc=8.235*k_gas/D;
                  124. F(2)=m_grate*Cp_grate*X(2)/dt-m_grate*Cp_grate*X_old(2)/dt-(Hr+Hc)*(pi/4*D^2-n*a)*X(1)+(Hr+Hc)*(pi/4*D^2-n*a)*X(2);
                  125. Z(2,1)=-(Hr+Hc)*(pi/4*D^2-n*a);
                  126. Z(2,2)=m_grate*Cp_grate/dt+(Hr+Hc)*(pi/4*D^2-n*a);
                  127. for i=2:N+1
                  128. k_gas=all(13)*X(5*(i-1)+1)+all(14);
                  129. Hsf=6*k_gas/D;
                  130. F(5*(i-1)+2)=-Hsf*Asf*X(5*(i-1)+1)+Hsf*Asf*X(5*(i-1)+2)+rho_pot*Cp_pot*(1-phi)*pi*D^2*dx/4/dt*X(5*(i-1)+2)-rho_pot*Cp_pot*(1-phi)*pi*D^2*dx/4/dt*X_old(5*(i-1)+2);
                  131. Z(5*(i-1)+2,5*(i-1)+1)=-Hsf*Asf;
                  132. Z(5*(i-1)+2,5*(i-1)+2)=Hsf*Asf+rho_pot*Cp_pot*(1-phi)*pi*D^2*dx/4/dt;
                  133. end
                  134. for i=1:N+1
                  135. Hr=eps*sig*(X(5*(i-1)+3)+X(5*(i-1)+1))*(X(5*(i-1)+3)^2+X(5*(i-1)+1)^2);
                  136. k_gas=all(13)*X(5*(i-1)+1)+all(14);
                  137. Hc=8.235*k_gas/D;
                  138. F(5*(i-1)+3)=-(Hr+Hc)*X(5*(i-1)+1)+(Hr+Hc)*X(5*(i-1)+3)+X(5*(i-1)+3)/(t/2*(1/k+1/k_dash))-X(5*(i-1)+4)/(t/2*(1/k+1/k_dash))+rho_brick*pi*D*h*t*Cp_brick/dt*(X(5*(i-1)+3)-X_old(5*(i-1)+3));
                  139. Z(5*(i-1)+3,5*(i-1)+1)=-Hr-Hc;
                  140. Z(5*(i-1)+3,5*(i-1)+3)=Hr+Hc+1/(t/2*(1/k+1/k_dash))+rho_brick*pi*D*h*t*Cp_brick/dt;
                  141. Z(5*(i-1)+3,5*(i-1)+4)=1/(t/2*(1/k+1/k_dash));
                  142. F(5*(i-1)+4)=-X(5*(i-1)+3)/(t/2*(1/k+1/k_dash))+2*X(5*(i-1)+4)/(t/2*(1/k+1/k_dash))-X(5*(i-1)+5)/(t/2*(1/k+1/k_dash))+rho_brick*pi*D*h*t*Cp_brick/dt*(X(5*(i-1)+4)-X_old(5*(i-1)+4));
                  143. Z(5*(i-1)+4,5*(i-1)+3)=-1/(t/2*(1/k+1/k_dash));
                  144. Z(5*(i-1)+4,5*(i-1)+4)=2/(t/2*(1/k+1/k_dash))+rho_brick*pi*D*h*t*Cp_brick/dt;
                  145. Z(5*(i-1)+4,5*(i-1)+5)=-1/(t/2*(1/k+1/k_dash));
                  146. F(5*(i-1)+5)=-X(5*(i-1)+4)/(t/2*(1/k+1/k_dash))+X(5*(i-1)+5)/(t/2*(1/k+1/k_dash))+Hext*X(5*(i-1)+5)-Hext*T_amb+rho_brick*pi*D*h*t*Cp_brick/dt*(X(5*(i-1)+5)-X_old(5*(i-1)+5));
                  147. Z(5*(i-1)+5,5*(i-1)+4)=-1/(t/2*(1/k+1/k_dash));
                  148. Z(5*(i-1)+5,5*(i-1)+5)=1/(t/2*(1/k+1/k_dash))+Hext+rho_brick*pi*D*h*t*Cp_brick/dt;
                  149. end
                  150. F(5*(N+1)+1)=X(5*(N+1)+1)-X(5*(N+1)+N+3)-X(5*(N+1)+N+4)-n_holes*X(5*(N+1)+N+5)-Mf;
                  151. Z(5*(N+1)+1,5*(N+1)+1)=1;
                  152. Z(5*(N+1)+1,5*(N+1)+N+3)=-1;
                  153. Z(5*(N+1)+1,5*(N+1)+N+4)=-1;
                  154. Z(5*(N+1)+1,5*(N+1)+N+5)=-n_holes;
                  155. F(5*(N+1)+2)=pi*D^2*h*gas_constant/4/dt/X(1)-X(5*(N+1)+1)+X(5*(N+1)+2)-pi*D^2*gas_constant*h/4/dt/X_old(1);
                  156. Z(5*(N+1)+2,1)=pi*D^2*h*gas_constant/4/dt/X(1)^2;
                  157. Z(5*(N+1)+2,5*(N+1)+1)=-1;
                  158. Z(5*(N+1)+2,5*(N+1)+2)=1;
                  159. for i=3:N+2
                  160. F(5*(N+1)+i)=phi*pi*D^2*dx*gas_constant/4/dt/X(5*(i-2)+1)+X(5*(N+1)+i)-X(5*(N+1)+i-1)-phi*pi*D^2*dx*gas_constant/4/dt/X_old(5*(i-2)+1);
                  161. Z(5*(N+1)+i,5*(i-2)+1)=phi*pi*D^2*dx*gas_constant/4/dt/X(5*(i-2)+1)^2;
                  162. Z(5*(N+1)+i,5*(N+1)+i-1)=-1;
                  163. Z(5*(N+1)+i,5*(N+1)+i)=1;
                  164. end
                  165. mu_gas=all(10)*X(1)+all(11);
                  166. F(5*(N+1)+N+3)=X(5*(N+1)+N+5+1)-X(5*(N+1)+N+5+2)+0.5*T_amb/air_constant*(X(5*(N+1)+N+3)/A1_ent)^2-0.5*(X(1)+X(6))/2/gas_constant*(X(5*(N+1)+2)/(pi/4*D^2))^2-X(5*(N+1)+3)^2*(-0.5*T_amb/air_constant/A1_ent^2+0.5*0.5*T_amb/air_constant/A1_ent^2+(1-A1_ent/(D*h))^2*0.5*T_amb/air_constant/A1_ent^2+(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(X(1)/(pi/4*D^2*gas_constant))^2+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet*X(1)/gas_constant/(pi/4*D^2)^2)-gas_constant/X(1)*g*(h-h_inlet)-(ke_grate+kc_grate)*0.5*X(1)/gas_constant*(X(5*(N+1)+2)/(n*a))^2;
                  167. Z(5*(N+1)+N+3,1)=-0.5/2/gas_constant*(X(5*(N+1)+2)/(pi/4*D^2))^2-X(5*(N+1)+N+3)^2*((ke_grate+kc_grate)*0.5/gas_constant/(n*a)^2+0.25*0.5/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(1/(pi/4*D^2*gas_constant))^2*2*X(1)+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet/gas_constant/(pi/4*D^2)^2)-gas_constant/X(1)^2*g*(h-h_inlet)-(ke_grate+kc_grate)*0.5/gas_constant*(X(5*(N+1)+2)/(n*a))^2;
                  168. Z(5*(N+1)+N+3,6)=-0.5/2/gas_constant*(X(5*(N+1)+2)/(pi/4*D^2))^2;
                  169. Z(5*(N+1)+N+3,5*(N+1)+2)=-2*X(5*(N+1)+2)*(0.5*(X(1)+X(6))/2/gas_constant/(pi/4*D^2)^2-(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2);
                  170. Z(5*(N+1)+N+3,5*(N+1)+N+3)=-2*X(5*(N+1)+N+3)*(-0.5*T_amb/air_constant*(1/A1_ent)^2-0.5*T_amb/air_constant/A1_ent^2+0.5*0.5*T_amb/air_constant/A1_ent^2+(1-A1_ent/(D*h))^2*0.5*T_amb/air_constant/A1_ent^2+(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(X(1)/(pi/4*D^2*gas_constant))^2+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet*X(1)/gas_constant/(pi/4*D^2)^2);
                  171. Z(5*(N+1)+N+3,5*(N+1)+N+5+1)=1;
                  172. Z(5*(N+1)+N+3,5*(N+1)+N+5+2)=-1;
                  173. F(5*(N+1)+N+4)=X(5*(N+1)+N+4)^2*(-0.5*T_amb/air_constant/A2_ent^2+0.5*0.5/air_constant*T_amb/A2_ent^2+(1-A2_ent/(D*h))^2*0.5/A2_ent^2/air_constant*T_amb+0.25*0.5/gas_constant*X(1)/(pi/4*D^2)^2+150*mu_gas*D/2*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*(X(1)/(A2_ent*gas_constant))^2+1.75*D/2*(1-phi_fuel)/phi_fuel^3/Dp_fuel*X(1)/gas_constant/A2_ent^2)-X(5*(N+1)+N+3)^2*(-0.5*T_amb/air_constant/A1_ent^2+0.5*0.5*T_amb/air_constant/A1_ent^2+(1-A1_ent/(D*h))^2*0.5*T_amb/air_constant/A1_ent^2+(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(X(1)/(pi/4*D^2*gas_constant))^2+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet*X(1)/gas_constant/(pi/4*D^2)^2);
                  174. Z(5*(N+1)+N+4,1)=X(5*(N+1)+N+4)^2*(0.25*0.5/gas_constant/(pi/4*D^2)^2+150*mu_gas*D/2*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*(1/(A2_ent*gas_constant))^2*2*X(1)+1.75*D/2*(1-phi_fuel)/phi_fuel^3/Dp_fuel/gas_constant/A2_ent^2)-X(5*(N+1)+N+3)^2*((ke_grate+kc_grate)*0.5/gas_constant/(n*a)^2+0.25*0.5/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(1/(pi/4*D^2*gas_constant))^2*2*X(1)+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet/gas_constant/(pi/4*D^2)^2);
                  175. Z(5*(N+1)+N+4,5*(N+1)+N+3)=-2*X(5+6*N+N+3)*(-0.5*T_amb/air_constant/A1_ent^2+0.5*0.5*T_amb/air_constant/A1_ent^2+(1-A1_ent/(D*h))^2*0.5*T_amb/air_constant/A1_ent^2+(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(X(1)/(pi/4*D^2*gas_constant))^2+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet*X(1)/gas_constant/(pi/4*D^2)^2);
                  176. Z(5*(N+1)+N+4,5*(N+1)+N+4)=2*X(5*(N+1)+N+4)*(-0.5*T_amb/air_constant/A2_ent^2+0.5*0.5/air_constant*T_amb/A2_ent^2+(1-A2_ent/(D*h))^2*0.5/A2_ent^2/air_constant*T_amb+0.25*0.5/gas_constant*X(1)/(pi/4*D^2)^2+150*mu_gas*D/2*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*(X(1)/(A2_ent*gas_constant))^2+1.75*D/2*(1-phi_fuel)/phi_fuel^3/Dp_fuel*X(1)/gas_constant/A2_ent^2);
                  177. mu_gas=all(10)*X(1)+all(11);
                  178. F(5*(N+1)+N+5)=-X(5*(N+1)+3)^2*(-0.5*T_amb/air_constant/A1_ent^2+0.5*0.5*T_amb/air_constant/A1_ent^2+(1-A1_ent/(D*h))^2*0.5*T_amb/air_constant/A1_ent^2+(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(X(1)/(pi/4*D^2*gas_constant))^2+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet*X(1)/gas_constant/(pi/4*D^2)^2)-gas_constant/X(1)*g*(h-h_inlet)+air_constant/T_amb*g*h+X(5*(N+1)+5)^2*(-0.5*T_amb/air_constant/A_holes^2+0.5*0.5*T_amb/air_constant/A_holes^2+(1-A_holes/(D*h))^2*0.5*T_amb/air_constant/A_holes^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2);
                  179. Z(5*(N+1)+N+5,1)=-X(5*(N+1)+3)^2*((ke_grate+kc_grate)*0.5/gas_constant/(n*a)^2+0.25*0.5/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(1/(pi/4*D^2*gas_constant))^2*2*X(1)+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet/gas_constant/(pi/4*D^2)^2)-gas_constant/X(1)^2*g*(h-h_inlet)+X(5*(N+1)+5)^2*0.25*0.5/gas_constant/(pi/4*D^2)^2;
                  180. Z(5*(N+1)+N+5,5*(N+1)+N+3)=-2*X(5+6*N+N+3)*(-0.5*T_amb/air_constant/A1_ent^2+0.5*0.5*T_amb/air_constant/A1_ent^2+(1-A1_ent/(D*h))^2*0.5*T_amb/air_constant/A1_ent^2+(ke_grate+kc_grate)*0.5*X(1)/gas_constant/(n*a)^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2+150*mu_gas*(1-phi_fuel)^2/phi_fuel^3/Dp_fuel^2*h_inlet*(X(1)/(pi/4*D^2*gas_constant))^2+1.75*(1-phi_fuel)/phi_fuel^3/Dp_fuel*h_inlet*X(1)/gas_constant/(pi/4*D^2)^2);
                  181. Z(5*(N+1)+N+5,5*(N+1)+N+5)=2*X(5+6*N+N+5)*(-0.5*T_amb/air_constant/A_holes^2+0.5*0.5*T_amb/air_constant/A_holes^2+(1-A_holes/(D*h))^2*0.5*T_amb/air_constant/A_holes^2+0.25*0.5*X(1)/gas_constant/(pi/4*D^2)^2);
                  182. F(5*(N+1)+N+5+1)=X(5*(N+1)+N+5+1)-air_constant/T_amb*g*H;
                  183. Z(5*(N+1)+N+5+1,5*(N+1)+N+5+1)=1;
                  184. for i=2:N+1
                  185. mu_gas=all(10)*X(5*(i-1)+1)+all(11);
                  186. F(5*(N+1)+N+5+i)=X(5*(N+1)+N+5+i)-X(5*(N+1)+N+5+i+1)+0.5*(X(5*(i-2)+1)+X(5*(i-1)+1))/2/gas_constant*(X(5*(N+1)+i)/(phi*pi/4*D^2))^2-0.5*(X(5*(i-1)+1)+X(5*i+1))/2/gas_constant*(X(5*(N+1)+i+1)/(phi*pi/4*D^2))^2-g*dx*gas_constant/X(5*(i-1)+1)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(5*(i-1)+1)*(X(5*(N+1)+i)+X(5*(N+1)+i+1))/(2*gas_constant*pi/4*D^2))^2-1.75*(1-phi)/phi*3/Dp*dx*X(5*(i-1)+1)/gas_constant*((X(5*(N+1)+i)+X(5*(N+1)+i+1))/(2*pi/4*D^2))^2;
                  187. Z(5*(N+1)+N+5+i,5*(i-2)+1)=0.5/2/gas_constant*(X(5*(N+1)+i)/(phi*pi/4*D^2))^2;
                  188. Z(5*(N+1)+N+5+i,5*(i-1)+1)=0.5/2/gas_constant*(X(5*(N+1)+i)/(phi*pi/4*D^2))^2-0.5/2/gas_constant*(X(5*(N+1)+i+1)/(phi*pi/4*D^2))^2-g*dx*gas_constant/X(5*(i-1)+1)^2-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*((X(5*(N+1)+i)+X(5*(N+1)+i+1))/(2*gas_constant*pi/4*D^2))^2*2*X(5*(i-1)+1)-1.75*(1-phi)/phi*3/Dp*dx/gas_constant*((X(5*(N+1)+i)+X(5*(N+1)+i+1))/(2*pi/4*D^2))^2;
                  189. Z(5*(N+1)+N+5+i,5*i+1)=-0.5/2/gas_constant*(X(5*(N+1)+i+1)/(phi*pi/4*D^2))^2;
                  190. Z(5*(N+1)+N+5+i,5*(N+1)+i)=0.5*(X(5*(i-2)+1)+X(5*(i-1)+1))/2/gas_constant*(1/(phi*pi/4*D^2))^2*2*X(5*(N+1)+i)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(5*(i-1)+1)/(2*gas_constant*pi/4*D^2))^2*2*(X(5*(N+1)+i)+X(5*(N+1)+i+1))-1.75*(1-phi)/phi*3/Dp*dx*X(5*(i-1)+1)/gas_constant*(1/(2*pi/4*D^2))^2*2*(X(5*(N+1)+i)+X(5*(N+1)+i+1));
                  191. Z(5*(N+1)+N+5+i,5*(N+1)+i+1)=-0.5*(X(5*(i-1)+1)+X(5*i+1))/2/gas_constant*(1/(phi*pi/4*D^2))^2*2*X(5*(N+1)+i+1)-150*mu_gas*(1-phi)^2*dx/phi^3/Dp^2*(X(5*(i-1)+1)/(2*gas_constant*pi/4*D^2))^2*2*(X(5*(N+1)+i)+X(5*(N+1)+i+1))-1.75*(1-phi)/phi*3/Dp*dx*X(5*(i-1)+1)/gas_constant*(1/(2*pi/4*D^2))^2*2*(X(5*(N+1)+i)+X(5*(N+1)+i+1));
                  192. Z(5*(N+1)+N+5+i,5*(N+1)+N+5+i)=1;
                  193. Z(5*(N+1)+N+5+i,5*(N+1)+N+5+i+1)=-1;
                  194. end
                  195. F(5*(N+1)+N+5+N+2)=X(5*(N+1)+N+5+N+2);
                  196. Z(5*(N+1)+N+5+N+2,5*(N+1)+N+5+N+2)=1;
                  197. X=X-Z\F;
                  198. X'
                  199. end

                  Written, reviewed, revised, proofed and published with