Loading...

Summer Research Fellowship Programme of India's Science Academies

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  
 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)
 ΦPorosity 
m Mass floe rate (kg/s) 
 kThermal conductivity (W/m/K) 
 ttime (s) 
 KCo-efficient of Pressure loss  
Subscripts  
 inInlet 
 out Outlet
 g / gasFlue gas 
 refReference 
 FFuel 
 airAir 
 amb / a Ambient 
 iNo of control volume 
  Superscripts
 tAt time t 
 t+ΔtAt time t+Δt 

 

INTRODUCTION

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.

Image_2.png
    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 = Atmospheric PressureRgas/air 

    Density of Gas/Air = Gas/Air constantTemperature 

    Concepts

    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) = 0.5 V2 ρ2  ...(1)

    Inlet_1.png
      Expansion and Contraction losses at inlet

      Exit loss (Expansion) = (1-AinAout)2 Vin2 ρ2    ...(2)

      Pressure loss due to bend = 14 V2 ρ2      ...(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

      P_Fuel = 150 μgas (1-ϕ)2 u2 (hgrate-hchannel)ϕ3 Dp2+1.75 ρ (1-ϕ) u2 (hgrate-hchannel)ϕ3 Dp ...(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.

      Dp = 6*Fuel particle volume Fuel particle surface area  ...(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:

      ρgt Cvair (ρgt+t-ρgt) Vt=mt Cpair (Tint+t-Tref) -mout Cpgas (Toutt+t-Tref) +mf CV +(Hr+Hc) Awall(Tw1-Tg)t+t  ...(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.jpg
        Grate in Pottery Furnance at Khurrampur 

        Pressure loss due to grate contraction = 0.42 (1-AopenD (hgrate-hchannel)) vgrate2 ρc2  ...(7)

          if AopenAF < 0.762  

        else = (1-AopenAF)2  Vgrate2 ρc2 

        Pressure loss due to grate expansion = (1-AopenAF)2 Vgrate2 ρc2    ...(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;

        ϕ = 1-Total volume of potsVolume of pottery chamber ...(9) 

        The height of each control volume is given by dx = H-hgrateN  ...(10)

        The first law equation for pottery in pottery chamber is;

        ρpot Vpot Cppot (TPit+t-TPitt) = Hsf Asf (Tgi-TPi)  ...(11)

        The energy equation for gas present in control volume is;

        ρgit Vg Cvg (Tgit+t-Tgitt) = minit Cp (Ti+Ti-12-Tamb)t+t - moutit Cp (Ti+Ti-12-Tamb)t+t +Hsf Asf (Tpi-Tgi) + (Hr+Hc) Awall (Twi-Tgi)  ...(12)

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

        as used by Oza et al.[1]

        Pdrop = 150 μgas (1-ϕ)2 v2 dxϕ3 Dp + 1.75 (1-ϕ) ρ v2 dxϕ3 Dp ​ ...(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.jpeg
          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.

          wall image_1.png
            Electrical anology of thermal resistance of rat trap bond construction

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

            K  Awall  (Tw1  Tw2)t  =  (Hr+Hc)  Awall  (TgTw1)\frac{K\;A_{wall\;}(T_{w1}\;-T_{w2})}t\;=\;(H_r+H_c)\;A_{wall}\;(T_g-T_{w1})  ...(14)

            K Awall(Tw1-Tw2)t=K' Awall(Tw2-Tw3)t ...(15)

               K'Awall(Tw2-Tw3)t=KAwall(Tw3-Tw4)t  ...(16)

            K Awall(Tw3-Tw4)t=Hext Awall(Tw4-Tamb)   ...(17)

            Pressure Evaluation

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

            Pi + 12ρi vi2 + ρi g zi = Pi+1 + 12 ρi+1 vi+12 +ρi+1 g zi+1 + Pressure losses ..(18) 

            Determination of heat transfer co-efficient

            H = Nu kgasD  ...(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 AND RECOMMENDATIONS

            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

            all=xlsread('main.xlsx','properties','F3:F42');

            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