A Study of a Genetic Algorithm for Polydisperse Spray Flames
TTechnion - Israel Institute of Technology
Faculty of Aerospace EngineeringA Study of a Genetic Algorithm for PolydisperseSpray Flames
Final Project Report Towards M.Eng in Aerospace EngineeringAdvisor : Prof. Barry GreenbergSubmitted by : Daniel Engelsman a r X i v : . [ c s . N E ] A ug ontents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
The Optimization Model 36
References 57Appendices 60
Appendix A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Appendix B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62Appendix C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 bstract
Modern technological advancements constantly push forward the human-machineinteraction. Nowadays, finding an application whose algorithm does not utilize a MachineLearning (ML) methods, is quite rare. The reason for that is their capability of solvingabstract problems that were so far not even expressible, and thus remained enigmatic.Evolutionary Algorithms (EA) are an ML subclass inspired by the process of naturalselection - ”
Survival of the Fittest ”, as stated by the Darwinian Theory of Evolution. Themost notable algorithm in that class is the Genetic Algorithm (GA) - a powerful heuristictool which enables the generation of a high-quality solutions to optimization problems. Inrecent decades the algorithm underwent remarkable improvement, which adapted it into awide range of engineering problems, by heuristically searching for the optimal solution.Despite being well-defined, many engineering problems may suffer from heavy analyticalentanglement when approaching the derivation process, as required in classic optimizationmethods. Therefore, the main motivation here, is to work around that obstacle.In this piece of work, I would like to harness the GA capabilities to examine optimalitywith respect to a unique combustion problem, in a way that was never performed before.To be more precise, I would like to utilize it to answer the question : ”What form of aninitial droplet size distribution (iDSD) will guarantee an optimal flame ?”To answer this question, I will first provide a general introduction to the GA method, thendevelop the combustion model, and eventually merge both into an optimization problem.
Acknowledgments
First and foremost, I would like to thank Ms. Debbie Warril.Without you I wouldn’t have made it this far. As simple as that. Thank you.I would also like to pay my gratitude to the project supervisor - Professor Barry Greenberg.For the professional academic guidance and kind hearted support all along the way.Last but not least, I would like to thank my beloved parents and family, just for being whoyou are. I believe we are lucky people in this world.1 omenclature
Latin symbols B Dimensionless reaction [ − ] B i,i +1 i -th integral coefficient [1 / sec] C i i -th integral coefficient [1 / sec] C P Specific heat at constant pressure [ J/ ( Kg · K )] c Half inner channel size (normalized) [ − ] d Droplet diameter [ m ] D g Mass diffusion coefficient [ m /sec ] E Evaporation rate [ m /sec ] F Drag acceleration [ m/sec ] K Thermal diffusion coefficient [ m /sec ]¯ E Normalized evaporation rate [ m ] L Half inner channel size [ m ] m Mass fraction [ − ]˙ m Mass flux [
Kg/ ( m · sec )] N Number of sections [ − ] n Droplet size probability function [Droplets] q Heat flux [ J/ ( m · sec )] R Half external channel size [ m ]˜ R Droplet volumetric change rate [ m /sec ] S Mass source / sink element [1 /sec ]¯ S Heat source / sink element [
K/sec ] t Time [ sec ]˜ T Temperature [ K ] T Dimensionless temperature [ − ] T max Normalized maximum tip flame temperature [ − ] U g Characteristic flow velocity [ m/sec ] V Oxidizer’s initial mass fraction [ − ] v Droplet volume [ m ] x Cartesian coordinate perpendicular to the flow [ m ] y Cartesian coordinate parallel to the flow [ m ]2 reek symbols α Coefficient of the integral property [ − ] γ Normalized mass fraction [ − ] γ d Normalized mass fraction (droplets) [ − ] γ F Normalized mass fraction (fuel) [ − ] γ O Normalized mass fraction (oxidizer) [ − ] γ T Normalized mass fraction (temperature @ S-Z) [ − ] δ Normalized initial mass fraction (fuel) [ − ]∆ Damk¨ohler number for evaporation [ − ]∆ i Integral coefficient of i -th section [ − ] η Normalized coordinate parallel to the flow [ − ] η max Normalized maximum flame height [ − ]Λ Normalized latent heat [ − ] ξ Normalized coordinate perpendicular to the flow [ − ] ρ Density [
Kg/m ] ν Stoichiometric coefficient [ − ] Shortcuts
Da Damk¨ohler number for reaction [ − ]DoF Degree of FreedomEA Evolutionary AlgorithmsF/O Fuel / Oxidizer ratio [ − ]GA Genetic AlgorithmiDSD initial Droplet Size DistributionLe Lewis number [ − ]LHS Left Hand SideML Machine LearningPe P´eclet number [ − ]rev. Reversal (point)RHS Right Hand SideSMD Sauter Mean Diameter [ m ]S-Z Schwab-Zeldovich transformationWe Webber number [ − ]3 Introduction
The GA consists of 3 mechanisms that reflect the natural selection process (survival of thefittest), where the fittest is selected for producing the next generation’s offspring [1].
Natural selection - The individual’s / parent probability to be selected to produce theoffspring of the next generation. Its genetic variation determines its survival chances, andthus points which chromosomes are to be preserved and multiplied between generations. Innature, it’s caused by forced competitive interaction between different populations andindividuals, that rewards the fittest among them (best intelligence, physique etc.).
Crossover - The recombination of the genetic information of two parents (in nature knownas mating). Better individuals will participate in the production of the next generation,such that the last generation will hold the best former genetic qualities. Each singlecrossover (out of hundreds) is subjected to a random partitioning, across every generation.
Mutation - Random genetic alteration during the recombination process. The bigger thesample space is, the wider the diversity becomes, and so do the chances for new improvedfeatures. A positive / negative feature caused by a mutation will be reflected in theindividual’s fitness quality, resulting in higher / lower chances to transfer its genes.4
Chromosome - the individual’s set of properties that represent a candidate solution.As such, the chromosome is subjected to modifications at every generation. • Population - set of n -random candidate solutions, given an optimization problem.The above mechanisms are applied in a loop, where each iteration (=generation) the fittestindividuals are extracted and go through genetic recombination. This process continuesuntil a termination criterion is met, and the best candidate solution is received. From leftto right is the evolution process from 1st random initialization until the 78th generation :Consider the above fitness function to be an (cid:96) -norm : f GA ( x ti ) = (cid:107) x ti (cid:107) = (cid:80) j | x i,j | .The optimal candidate solution at generation t means : arg max x t f GA ( X t ) .Note the fitness (score) evolution from t = 0 until t = 78 across all of the population.As seen, the chromosomes are consisted of atomic sequences named Genes. Assuming anarbitrary chromosome that’s composed of n genes we get x i ∈ R n , whereas the differentchromosomes can be seen as a points in the R n space, whose genes are coordinates.The GA then, is responsible for finding the best performance among all candidates, as theyare measured by a fitness function, or equivalently as they are projected onto a metric axis.The optimal solution (minimum / maximum) is actually the candidate whose scoreperforms best in the R n +1 , namely the closest to the global extremum (see Appendix A).Classic optimization techniques utilize a closed form objective (=fitness) functions that areconveniently differentiable. By calculating the roots of their first and second derivatives,one can extract solutions in the form of minima, maxima or a saddle point [2].However, in complex analytical cases (as in ours), one would rather work around thattiring derivation process which can impose significant challenges, and implement instead5olution oriented heuristic methods. For illustration, consider the following 3D function : f ( x, y, z ) = sin ( x − x ) x − x · sin ( x − y ) x − y · (cid:16) − z + z (cid:17) ; ( x , y , z ) = (10 , , { ≤ x ≤ , ≤ y ≤ , − ≤ z ≤ } The GA manages to find a global solution, being slightly dependent on the mesh resolution.The individual’s score at a given generation is measured by its performance to a desiredfitness function. That function acts as a comparison measure, where the best candidate isthe one whose score is optimal. That optimality can be either minimum or maximum,depending on the problem’s nature - concave, convex or non-convex (see Appendix B).One of the main advantages of the GA is its indifference to the internal workings of thefitness function, namely it refers to it as a ”black box”, evaluating different points f GA ( x )due to its coding methodology. As introduced above, and detailed thoroughly here [3] , the GA are powerful heuristic searchmethods, successfully used to find optimal or near-optimal solutions in many complexdesign spaces. Early implementations of the GA in context of combustion were made back6n the 90’s, as Runhe Huang (1995) [4] showed an implementation on a combustion controlproblem. Instead of learning a control action for every point encountered, a GA was usedto learn control actions for a set of limited number of prototype states, and afterwardsapplying nearest neighbour matching to extract the optimal rule.Danielson et al. (1998) [5], presented the GABSys (GA Bond Graph System), anoptimization tool that utilized a 2-stroke combustion engine model :They parametrized several related factors (geometric, thermodynamic etc.), and expressedthe objective (fitness) function in terms of fuel consumption, power etc. (along cycles),such that eventually a performance space could be spanned and satisfy local extrema :”Although there is no guarantee that the result is optimal, the resulting engine is still veryimpressive considering that after fewer than 3,000,000 designs from a domain of 1.329 · ,the GA selects a constructable fuel-efficient engine design”.Polifke et al. (1998) [6] used the GA for combustion reactive mechanism to carry out thesubtle optimization process, with a minimum human effort. Harris et al. (1999) [7] usedthe GA for determining the optimal reaction rate parameters of the O/F mixture.Vossoughi & Rezazadeh (2005) [8] introduced a multi-objective GA for an engine controlunit, where the objective functions were tailored to the calibration parameters in sought ofoptimal configuration. Quite similarly, Rose et al. (2009) [9] implemented the GA on agas-exchange system of combustion engine, producing a significantly higher power outputthan was achieved through a basic manual optimization procedure.Shtauber & Greenberg (2010) [10] conducted a wide study of polydisperse spray diffusionflames. By analytical and numerical investigation, they have shown the iDSD influence7pon the flame properties and its sensitivity to extinction. It is worth mentioning that thisproject is considerably a continued work on Shtauber’s thesis, but focuses primarily onfinding the optimal flame properties, using the GA.Sikalo et al. (2015) [11] described an automatic method for the optimization of reactionrate constants of reduced reaction mechanisms. Based on GA, the technique aimed atfinding new reaction rate coefficients that minimize the error introduced by the precedingreduction process. The error was defined by an objective function that covers regions ofinterest where the reduced mechanism may deviate from the original mechanism.Kaplan et al. (2015) [12] presented a general approach for developing an automatedprocedure to determine optimal reaction parameters for a simplified model to simulateflame acceleration and deflagration-to-detonation (DDT) in a methane-air mixture. Thelaminar flame profile was computed using reaction parameters in a 1D Navier-Stokes code,and matched the profile obtained by a detailed chemical reaction mechanism.Pan et al. (2018) [13] utilized the GA in a boiler combustion control system, to optimizethe bias coefficients that maintain the excess air ratio at the optimal combustion intervalunder variable load conditions (the blue plot) :Liu et al. (2019) [14] presented a GA implementation on dual diesel/natural injectionparameters, where the indicated specific fuel consumption (NO x and CH ) emissions areselected as the optimization objectives. Similarly, Zhao et al. (2019) [15] focused onreducing unburned carbon by optimizing operating parameters via a novel high-efficientGA, which was experimentally validated. 8 ummary Many of the above researchers had no relation to the combustion physical aspects. Instead,they used it as a convenient optimization framework, for its convenient modelability andbeing experimentally validable. The main implementations were :( ◦ ) Holistic analysis of an engineering systems (mainly combustion configuration) e.ginternal combustion or spark ignition engines, heat exchangers, chemical reactors etc.( ◦ ) State space representation of control systems and attempt to optimize a desiredvariable (power, emission, efficiency, fuel consumption etc.)Applying the GA on big frameworks may provide high-level understanding of the engineefficiency, emission aspects or different dynamic profiles. However, smaller focus areas thatactually comprise the problem’s inner core, might be lost. More precisely, they are not evenexpressed in the cost function, and are thus overridden by macroscopic interests. In the absence of any research that applied the GA with a well-defined combustion model,I aim to focus on smaller scopes of interest, primarily on the flame characteristics.At first, by being able to express the temperature field and investigating its reactions to awide range of parameters. To that end, factors like the maximum flame height and themaximum tip flame temperature will be serve as indicators.Afterwards, I would like to gain control on the GA model by being able to executeoptimization schemes in a growing complexity order, by either extreme chemical scenariosor by maximizing the degree of freedom (DoF).Finally, when full integration is achieved, the GA will be harnessed for the sake ofoptimization scenarios in order to shed light on the principal factors that may optimize thecurrent combustion model. These steps comprise the current piece of work, in a way thatwas not conducted before, especially on the seam between GA and combustion. I aim toinnovate by investigating the cause and effects evoked as a result of the optimality, andvalidate them according to the literature. 9
The Combustion Model
In this section I will develop the governing equations describing the mathematical model ofthe polydisperse spray diffusion flame. That is by presenting the underpinningassumptions, equations, normalization, boundary conditions and full solution. Afterwards Iwill validate them with a set of results, which will be followed by discussions.
The big picture
Based on the classic flame model of Burke & Schumann (1928) [16], the F/O interaction isseparated by a steady state diffusion flame, in a laminar parallel co-flow :This model was further elaborated by Greenberg (1989) [17] by assuming that the liquidfuel droplets were homogeneously suspended in an inert gas stream.Using Tambour (1985) [18] sectional approach for describing the spray polydispersity, inaddition to the droplets evaporation rate, it is necessary to refer to the droplets differentsizes, as a result of the iDSD and the evaporation rate mechanism ( d law ).10 .1 The Spray Equation Tambour’s sectional approach assumes discrete distribution of the droplets, whose mostfundamental size of is called Monomer. Its size will dictate the field’s resolution, and itexpresses the spray as a probabilistic function of the droplets : n = n ( t, x, y, v ) → continuous n ( t, x, y, v ) · dx · dy · dv (2.1)Each droplet of the spray is indicated by j ∈ N that denotes the number of monomerscarried inside (e.g number of molecules in a droplet) : dn j dt = − E j n j + E j +1 n j +1 , j = 1 , , ... (2.2) E j indicates the evaporation rate of the single j -th droplet, such that it is not dependenton the environment’s temperature, but on the droplet’s size, based on d law . We’ll dividethe spray into N sections and define an integral property (IP) for the i -th section : Q i ( t, x, y, v ) (cid:44) (cid:90) v Hi v Li αv ˆ γ n ( t, x, y, v ) dv (2.3)The volume defines the i -th section and is bounded within v ∈ [ v L i , v H i ) . The ˆ γ coefficientdefines the property whereas ˆ γ = { , , } based on the IP’s dependence on the number ofdroplets, their volume or their surface area. Then α is set into a desirable property (e.gdensity) such that the integration would yield the i -th mass section. Further he showed : dQ i dt = − C i Q i + B i,i +1 Q i +1 , i = 1 , . . . , N (2.4)Where the general integral coefficients are : B i,i +1 = (cid:16) v H i v L i +1 (cid:17) ˆ γ I E ( v H i ) v H i +1 − v L i +1 , B N,N +1 = 0 (2.5) C i = (cid:16) v H i − v L i (cid:17) ˆ γ II E ( v L i ) v H i − v L i + 1 v H i − v L i (cid:90) v Hi v Li v ˆ γ E ( v ) dv ˆ γ III (2.6)I. Expresses Q i ’s growth by droplets addition from i + 1 → i section.II. Expresses Q i ’s diminution by droplets downgrade from the i -th section.III. Expresses Q i ’s diminution by droplets evaporation in the i -th section.11sing the general conservation equation given by Williams (1985) [19] : ∂n∂t + ∂∂v ( ˜ R (cid:124)(cid:123)(cid:122)(cid:125) ∂v∂t n ) + ∇ · ( U d n ) + ∇ U d · ( F (cid:124)(cid:123)(cid:122)(cid:125) Drag accel. n ) = Γ (2.7)Where Γ is a source of droplets resulted from the collision rate. However, we should recallsome of the prior assumptions made in the thesis [10] : ∇ U d · ( F n ) (cid:124) (cid:123)(cid:122) (cid:125) v drop = v gas = 0 ∂n∂t (cid:124)(cid:123)(cid:122)(cid:125) steady state = 0 Γ (cid:124)(cid:123)(cid:122)(cid:125) negligible collision = 0 U d = U g (cid:124) (cid:123)(cid:122) (cid:125) negligible collision = 0 (2.8)Using the sectional approach based on these assumptions, we get the following equation : ∇ · (cid:0) U g Q i (cid:1) = − C i Q i + B i,i +1 Q i +1 , i = 1 , . . . , N (2.9)Plugging (ˆ γ = 1 , α = ρ d /ρ T ot ) in (2.3) we get the mass fraction equation ( Q i = m d i ) U g ∂m d i ∂y = − C i m d i + B i,i +1 m d i +1 , i = 1 , . . . , N (2.10)The continuous sectioning (cid:32) d H i − → d L i d H i → d L i +1 (cid:33) of the integral coefficients (2.5, 2.6) becomes : B i,i +1 = 32 E (cid:34) d L i +1 d H i +1 − d L i +1 (cid:35) , i = 1 , . . . , N (2.11) C i = 32 E (cid:20) d H i − d L i d H i − d L i (cid:21) , i = 1 , . . . , N (2.12)12 .1.1 Boundary conditions The droplets are described in a 1st order equation whose BC are the liquid fuel massfraction at the channel’s exit. Away from the nozzle (”far field”) the polydispersity isassumed to be homogeneous: m d i = m Tot, fuel · δ i , ≤ x ≤ L , L < x ≤ R , i = 1 , . . . , N (2.13)Applying the following normalizations :( ξ, η, c ) (cid:44) (cid:18) xR , yD g U g R , LR (cid:19) (2.14)( γ d i ) (cid:44) ( m d i /m Tot, fuel ) (2.15)( ψ i , ∆ i ) = R D g ( B i , C i ) (2.16)And we get the dimensionless spray equation : ∂γ d i ∂η = − ∆ i γ d i + ψ i γ d i +1 (2.17)Whereas the integral coefficients are defined as :¯ E = E (cid:18) R D g (cid:19) ∆ i ≤ i ≤ N = 3 E (cid:18) d H i − d L i d H i − d L i (cid:19) ψ ii ≤ i ≤ N − = 3 ¯ E (cid:32) d i +1 d H i +1 − d L i +1 (cid:33) (2.18)And the dimensionless BC are : γ d i = m d i m T ot,fuel = δ i , ≤ ξ ≤ c , c < ξ ≤ , i = 1 , . . . , N (2.19)13 .1.2 Analytical solution The spray equation (2.17) contains the coefficients (2.18) and the compatible BC (2.19).Since subsequent equations are mutually dependent, we’ll propose an iterative approach : γ d j = N (cid:88) i = j Ω ij e − ∆ i η ; ∂γ d j ∂η = − N (cid:88) i = j ∆ i Ω ij e − ∆ i η ; Ω ij ( i>j ) = ψ j ∆ j − ∆ i Ω i,j +1 (2.20)Ω ij is an influence coefficient. By plugging inside the dimensionless spray equation (2.17) : − N (cid:88) i = j ∆ i Ω ij e − ∆ i η = − ∆ j N (cid:88) i = j Ω ij e − ∆ i η + ψ j N (cid:88) i = j +1 Ω i,j +1 e − ∆ i η (2.21)Applying solution on the BC : γ d j (0) = N (cid:88) i = j Ω ij = δ j (2.22)Ω jj = N (cid:88) i = j δ i − N (cid:88) i = j +1 Ω ij = γ d j (0) − N (cid:88) i = j +1 Ω ij , ( j < N ) (2.23)Ω NN = N (cid:88) i = j δ i = γ d N (0) (2.24)The above solution is approximated as continuous, despite the mass fraction beingdiscontinuous as droplets may join or leave the i -th section (=discrete phenomenon).This approximation allows the analytical solution as it presumes that the average sprayinjection may contain up to hundred thousands of droplets. That way, joining of a singledroplet from larger section, is negligible, as it is smaller by several order of magnitudes.To sum up the droplets solution : γ d j = N (cid:88) i = j Ω ij e − ∆ i η (2.25)Ω ij ( i>j ) = ψ j ∆ j − ∆ i Ω i,j +1 (2.26)Ω jj = γ d j (0) − N (cid:88) i = j +1 Ω ij = δ j − N (cid:88) i = j +1 Ω ij (2.27)14 .2 The Gaseous Phase Equation According to Fick’s 1st law for diffusion, the mass flux is linear with its spatial gradient :˙ m A = − D A ∇ ρ A (2.28)Combining it with the continuity equation provides its 2nd law, AKA the diffusionequation and refers to the concentration change as a function of time. Using that, thegaseous phase equation in terms of mass fraction will be expressed as : U A ∂m A ∂y = D A (cid:18) ∂ m A ∂x + ∂ m A ∂y (cid:19) + S A (2.29)This equation is valid for both gaseous fuel and oxidizer where S A is a sink / source ofelement A in terms of rate. Assuming equal diffusion coefficients and equal velocities forboth gaseous fuel and oxidizer :Fuel : U g ∂m g, fuel ∂y = D g (cid:16) ∂ m g, fuel ∂x + ∂ m g, fuel ∂y (cid:17) + S g, fuel-reac + S d, fuel (2.30)Oxidizer : U g ∂m O ∂y = D g (cid:16) ∂ m O ∂x + ∂ m O ∂y (cid:17) + S O − reac (2.31)Using Schwab-Zeldovich (S-Z) transformation will help us uniting both equations : m (cid:44) m g, fuel − m O /ν (2.32)Where the stoichiometric ratio ( ν ) is originated at :Fuel + ν Oxygen → Heat + ProductsAnd the reactant elements are active only at the reaction zone such that : S g, fuel-reac = S O , fuel-reac /ν (2.33)So by implementing the transform : (2.30) - (2.31) / v we can get rid of the reactants andthe nonlinear reaction rate term : U g ∂m∂y = D g (cid:16) ∂ m∂x + ∂ m∂y (cid:17) + S d (2.34) S d is the source element expressing the droplets evaporation rate and their contribution tothe gaseous phase. Using this we can find the flame shape by solving for m = 0 .15 .2.1 Boundary conditions The gaseous phase equation requires 2 BC for each axis. The BC in the channels’ exitcontain diffusive flux elements resulting from the mass fraction gradient : y = 0 : m g, fuel − D g U g, fuel ∂m g, fuel ∂y = m Tot, fuel (cid:16) − (cid:80) Ni =1 δ i (cid:17) , ≤ x ≤ L , L ≤ x ≤ R (2.35)Similarly, we get in the oxidizer equation : y = 0 : m O − D g U , fuel ∂m O ∂y = , ≤ x ≤ Lm O ( y = 0) , L ≤ x ≤ R (2.36)And by using S-Z transformation : m − D g U g, fuel ∂m∂y = m Tot, fuel (cid:16) − (cid:80) Ni =1 δ i (cid:17) , ≤ x ≤ L − m O ( y = 0) /ν , L ≤ x ≤ R (2.37)Applying the following assumptions : ∂m∂x (cid:124)(cid:123)(cid:122)(cid:125) Symmetry (cid:12)(cid:12)(cid:12)(cid:12) y ≥ x =0 = 0 ∂m∂x (cid:124)(cid:123)(cid:122)(cid:125) Impenetrable channel (cid:12)(cid:12)(cid:12)(cid:12) y ≥ x = R = 0 ∂m∂y (cid:124)(cid:123)(cid:122)(cid:125) Thermodynmic equilibrium (cid:12)(cid:12)(cid:12)(cid:12) y →∞ ≤ x ≤ R = 0 (2.38)Axes normalization is similar as before (2 . , but we’ll add :( γ, V ) = ( m, m O ( y = 0) /ν ) m Tot, fuel (2.39)And we get the dimensionless gaseous phase equation : ∂γ∂η = ∂ γ∂ξ + 1 P e · ∂ γ∂η + ¯ S d (2.40)Where (cid:0) P e, ¯ S d (cid:1) stands for Peclet number and the normalized source element : P e = U d RD g ¯ S d = S d R D g m Tot, fuel (2.41)16he source element in the gaseous fuel equation equals to the droplets evaporation rateand stems from the overall rates over N sections :¯ S d (cid:44) N (cid:88) j =1 ∆ i γ d i − ψ i γ d i +1 (2.42)Such that the dimensionless gaseous phase equation : γ − P e ∂γ∂η = − (cid:80) Ni =1 δ i , ≤ ξ ≤ c − V , c ≤ ξ ≤ ∂γ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) η ≥ ξ =0 , = 0 ∂γ∂η (cid:12)(cid:12)(cid:12)(cid:12) η →∞ ≤ ξ ≤ = 0 (2.44) Using S-Z transform, the normalized mass fraction fulfils : γ (cid:44) γ F − γ O .We can solve it as a sum of the following equations :Homogenous : ∂γ h ∂η = ∂ γ h ∂ξ + 1 P e · ∂ γ h ∂η (2.45)Particular : ∂γ p ∂η = ∂ γ p ∂ξ + 1 P e · ∂ γ p ∂η + ¯ S d (2.46) Homogeneous solution
Both equations will be solved using the separation of variables method - γ h = f h ( η ) · g h ( ξ ) df h dη · g h = f h · d g h dξ + 1 P e · d f h dη · g h (2.47) (cid:18) df h dη − P e · d f h dη (cid:19) /f h = d g h dξ /g h = Const. (cid:44) − α (2.48)17olve g h ( ξ ) d g h dξ = − αg h → g h = C sin( √ αξ ) + C cos( √ αξ ) (2.49)Solve f h ( ξ ) : d f h dη − ( P e ) df h dη − ( αP e ) f h = 0 (2.50)Guess the exponential solution of the form of f h = C e C η and plug it above :Plugging ⇒ · ( 1 C e C η ) ⇒ C − ( P e ) C − ( αP e ) = 0 (2.51) C , = 12 [( P e ) ± (cid:112) ( P e ) + 4 α ( P e ) ] < (cid:124)(cid:123)(cid:122)(cid:125) fulfill BC C = ( P e )2 (cid:20) − (cid:115) α ( P e ) (cid:21) (cid:44) q n (2.53)And finally we get the homogeneous equation solution : γ h = f h · g h = e q n η ( C sin( √ αξ ) + C cos( √ αξ )) (2.54) Particular solution
Also here, using separation of variables − γ p = f p ( η ) · g p ( ξ ) (cid:18) df p dη − P e · d f p dη (cid:19) · g p = d g p dξ · f p + ¯ S d ; ¯ S d = − N (cid:88) j =1 ∂γ d j ∂η (2.55)Based on the dimensionless spray BC (2.20), we’ll substitute ∂γ dj ∂η such that :¯ S d = − N (cid:88) j =1 ∂∂η N (cid:88) i = j Ω ij e − ∆ i η = N (cid:88) j =1 N (cid:88) i = j ∆ i Ω ij e − ∆ i η (2.56)Ω ij ( ξ ) ⇒ ¯ S d ( ξ, η ) , such that Ω ij can be seen as an influence coefficient on E, whereas i indicates the droplets influential section and j indicates the influenced section. An18lternative summation to these developments proposes shifting the indices such that : N (cid:88) j =1 N (cid:88) i = j ∆ i Ω ij e − ∆ i η = N (cid:88) i =1 i (cid:88) j =1 ∆ i Ω ij e − ∆ i η (2.57)Such that ¯ S d can be represented as : S d = N (cid:88) i =1 ∆ i e − ∆ i η (cid:124) (cid:123)(cid:122) (cid:125) f ( η ) i (cid:88) j =1 Ω ij (cid:124) (cid:123)(cid:122) (cid:125) g ( ξ ) ; i (cid:88) j =1 Ω ij (cid:44) H i ( ξ ) (cid:124) (cid:123)(cid:122) (cid:125) Heavisidefunction = , ξ > , else (2.58)The latter term is dependent on the iDSD and is nullified at H ( ξ > c ) = 0 . By usingFourier series we can present the Heaviside function discretely : H i ( ξ ) = 12 k i + ∞ (cid:88) m =1 k m i cos( mπξ ) (2.59)Equivalently with ¯ S d particular solution : γ p = N (cid:88) i =1 (cid:32) b i + ∞ (cid:88) n =1 b n i cos( nπξ ) (cid:33) e − ∆ i n (2.60)After plugging inside (2.55) and grouping the elements : N (cid:88) i =1 (cid:34) − (cid:18) ∆ i + ∆ i P e (cid:19) (cid:32) b i + ∞ (cid:88) n =1 b n i cos( nπξ ) (cid:33) + ∞ (cid:88) n =1 ( nπ ) b n i cos( nπξ ) (cid:35) e − ∆ i η . . . (2.61)= N (cid:88) i =1 ∆ i e − ∆ i η (cid:32) k i + ∞ (cid:88) m =1 k m i cos( mπξ ) (cid:33) Equating the coefficients to extract b i , b n i :12 N (cid:88) i =1 (cid:18) ∆ i − ∆ i P e (cid:19) b i e − ∆ i η = 12 N (cid:88) i =1 ∆ i k i e − ∆ i η (2.62) ⇒ b i = − ∆ i ∆ i + (cid:0) ∆ i P e (cid:1) · k ,i (2.63)Similarly with b n i ⇒ b n i = − ∆ i ∆ i + (cid:0) ∆ i P e (cid:1) − ( nπ ) · k n i (2.64)19ecall the γ = γ h + γ p solution such that finding the coefficients of C , C , α can be doneby matching the solution to the BC : ∂γ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 = √ αe q n η C = 0 ⇒ C = 0 (2.65) ∂γ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =1 = −√ αe q n η C sin( √ α ) = 0 ⇒ C n (cid:54) = 0 , α = ( nπ ) n ∈ N (2.66)Therefore the homogeneous solution is : γ h = ∞ (cid:88) n =1 e q a η C n cos( nπξ ) (2.67)where q n (cid:44) ( P e )2 (cid:20) − (cid:115) nπ ) ( P e ) (cid:21) (2.68)Applying the BC in gaseous phase equation (2.43) : (cid:18) γ − P e ∂γ∂η (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) η =0 = − (cid:80) Ni =1 δ i , ≤ ξ ≤ c − V , c ≤ ξ ≤ RHS (cid:44) d ∞ (cid:88) n =1 d n i cos( nπξ ) (2.70) d = 2 (cid:90) c (cid:16) − N (cid:88) i =1 δ i (cid:17) dξ + 2 (cid:90) c ( − V ) dξ = 2 c (cid:16) − N (cid:88) i =1 δ i (cid:17) + 2 V ( c −
1) (2.71)Such that : d = 2 c (1 − N (cid:88) i =1 δ i + V ) − V (2.72)Similarly with d n : d n = 2 (cid:90) c (cid:16) − N (cid:88) i =1 δ i (cid:17) cos( nπξ ) dξ − (cid:90) c V cos( nπξ ) dξ (2.73) d n = 2 (cid:16) − (cid:80) Ni =1 δ i (cid:17) nπ sin( nπc ) + 2 Vnπ sin( nπc ) (2.74)20pplying η = 0 on the RHS : RHS = (cid:32) − N (cid:88) i =1 δ i + V (cid:33) (cid:34) c + 2 π ∞ (cid:88) n =1 sin( nπc ) n cos( nπξ ) (cid:35) − V (2.75)Developing the LHS expression : γ − P e ∂γ∂η = C + N (cid:88) i =1 e q a η C n cos( nπξ ) (cid:16) − q n P e (cid:17) · · · + N (cid:88) i =1 e − ∆ iη (cid:32) b i + ∞ (cid:88) n =1 b n i cos( nπξ ) (cid:33) (cid:18) i P e (cid:19) (2.76)Applying LHS (cid:12)(cid:12)(cid:12) η =0 = RHS and sorting the elements : (cid:18) γ − P e ∂γ∂η (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) η =0 = C + N (cid:88) i =1 b (cid:18) i P e (cid:19) . . . (2.77)+ ∞ (cid:88) n =1 (cid:34) C n (cid:16) − q n P e (cid:17) + N (cid:88) i =1 b n i (cid:18) i P e (cid:19)(cid:35) cos( nπξ ) = (cid:32) − N (cid:88) i =1 δ i + V (cid:33) (cid:34) c + 2 π ∞ (cid:88) n =1 sin( nπc ) n cos( nπξ ) (cid:35) − V Equating the coefficients to extract C , C n : C = c (1 + V ) − V − N (cid:88) i =1 cδ i + 12 b i (cid:18) i P e (cid:19) (2.78) C n = 2(1 + V ) sin( nπc ) nπ (cid:0) − q n P e (cid:1) − N (cid:88) i =1 b n i (cid:0) ∆ i P e (cid:1) nπ + 2 δ i sin( nπc ) (cid:0) − q n P e (cid:1) nπ (2.79)Finally, after resorting the elements we can write down the full expression for γ : γ = c (1 + V ) − V + 12 N (cid:88) i =1 ( k i − cδ i + b i e − ∆ i η (cid:1) + 2 π (1 + V ) ∞ (cid:88) n =1 sin( nπc ) n (cid:16) − q n P e (cid:17) e q n η cos( nπξ ) . . . (2.80)+ ∞ (cid:88) n =1 N (cid:88) i =1 b n i (cid:32) e − ∆ i η − (cid:0) ∆ i P e (cid:1) nπ + 2 δ i sin( nπc ) (cid:0) − q n P e (cid:1) nπ e g n η (cid:33) cos( nπξ )Whereas the following terms ( k , k n i , b i , b n i , q n )= (2.58, 2.59, 2.63, 2.64, 2.68) .21 .3 The Temperature Equation Similarly to the gaseous phase development from the conservation of mass law, we canwrite the temperature equation from the conservation of energy law. According to Fourierlaw of thermal conduction, the heat flux is linear to the spatial gradient of the temperature: q = − λ (cid:124)(cid:123)(cid:122)(cid:125) Thermal conduction · ( ∇ ˆ T ) (cid:124) (cid:123)(cid:122) (cid:125) Dimensionless temperature (2.81)Using the continuity equation for the energy yields the heat equation: U g ∂ ˆ T∂y = K (cid:32) ∂ ˆ T∂x + ∂ ˆ T∂y (cid:33) + ¯ S reaction (cid:124) (cid:123)(cid:122) (cid:125) Energy source + ¯ S d, vapor. (cid:124) (cid:123)(cid:122) (cid:125) Energy sink + ¯ S d, burning (cid:124) (cid:123)(cid:122) (cid:125) Energy source ; K = λρC P (2.82)The ¯ S terms denote the energy transferred along the process. However, it is customary toassume that no energy transferred in the ∆ ˆ T between the droplets and its carrier gas. The BC at the channel’s exit includes a thermal diffusive flux element. The temperature ofboth gaseous fuel and oxidizer (= ˆ T ) and equals to the droplets’ evaporation temperature :ˆ T − KU , fuel ∂ ˆ t∂y = ˆ T (2.83)Applying the following assumptions : ∂ ˆ T∂x (cid:124)(cid:123)(cid:122)(cid:125)
Symmetry (cid:12)(cid:12)(cid:12)(cid:12) y ≥ x =0 = 0 ∂ ˆ T∂x (cid:124)(cid:123)(cid:122)(cid:125)
Insulated walls (cid:12)(cid:12)(cid:12)(cid:12) y ≥ x = R = 0 ∂ ˆ T∂y (cid:124)(cid:123)(cid:122)(cid:125)
Thermodynamic equilibrium (cid:12)(cid:12)(cid:12)(cid:12) y →∞ ≤ x ≤ R = 0 (2.84)Axes normalization is according to (2.14) and we’ll normalize the reference temperature (cid:16) ˆ T ref (cid:17) and the injected fluid (cid:16) ˆ T (cid:17) temperature :( T ) = (cid:32) ˆ T − ˆ T ˆ T ref (cid:33) ; ˆ T ref = ˜ q reac m Tot, fuel C P (2.85)22he energy sinks and sources elements will be normalized as such : (cid:0) ¯ S reac , ¯ S d,v , ¯ S d, b (cid:1) = R D g ˆ T ref (cid:16) ˜ S reaction , ˜ S d , vapor . , ˜ S d , burning (cid:17) (2.86)Such that the dimensionless temperature equation : ∂T∂η = Le · (cid:18) ∂ T∂ξ + 1 P e ∂ T∂η (cid:19) + ¯ S reac + ¯ S d,v + ¯ S d,b (2.87)We shall work with Lewis number of : Le (cid:44) KD g = λρC P D g = 1 (2.88)Meaning that the rate of the reactants’ mass diffusion towards the reaction zone, is equalto the heat’s thermal diffusion from the reaction zone. Now, let us elaborate on thedimensionless sources elements. ¯ S d,v expresses the heat absorption resulted from thedroplets’ evaporation and is equal to the evaporation rate multiplied by latent heat :¯ S d,v = − Λ S d H ( c − ξ ) ; H ( ξ ) (cid:124) (cid:123)(cid:122) (cid:125) Heavisidefunction = , ξ > , else (2.89)Note that the Heaviside step function acts as an ”on / off” switch controlling the source.¯ S d,b expresses the external heat emitted from the moving burned droplets, after beingignited. Using former assumptions we get the equality E burned = E pre-ignited such that :¯ S d,b = ¯ S d H ( c − ξ ) H ( η − η f ) (2.90) η f denotes the flame height with respect to the ξ axis and H ( ξ ) indicates the obtainedheat regions. Using S-Z transform we’ll define : γ T (cid:44) T + γ F (2.91)Such that by utilizing Eqs. (2.14, 2.87) we get : ∂γ T ∂η = ∂ γ T ∂ξ + 1 P e · ∂ γ T ∂η + (1 − Λ) ¯ S d H ( c − ξ ) (2.92)23he dimensionless BC for the temperature equation : T = T + (cid:18) λρC P D g (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) Le =1 (cid:32) D g U g R (cid:33) ∂T∂η ⇒ η =0 T − P e ∂T∂η = T (2.93)And on the normalized axes ∂T∂ξ (cid:12)(cid:12)(cid:12)(cid:12) η ≥ ξ =0 , = 0 ∂T∂η (cid:12)(cid:12)(cid:12)(cid:12) η →∞ ≤ ξ ≤ = 0 (2.94)The dimensionless BC for the gaseous fuel : γ F − P e ∂γ F ∂η = (cid:40) − (cid:80) Ni =1 δ i , ≤ ξ ≤ c , c ≤ ξ ≤ ∂γ F ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) η ≥ ξ =0 , = 0 ∂γ F ∂η (cid:12)(cid:12)(cid:12)(cid:12) η →∞ ≤ ξ ≤ = 0 (2.96)Finally we get the full expression for S-Z transform (2.91) such that dimensionless γ T : γ T − P e ∂γ T ∂η = (cid:40) − (cid:80) Ni =1 δ i + T , ≤ ξ ≤ cT , c ≤ ξ ≤ ∂γ T ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) η ≥ ξ =0 , = 0 ∂γ T ∂η (cid:12)(cid:12)(cid:12)(cid:12) η →∞ ≤ ξ ≤ = 0 (2.98) Recall γ T (cid:44) T + γ F (4 .
11) and the dimensionless temperature equation (4.12) ∂γ T ∂η = ∂ γ T ∂ξ + 1 P e · ∂ γ T ∂η + (1 − Λ) S d H ( c − ξ ) (2.99)24sing the BC elaborated at (2.97, 2.98) we can solve γ T as a sum of γ T = γ T h + γ T p Homogenous : ∂γ T h ∂η = ∂ γ T h ∂ξ + 1 P e · ∂ γ T h ∂η (2.100)Particular : ∂γ T p ∂η = ∂ γ T p ∂ξ + 1 P e · ∂ γ T p ∂η + (1 − Λ) ¯ S d (2.101) Homogeneous solution
Using the same solution method as in the gaseous phase section γ T h = f T h · g T h = e q n η (cid:0) C sin( √ αξ ) + C cos( √ αξ ) (cid:1) (2.102) Particular solution
We’ll define the following variable, and plug terms from (2.58) :¯ S (cid:48) d = (1 − Λ) ¯ S d = (1 − Λ) N (cid:88) i =1 ∆ i e − ∆ i η (cid:124) (cid:123)(cid:122) (cid:125) f ( η ) i (cid:88) j =1 Ω ij (cid:124) (cid:123)(cid:122) (cid:125) g ( ξ ) (2.103)Also here, we’ll define: H (cid:48) i ( ξ ) (cid:44) (1 − Λ) H i ( ξ ) = (1 − Λ) i (cid:88) j =1 Ω ij (2.104)Using Fourier ⇒ H (cid:48) i ( ξ ) = 12 k (cid:48) i + ∞ (cid:88) m =1 k (cid:48) m i cos( mπξ ) (2.105)Its coefficients are dependent on the following relations (2.58, 2.59) : (cid:0) k (cid:48) i , k (cid:48) m i (cid:1) = (1 − Λ) ( k i , k m i ) (2.106)Such that the homogeneous equation solution is ( q n @ 2.68 ) : γ T p = N (cid:88) i =1 (cid:32) b (cid:48) i + ∞ (cid:88) n =1 b (cid:48) n i cos( nπξ ) (cid:33) e q n η (2.107)25sing similar relations as before (s.t. η = 0 ): ⇒ b (cid:48) i = − ∆ i ∆ i + (cid:0) ∆ i P e (cid:1) · k (cid:48) ,i (2.108)Similarly with b (cid:48) n i ⇒ b (cid:48) n i = − ∆ i ∆ i + (cid:16) ∆ i P e (cid:17) − ( nπ ) · k (cid:48) n i (2.109)Applying ( V = 0 , η = 0) on the RHS : RHS (cid:12)(cid:12)(cid:12) η =0 = c (cid:32) − N (cid:88) i =1 δ i + T (cid:33) + 2 π (cid:32) − N (cid:88) i =1 δ i (cid:33) ∞ (cid:88) n =1 sin( nπc ) n cos( nπξ ) (2.110)Based on (2.77) development we’ll sort the elements: (cid:18) γ − P e ∂γ∂η (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) η =0 = C (cid:48) + N (cid:88) i =1 b i (cid:18) i P e (cid:19) . . . (2.111)+ ∞ (cid:88) n =1 (cid:34) C (cid:48) n (cid:16) − q n P e (cid:17) + N (cid:88) i =1 b (cid:48) n i (cid:18) i P e (cid:19)(cid:35) cos( nπξ ) (2.112)Equating the coefficients to extract C (cid:48) , C (cid:48) n C = c + T − N (cid:88) i =1 cδ i + 12 b (cid:48) i (cid:18) i P e (cid:19) (2.113) C (cid:48) = 2 sin( nπc ) nπ (cid:0) − q n P e (cid:1) − N (cid:88) i =1 b (cid:48) n i (cid:0) ∆ i P e (cid:1) nπ + 2 δ i sin( nπc ) (cid:0) − q n P e (cid:1) nπ (2.114)Finally, after resorting the elements we can write down the full expression for γ T : γ T = c + T + 12 N (cid:88) i =1 (cid:0) k (cid:48) i − cδ i + b (cid:48) i e − ∆ i η (cid:1) + 2 π ∞ (cid:88) n =1 sin( nπc ) n (cid:0) − q n P e (cid:1) e q n η cos( nπξ ) . . . (2.115)+ ∞ (cid:88) n =1 N (cid:88) i =1 b (cid:48) n i (cid:32) e − ∆ i η − (cid:0) ∆ i P e (cid:1) nπ + 2 δ i sin( nπc ) (cid:0) − q n P e (cid:1) nπ e q n η (cid:33) cos( nπξ ) (2.116)Whereas the following terms (cid:0) k (cid:48) , k (cid:48) n i , b (cid:48) i , b (cid:48) n i , q n (cid:1) = (2.106, 2.106, 2.108, 2.109, 2.68) .26 .4 Results Based on (2.91) the dimensionless temperature field can be expressed as : T ( ξ, η ) = γ T − γ if ξ ≤ ξ ( γ ≈ γ T if ξ ( γ ≈ < ξ ≤ R (2.117)The flame front contains the set of points that fulfill the stoichiometric ratio between thereactants and thus lies on the black dashed line ξ ( γ ≈
0) . In the gaseous fuel zone themass fraction satisfies - γ = γ F , where in the oxidizer zone it’s nullified - γ F = 0 :We can see that the flame base is formed at half the inner channel c = LR , while itsmaximum height is obtained at the middle of the symmetric flame. Trivially, differentinput setups would result in different values of the temperature field.However, of all existing points, we are interested in two specific points that might shedlight on the performance of a given execution, obtained via the following functions : ◦ η max = f η ( ¯ E, d, δ ) ⇒ Returns the highest value of the flame front. ◦ T max = f T ( ¯ E, d, δ ) ⇒ Returns the highest value of the temperature field.27he temperature field can be also seen as a set of 3D points - (cid:0) ξ, η, T ( ξ, η ) (cid:1) :Both variables at the edge of the front flame will serve as the main indicators of a givenexecution. In the following subsections we will further inspect them, in order to establish ageneral intuition with respect to several different factors.28 .4.1 Maximum flame height The following simulations examine the model’s sensitivity in terms of maximum flameheight, namely η max reaction to several initial parameters : (i) liquid fuel fraction - δ i (ii)the droplets section - d and (iii) the evaporation rate - ¯ E . Different reactions to the iDSDcomposition were observed, and thus the following scenario will present them. Monosectional iDSD
An iDSD is said to be monosectional when the liquid fuel fraction occupies only onesection, namely the initial droplet size is uniform. Consider the following set ofmonosectionals, executed at a wide range of sections and varying amounts of liquid fuel :Given five section sizes - d = { } , each contains five executions of different amountsof liquid fuel - δ d = { . . . } . Note the overall tendency of η max to increase alongthe section sizes, in a typical tradeoff with the evaporation rate ¯ E .The red dashed line denotes the gaseous flame height towards highest ¯ E →∞ , see detaileddiscussion next. 29 olysectional iDSD The polysectional scenario is an iDSD characterized by a multimodal distribution :In this scenario, executions i = [1 ,
9] are polysectionals composed of random iDSDinitialization, whereas each total sum equals to one. Execution i = 10 acts as a ”controlgroup” as it has only one section, namely a monosectional of δ d =4 = 1 . E zones (vertical green separation) :¯ E low := η max is governed mainly by the amount of liquid fraction ( δ ) (2.118)¯ E mid := η max is dependent significantly by the evaporation rate ( ¯ E ) (2.119)¯ E high := η max approaches to the gaseous flame height (dashed red) : (2.120)lim ¯ E →∞ f ( ¯ E ) ≈ . ∀ { d, δ } Additionally, let us define the intersection between ¯ E low and ¯ E mid as a reversal point - ¯ E rev. from which bigger fuel fractions yield higher flame, as opposed to before ( ¯ E < ¯ E rev. ) .Note the { d − ¯ E } correlation, and its impact on the flame height : ◦ Small particles ( d = 1) yield optimal flame at relatively low ¯ E mid ≈ ∀ δ . ◦ Larger particles ( d = 5) yield optimal flame relatively higher ¯ E mid ≈ ∀ δ . In this scenario I would like to go through the same process regarding T max .31 onosectional iDSD Starting off with the same setup of monosectional iDSD executed previously :In a closer look around the optimal zone ¯ E mid :Either here, in the absence of any constraints, we get T max ↑ for any { δ ↑ , d ↑ } .32 olysectional iDSD Using a random initialization for all executions except i = 6, which is a monosectional :Mind the graphs curvatures that later will be discussed. At a closer look around ¯ E mid :Similarly to what we saw at η max , the ”hottest” T max is obtained for i = 6, which is the33nly sample that is not a polysectional iDSD. Both η max and T max exhibited a typical pattern as elaborated on (2.118 - 2.120). At bothcases, the polysectional iDSD performed poorer in comparison with the monosectional.Towards higher evaporation rates ( ¯ E high ), it seemed that all of the graphs approach thegaseous flame height, regardless their initial conditions. Here - η max ≈ . ∀ { d , δ } .Quite similarly, the temperature graphs show an asymptotic behavior towards ¯ E high . Thegaseous flame temperature satisfies T f ≈ . E , where T max reaction’s benefits with smaller amounts of the liquid fuel :34ooking closely, we can tell that slow evaporation rates benefit with small δ amountscontrarily to bigger δ that perform significantly poorer. That mechanism can be explainedby the liquid ”overload” that absorbs more heat during the evaporation process :The qualitative nature of the graphs shows how the ¯ E rev. changes with respect to d sectionand acts as a ”role reversal” between δ and ¯ E . Note that the evaporation rate axis islogarithmic, and hence the differences turn more dramatic as we move forward.The polysectional iDSD have shown poorer performances in comparison with themonosectional, as it affects directly the droplets volume change :(i) It downsizes the droplets size inside a given section.(ii) It is responsible for the joining of new droplets from higher sections.It is therefore no surprise that either η max or T max curvatures, are appearing in accordancewith occupied sections. 35 The Optimization Model
In this section the combustion model will be adjusted to the optimization model. At first, aformal development of the GA will be presented, then the combustion model will be mergedwithin. Consider the following informal description that demonstrates that process :
Algorithm 1:
Pseudocode of the Genetic AlgorithmInitialize random population (cid:16) X t ∈ R n × (cid:17) : t ← for x j ∈ X do x j ← random (cid:0) ¯ E, d, δ (cid:1) ∈ R end Commence cost minimization until termination condition is met : while t < GA
MaxIteration or GA Cost > GA
T ermination do Evaluate Fitness(): Y t ← f ( X t )Selection():[ x It , x IIt ] ← X t (cid:16) max( Y t ) k =2 (cid:17) Crossover( X t ):[ x (cid:48) t , x (cid:48)(cid:48) t ] ← recombine (cid:16) x It , x IIt (cid:17)
Mutation(): X t +1 ← M utate (cid:16) x (cid:48) t , x (cid:48)(cid:48) t (cid:17) t ← t + 1 end Extract best candidate from last generation :
Return x ∗ j ← X t (cid:16) max( Y t ) k =1 (cid:17) Legend f ( · ) Fitness function ( f η or f T ) x t Chromosome at time ty t Chromosome’s fitness at time tX t Population at time tY t Population’s fitness at time t max( Y t ) k Find k largest elements X t (cid:0) g ( · ) (cid:1) Extract x t ∈ X t that satisfies g ( · )36 ote : The GA text denotes the GA’s internal functions which are problem-independent,and are decoded in accordance with the algorithmic design. Moreover, the GA uses abit-string representation to encode the solutions, in the form of chromosomes. The canonical form of an optimization problem [20] is commonly written as :minimize x f ( x ) (3.1)subject to g i ( x ) ≤ , i = 1 , . . . , m (3.2) h i ( x ) = 0 , i = 1 , . . . , p (3.3)where f : R n → R The objective function g i ( x ) ≤ h j ( x ) = 0 Equality constraintAs mentioned in Definition [ B.3 ], a maximization problem can be solved analogously bysimply negating the objective function : max x f ( x ) ⇔ min x − f ( x ) .The first technical obstacle I encountered, was Matlab’s inability to execute equalconstraints while one of the search variables is discontinuous. The GA becomes inoperable,as the section variable ( d ) is restricted to be an integer N > , and thus breaks thecontinuity assumption of the search space [21]. Commonly, these class of problems can beaddressed as ”Mixed-integer linear programming” (MILP) [22],[23].However, after enough search the answer was found in one of Matlab’s forums [24], thatproposed to achieve the equality constraint by setting two inequality constraint :desired : ax + bx = c (impossible) (3.4)workaround : ax + bx ≤ c & − ax − bx ≤ − c (3.5)Practically, the executions were conducted as in Eq. (3.5), but for simplicity reasons I willwrite them regularly as an equality constraint, as in Eq. (3.4).37s mentioned above, there are two types of problem’s constraints :(3.2) - Inequality constraint (3.3) - Equality constraintWhile inequality imposes a half-constraint on a variable between a certain point onwards,the equality constraint imposes a full-constraint. In other words, an equality constraint ona single variable, degenerates its participation in the optimization process, and thusreduces the overall degree of freedom (DoF) by one. Contrarily, a variable that is notsubjected to any constraint is said to be unconstrained, as it is completely free to besampled within the domain x k ∈ ( −∞ , ∞ ), and thus contributing 1 DoF to the system. Here, the objective functions can return either η max or T max . The inequality constraint isthe liquid fraction that’s constrained to be smaller than a desired value, and the equalityconstraint is the constant evaporation rate. On Matlab :Thereby, the program is free to find the optimal iDSD as d is an unconstrained variable :max ¯ d f ( ¯ E, d, δ )subject to 1 ≤ d i ≤ ∀ i = 1 , ..., n sections (cid:88) i δ i ≤ δ const ≤ E = ¯ E const The i index denotes the serial section number, which cannot be greater than the totalnumber of sections ( i ≤ d i (cid:98) = δ i ).38onsider the following 3 DoF toy-example, which optimizes the iDSD (cid:12)(cid:12) i ≤ :We can see that out of all feasible combinations, the optimal distribution consists of the2nd, 6th and 9th sections, each correlates to a different index of liquid fraction.The following figure shows the on-line optimization process as presented by Matlab :The x-axis denotes the iterative process of forming better generations that can betterperform the fitness function on the y-axis (penalty value). Note that the values here areexpressed as negative, as by default the algorithmic is designed to minimize a penalty.Therefore the problem was equivalently adjusted, by changing the output sign [ B.3 ] .39
Results
This section presents a wide set of examinations of the optimal flame height and theoptimal tip flame temperature. Each examination is conducted at growing levels of DoF,over several evaporation rates, starting from a single DoF towards 9 DoF (total sections).However, in order to fit with the algorithm requirements, it must fulfill two constraints :(i) (cid:80) i δ i ≤ . E = const. - the evaporation rate must be constant with timeThereby, the algorithm will be free to compute the optimal iDSD in return. Note : d (number of sections) and DoF (degree of freedom) are used interchangeably. Starting from a single section optimization, we’ll utilize (cid:80) i δ i ≤ . turquoise column should be less than or equal to 0 .
7, at some calculations itturns to be greater due to small numerical noise, stemmed by residues from the GAinitialization. Described in details at the following document [25].40rom here on, we’ll rise the number of sections to DoF = { , , , } :The red marked cells are designated to emphasize the repeated pattern of a given sectionwith respect to certain evaporation rate, namely optimal ¯ E .41onsider the following executions for DoF = { , , , } : Overall, we can see that across all of the executions, the flame height exhibits values withina bounded range of η max ∈ [0 . , . d ∗ = { , , , , , } patternappears across all executions, regardless the number of DoF / sections.However, the optimization ”quality” shows growing signs of decay given more sections tooptimize, as more liquid fuel occupies the neighboring sections, on account of a singlesection (which proved to guarantee optimality).42um it in a scatter plot :Overall, given more sections to optimize we can say that :(i) Maximum flame height decreases - d ↑ ⇒ η max ↓ (ii) Maximum liquid fuel decreases - d ↑ ⇒ max(¯ δ ) ↓ Similarly to η max scheme, we’ll do the same with the maximum tip flame temperature :Note that the same pattern obtained before, is achieved either here. The reason for that isthe close relation of η max and T max , as mentioned here, where η max location tends tocoincide with T max , and vice versa. 43ising the number of sections - DoF = { , , , } :44onsider the following executions for DoF = { , , , } Same optimal patterns d ∗ = { , , , , , } as obtained before, seem to appear here aswell. The main difference however, is the T max values whose range is slightly above that of η max , and the relative differences between the values of T max are smaller than those of η max .45 .2.1 Discussion Let us sum the above tables in a scatter plot :Either here the same behavior occurs, when more sections are optimizable :(i) Maximum tip flame temperature decreases - d ↑ ⇒ T max ↓ (ii) Maximum liquid fuel decreases - d ↑ ⇒ max(¯ δ ) ↓ However, the range of T max is much smaller, and relative differences are aggravated as d ↑ . An interesting phenomenon that can be seen at both η max and T max plots, is the lowerperformances of : f ( ¯ E =2 · ) < f ( ¯ E =10 ) . Having no other constraint, shouldn’t the readerexpect that higher evaporation rates will guarantee better performance ? He definitelyshould. But one must keep in mind that the discrete sectioning of the droplet size intomax ( d ) ≤
9, allows a narrow ”optimal evaporation zone” for each section size.Therefore, the local decrease of T max ( ¯ E =20 , ) is because of the predetermined ¯ E = 20 , .3 Optimization summary In this section I executed two optimization scenarios, each composed of a varied DoF,constrained only by (cid:80) i δ i ≤ . E const . We saw that the GA was able to find anoptimal iDSD, across different parametric configurations. That is to say, that by properlydelivering a complex problem to the GA, a near-optimal solution will (eventually) be found.Moreover, it was empirically shown that optimality of both η max and T max is guaranteedmostly when the liquid fuel occupies only one single section. Namely, optimal performanceis guaranteed when the iDSD is monosectional. But why is that ? The answer → [4.5]The DoF number was found as an important player in the optimization game, as higherDoF decay the results’ quality. For that reason, the following scheme will settle for 6 DoF. In this part I would like to examine the P´eclet number influence on a given optimal iDSD,which so far was determined by default as
P e = 10 .What is the P´eclet Number ?A class of dimensionless numbers relevant in the study of transport phenomena thatdenotes the ratio of the advection rate by the diffusion rate driven by an appropriategradient [26]. In the context of mass transfer, the P´eclet number is the product of the
Reynolds number and the
Schmidt number :
P e R = Re R Sc = U g RD g ◦ D g - Mass diffusion coefficient ◦ U g - Local flow velocity ◦ R - Characteristic length (half external channel)The next section consist of two scenarios of both η max and T max .Each scenario comprises of two different subsections :(i) A general sensitivity test of a monosectional iDSD to P e number (for perspective).(ii) A specific sensitivity test of the polysectional iDSD obtained at 6 DoF optimization.47 .4.1 Pe vs. flame height
Consider the following η ( (cid:80) δ =0 . max space presenting 3 different monosectionals - d = { , , } :Notice the ambivalent relation of P e depending on the evaporation rates (rev. ≡ reversal) : ¯ E < ¯ E rev. , P e ↑ ⇒ η max ↑ ¯ E > ¯ E rev. , P e ↑ ⇒ η max ↓ Now, consider the polysectional iDSD obtained at η max optimization in 6 DoF :Now we shall examine their sensitivity to several P e numbers, where the rounded rectanglein the bottom left side denotes the gaseous flame height sensitivity to
P e .48 .4.2 Discussion
In the first figure we can see a general solution space where some ¯ E regions benefit withsome P e while others are harmed. Small
P e numbers at low evaporation rates ( ¯ E P e = 1000 seems almost negligible, for any evaporation rate. In the second(above) figure, the ambivalent policy no longer exists, and instead we get : P e ↓ ⇒ η max ↑ ∀ ¯ E The reader may wonder how is that possible ?The answer goes back to the optimal iDSD, as stated here. Although being optimized as apolysectional, about ≈ 96 % of the liquid fuel is concentrated in one section , making it asif it was a monosectional iDSD. Knowing that, we can see that all these optimal iDSDstake place after the reversal points, where smaller P e numbers benefits with η max .49 .4.3 Pe vs. flame temperature Consider the following T ( (cid:80) δ =0 . max space presenting 3 different monosectionals - d = { , , } :The ambivalent relation seen here seems to be completely opposite to the previous one : ¯ E < ¯ E rev. , P e ↑ ⇒ T max ↓ ¯ E > ¯ E rev. , P e ↑ ⇒ T max ↑ Looking closer we can see that the reversal point ( ¯ E rev. ) is closer to the global maximum :50onsider the polysectional iDSD obtained at T max optimization at 6 DoF :Now we shall examine this iDSD sensitivity to several P e numbers, where the roundedrectangle in the bottom left side denotes the gaseous flame temperature sensitivity to P e .As mentioned above the overall policy of T max is completely opposite to that of η max .However, the relative differences are less dramatic : η ( P e =3) max η ( P e =1000) max ≈ . > T ( P e =1000) max T ( P e =3) max ≈ . η ( E →∞ ) max , here the gaseous T ( E →∞ ) max shows a weak response to P e .51 .4.4 Discussion The first figure exhibited a general solution space of T max where some ¯ E regions benefitedwith P e while others were harmed. In the ¯ E E > ¯ E rev. , where higher P e number benefits with T max .As an auxiliary argument, consider the following reference from the thesis (p. 117) [10],presenting the extinction maps that describe the flame’s sensitivity to the iDSD and P e :Without deep diving into Li˜n´an’s diffusion flame theory [27], and without loss of generality,we can carefully say that higher P e numbers reduce the extinction regions.Equivalently, that is to say that higher P e numbers expand the existence regions of theflame itself. Note that this example’s iDSD is arbitrary but it still provides a usefulqualitative information that supports my findings.52 .5 Empirical validation The above executions have shown the monosectional superiority in achieving optimalperformances. In attempt to explain that, I would like to provide another point of view foranalysis, that might shed some light on that behavior and utilize as a reliable sanity check.Consider the following image, presenting seven different initial distributions containing thesame total liquid fraction. The first iDSD is monosectional, but as we go forth thestandard deviation increases and spreads to more sections :The 1st iDSD (:=’data1’) contains only one section δ d =5 = 0 . δ d =4 + δ d =5 + δ d =6 = 0 . 05 + 0 . . 05 = 0 . black ), and themore the distribution is spread (more sections are ”occupied”), so η max turns lower.However, one may come up and claim that at a certain closed interval, the monosectionaliDSD actually performs the poorest, while others are optimal :Indeed, the polysectional iDSDs do perform better on lower ¯ E , as the distribution spreadwider. However, that argument is only half true. Let us add two more monosectionaldistributions to the current figure. 54he dashed line from left is a monosectional of δ d =2 , and the right one is of δ d =7 :Note that any non-monosectional distribution is situated beneath some of themonosectional. Adding more monosectionals and we get the full picture :The conclusion is clear, any polysectional iDSD is always bounded between neighboringmonosectionals. Therefore the latter performs better at any evaporation rate. As atouchstone, I added a uniform distribution upon all sections - δ d =[1:9] = . ( dotted pink ) .55 Conclusions In this research project I integrated a combustion model within a Genetic Algorithm in aninnovative manner that was not performed before. The main research question :”What form of an iDSD will guarantee the optimal flame performances ?”Was reviewed from different point of views, at first to gain a bird’s-eye view and afterwardsby a closer inspection. The optimization scheme left no room for doubt about the optimaldistribution for both η max or T max , which was found to be the monosectional iDSD.Moreover, a sensitivity test of the model with respect to P´eclet number was provided,followed by an empirical validation that supported the research findings.To conclude, the GA was successfully harnessed for the sake of the developed engineeringproblem, thereby reinforcing the idea of problems that were so far non-optimizable, cannow be optimized and provide an heuristic optimal solution. Along my working process I happened to think through several directions for future work : ◦ Experimentally - Given appropriate equipment, do the optimal values found convergewith real laboratory experiments ? ◦ Technically - Do the number of sections (9 by default) have any influence upon theresults ? Would more sections necessarily lead to a more continuous results space ? ◦ Complexity Analysis - Long and costly computation process can be analyzed, in orderto point to bottlenecks and to make the algorithm more efficient. ◦ More approaches - Nowadays, the GA main competitor is the AI’s top notch approach,the Reinforcement Learning. In this approach, an agent in a given environment (inputsetup), reflects the user interest to optimize a certain value function. The agentimproves his actions due to external rewards with respect to his actions. Over time heattempts to arrive at an optimal state, namely a set of actions that reward him most,thereby guaranteeing optimality. 56 eferences [1] Beyer, H. G. (2001). The theory of evolution strategies. Springer Science[2] Paul mcgregor (2006). Relative Minimums and Maximums, Calculus III course, LamarUniversity, Texas[3] Darrell Whitley (1994). A genetic algorithm tutorial. Computer Science Department,Colorado State University. Fort Collins, CO 80523, USA[4] Runhe Huang (1995). Evolving Prototype Rules and Genetic algorithm in aCombustion Control. 1995 IEEE-IAS, International Conference on IndustrialAutomation and Control Conference.[5] B. Danielson J. & Foster D. Frincke (1998). Using Genetic Algorithms to Breed aCombustion Engine. IEEE World Congress on Computational Intelligence, 1998,Anchorage, Alaska, USA[6] Wolfgana Polifke, Weiqun Geng & Klaus Dobbeling (1998). Optimization of RateCoefficients for Simplified Reaction Mechanisms with Genetic Algorithms. Combustionand Flame, 113(1/2), 119–134.[7] S.D. Harris, Elliott, L., Ingham, D. B., M. Pourkashanian & C. W. Wilson, (2000). Theoptimisation of reaction rate parameters for chemical kinetic modelling using geneticalgorithms. In ASME Turbo Expo 2002: Power for Land, Sea, and Air (pp. 563-572).[8] G. R. Vossoughi & Siavash Rezazadeh (2005). Optimization of the Calibration for anInternal Combustion Engine Management System Using Multi-Objective GeneticAlgorithms. Evolutionary Computation, 2005. The 2005 IEEE Congress on, Volume: 2[9] C. D. Rose, S. R. Marsland & D. Law (2009). Optimisation of the Gas-ExchangeSystem of Combustion Engines by Genetic Algorithm. 2009 4th InternationalConference on Autonomous Robots and Agents.[10] Shtauber, I. & Greenberg, J.B. (2010), A study of Polydisperse Spray Diffusion Flamesand their Extinction in Co-flow. Final Paper towards M.Sc in Aerospace Engineering[11] Nejra Sikalo, Olaf Hasemann, Christof Schulz, Andreas Kempf & Irenaus Wlokas(2015). A Genetic Algorithm-Based Method for the Optimization of Reduced KineticsMechanisms. International Journal of Chemical Kinetics 475712] Carolyn R. Kaplan, Alp Ozgen & Elaine S. Oran (2017).Chemical-diffusive models forflame acceleration and transition-to-detonation: genetic algorithm and optimisationprocedure. Combustion Theory and Modelling, 2019[13] Hongguang Pan, Weimin Zhong, Zaiying Wanga & Guoxin Wanga (2017).Optimization of industrial boiler combustion control system based on genetic algorithm.Computers and Electrical Engineering 70 (2018) 987–997[14] Jie Liua,b, Biao Maa & Hongbo Zhaoa (2019). Combustion parameters optimizationof a diesel/natural gas dual fuel engine using genetic algorithm. Fuel 260 (2020) 116365[15] Yiding Zhao, Qinghe Wu, Heng Li, Shuhua Ma & Ping He (2019). Optimization ofThermal Efficiency and Unburned Carbon in Fly Ash of Coal-Fired Utility Boiler viaGrey Wolf Optimizer Algorithm. 2010 International Conference on Electrical andControl Engineering[16] Burke, S. P., and T. E. W. Schumann. ”Diffusion flames.” Industrial & EngineeringChemistry 20.10 (1928): 998-1004.[17] Greenberg, J.B., ”The Burke-Schumann Diffusion Flame Revisited-With Fuel SprayInjection”, Combustion and Flame 77, pp. 229-240, (1989).[18] Tambour, Y., A Lagrangian Sectional Approach for Simulating Droplet SizeDistribution of Vaporizing Fuel Sprays in a Turbulent Jet, Combustion and Flame 61Issue 1, pp. 15-28, (1985).[19] Williams, F.A., Phys. Fluids 1, pp. 541-545, (1958). See also ”Combustion Theory”,2nd Edition, The Benjamin/Cummings Publishing: Menlo Park, CA, (1985)[20] Boyd, Stephen P.; Vandenberghe, Lieven (2004). Convex Optimization page 143 (pdf).Cambridge University Press. p. 129. ISBN 978-0-521-83378-3.[21] V. Jeyakumar; Alexander M. Rubinov (9 March 2006). Continuous Optimization:Current Trends and Modern Applications. Springer Science & Business Media. ISBN978-0-387-26771-5. Continuous Optimization[22] Javier Larrosa, Albert Oliveras, Enric Rodrıguez-Carbonell (2019), CombinatorialProblem Solving (CPS). Mixed Integer Linear Programming[23] Sakawa M. (2002) Genetic Algorithms for Integer Programming. In: GeneticAlgorithms and Fuzzy Multiobjective Optimization. Operations Research / ComputerScience Interfaces Series, vol 14. Springer, Boston, MA.5824] Matlab Help Center, Solving Mixed Integer GA Optimization Problems [link]. Basedon Deb, K. (2000). An efficient constraint handling method for genetic algorithms.Computer methods in applied mechanics and engineering, 186(2-4), 311-338 [link].[25] Darrel Whitley (1997), A Genetic Algorithm Tutorial, Computer Science Department,Colorado State University, Fort Collins.[26] Patankar, S. (2018). Numerical heat transfer and fluid flow. Taylor & Francis[27] Linan, A. (1974). The asymptotic structure of counterflow diffusion flames for largeactivation energies. Acta Astronautica, 1(7-8), 1007-1039.59 ppendices Appendix A - Extremum definition Consider the following 2D [ A ] differentiable function f ( x, y ) , that satisfies f : X → R :Any point within the domain ( X ∈ R ) can be calculated by a desired function, and thenprojected onto a new space (cid:0) x, y, f ( x, y ) (cid:1) ∈ R . The global maximum of such function f is the point (cid:0) c, d, f ( c, d ) (cid:1) , if there exists some region surrounding ( c, d ) for which [ A ] : f ( x, y ) ≤ f ( c, d ) ∀ ( x, y )A function of two variables f has a critical point at ¯ x = ( c, d ) if : f ( c, d ) x = 0 and f ( c, d ) y = 0Such that any critical point ( ¯ x ∈ ¯ x ) within the domain will satisfy : ∇ f (¯ x ) = 0The global extremum is said to be point - ¯ x ∗ , which satisfies one of the following :min x f ( X ) : f (¯ x ∗ ) ≤ f (¯ x )max x f ( X ) : f (¯ x ∗ ) ≥ f (¯ x )60 -dimensional example Given a more general case of a several real variables function that associates an arbitraryn-dimensional point of ¯ x ∈ R n to f : X → R . For illustration, a 3D function [ A ] : red - domain purple - imageThese points can be calculated by an appropriate function such as f ( ¯ x (cid:122) (cid:125)(cid:124) (cid:123) x , . . . , x n ), and thenprojected onto (cid:0) ¯ x, f (¯ x (cid:1) ∈ R n +1 space. Derivation of f provides a system of n equations : ∇ f = (cid:16) ∂f∂x , . . . , ∂f∂x n (cid:17) = ¯0Whereas either here, ¯ x is the set of critical points that nullify the gradient : ∇ f (¯ x ) = 0Similarly to before, the global extrema would yield :min x : f (¯ x ∗ ) ≤ f (¯ x ) ⇔ max x : f (¯ x ∗ ) ≥ f (¯ x ) References • Boris P., Baranenkov, G. Moscow University (1964). Problems in mathematical analysis. • Stewart, James (2008). Calculus: Early Transcendentals (6th ed.). Brooks/Cole. • Craig A. Tovery, Georgia Institute of Technology (2010). Multidimensional Optimization61 ppendix B - Convexity vs. concavity Definition B.1 - A real-valued function f is considered convex if the line segmentbetween any two points on graph of the function f ( x ) lies above or on the graph. Definition B.2 - A subset C is considered convex set if, with any two points, it containsthe whole line segment that joins them.Let X be a convex set in a vector space and let f : X → R be a function. f is convex if : ∀ x , x ∈ X ∀ t ∈ [0 , 1] : f ( tx + (1 − t ) x ) ≤ tf ( x ) + (1 − t ) f ( x ) Definition B.3 - f is said to be concave if ( − f ) is convex (negative of a convex function).Knowing the function’s type dictates the objective function’s type (minimize / maximize). Optimization type An optimization problem is said to be convex if its objective function is a convexfunction [ ] and the constraint set is a convex set [ ] . The optimization task is to find theglobal extremum (convex = minimum, concave = maximum), and is expressed typically as :minimize x f ( x )s . t x ∈ C g i ( x ) ≤ , i = 1 , . . . , mh i ( x ) = 0 , i = 1 , . . . , p, Where f and C := { g , . . . , g m , f , . . . , f p } are convex.The important characteristic we get is that the local minimum is a global minimum. Definition B.4 - C is said to be a non-convex set in the presence of at least one singleuncontained convex combination. Analogously, non-concave set are with uncontainedconcave combination. Conventionally, both terms are referred as non-convex.An optimization problem that violates either one of these conditions [ B.1 ] , [ B.2 ] , i.e. utilizes anon-convex objective function, or a non-convex constraint set [ B.4 ] , is considered as a non-convex optimization problem. 62 non-convex optimization may have multiple locally optimal points such that identifyingwhether the problem has no solution or if the solution is global, is a challenge per se. References • Meyer, Robert (1970), The validity of a family of optimization methods [link] • Z˘alinescu, C. (2002). Convex analysis in general vector spaces [link] • Kjeldsen, Tinne Hoff. (2006), Convexity and Mathematical Programming [link] • Shashi Kant Mishra (2011), Topics and applications in Nonconvex Optimization [link] Appendix C - Code Access As part of my worldview, not only ideas should be accessible to everyone who desires, butalso their implementation means. Optimally, alongside a clear installation instructions andan execution scheme. The world today is flat, very much thanks to the open source cultureat both academy and industry.Aside from global access, it promotes a transparency climates where the laymen is capableof challenging the biggest researchers by validating their results, and check whether itcorresponds to their publications. Using the most popular platform, I would like to referthe reader to GitHub.com website where the project’s contents can be found : https://github.com/Daniboy370/Masters-Project - f inf in