MATHEMATICAL MODELING AND IDENTIFICATION OF THE MATHEMATICAL MODEL PARAMETERS OF DIESEL FUEL INJECTION SYSTEMS

The aim of the paper is to highlight, if not all, then at least the most important benefits of the application of mathematical modeling in the development of a new solution or in the process of improvement and the optimization of an existing product. The thesis presents the basic principles of the model creation without providing a detailed mathematical interpretation of the process; therefore, the model includes a substantial number of unknown or approximately known parameters and it is not possible to create a model of the transparent cybernetic box type. In order to achieve a synthesis of the model and the experimental results, it is necessary to apply the optimization methods for the identification of unknown model parameters. The thesis presents the experimental results along with the results of modeling and identification as well as the analysis of the process of identification of unknown model parameters.


Introduction
Te fuel injection system for diesel engines is a mechanical-hydraulic system in which the processes happen very fast and in which very high fuel pressure develops -with motors of the latest generation even over 2,000 bar.
The fuel injection process itself lasts for less than 2 milliseconds and it can be compared only with the duration of a gunshot.If we add to it the fact that the system should ensure the repeatability of the process during the working life of the engine of 5,000 and even 10,000 hours, with about 1,000 cycles per each cylinder (i.e.4,000 cycles for the entire motor) during each minute in a high-speed four-cylinder four-stroke diesel engine, it becomes clear how much inventiveness and knowledge is necessary for the design and creation of such a system.
In order to critically analyze and evaluate various impacts, the structural, functional, as well as various external disturbances, a good knowledge of the physics of the processes within the system is necessary.The conditions for defining the relations between the input and output parameters and the impacts of the disruption on the operation of the system are thus created.There are two methods available -the empirical one and the theoretical one.
The first method provides reliable information in terms of perceiving the integral, macro characteristics of the process, since the experimental methods do not allow access to all physical values of the process.Evaluation of the physical quantities that cannot be determined experimentally requires the introduction of assumptions, which raises a question of reliability of conclusions obtained in such a way.It is even possible to make wrong judgments.And when it comes to the experimental method, if we take into consideration the subjective characteristics of researchers such as experience, intuition and talent, the process of identifying the most optimal solution from the set of possible solutions is expensive and time-consuming.
Another procedure, analytical one, implies the creation of an idealized mathematical model for the simulation of the mechanical and hydrodynamic processes in the system.If the mathematical model fully imitated the real system, preconditions would be created for the analysis of the process both on the macro and the micro level and in the aspects that are unreachable by experimental methods.
Experience shows that, due to the complexity of the process in the system, it is most rational to combine both methods.
In the following section, an example will be presented of mathematical modeling of the pump -pipe -injector system, with an emphasis on the process of identification of unknown model parameters by applying the optimization methods.This system, in different varieties, has been long in use with diesel engines and there are a large number of publications devoted to the mathematical interpretation of the process within the system and the methods of model solving.The process proposed here has been im-plemented for practical purposes, applied in several Research & Development projects, and the results have been implemented in materialized systems.
The whole "package" also includes the process of the system optimization; however, it will not be discussed here due to the limited space.

Mathematical modeling of the process
There is a certain disagreement about the definition of the model and the role of the modeling of physical processes and systems.Therefore, due toits comprehensiveness, the often quoted definition of M. Pešelis given here: "The term 'model' is used in the broadest sense as a structure for storing knowledge, in which the leading role is played by the idea of reduction of information.The reduction of information is in the closest possible way related to the establishment of relations between different groups of information, and therefore to the issues of acquainting with the structure and function of the system.Therefore, the development of the models represents an important aid to obtain new information on the basis of already existing special knowledge" (Jankov, 1984, p.20).
The end result of the development and implementation of the mathematical model is to obtain reliable information on the basis of which researchers will be able to choose the best solution in the given circumstances, having in mind certain limitations.The very procedure of the construction and application of the model consists of several steps: 1.The formulation of the problem.
2. The collection of the data on the system and the identification of external influences in case the system interacts with other systems.
3. The development of a mathematical model using mathematical methods.
4. The analysis of the impact of certain (according to researchersimportant) parameters of the model parameters on the system processes.
5. The identification of unknown model parameters, in case the model is in the form of a cybernetic semi-transparent box.
6.The optimization of the system -application of the optimization method on the identified model.
7. The evaluation of the obtained results, in which available experimental results are of particular importance as well as the experience and knowledge of researchers.The role of researchers is indispensable because the model cannot provide all the answers.8.The last phase is the implementation of the results into the process of designing a new or modifying the existing solution.
When creating a mathematical model of a complex system, it may happen that the mathematical interpretation of physical processes is superfluously cumbersome or that certain processes cannot be described analytically.In this case, simplifications must be introduced under certain assumptions which should be taken into account when making decisions at the stage of evaluation and implementation of results.On the other hand, the model provides information complex to the extent of how complexly the physical process is mathematically described.
It has already been pointed out that there are a number of published papers devoted to the problems of mathematical modeling of fuel injection systems, and it is up to researchers to apply the already developed model with or without alterations or by respecting the physical laws to build their own solution.Therefore, the thesis will expose only the basic principles common to all the models analyzed by the author, without introducing details.
In line with the form, the system can be divided into components.Their number depends on the configuration of the system.The elementary system consists of the piston mechanism, a working space in the pump's cylinder, a delivery valve, a pipe, a chamber in the injector and the injector needle with all the moving parts.Control (regulatory) elements on the pump or the nozzle can also be found.Fuel is a separate component, the physical characteristics of which, in real conditions and according to the experience of the author based on experimental results, deviate significantly from the values obtained in stationary conditions.
In the mathematical model used, a separate software module is formed for each component and mutual influences are defined.There are a number of reasons for dividing the system into components, and only two will be mentioned.First, it is easier to monitor and understand the semantics, i.e. the contents and the meaning of the model, as well as the limitations and conditions for the similarity of the modeling results and the characteristics of the real system.Secondly, it allows an easy application of the model to other systems of various forms by replacing the existing modules or adding new ones.
Further on, only the basics of the creation of the mathematical interpretation of physical processes in the system will be displayed.Displaying all the details of only one of possible options requires significantly more space.
The first assumption -the system can be extracted from the engine, i.e. a higher-level system, and considered separately while defining and respecting the impact of the engine, but also that of the environment.
For all the chambers in certain parts of the system, the expressions for determining pressure and other fuel physical properties are obtained from mass balance expressed though mass flows, which, assuming that fuel is compressible, homogeneous and isentropic and that pressure is equal in the whole area when mass is equal to the result of multiplying the density and volume m = ρV, reads as follows: And from the equation of state in the form: is the balance of the mass flow through the inlet and outlet orifices for communication of the control space with the adjacent parts of the system, аnd / g dm dt represents undesired fuel loss from the monitored space.Losses occur with fuel flowing through technological gaps in the pump and the injector.It should be emphasized that the fuel density and the control volume value are not constant.Density depends on pressure, temperature and the concentration of the gaseous phase, the presence of which, according to empirical data, cannot be ignored.Volume is variable for most of the mentioned spaces and depends on the current position of the system moving parts.The change in volume caused by the elastic deformation of the parts is practically negligible.
If we define the fuel compressibility via the so-called compressibility coefficient , the expression for determining the pressure in one of the system chambers will be formed as: where the first Sum in the parentheses is the balance of the fuel flow at the boundaries of the control space with the volume V, and the second Sum is the balance of the shifts of space boundaries (the piston, the valve, etc.).
It should be particularly emphasized that the engine, directly linked to the analyzed system, is a cyclic machine by the nature of the work process.In the time period that includes a multiple number of cycles, it is an imperative to have the repeatability of the cycle.This means the repeatability of the angular velocity from cycle to cycle, i.e. the equality of the mean angular velocity of the engine crankshaft, i.e. the angular velocity of the pump's camshaft of the observed system from cycle to cycle.During one working cycle, angular speed is variable; therefore, in order to understand the process, it is more convenient to regard the fuel injection system as a function of the camshaft position (as a reference shaft), and not as a function of time.This is achieved by introducing the following expression: where φ is the angular position of the reference shaft in relation to a certain reference position, and ω is a current value of the angular velocity of the shaft in the position φ.
The modeling of the fluid flow through the inlet and outlet orifices, generally irregular in shape, would significantly burden the model.Practice shows that, with sufficient accuracy, Bernoulli's expression for the flow ratecan be used with the introduction of a correction, the flow coefficient μ of a constant value, with which the theoretically idealized results approach the real process.Thus, the volume flow Qthrough the geometric flow cross section Acan be determined from the following expression, with a difference in pressure Δp: The form of a gap between the cylindrical surfaces through which losses occur is unknown because the straightness and roundness of the parts' surfaces are of the same order as the nominal gap.Neither is known any eccentricity of the parts, so it is practically impossible to make a mathematical model to calculate the loss of working fluid.In literature, there are a lot of semi-empirical expressions, out of which, based on the experimental results, the expressions of Astakhov (Astakhov, 1972) simulate the real process the best.
The position, i.e. the movement of freely moving parts of the system of mass m (the valve, the absorber, the injector needle, etc.) is determined from the equality of inertial forces and forces F i which act on the movable part: Friction forces as a constant value and damping forces linearly dependent on the movement speed should also be included into the sum of the forces.Also, from the standpoint of the stability of the model integration process, the following has been proven useful to include in the model: deformation of the system elements restricting the movement of the parts in the boundary positions, as well as the friction forces during the movement in the deformation zone.This eliminates the discontinuity of the calculated values and increases the stability of the integration process.
The flow in the pipe is three-dimensional and unsteady, with variable physical characteristics of the fuel.Mathematical expressions to describe such flow are very complex and unsuitable for practical use.It has turned out that, with certain simplifications, the expressions of the mass and impulse conservation law can be transformed into a form suitable for solving while preserving the process physical aspect.If we introduce the following assumptions: -The pipecross section is constant and multiple times smaller than the length, so it can be adopted that the flow is one-dimensional and linear; -Disturbances are transmitted at the speed of sound by waves perpendicular to the movement direction -the pipe; -The change of the working fluidstate is isentropic and the speed of sound can be calculated from Laplace's expression 2 / a p     -The fuel characteristics depend only on pressure, and -Gravitation is negligible, the expressions acquire a simpler form: According to Darsy and Weisbach, the hydraulic friction force at a unit volume (the last item in the second expression) is a non-linear function of the fluid velocity, which is not convenient from the point of integration.The expression can be linearized by introducing the coefficient of hydraulic resistance K that linearly depends on the flow velocity; however, it is adopted that it has a constant value in a single integration step and is calculated for each step on the basis of the mean flow velocity.Thus, the previous expressions obtain the following form: And finally, fuel is regarded as a separate component; its physical characteristics in stationary conditions depend on pressure and temperature and can be determined with the help of published empirical expressions (Huber, Schafitz, 1966).However, the experimental results obtained in real systems have significantly different values.
This is most apparent with the speed of propagation of disturbances in the system pipes, the value of which, according to experimental results, reaches a significantly lower value than the value obtained in stationary conditions.The reason for this is the presence of the gaseous phase in the fuel which is uniformly dispersed in the liquid phase in the form of tiny bubbles (Murzin, 1990).On the basis of a specific or adopted (presumed) gaseous phase concentration value, it is possible to determine the fuel characteristics as a function of pressure and temperature (Fox, 1977), (Thiemann, 1988).Finally, we can conclude that the fuel characteristics, density, compressibility and also the speed of propagation of disturbances depend on the pressure, temperature and concentration of the gaseous phase.

Solving the mathematical model
The mathematical model consists of a set of two partial nonlinear hyperbolic differential equations and a set of non-homogeneous nonlinear differential Cauchy -type equations of the first or second order, the number of which depends on the number of components.There is also a set of algebraic expressions and a set of restrictions and logical conditions that alter the system structure depending on the process phases.The first set of differential equations describes the fluid movement in the flow plane/coordinate system consisting of the pipe longitudinal axis and the time coordinate.The other set describes boundary conditions at the pipe end.
The fluid flow can be viewed as a wave phenomenon -as a flow induced by the propagation of waves of finite intensity, which at each point of the coordinate system determine the fuel condition.In this sense, the differential equations that describe the flow in the pipe are transformed intoa so-called "characteristic" form (Fomin, 1973) which determines the change in fluid pressure and speed along the characteristics -the trajectories of waves at the coordinate system.
They are solved using the so-called method of characteristics, but in order to apply numerical methods in the expressions, finite differences should be used instead of differentials.Consequently, the coordinate system should be subdivided by the rectangular network Δх-Δt, where Δх is a step along the pipe, and Δta step along the timeline axis, provided that the characteristic does not cross over from one field to another, i.e. that the requirement (11) in which the right members of the equations represent the average value of the sum, i.e. the differences of flow velocity and the wave propagation speed within one integration step Δt.Along the characteristics, the changes in pressure and flow velocity are calculated, thus determining their values at all nodes at the "level"   1 t j t    , based on the already calculated values in all the nodes along the pipe for the "level" t j t   .
This procedure must be synchronized with the integration of the equations of the boundary conditions which determine the pressure at the beginning and at the end of the pipe.
The second method consists of calculating Riemann's invariants (Hanjalić, 1978) along the characteristics out of which the velocity and pressure in the flow field are calculated.Both methods provide the same results with approximately the same computing time.
When selecting numerical methods for solving differential equations of boundary conditions, it should be kept in mind that the model changes the structure and that the gradients of physical values of the fuel, primarily pressure and speed, have extremely high values.Experience has shown that, due to the integration process stability, it is more convenient to use only the methods that do not require the calculations of higher order differentials.Such is, for example, the Euler-Cauchy's method, a so-called predictor-corrector method with a variable step and an iterative process of correction.The integration step varies depending on the calculated error as a difference between the results of two successive iterations in the same integration step.If, for a limited number of iterations (not more than fi-ve) the error is greater than the adopted maximum value, the step should be reduced, and if it is lower than the lowest adopted limit value, the next step should be increased.The limit values are determined for each of the calculated values so as to obtain sufficient accuracy within a reasonable computing time.
What is provided by the mathematical model.If we interpret it as a set of rules with which the set of input parameters is mapped onto a set of output parameters, then a change of one of the input data reflects onto some or all of the output data.This allows the implementation of systematic research -a systematic changing of one by one input data and determining the degree of their impact on the process in the system, by comparison and qualitative and quantitative evaluation of representative, output parameters.The scope of input parameters changes is determined in accordance with the model structure, physical meaning of a value and researcher's experience.
Although the mathematical model is exposed only in its basic form, without detailed mathematical interpretations of the process, it can be concluded that it contains a variety of values that are unknown or just approximately known.The model, therefore, has a form of a semitransparent cybernetic box with a set of unknown parameters that significantly affect the process.Also, there is a second group of parameters introduced into the model primarily to improve the stability of the numerical integration procedure.Their influence on the results is not significant.It is highly unlikely that researchers, solely based on their knowledge and experience, can determine the actual values of quantities such as friction and damping forces during the motion of freely moving parts, flow coefficients, gaseous phase concentration in the fuel, etc..The unknown parameters may also be some constructive data that cannot be reached without the destruction of the analyzed system.A model with experimentally determined unknown parameters can be used for a systematic analysis in order to acquire knowledge about the impacts of the input parameters on the process.If experimental results from similar systems are also available, it is possible to define quite reliably a relatively narrow range of unknown parameters so that the model can be used in a design phase as well.
Mathematical modeling is an essential tool for researchers in modifying, improving and optimizing realized systems or in a design of a similar solution where the results of experimental studies of the basic system are available.Experimental results provide an opportunity to identify unknown parameters and apply the modeling results with sufficient reliability in designing a desired solution.

Identification of unknown parameters of the mathematical model
In order for a computational model to fully simulate the process in the system, it is necessary that the difference between a series of modeling results and the corresponding series of experimental results for one or more characteristic process quantities is equal to zero in an ideal case.Let us denote a series of modeling results as a function of the reference shaft position and the vector of model parameters with

 
, is a new unknown parameter with which the measuring results are coordinated with the modeling results regarding the angular position of the reference shaft, and 0 e V is a parameter that should also be determined in the case when the absolute level of the measured value must be corrected.The difference is the total error as a difference between the combination of random and systematic errors in measurement and modeling errors.It cannot be equal to zero either, so the task of identifying unknown parameters comes down to a search for the minimum of this difference.The mathematical model is complex, the number of unknown parameters is large, they are mutually independent and the task cannot be solved by deterministic methods.It remains to solve the problem with statistical methods.We should formulate the error functional in the form of a sum of the squares of differences of modeled and measured (n) values of the characteristic value.
In the expression, x denotes a vector of unknown model parameters, and only one estimation value is included in the sum of squares, although there can be more of them.The problem comes down to searching for the error functional minimum, for which optimization methods are used.At the same time, unknown parameters cannot reach an arbitrary value, but in order to preserve the process physicality, we must define the boundaries within which an optimal value of the parameters can be sought.Any combination of parameters within the established boundaries represents a valid solution, and a set of all permissible solutions represents valid space.The mathematical model is complex, therefore it is expected that, in a set of allowed solutions, a greater number of local functional minima are obtained, and those are obtained within the search for the minimum with different initial values of unknown parameters.The search for an absolute minimum is often cumbersome, and often not even necessary in practice, so finally, the solution that corresponds to one of the local minima of the functional is accepted.
The task would be solved by analytical methods if the model had a form that allows calculating the gradients.This is not possible in this case, so the use of numerical methods is necessary.According to the experience of other researchers, out of a whole range of developed optimization methods for determining the minimum sum of squares (Tomić, 1987), (Vengerskij, et al, 1982), the highest efficiency was demonstrated by Marqurdt's method (Gill, et al, 1981) which will not be discussed further here due to space limitations.
An example of the application of the mathematical model with the identification of unknown model parameters will be shown here together with the experimental results obtained within the project of defining the fuel injection system that would be delivered by an alternative manufacturer, and that would be applied to the existing engine in a serial production as a replacement solution.The project was implemented by the Factory of Engines for Special Purposes in Pale, where the research was carried out.It should be particularly emphasized that the results presented are of a preliminary nature, they were obtained on a system that does not meet the requirements to be applied on the specific engine and they do not contain potentially confidential information for the engine manufacturer.
The experimental studies were carried outside the engine, on a test stand for the purpose of testing the system using all the elements intended to be installed in the engine.Two characteristic process quantities were measured -the fuel pressure at the pipe end and the injector needle position.The following quantities were also measured: the position of the pump camshaft, the fuel pressure and temperature and the mean value of the volumetric flow rate of fluid through the injector, i.e. the cumulative quantity for a selected number of operating cycles.The following measuring devices were used: -Piezoelectric pressure sensor AVL 5QP, -Sensor AVL 3056-A01, -Needle lift inverter -hand-made inductive half bridge, -Sensor AVL 3076-A01, and -Encoder of camshaft angular position: optical incremental The obtained electrical quantities are digitalized using the National Instruments data acquisition system USB 621 and stored in a personal computer.
The obtained results include systematic and random measurement errors.Out of systematic errors, two have been identified and previously discussed -the angular position of the cam profile in relation to the flow of the measured quantities and the absolute level of the measured signal of fluid pressure as a result of operating principles of piezoelectric sensors.Both errors are effectively eliminated by the procedure of identification of unknown parameters.On the other hand, random errors are of stochastic nature and they are neither possible to identify nor to eliminate.However, when random errors are subjected to normal distribution with the zero mean value, then their level becomes significantly reduced by calculating the average cycle obtained by averaging the characteristic quantity values (measured over a number of consecutive cycles, for each angular shaft position).In this case, the influence of random errors can be disregarded (Tomić, 1987).Figs. 1 and 2 display the pressure at the pipe endaver aged for 50 consecutive cycles as well as the statistical parameters -standard deviation and Pearson's coefficients.It should be noted that the measured pressure value between two consecutive injections is equal to zero, which is a result of the fact that the measurement was performed with a piezoelectric sensor.That is why this quantity is of one of the unknown parameters that should be determined.
The value of the asymmetry coefficient is close to zero and the shape coefficient varies around a value of 3 (from 2 to 4) except in the process phases where the gradient values of the measured quantity are extremely high.This can be explained by anerror in measuring the camshaft angular position.in the process phases where the gradient of the measured quantity ishigh (up to 350bar per angle), a measuring error of one angle increment (0.1 degree) results in aconsiderable error.The value of Pearson's coefficients shows that random errors are subject to normal distribution; consequently, their impact on the measured quantity can be neglected if the values measured over a number of consecutive cycles are averaged.It has already been emphasized that the identification process must be repeated multiple times with different initial values of the unknown parameters, i.e. more local minima of the functional should be found, and then one chosen as a final solution.The following illustrations display in parallel: the results of the measurements and the results of modeling with the initial and the "exact" (identified) values of the unknown parameters for two operation modes determined by the total flow through the nozzle Vc during one cycle and by the camshaft rotational speed n.In order to illustrate the convergence of the suggested method, in both cases a large number of parameters are identified -10.For the same reason, the illustration shows the results of "computing" processes, with which the initial values of certain model parameters are significantly different from the "true" ones, which resulted in a large difference between the results of the identified model and the results obtained with the initial parameter values.The following diagrams show the course of the change of the functional and several unknown model parameters during the identification process.In the first case (Fig. 5), which corresponds to the results of Fig. 3, the allowed space in which we seek a solution is narrowed.The result is a high value of the functional and reaching the threshold values of model parameters.At first glance (Fig. 3) the results of measurements and modeling match relatively well, so a question arises whether the threshold values of other model parameters are established correctly, i.e. whether the identified model represents a real system.This example shows that, when assessing the modeling results, researchers must analyze in detail the whole process, not only the modeled and experimental values of the selected characteristic process quantities.In the second example (Figs. 4 and 6), the threshold values of the parameters are determined taking into account all the known system characteristics; therefore, the final value of the functional is much lower, and the model parameters are within the permitted area.
The comparison ofthe experimental and calculated values (Fig. 4) shows that the model appropriately simulates the operational process and that the results can be used for system design and optimization.In addition, the differences observed during the change of the analyzed value indicate that the model does not fully describe the real physical process in the system, but that it can and should be improved in further research.

Conclusion
Mathematical modeling, with the application of optimization methods for identifying model parameters, enables connecting theoretical and empirical relations, on the one hand, with experimental results, on the other.This provides a number of benefits to researchers.
-It is possible to also gain insight into the process phases in accessible by experimental methods.-The relationship between input and output values can be established at minimum cost -without cumbersome and detailed experimental research.-It is possible to carry out a systematic analysis of all factors influencing the course of hydro-mechanical processes within the system and therefore to acquire relevant new knowledge.-The mathematical interpretation of the process allows the use of optimization methods and effective identification of possible solutions leading to the choice of the best one.Since the processes within the system are complex, a detailed mathematical interpretation leads to the creation of an unnecessarily cumbersome model; therefore, some necessary simplifications that do not jeopardize the physicality of the process have proved to be justified in practice.This approach provides a model that appropriately simulates a real physical process.However, the model cannot be integrated by analytical methods and it is necessary to create computer programs along with the application of numerical methods.
Finally, it is possible to realize a theoretical-experimental procedure that takes the shortest route to the objective -to find the optimal solution and implement it in the development project.
and a corresponding series of the experimental results of the same characteristic values of the process with the index i represents the position of values of the results in a corresponding sequence.The phase shift 0 cost  