Aviation by pengxuebo


Session 5. Aviation


                                                Jonas Butkevičius

                                     Vilnius Gediminas Technical University
                                   Plytinės g. 27, Vilnius, LT-10105, Lithuania
                            Phone: (850)2745099. E-mail:Jovas.Butkevicius@ti.vgtu.lt

        The entry of Lithuania into European Union has changed the situation in macroeconomics. It has
influenced the changing situation of competition and trade including air transport.
        In this article the liberation of air transport market, development of air fleet, increasing coverage of flight,
enlargement the geography of routes, development of charter flight, change of airports and air navigation
system‘s work, consolidation in deference of unsure aviacompanies, change in opinion of specialists, influence
to attract new aviacompanies, and consolidation competition is analysed.
Keywords: air transport, market, liberation, flight park, passenger transportation


      After May 1, 2004, when Lithuania became the Member State of the European Union (EU), the
macroeconomic environment of the country has changed [1]. This caused changes in competition and
business situation of the entire national economy, including transport sector as well [2]. Therefore it is
important to identify the effect of joining the EU on Lithuanian transport system.
      This article analyses how Lithuania‘s entry to the EU has affected the national air transport.


       On April 1, 1997 all restrictions on entering the local market (cabotage) were repealed in the
European Union. From that time the air transport markets of EU Member States are considered to be
totally liberalised. Principal features of the liberalisation are as follows:
          − open market enabling flights between any airport within the European Union;
          − common licensing standards for carriers in all EU countries;
          − free rating of tariffs;
          − equal rules applied for regulation of operation of regular and charter airlines.
       After liberalisation of air transport market airlines received equal conditions for fair competition
in the EU market, without restrictions of frequency of flights, costs and route network. Consequently,
in the struggle for attracting the clients it is not only necessary to be competent in the field of tariffs,
but also to regard such parameters as the frequency of flights, convenience of arrival and departure
slots, quality of services supplied to passengers and the network of routes.
       After joining the EU, Lithuanian airlines are also directly influenced by the liberalisation of
aviation market. The impact has produced positive and negative effects.
       Among positive effects the fact that national airlines can open their flights to any airport of the
EU without bilateral agreements, without parity evidence, and without other bureaucratic restrictions
can be mentioned. The geography of flights is constantly expanding thus influencing the growth of
passenger flows and the increase of profits.
       Among other important effects of market liberalisation on national airlines are the following:
          – EU carriers are also able to choose airports of our country without limits, therefore new
             companies, such as “Air Baltic”, entered our market, and this boosted the competition for
             Lithuanian companies;
          – the increased competition reduced the costs of services, particularly because of the entry of
             “cheap flights” airlines, thus causing the decrease of airlines‘ profits.

         The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

        Liberalisation of airline companies has a positive influence on Lithuanian citizens’ travelling by
air, as the ticket prices go down, and the geography of flights from national airports has expanded, and
consequently airlines are forced to improve the quality of services, etc.
        Liberalisation of aviation market also has a positive influence on national airports because the
handling of the airports and the increasing flows of passengers bring additional profits to these airports.


       A positive result of Lithuania’s entry to the EU is the fact that airlines has to renew their aircraft
according to EU standards, and this influences directly the quality of services delivered to passengers,
improves the environmental and flight safety. For airlines the purchase (or leasing) of new aircraft is
related to considerable additional expenses.
       Lithuanian aircraft fleet has most intensely changed by the airline companies “Lietuvos
avialinijos” (“Lithuanian Airlines”, LAL, present name: AB “flyLAL-Lithuanian Airlines”) and
“Aurela”. “Lithuanian Airlines” have disposed of 2 airplanes Boeing 737-200, which do not meet the
EU standards of noise emission. “Aurela” has disposed of 2 airplanes JAK-42 because of the same
reasons. In 2006 “Lithuanian Airlines ” operated 5 airplanes of Boeing 737-500 type and 4 airplanes
of SAAB-2000 type (in 2003 it owned 3 of the type) the noise emission of which meets the EU
standards of noise emission. “Aurela” disposed of Soviet-type airplanes JAK-42 and purchased 1
airplane of Boeing-737-300 and 1 airplane of HS 125-700 B type. The airplane Boeing 737-700
owned by “Aurela” makes charter flights to non-EU countries. The airline “Aviavilsa” will have to
change 2 airplanes of AN-26B type for the shipment of DHL because of the lack of certificate in line
with EASA PART 145.


       The amount of serviced passengers in national airports have been increasing constantly – if in
2000 services are delivered to 580.043 thousand passengers, so in 2001 the number grows up to
649.968 thousand passengers (growth rate 12.1%), and in 2002 – 700.853 thousand passengers
(growth rate 7.8%), in 2003 – 788.248 thousand passengers (growth rate 12.5%). Particularly
extensive growth starts after Lithuania’s entry to the EU: in the year 2004 the amount of passengers
reached 1112.916 thousand passengers (growth rate 41.2%), and in 2005 the amount even reached
1453.222 thousand (growth rate 30.6%). If compared with the years 2003 and 2005, the amounts of
passengers increase actually by 84.4%. The main reason of such a growth of passenger amounts is the
fact that Lithuania joined the EU, and this enabled travelling to EU countries without visas, also the
extended geography and amount of flights, etc.
       In 2005 the largest amount of passengers is as follows to and fro the United Kingdom – 15.3%,
Germany – 12.8%, Denmark – 12.3% and Ireland – 6.6%.


       In 2005 from Lithuanian airports 39282 flights has been made by 16% more than in 2004. In
2005 from Vilnius airport 29193 flights has been made (i. e. by 23% more than in 2004), from Kaunas
airport 4611 flights has been made (decreased by 5%) and from Palanga airport 5478 flights has been
made (growth rate 4%). The growth rate of flights is caused by the fact that Lithuania has joined the EU.
       From 2003 to 2005 total quantity of flights from Lithuanian airports increases by 41.1%, of
which from Vilnius airport grew by 59.2%, from Kaunas airport increases by 13.1%, from Palanga
airport increases by 3.7 %.
       During the year 2005 Lithuanian airlines have carried 674.3 thousand passengers – it is by 14%
more than during the year 2004. This is also the result of Lithuania‘s entry to the European Union.


      The development of new flight routes in air transport is an important fact of Lithuania‘s entry to
the EU.

Session 5. Aviation

       When Lithuania joined the EU, the market became liberalised, and airlines received an
opportunity to choose flight routes by themselves. In 2004, from Vilnius airport new destination
flights have been opened to six towns: Milan (LAL), Dublin (LAL and “Air Baltic”), Cologne (“Air
Baltic”), Hamburg (“Air Baltic”), Munich (“Air Baltic”) and Oslo (“Air Baltic”).
       Competing with other airlines, the company “Lietuvos avialinijos” (LAL, present name: AB
“flyLAL-Lithuanian Airlines”) is prepared to open 11 new flight routes through the airport of Palanga:
to Kaliningrad, Latvia, Germany and Scandinavia. Furthermore, the LAL airlines are prepared to open
new flight routes from Vilnius to Dubrovnik in Croatia, Thessaloniki in Greece, Istanbul and Antalya
in Turkey.
       “Air Baltic” is also planning the development of flights from Vilnius airport.


       Former system of bilateral agreements prevents airlines from free development of transportation
for meeting the growing demands for travelling by air transport. It is particularly obvious in the market
of charter flights where the demand is rapidly growing, and the capacities of foreign airlines to make
flights from Lithuania are narrow.
       The market of charter flights is liberalised by Lithuania joining the EU. This stimulated the
opening of new charter routes and entry of new airlines to the market of charter flights.
       In 1995 the passenger flows by charter flights make 11.2 % of the total amount of air transport
passengers. In 2000 passengers of charter flights make up to 12.2%, and in 2004 – 16.6%. In 2004
Lithuanian airlines, as well as those from abroad (Egypt, Spain, Tunisia, Ukraine), make 2998 charter
flights and carry 166983 passengers. Besides, 1986 passengers are carried by local flights. In 2004,
from Vilnius airport two new airlines started their charter flights – “Joanos avialinijos” and
“Gintarinės avialinijos”. New direct charter flights to Spanish resorts in Gran Canaria, Tenerife and
Palma de Mallorca were started in 2004.


       Lithuania‘s entry to the EU and liberalisation of aviation market make favourable conditions for
airlines with cheap flights to come into the national market. Thus it causes a strong competition to
national airlines.
       In 2005 “cheap flights” airlines make 24% of inland regular flights within Europe. They deliver
services for 900 flight routes linking more than 200 towns. They own over 250 airplanes and carriy
almost 60 million passengers per year (i.e. approximately 16% of the inland passengers of Europe).
       As mentioned before, the entry of “cheap flights” airlines into the national market caused a
strong competition for national airline companies. In 2004 “Air Baltic” carried 96817 passengers that
made 9.7% of passengers of Vilnius airport, and LAL carried 430417 passengers that made 46.3%.
However forecasts of Vilnius International Airport show that in 2008 “Air Baltic” will already carry
535.2 thousand passengers, and LAL (present name: AB “flyLAL-Lithuanian Airlines”) will carry
696.6 thousand passengers. Thus “Air Baltic” will obtain even bigger share of the market and will
pose even stronger competition to LAL.


       Air navigation services in Lithuania are delivered by State Enterprise “Oro navigacija”. It
should be noted that after the entry to the EU the amount of air navigation services grew considerably
in Lithuania: if in 2003 – 86805 flight services were delivered, so in 2004 their amounts grew by
25.4% and made 108875 flights. Furthermore, in 2005 the amount of flight services delivered
increased yet by 12.2% and made 122209 deliveries of flights. In the period of 2003-2005 the amount
of serviced flights grew even by 40.8%. In 2005 – 68% airplanes flew by transit, during the year the
amount of transit flights grew by 9.2%.
       Regarding the SE “Oro navigacija”, principal positive effects of Lithuania‘s entry to the EU are
the following:
        – disappearance of customs duties in purchasing of equipment from EU countries;

         The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

       – rise of factors enabling lower prices for the delivery of services;
       – emerging opportunities for using financial support from EU Structural funds.
      An important change concerning the SE “Oro navigacija” is the fact that the implementation of
the “Single European Sky” concept foresees the implementation of functional air space components
regarding the procedural airplane demands, regardless of national borders of countries.


       An important result of Lithuania‘s entry to the EU is the factor that business stakeholders have
to consider the level of labour force – this is underlined in the questionnaires and interviews of all
specialists employed in civil aviation.
       The problem of pilots is one of the most important subjects. In Lithuania pilots are qualified in
line with EU standards. Licences issued by CAA are recognised in all EU countries. Before the
Lithuania‘s EU membership here was a surplus in pilots. However, after joining the EU, a certain
proportion of pilots have been lured away by foreign (especially by Latvian) airlines, where the wages
are more considerable. So at present a lack of pilots is already palpable. The same can be said about
the lack of aviation engineers.
       On January 1, 2006 there are 2302 specialists of civil aviation in Lithuania, they are holders of
2991 issued licences, from which: 201 licences of airlines’ transport pilot, 423 licences – commercial
aviation pilot, 70 licences – flight engineer, 20 licences – navigator, 134 licences – air traffic manager.
       Concerning the migration of pilots national airlines are forced to raise the wages for pilots and
high rank specialists. This, in turn, raises the costs of services of airlines and reduces their profits.
From the viewpoint of the staff this is a positive result as the wages are growing, thus corresponding to
the professional standards of the staff.


       Great advantage of EU membership for Lithuania is a changed attitude to Lithuania as a
country, whereas the EU membership is a factor guaranteeing that each EU Member State is a reliable
country. Therefore the fact of joining the EU and the improvement of the national image brings the
following positive results to the air transport sector:
       – rating of national airlines has improved, thus making easier negotiations on slots, aircraft
handling in airports of other countries, supply of spare parts for aircraft, etc.
       – pilots trained in Lithuania have more opportunities of working in airlines of other countries;
       – EU membership guarantees a higher standard of civil aviation, therefore the flows of tourists
increase in the national airports, and foreign airlines are more encouraged to select flights through
Lithuanian airports, etc.;
       – different contracts with foreign countries in the field of civil aviation are easier to make;
       – no customs duties required for purchasing equipment in the EU countries;
       – other EU countries get more information about Lithuania, our country has become more
attractive and more reliable to other EU countries, etc.
       All the above mentioned factors are basing on the improved national image of Lithuania. This
applies not only to the air transport, but to other modes of transport as well.


       1. EU membership enables Lithuanian airlines to open flights without any restrictions to all
eligible airports of the EU. On the other hand, new airline companies have entered Lithuanian market,
thus boosting the competition for national airlines.
       2. A positive factor is that national airlines has to renovate their aircraft fleet in line with EU
       3. EU membership enabled non-visa entrance to EU countries, therefore the amount of
passengers in Lithuanian airports has grown considerably; if compared with the period of 2003 and
2005, the growth rate reached even 84.4%.

Session 5. Aviation

       4. The proportions of flights from Lithuanian airports have increased significantly, and the
geography of flights has expanded distinctly.
       5. Proportions of charter flights from Lithuanian airports have increased as well.
       6. After Lithuania’s entry to the EU “cheap-flights” airlines come to the national market and
cause strong competition to Lithuanian airline companies.
       7. Joining the EU resulted in the increased proportions of air navigation services in Lithuania.
       8. An important effect of joining the EU is the necessity for business stakeholders to consider
the level of labour force.
       9. EU membership improves the positive image of Lithuania as a State.


1. Butkevičius, J. The effect of Lithuania‘s entry into the European Union on the national transport
   system and the transport system development: manuscript of monograph. Vilnius: Technika,
   pp. 35-118.
2. Butkevičius, J. Development of passenger transportation by railway from Lithuania to European
   States, Transport – 2007, Vol. XXII, No 2, 2007, pp. 73-79.

          The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”


                                 Eugene Kopytov, Vladimir Labendik,
                                  Sergey Yunusov, Alexey Tarasov

                                 Transport and Telecommunication Institute
                                   Lomonosova 1, Riga, LV-1019, Latvia
                        Phone: (+371)7100590. Fax: (+371)7100660. E-mail: tsi@tsi.lv

      The tasks of neural networks’ application in the airplane power plants automatic diagnostic systems are
considered as well as their peculiarities and advantages.
Keywords: aircraft power plant, technical maintenance, status condition, diagnostics system, neural network


       The modern air engine is a complicated technical object that embodies the most progressive
technologies of science and engineering. The provision of needed economic characteristics of the air
engine power plant is possible by the rising the level of thermodynamic cycle operation parameters,
that mostly have an influence on the air engine’s resource and reliability. During the gas turbine
engines operation the various tasks of aircraft power plant vital activities guaranteeing and airplane
systems and energy plants functioning are solved in the onboard computer. Opportune and qualitative
aircraft engine parameters’ monitoring and diagnostics make it possible to realize the exploitation of
the engine based on its status condition. In spite of the big number of various gas turbine engines
diagnostics and monitoring methods, none of them is completely universal one. It can be explained
mostly by a very high complication of the air engine: there are many parameters, links, the processes
are non-linear, there are many modes of operation, etc. All it assumes the need of complex methods’
application for the solving the air-engine’s parameters’ control and diagnostics tasks. The promising
method of air engine managing and conditions’ control is the neural networks’ application [1-3] that
grants the solving of the wide spectrum of tasks.


       The main tasks that can be solved by neural networks application are the following:
              the identification of gas turbine engine mathematical models;
              the classification of aircraft power plant modes of operation;
              the control and diagnostics of aircraft power-plant status condition;
              the analysis of trends of aircraft power-plant status condition parameters, measured
              onboard, aimed for forecast of these parameters variation depending on the life length.
       The usage of gas turbine engines mathematical models is needed for many research and
practical tasks solving, for example, for the synthesis and analysis of the power-plant managing during
the various regimes of flight.
       The ideal mathematical model (MM) must satisfy the following mutually contradictory
              to describe adequately the connections between parameters and the processes in gas
              turbine engine;
              to provide the given precision of the parameters’ calculation;
              to be convenient for usage in calculations and modelling;
              to be adaptable (learnable) for the individual copy of the engine, etc.
        In practice usually the set of mathematical models of different complexity is used. Each model
satisfies the part of requirements given above and has different areas of application.

Session 5. Aviation

       The most precise and complicated models are non-linear gas turbine engine components’
models, but they are used for the solving of research tasks and optimization of the characteristics of
power-plants control system. Simplified non-linear models are widely spread, for example, regression
models as well as the models with variable coefficients. But even using this approach it’s not always
possible to receive the universal models because of their different representation forms and the need of
models’ coefficients correction. Classical multidimensional functions approximation methods do not
allow to realize simple mechanisms of the mathematical models structure choosing. The realization of
classical interpolation methods on the base of spline-functions requires the considerable computing
resources. In such a situation as a rule the providing of the calculations in the real-time is problematic.
       For the air engine throttle control model the altitude-velocity and throttle performances are
described by the authors using the similar modes [4]. Also in the work [4] the simple and precious
mathematical model is developed.
       Today in the works about the power plants management systems designing the approach of
neural network usage for the mathematical model identification is proposed. The main idea of such
approach usage is in the process of network learning that means the adjustment of the big number of
coefficients (synaptic weights between neurons). The building of the neural network models is based
on the standard procedures of the neural network structure and their learning methods. Multilayered
neural network organization makes it possible to carry out the parallel calculations that support the
solving of the characteristics approximation task in the real time [3].
       Thermo-gas parameters of the aircraft engine (temperature, pressure, air consumption, etc.) in
the different sections of the gas-air flow duct as well as mode and exit parameters (the rotor rotating
frequencies, throttle, fuel consumption, etc.) are the carriers of information about the aircraft power
plant condition. That’s why they can be used for the definite classes of the aircraft engine conditions
(properly functioning, non-properly functioning) recognizing with the aid of different mathematical
models of gas turbine engines.
       Using the power plant parameters trends the classes of the air-engine can be recognized by the
uncovering the correlations between measured and calculated using mathematical model parameters.
Statistical characteristics of the controlled parameters registration results (also caused by the appearing
and developing of the defects in the aircraft power plant) make it possible to forecast the changes of
the aircraft power plant condition changes during the exploitation process [5].
       In the case when the number of the measurable parameters is not sufficient for the linear
mathematical model development the neural network can be used as the model for the standard and
defected engine’s model. The analysis of the given parameters’ deviations in the time is carried out by
the calculation of the metrical distance between standard engine’s data basing on the neural network
and data received during the exploitation of the engine. The results of the quantitative modelling are
evidence of complex monitoring and exploitation management possibility of the aircraft power plant
using new neural network technologies. Such methods extend and expand “classical” methods that can
raise the reliability of management and trustworthiness of the parameters control and aircraft engine
diagnostics as well as the decision making processes efficiency on-board while “critical” conditions
are detected.


       The mathematical model of the aircraft engine power plant is the model that generalizes many
local models, for example, models of the flow path of the air engine. The relationship between
elementary models of the physical processes can be described by the graph that makes possible to
formalize the model integrity and coherence research.
       Direct application of the neural network models, for example, for the diagnostic model of the
aircraft power plant’s flow path can be represented as multilayered perceptron with two hidden layers.
Perceptron’s inputs are the controlled parameters of the gas turbine engine (y1, ..., yn), its outputs are
the controlled parameters’ deviations from their nominal values (Δy1 ,.., Δyn ) . The distinctive feature of
the neural network as the diagnostic model of the aircraft power plant is that it’s weight matrix is
formed in the process of two-staged learning procedure basing on the models that are linear on the
local aircraft engine’s parameters.

         The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

       The most simple example of the aircraft power plant mathematical model is the throttle
characteristic of the aircraft engine in two coordinates (GF_coer, nc_coer), where GF_coer is an adjusted
value of fuel consumption in the combustion chamber, nc_coer is an adjusted value of the gas turbine
engine turbo compressor’s rotation frequency that is represented in the relative (non-dimensional)
units. As the rule the powerful deviation on the throttle characteristics’ relations can be observed
because of some disturbing factors (changes of the air engine’s flow path geometrical parameters,
errors of the main engine’s parameters’ and external conditions used for engine’s parameters’
conversion to standard atmospheric conditions measurements errors, etc.). Recent investigations show,
that in the case of modes operation near maximum the measured fuel consumption’s deviation from
the nominal value can be up to 4,5% in the case when such an error can be caused by uniformly
distributed errors of the direct measurements of the temperature and pressure as well as fuel
characteristics’ changes. If tolerance to the parameters’ deviations is less than range explainable by
random factors, then while engine parameters are regulated systematic deviations must be
compensated at the expense of engine parameters’ deviations from the their standard values.
       Investigations carried out show that the most appropriate neural network architecture that can
satisfy prescribed requirements for the aircraft power plant control quality is ensemble neural network.
Such architecture in the given case is represented as the following chain: radial basis function
networks (RBF) → perceptron → Kohonen’s neural network (KN) (see Figure 1).


E                                  RBF2                    PERCEPTRON


                               Figure 1. Ensemble neural network architecture

       In ensemble architecture the first (input) layer filters the measured data, the second
(intermediate) layer represents the thickener of the aircraft power plant status characteristics, and the
output layer is the status conditions classifier.
       The modelling of the given task using PC in the case of aircraft power plant’s work in different
modes and with different combinations of parameters let to denote 3-5 most informative parameters
that maximally affect the normal functioning process of the aircraft engine. In this case the usage of
five RBF networks each with five input parameters (state vector) and two output parameters (state “1”

Session 5. Aviation

– normal, state “0” – failure) is considered as the optimal one. Perceptron in the ensemble is the field
that concentrates five RBF’s outputs. Kohonen’s network (classifier) has two inputs and one output.
This network with a high degree of precision carries out the classification (recognition) of the aircraft
power plant status condition, also in the case of partial or full uncertainty about its parameters.
Computer experiment has shown that ensemble neural network recovered the absent data steady and
determined the condition of the power plant. The estimation of the power plant according to its
exploitation time based on neural network is carried out in the following way. The standard model of
the power plant (parameters of engine) received in the process of factory development testing, which
is stored in the neural network basis as the individual informational “portrait”, is compared to the
power plant parameters during the exploitation [6]. In the process of the comparative analysis the
special metrics is calculated. Its value can be used to estimate air engine status condition and to build
the separating plane. For such task solution the 3-layered perceptron (2 inputs, 14 hidden neurons, 4
outputs) is used. It is trained using the error back-propagation algorithm.


       The application of the neural network apparatus is effective in the solution of many tasks:
aircraft engine power plant’s “image” identification, control, mode classification, diagnostics, trend
analysis, forecasting, etc.
       The problem of aircraft power plant and its subsystems control and diagnostics is a complicated
task that is concerned with the need of taking into account many factors also the uncertainty factors.
The application of the artificial intelligence methods based on the neural networks makes it possible to
find new ways of this problem solving that are based on the using the knowledge and experience of the
experts. These are image recognition theory, learning theory, theory of adaptation to the changing
external conditions, theory of the decision making in the case of the information deficiency, etc.
       Investigations considered in [7, 8] show that the neural networks using for analysis of the
aircraft power plants is effective and promising.


1. Yefimov, V.V., Yakovkin, V.A. Method of technical diagnostics based on Neural Networks /
   Priborostroenie, Vol. 9, 1999. (In Russian)
2. Vasilyev, V.I., Zhernakov, S.V., Urazbakhtina, L.B. Neural Network Monitoring of Gas Turbine
   Engine. Neurocomputers: design and application, Radiotehnika, Vol. 1, 2001, pp. 37-43. (In
3. Neurocomputers in Aviation / Ed. by Vasilyev V.I., Ilyasov B.G., Kusimov. In: Book 14: study
   book for high school establishments. Moscow: Radiotehnika, 2004. 496 p. (In Russian)
4. Labendik, V., Ozolinsh, I., Zvanchuk, P. Method for turbofan thrust control in flight using similar
   engine mode of operation work. In: Abstract collection of II International Scientific and Technical
   Conference «Aero Engines of XXI Century», Vol. II, 06-09.12.2005, Moscow, Russia. Moscow:
   CIAM, 2005, pp. 244-246.
5. Zhernakov, S.V. Definition of aviation engine parameters using active expert system,
   Aviacionnaja promishlennost, Vol. 4, 2001, pp. 24-28. (In Russian)
6. Artificial tools of aviation engines reliability diagnostics and forecasting / V.I. Dubrovin,
   S.A. Subbotin, A.V. Boguslaev and others. Zaporozhje: OAO “Motor-Sich, 2003. 279 p. (In
7. Kopytov, E., Labendik, V., Osis, A., Tarasov. A. Method of aviation engine diagnostics in the case
   of partial loss of information. In: Abstract collection of II International Scientific and Technical
   Conference «Aero Engines of XXI Century», Vol. II, 06-09.12.2005, Moscow, Russia. Moscow:
   CIAM, 2005, pp. 246-247.
8. Kopytov, E., Labendik, V., Osis, A., Tarasov, A. Neural Networks Application for Analysis of
   Flight Information in Aircraft Engine Diagnostic System. In: Transport and Telecommunication,
   2006, vol. 7(2). Riga: TTI, 2006, pp. 287-294.

            The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”


    Konstantin N. Nechval1, Nicholas A. Nechval2, Gundars Berzins2, Maris Purgailis2,
  Juris Krasts2, Uldis Rozevskis2, Kristine Rozite2, Vladimir Strelchonok3, Natalie Zolova2
         Mathematical Methods and Modelling Department, Transport and Telecommunication Institute
                                  Lomonosova 1, Riga, LV-1019, Latvia
                                         E-mail: konstan@tsi.lv
                              Mathematical Statistics Department, University of Latvia
                                     Raina Blvd 19, Riga, LV-1050, Latvia
                                            E-mail: nechval@junik.lv
                                     Informatics Department, Baltic International Academy
                                            Lomonosova 1, Riga, LV-1019, Latvia
                                                    E-mail: str@apollo.lv

        This paper introduces a new technique for early identification of fatigue cracks, namely the constant false
alarm rate (CFAR) test. This test works on the null hypotheses that a target vibration signal is statistically similar
to a reference vibration signal. In effect, this is a time-domain signal processing technique that compares two
signals, and returns the likelihood whether the two signals are similar or not. The system monitors the vibration
signal of the rotor as it cycles, and compares that vibration signal with, say, the original vibration signal. The
difference vector reflects the change in vibration over time. By developing of a crack, the vector changes in a
characteristic way. Thus, it is possible, during CFAR test, to determine whether the two signals are similar or
not. Therefore, by comparing a given vibration signal to a number of reference vibration signals (for several
crack scenarios) it is possible to state which is the most likely condition of the rotor under analysis. The CFAR
test not only successfully identifies the presence of the fatigue cracks but also gives an indication related to the
advancement of the crack. This test, despite its simplicity, is an extremely powerful method that effectively
classifies different vibration signals, allowing for its safe use as another condition monitoring technique.
Keywords: vibration signal, changes, recognition, CFAR test


       The machines and structural components require continuous monitoring for the detection of
cracks and crack growth for ensuring an uninterrupted service. Non-destructive testing methods like
ultrasonic testing, X-ray, etc., are generally useful for the purpose. These methods are costly and time
consuming for long components, e.g., railway tracks, long pipelines, etc. Vibration-based methods can
offer advantages in such cases [1]. This is because measurement of vibration parameters like natural
frequencies is easy. Further, this type of data can be easily collected from a single point of the
component. This factor lends some advantages for components, which are not fully accessible. This
also helps to do away with the collection of experimental data from a number of data points on a
component, which is involved in a prediction based on, for example, mode shapes.
       Nondestructive evaluation (NDE) of structures using vibration for early detection of cracks has
gained popularity over the years and, in the last decade in particular, substantial progress has been
made in that direction. Almost all crack diagnosis algorithms based on dynamic behaviour call for a
reference signature. The latter is measured on an identical un-cracked structure or on the same
structure at an earlier stage.
       Dynamics of cracked rotors has been a subject of great interest for the last three decades and
detection and monitoring have gained increasing importance, recently. Failures of any high speed
rotating components (jet engine rotors, centrifuges, high speed fans, etc.) can be very dangerous to
surrounding equipment and personnel (see Figure 1), and must always be avoided. Jet engine disks
operate under high centrifugal and thermal stresses. These stresses cause microscopic damage as a
result of each flight cycle as the engine starts from the cold state, accelerates to maximum speed for

Session 5. Aviation

take-off, remains at speed for cruise, then spools down after landing and taxi. The cumulative effect of
this damage over time creates a crack at a location where high stress and a minor defect combine to
create a failure initiation point. As each flight operation occurs, the crack is enlarged by an
incremental distance. If allowed to continue to a critical dimension, the crack would eventually cause
the burst of the disk and lead to catastrophic failure (burst) of the engine. Engine burst in flight is
rarely survivable.

                                   Figure 1. Jet engine fan section failure

       In Figure 2 is given micrograph showing one of the cracks detected in the bladed disk assembly
of the High Pressure Turbine for over 17,000 cycles.

            Figure 2. Micrograph showing one of the cracks detected in the High Pressure Turbine
                                         (bladed disk assembly)

      Schematic showing the orientation of the columnar grains at the leading and trailing edges in the
blade root region of the failed turbine rotor blisk is presented in Figure 3.

            The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

       Figure 3. Schematic showing the orientation of the columnar grains at the leading and trailing edges
                             in the blade root region of the failed turbine rotor blisk


       Suppose that we desire to compare a target vibration signal and a kth reference vibration signal,
which have p response variables. Let xij(k) and yij be the ith observation of the jth response variable of
the kth reference signal and the target signal, respectively. It is assumed that all observation vectors,
xi(k)=(xi1(k), ..., xip(k))′, yi=(yi1, ..., yip)′, i=1(1)n, are independent of each other, where n is a number of
paired observations. Let zi(k) = xi(k)−yi, i=1(1)n, be paired comparisons leading to a series of vector
differences. Thus, in order to compare the above signals, and return the likelihood whether the two
signals are similar or not, it can be obtained and used a sample of n independent observation vectors
Z(k)=(z1(k), ..., zn(k)). Each sample Z(k), k∈{1, …, m}, is declared to be realization of a specific
stochastic process with unknown parameters. It is assumed here that zi(k), i=1(1)n, are independent
p-multivariate normal random variables (n≥p+2) with common mean a(k) and covariance matrix
(positive definite) Q(k). A goodness-of-fit testing for the multivariate normality is based on the
following theorem.
       Theorem 1 (Characterization of the multivariate normality). Let zi(k), i=1(1)n, be n
independent p-multivariate random variables (n≥p+2) with common mean a(k) and covariance matrix
(positive definite) Q(k). Let wr(k), r = p+2, …, n, be defined by

            r − ( p + 1) r − 1                                                                               ⎛ S (k )    ⎞
wr (k ) =                      (z r (k ) − z r −1 (k ) )′ S r−11 (k )(z r (k ) − z r −1 (k )) = r − ( p + 1) ⎜ r
                                                              −                                              ⎜ S (k )
                                                                                                                      − 1⎟,
                  p        r                                                                          p      ⎝ r −1      ⎠

r= p+2, …, n,                                                                                                                 (1)


                r −1
z r −1 (k ) =   ∑ z (k ) /(r − 1),
                i =1
                       i                                                                                                      (2)

                r −1
S r −1 (k ) =   ∑ (z (k ) − z
                i =1
                           i    r −1 ( k ))(z i ( k )   − z r −1 (k ))′,                                                      (3)

Session 5. Aviation

then the zi(k) (i=1, …, n) are Np(a(k),Q(k)) if and only if wp+2(k), …, wn(k) are independently
distributed according to the central F-distribution with p and 1, 2, . . . , n−(p+1) degrees of freedom,
        Proof. The proof is similar to that of [2] and so it is omitted here. ˛
        Goodness-of-fit testing for the multivariate normality. The results of Theorem 1 can be used
to obtain test for the hypothesis of the form H0: zi(k) follows Np(a(k),Q(k)) versus Ha: zi(k) does not
follow Np(a(k),Q(k)), ∀i = 1(1)n. The general strategy is to apply the probability integral transforms of
wk, ∀k = p+2(1)n, to obtain a set of i.i.d. U(0,1) random variables under H0 [2]. Under Ha this set of
random variables will, in general, not be i.i.d. U(0,1). Any statistic, which measures a distance from
uniformity in the transformed sample (for instance, Kolmogorov-Smirnov statistic) can be used as a
test statistic.
      Testing for similarity of the two signals. In this paper, for testing that the two signals (target
signal and reference signal) are similar, we propose a statistical approach that is based on the
generalized maximum likelihood ratio. We have the following hypotheses:
      H0(k): Similarity is valid for the acceptable range of accuracy under a given experimental
      H1(k): Similarity is invalid for the acceptable range of accuracy under a given experimental
      Thus, for fixed n, the problem is to construct a test, which consists of testing the null hypothesis

H0(k): zi(k) ∼ Np(0,Q(k)), ∀i = 1(1)n,                                                                       (4)

where Q(k) is a positive definite covariance matrix, versus the alternative

H1(k): zi(k) ∼ Np(a(k),Q(k)), ∀i = 1(1)n,                                                                    (5)

where a(k)=(a1(k), ..., ap(k))′ ≠ (0, ..., 0)′ is a mean vector. The parameters Q(k) and a(k) are unknown.


       In order to distinguish the two hypotheses (H0(k) and H1(k)), a generalized maximum likelihood
ratio (GMLR) statistic is used. The GMLR principle is best described by a likelihood ratio defined on
a sample space Z with a parameter set Θ, where the probability density function of the sample data is
maximized over all unknown parameters, separately for each of the two hypotheses. The maximizing
parameter values are, by definition, the maximum likelihood estimators of these parameters; hence the
maximized probability functions are obtained by replacing the unknown parameters by their maximum
likelihood estimators. Under H0(k), the ratio of these maxima is a Q(k)-free statistic. This is shown in
the following. Let the complete parameter space for θ(k)=(a(k),Q(k)) be Θ={(a(k), Q(k)): a(k)∈Rp,
Q(k)∈Qp}, where Qp is a set of positive definite covariance matrices, and let the restricted parameter
space for θ, specified by the H0(k) hypothesis, be Θ0={(a(k),Q(k)): a(k)=0, Q(k)∈Qp}. Then one
possible statistic for testing H0(k): θ(k)∈Θ0 versus H1(k): θ(k)∈Θ1, where Θ1=Θ−Θ0, is given by the
generalized maximum likelihood ratio

LR = max LH1 ( k ) (Z(k ); θ(k ) )              max LH 0 ( k ) (Z(k ); θ(k ) ).                              (6)
       θ ( k )∈Θ1                              θ ( k )∈Θ 0

Under H0(k), the joint likelihood for Z(k) is given by

                                                                 ⎛      n
                                                                       ∑ z′ (k )[Q(k )]
                                                      − n/2
LH 0 ( k ) (Z(k ); θ(k )) = (2π )   − np/2
                                             Q( k )           exp⎜ −
                                                                 ⎜            i
                                                                                           z i (k ) / 2 ⎟.
                                                                                                        ⎟    (7)
                                                                 ⎝     i =1                             ⎠

           The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

Under H1(k), the joint likelihood for Z(k) is given by

                                                             ⎛         n
                                                                   ∑ (z (k ) − a(k ))′[Q(k )]
                                                  − n/2
LH1 ( k ) (Z(k ); θ(k )) = (2π ) − np/2 Q(k )             exp⎜ −
                                                             ⎜                i
                                                                                                      ( z i (k ) − a(k ))/2 ⎟.
                                                                                                                            ⎟     (8)
                                                             ⎝     i =1                                                     ⎠

It can be shown that

                                                               − n/2
  max LH j ( k ) (Z(k ); θ(k ) ) = (2π ) − np/2 Q (k )                     exp(−np / 2), j = 0, 1,                                (9)
θ ( k )∈Θ j                                       j

        ˆ                      ˆ
where Q 0 ( k ) =Z(k)Z′(k)/n, Q1 (k ) =(Z(k) − a (k)u′)(Z(k) − a (k)u′)′/n, and a (k)=Z(k)u/u′u are the well-
                                               ˆ               ˆ                ˆ
known maximum likelihood estimators of the unknown parameters Q(k) and a(k) under the
hypotheses H0(k) and H1(k), respectively, u=(1, ..., 1)′ is the n-dimensional column vector of units. A
substitution of (9) into (6) yields
                n/2             −n / 2
LR = Q 0 (k )         ˆ
                      Q1 (k )            .                                                                                       (10)

Taking the (n/2)th root, this likelihood ratio is evidently equivalent to
       ˆ        ˆ
LR • = Q 0 (k ) Q1 (k )         = Z(k )Z′(k ) / Z( k )Z′(k ) − (Z(k )u)(Z(k )u)′ / u′u .                                         (11)

Now the likelihood ratio in (11) can be considerably simplified by factoring out the determinant of the
p × p matrix Z(k)Z′(k) in the denominator to obtain this ratio in the form

LR • = 1 1 − (Z(k )u)′[Z(k )Z′(k )]−1 (Z( k )u) / n .              )                                                             (12)

This equation follows from a well-known determinant identity. Clearly (12) is equivalent finally to the

          ⎛n− p⎞            ⎛n− p⎞             −1
          ⎜ p ⎟(LR • − 1) = ⎜ p ⎟na′(k )[T(k )] a(k ),
Vn (k ) = ⎜    ⎟            ⎜    ⎟                                                                                               (13)
          ⎝    ⎠            ⎝    ⎠

where T(k ) = nQ1 (k ) . It is known that (a(k ), T( k )) is a complete sufficient statistic for the parameter
θ(k)=(a(k),Q(k)). Thus, the problem has been reduced to consideration of the sufficient statistic
 (a(k ), T( k )) . It can be shown that under H0, Vn is a Q(k)-free statistic which has the property that its
distribution does not depend on the actual covariance matrix Q(k). This is given by the following

      Theorem 2 (PDF of the statistic Vn(k)). Under H1(k), the statistic Vn(k) is subject to a non-
central F-distribution with p and n−p degrees of freedom, the probability density function of which is
                                                  −1         p
                              ⎡ ⎛ p n − p ⎞⎤ ⎡ p ⎤ 2           −1
f H1 ( k ) (vn (k ); n, q ) = ⎢Β⎜ ,       ⎟⎥ ⎢     ⎥ vn (k ) 2
                              ⎣ ⎝2    2 ⎠⎦ ⎣ n − p ⎦

                                  −n                     ⎛ p nq (k ) ⎡   n− p ⎤ ⎞
          ⎡              ⎤
             p                     2
       × 1+
                 vn (k ) ⎥               e − nq / 2 1 F1 ⎜ n ; ;      1+              ⎟, 0<vn(k)<∞.                              (14)
                                                                 2 ⎢     pvn (k ) ⎥ ⎟
          ⎢              ⎥
          ⎢ n− p         ⎥                               ⎜2 2        ⎣            ⎦ ⎠
          ⎣              ⎦

Session 5. Aviation

where 1F1(b;c;x) is the confluent hyper-geometric function, q(k)=a′(k)[Q(k)]−1a(k) is a non-centrality
parameter. Under H0(k), when q(k) = 0, (14) reduces to a standard F-distribution with p and n−p
degrees of freedom,
                                    −1                    p                      −n
                           ⎡ ⎛ p n − p ⎞⎤ ⎡ p ⎤ 2           −1 ⎡  p          ⎤    2
f H 0 ( k ) (vn (k ); n) = ⎢Β⎜ ,       ⎟⎥ ⎢     ⎥ vn (k ) 2 ⎢1 +      vn (k )⎥        ,   0<vn(k)<∞.   (15)
                           ⎣ ⎝2    2 ⎠⎦ ⎣ n − p ⎦              ⎣ n− p        ⎦

       Proof. The proof follows by applying Theorem 1 [3] and being straightforward is omitted here. ˛


       The CFAR test of H0(k) versus H1(k), based on Vn(k), is given by

       ⎧ ≥ h(k ), then H1 (k )
Vn (k )⎨                                                                                               (16)
       ⎩< h(k ), then H 0 (k ),

where h(k)>0 is a threshold of the test which is uniquely determined for a prescribed level of
significance α(k). It follows from (15) that this test achieves a fixed probability of a false alarm.
       If Vn(k)>h(k) then the kth reference vibration signal is eliminated from further consideration.
       If (m −1) reference vibration signals are so eliminated, then the remaining reference vibration
signal (say, kth) is the one with which the target vibration signal may be identified.
       If all reference vibration signals are eliminated from further consideration, we decide that the
target vibration signal cannot be identified with one of the m specified reference vibration signals.
       If the set of reference vibration signals not yet eliminated has more than one element, then we
declare that the target vibration signal may be identified with the k*th reference vibration signal,

k * = arg max (h(k ) − Vn (k )),                                                                       (17)
          k ∈D

where D is the set of simulation models not yet eliminated by the above test.


      The main idea of this paper is to find a test statistic whose distribution, under the null
hypothesis, does not depend on unknown (nuisance) parameters. This allows one to eliminate the
unknown parameters from the problem.


1.   Dimarogonas, A.D. Vibration of Cracked Structures: a State of the Art Review, Engineering and
     Fracture Mechanics, Vol. 55, 1996, pp. 831-857.
2.   Nechval, N.A., Nechval, K.N. Characterization Theorems for Selecting the Type of Underlying
     Distribution. In: Abstract Book of Communications of the 7th Vilnius Conference on Probability
     Theory and Mathematical Statistics & the 22nd European Meeting of Statisticians. Vilnius,
     Lithuania: TEV, 1998, pp. 352-353.
3.   Nechval, N.A. Radar CFAR Thresholding in Clutter under Detection of Airborne Birds. In:
     Proceedings of the 21st Meeting of Bird Strike Committee Europe. Jerusalem, Israel: BSCE,
     1992, pp. 127-140.

          The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”


                                           Konstantin N. Nechval

         Mathematical Methods and Modelling Department, Transport and Telecommunication Institute
                                  Lomonosova 1, Riga, LV-1019, Latvia
                                         e-mail: konstan@tsi.lv

        For important fatigue-sensitive structures of aircraft whose breakdowns cause serious accidents, it is
required to keep their reliability extremely high. In this paper, we discuss inspection strategies for such important
structures against fatigue failure. The focus is on the case when there are fatigue-cracks unexpectedly detected in
a fleet of aircraft within a warranty period (prior to the first inspection). The paper examines this case and
proposes stochastic models for prediction of fatigue-crack growth to determine appropriate intervals of the
inspections. We also do not assume known parameters of the underlying distributions, and the estimation of that
is incorporated into the analysis and decision-making. Numerical example is provided to illustrate the procedure.
Keywords: aircraft, fatigue crack, inspection interval


       Fatigue is one of the most important problems of aircraft arising from their nature as multiple-
component structures, subjected to random dynamic loads. The analysis of fatigue crack growth is one
of the most important tasks in the design and life prediction of aircraft fatigue-sensitive structures (for
instance, wing, fuselage) and their components (for instance, aileron or balancing flap as part of the
wing panel, stringer, etc.). An example of in-service cracking from B727 aircraft (year of manufacture
1981; flight hours not available; flight cycles 39,523) [1] is given on Figure 1.

                           Figure 1. Example of in-service cracking from B727 aircraft

       Several probabilistic or stochastic models have been employed to fit the data from various
fatigue crack growth experiments. Among them, the Markov’s chain model [2], the second-order
approximation model [3], and the modified second-order polynomial model [4]. Each of the models
may be the most appropriate one to depict a particular set of fatigue growth data but not necessarily
the others. All models can be improved to depict very accurately the growth data but, of course, it has

Session 5. Aviation

to be at the cost of increasing computational complexity. Yang’s model [3] and the polynomial model
[4] are considered more appropriate than the Markov’s chain model [2] by some researchers through
the introduction of a differential equation which indicates that fatigue crack growth rate is a function
of crack size and other parameters. The parameters, however, can only be determined through the
observation and measurement of many crack growth samples. If fatigue crack growth samples are
observed and measured, descriptive statistics can then be applied directly to the data to find the
distributions of the desired random quantities. Thus, these models still lack prediction algorithms.
Moreover, they are mathematically too complicated for fatigue researchers as well as design
       A large gap still needs to be bridged between the fatigue experimentalists and researchers who
use probabilistic methods to study the fatigue crack growth problems.


       Let us assume that a fatigue-sensitive component has been found cracked on n aircraft within a
warranty period. The cracking has not yet caused an accident, but the safety experts have told the
manager that if this item fails, an accident will be possible. It is clear that the part will have to be
redesigned and replaced. The manager’s dilemma is that redesigning the part, manufacturing the new
design, and installing it in the fleet will take, say, at least two years. The manager must decide how to
manage risk for the next two years. The alternatives include doing nothing and accepting the risk of
continued cracking and the possibility of an accident. An inspection program is usually instigated,
which should reduce the risk of failure, but due to uncertainties in aircraft loading histories, provides
no direct measurement of the criticality of the detected cracks. Generally, such a program will lead to
some aircraft being grounded, eliminating risk for those aircraft and reducing overall risk, but reducing
operational capability. This would leave precious few aircraft to spare before the service’s ability to
accomplish its mission became impaired. In such a scenario, the decision process involves a complex
probability problem concerning the likelihood of additional failures and acceptable risk. To compound
the difficulty little guidance is provided in aircraft design specifications for this situation. The situation
presented is not uncommon.
       The purpose of this paper is to present a more accurate stochastic crack growth analysis method,
while maintaining the simplicity of the proposed stochastic fatigue models, for the above problem. We
discuss the optimal relationship between the inspection time and the pre-specified minimum level of
reliability. To illustrate the proposed technique, a numerical example is given.


      The basis of most of the fatigue models is Paris-Erdogan’s law [5] relating the rate of growth of
crack size a to N cycles:

da ( N )
         = q[a ( N )]b                                                                                    (1)

in which q and b are parameters depending on loading spectra, structural/material properties, etc. We
fit da/dN vs a(N) with a function that we can integrate between limits (initial crack size, a0, and any
given crack size, a) to get a life prediction.
       In the linear region (see Fig. 2) we use Paris-Erdogan’s Equation (1) as follows. Integrating

 N          a( N )
     dN =    ∫
                     qa b
                          ,                                                                               (2)

we have

            The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

N − N0 =
             q( −b + 1)
                         [               −
                        a( N ) − b +1 − a0 b +1 .                                                               (3)

Thus, the crack growth equation representing the solution of the differential equation for Paris-
Erdogan’s law is given by

                1 ⎛ 1              1      ⎞
N − N0 =              ⎜ b −1 −            ⎟.                                                                    (4)
             q(b − 1) ⎝ 0
                      ⎜a             b −1 ⎟
                               a( N ) ⎠

                                 Figure 2. Crack growth rate versus crack size curve
                        (I = near-threshold region; II = linear region; III = instability region)

3.1. Sensitivity Analysis

          Consider the solution of the differential equation for Paris-Erdogan’s law written in the form of
(4) as:

              1 ⎛ 1          1 ⎞
N (a) =             ⎜ b −1 − b −1 ⎟,                                                                            (5)
           q(b − 1) ⎝ 0     a ⎟   ⎠

where a0 is the initial crack size at N0=0. The derivatives of the number of load cycles with respect to
the parameters q and b read:

dN ( a )    N (a)
         =−                                                                                                     (6)
 dq          q


dN ( a )     1 ⎡ 1 ⎛ ln a ln a0 ⎞         ⎤
         =      ⎢ ⎜ b −1 − b −1 ⎟ − N ( a)⎥.
                    ⎜a          ⎟                                                                               (7)
 db        b −1 ⎢ q ⎝
                ⎣         a0 ⎠            ⎥

Session 5. Aviation

From this one can see that the number of cycles to reach a certain crack size is very sensitive to
changes of the parameter q.


       The traditional analytical method of engineering fracture mechanics (EFM) usually assumes that
crack size, stress level, material property and crack growth rate, etc. are all deterministic values which
will lead to conservative or very conservative outcomes. However, according to many experimental
results and field data, even in well-controlled laboratory conditions, crack growth results usually show
a considerable statistical variability (as it is shown on Fig. 3).

                                    Figure 3. Constant amplitude loading fatigue test data curves

       Yet more considerable statistical variability is the case under variable amplitude loading (as
shown on Fig. 4).
       The basis of most data analyses seems to be to take logarithms in (1) and estimate b and q by
least squares in the equation

ln(da ( N ) / dN ) = ln q + b ln a ( N ).                                                             (8)

Unfortunately to use this equation estimates of da(N)/dN are required. Estimates of derivatives are
notoriously unreliable. If several repetitions of an experiment under the same conditions are made it is
not always clear how to combine the results. Moreover, as regressions model the properties of the
estimates of the coefficients in (8) are not the same as those of estimates of the coefficients in (4).
Thus it is sensible to ask why the estimation does not proceed directly from the data on crack size and
cycles through equation (4).
      It is interesting to note that if b were known q could be estimated from a straight line

  1             1
        −              = (b − 1)q ( N − N 0 )                                                         (9)
a0 −1
            a   b −1

and indeed such a plot for a few values of q is indicative of the nature of Paris-Erdogan’s equation in a
particular case.

         The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

                        Figure 4. Variable amplitude loading fatigue test data curves

       During the service of the components being assessed, there may be uncertainties in the applied
loading conditions, extrapolation of the material data to service conditions, component dimensions,
and nature, size and location of detected (postulated) defects, etc. These uncertainties/variations are
critical inputs to the crack growth assessment and can be taken into account using probabilistic
methodologies. There is now an extensive literature on the subject of the statistical nature of crack
growth. Most of the literature is concerned with model building and the agreement between the
general features of the model and the observed behaviour of the crack. However, little use has been
made of the statistical nature of the models to analyse experimental results.
       While most industrial failures involve fatigue, the assessment of the fatigue reliability of
industrial components subjected to various dynamic loading situations is one of the most difficult
engineering problems that remain. Material degradation processes due to fatigue depend upon material
characteristics, component geometry, loading history and environmental conditions. As a result,
stochastic models for crack growth have been suggested by many investigators in the last 15 years.
These include evolutionary probabilistic models, cumulative jump models and differential equation
(DE) models. DE models are the most widely used models for predicting stochastic crack growth
accumulation in the reliability and durability analyses of fatigue critical components.
       In practical applications of the stochastic crack growth analysis, either one of the following two
distribution functions is needed: the distribution of the crack size at any service time or the distribution
of the service time to reach any given crack size. Unfortunately, when the crack growth rate is
modelled as a random process, these two distribution functions are not amenable to analytical
solutions. As a result, numerical simulation procedures have been used to obtain accurate results. The
simulation approach is a very powerful tool; in particular, with modern high-speed computers.
However, it is a very time consuming procedure and therefore simple approximate analytical solutions
are very useful in engineering.
       The purpose of this paper is to present a more useful stochastic fatigue crack growth models by
using the solution of Paris-Erdogan’s law equation, which result in a simple analytical solution for
either the distribution of the service time to reach any given crack size or the distribution of the crack
size at any service time.
       The probability that crack size a(N) will exceed any given crack size a• in the service interval
(N0,N·), Pr{a(N·)>a•}, is frequently referred to as crack exceedance probability and can be found based
on the stochastic fatigue crack growth model. In addition to this probability distribution of crack size,
the probability distribution of cycles (or time) for a crack to grow from size a0 to a•, Pr{N(a•)≤ N •},
can also be found based on the above model. In fact, the probability that service time N(a•) will be

Session 5. Aviation

within the interval (N0, N •) for crack size to reach a• is identical to Pr{a(N·)>a•}.That is Pr{N(a•)≤ N •}
= Pr{a(N·) > a•}. To summarize the concept of the above derivation, the readers can refer to Fig. 5.

             Figure 5. Schematic diagram of crack size distribution and random time distribution


       These models allow one to describe the uncertainties in one or two parameters of the solution
(4) of the differential equation (1) for Paris-Erdogan’s law via parameters modelled as random
variables in order to characterize the random properties, which seem to vary from specimen to
specimen (see Fig. 3). In other words, the stochastic fatigue-crack growth parameter variability models
(with respect to the parameters b and q modelled as random variables) are given by

             1      ⎛ 1      1 ⎞
N − N0 =            ⎜ B −1 − B −1 ⎟ ,                                                                  (10)
           Q(B − 1) ⎜ a0
                    ⎝       a ⎟   ⎠

where (N−N0) is a joint random variable of B and Q. In fact, these models are suited to account for this
type of variability. The ones however cannot explain the variability of the crack growth rate during the
crack growth process. In particular, crack growth data (crack size versus service time and vice versa)
may be analysed using Eq. (10) by considering, for instance, two different approaches:
(i) B is identical for each specimen and Q varies from specimen to specimen, referred to as Case 1;
(ii) both B and Q vary from specimen to specimen, referred as Case 2.
        For Case 1, with B=1, the crack growth data for each specimen are best fitted by equation

           ⎡ a( N ) ⎤
        ln ⎢          ⎥
             a( N 0 ) ⎦
N − N0 = ⎣                                                                                             (11)

to obtain a sample value of Q, where a(N) ≡ a, a(N0) ≡ a0. For Case 2 equation (10) is used to best fit
the crack growth data for each specimen to obtain a set of sample values of B and Q. From the
statistical standpoint, B is considered to be a deterministic value and Q to be a statistical (random)
variable in Case 1, while both B and Q are considered to be statistical variables in Case 2. It is found
that the lognormal or Weibull’s distribution provides a reasonable fit for B and Q in both cases.

           The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

5.1. Weibull’s Crack Growth Parameter Variability Model

       Consider Case 1. The Weibull’s probability distribution function, F(q;σ,δ) , of Q is expressed as

               ⎧1 − exp[−(q / σ )δ ],     q ≥ 0,
F (q;σ , δ ) = ⎨                                                                                                   (12)
               ⎩0,                    otherwise,

in which F(q;σ,δ) is the probability that Q is smaller than or equal to an arbitrary value q; σ and δ are
distribution parameters representing the scale parameter and the shape parameter, respectively.

5.2. Crack Exceedance Probability

      For B=1, the probability that crack size a(N) will exceed any given (say, maximum allowable)
crack size a• can be derived and expressed as

                      ⎧    ln[a • / a ( N 0 )] ⎫       ⎡ ⎛ ln[a • / a ( N )] ⎞δ ⎤
                                               ⎬ = exp ⎢− ⎜
               •                                                         0 ⎟ ⎥
Pr{a ( N ) > a } = Pr ⎨Q >                                                        .                                (13)
                      ⎩       N − N0 ⎭                 ⎢ ⎜ σ ( N − N0 ) ⎟ ⎥
                                                          ⎝                  ⎠ ⎦

For B=b≠1 the maximum allowable crack size exceedance probability for a single crack is given by

                     ⎧    [a( N 0 )]− (b −1) − [a • ]− (b −1) ⎫       ⎡ ⎛ [a( N )]− (b −1) − [a • ]− (b −1) ⎞δ ⎤
                                                              ⎬ = exp ⎢− ⎜                                  ⎟ ⎥
Pr{a( N ) > a } = Pr ⎨Q >                                                      0
                                                                         ⎜ σ (b − 1)( N − N ) ⎟ ⎥ .                (14)
                     ⎩         (b − 1)( N − N 0 )             ⎭       ⎢ ⎝                                   ⎠ ⎦
                                                                      ⎣                            0

It will be noted that the crack exceedance probability can be used for assigning sequential in-service
inspections [6].


      These models allow one to characterize the random properties, which seem to vary during crack
growth (see Fig. 4), via crack growth equation with a stochastic noise V dependent, in general, on the
crack size a:

   1 ⎛ 1          1 ⎞
         ⎜ b −1 − b −1 ⎟ = N − N 0 + V                                                                             (15)
(b − 1)q ⎝ 0 a ⎟

  ⎛ 1 ⎛ 1             1 ⎞⎞
ln⎜          ⎜ b −1 − b −1 ⎟ ⎟ = ln( N − N 0 ) + V                                                                 (16)
  ⎜ (b − 1)q ⎜a      a ⎟⎟
  ⎝          ⎝ 0           ⎠⎠

 1      1
     −       = (b − 1)q( N − N 0 ) + V ,                                                                           (17)
a0 −1 a b −1

  ⎛ 1      1 ⎞
ln⎜ b −1 − b −1 ⎟ = ln[(b − 1)q( N − N 0 )] + V ,
  ⎜a                                                                                                               (18)
  ⎝ 0     a ⎟   ⎠

and so on, where V~N(0,σ2(b,q,N)), a0 ≡ a(N0), a ≡ a(N). They are suited to account for this type of
variability. The ones however cannot explain the variability of the crack growth rate from specimen to

Session 5. Aviation

6.1. Crack Exceedance Probability

     If V~N(0,σ2) in (15), then the probability that crack size a(N) will exceed any given (say,
maximum allowable) crack size a• can be derived and expressed as

                      ⎛⎡           [a ]− (b −1) − [a • ]− (b −1) ⎤  ⎞
Pr{a ( N ) > a • } = Φ⎜ ⎢ N − N 0 − 0                            ⎥ σ⎟,                                           (19)
                      ⎜                    (b − 1)q                 ⎟
                      ⎝⎣                                         ⎦  ⎠

where Φ(.) is the standard normal distribution function,

Φ(η ) =           ∫ exp(− x
                                  / 2) dx.                                                                       (20)
            2π   −∞

      If V~N(0,[(b-1)σ(N-N0)1/2]2) in (17), then the probability that crack size a(N) will exceed any
given (say, maximum allowable) crack size a• can be derived and expressed as

                      ⎛ ⎡ (b − 1)q( N − N 0 ) − ([a0 ]− (b −1) − [a • ]− (b −1) ) ⎤ ⎞
Pr{a ( N ) > a • } = Φ⎜ ⎢                                                         ⎥⎟ .                           (21)
                      ⎜              (b − 1)σ ( N − N 0 )1 / 2                      ⎟
                      ⎝ ⎣                                                         ⎦⎠

In this case, the conditional probability density function of a is given by

                                a −b               ⎛ 1 ⎡ (a − (b −1) − a − (b −1) ) − (b − 1)q ( N − N ) ⎤ 2 ⎞
f ( a; a0 , N 0 , N ) =                         exp⎜ − ⎢ 0                                            0
                                                                                                         ⎥ ⎟.
                                                                                                             ⎟   (22)
                        σ [2π ( N − N 0 )]1 / 2    ⎜ 2               (b − 1)σ ( N − N 0 )1 / 2
                                                   ⎝   ⎣                                                 ⎦ ⎠

6.2. Data Analysis for a Single Crack

       Consider the regression model corresponding to (17). Because the variance is non-constant (17)
is a non-standard model; however, on dividing by (N−N0)1/2 the model becomes

 a1− b − a1− b
           1/ 2
                = (b − 1)q( N − N 0 )1 / 2 + W ,                                                                 (23)
( N − N0 )

where W is normally distributed with mean zero and standard deviation (b−1)σ independent of N.
Thus if b is known the estimator of (b−1)q is just the least-squares estimator of the coefficient in
Equation (23) and the estimate of (b−1)σ is just the estimate of the variance of the regression. It
remains to determine what to do about b.
       Given the data describing a single crack, say, a sequence {(ai , N i )}in=1 , it is easy to construct a
log-likelihood using the density given by (22) and estimate the parameters b, q and σ by maximum
likelihood. The log-likelihood is

                                               1           n
                                                               ⎛ a1− b − ai1− b − (b − 1)q( N i − N 0 ) ⎞
                             i =1
L(b, q,σ ;{( ai , N i )} = −b ln ai − n ln σ −
                                               2         ∑     ⎜ 0
                                                          i =1 ⎝       (b − 1)σ ( N i − N 0 )1 / 2
                                                                                                        ⎟ .

Inspection shows that this differs from the standard least-squares equation only in the term –b∑lna,
where the subscript i has been dropped. The likelihood estimators are obtained by solving the

             The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

dL/db =0; dL/dq =0; dL/dσ =0.                                                                                                                                      (25)

In this case the equations have no closed solution. However, it is easy to see that the estimators for q
and σ given b are the usual least-squares estimators for the coefficients in (23) conditioned on b,

          1 ⎛ 1− b                                       ⎞⎛                  ⎞
                                             n                n
q (b) =      ⎜ na0 −
        b −1 ⎜
                                       ∑  i =1
                                                 ai1− b ⎟⎜
                                                        ⎟⎜    ∑
                                                         ⎠⎝ i =1
                                                                 ( Ni − N0 ) ⎟ ,

                1                 n
                                       [a1− b − ai1− b − q (b)(b − 1)( N i − N 0 )]2
σ 2 (b) =
             n(b − 1) 2           ∑
                                  i =1
                                                         Ni − N0
                                                                                     ,                                                                             (27)

and on substituting these back in the log-likelihood gives a function of b alone,
L(b) = −b        ∑ ln a − n ln[σ (b)] − n / 2.
                 i =1
                              i                                                                                                                                    (28)

Thus the technique is to search for the value of b that maximizes L(b) by estimating q and σ as
functions of b and substituting in L(b). In this study a simple golden-section search worked very

6.3. Pooling Data

       When several experiments have been performed it is possible to combine the log-likelihood
from each experiment to give estimators of the parameters of interest. Suppose that several
experiments have been performed. Each experiment is labelled with j, j runs from 1 to m, and yields nj
observations. The data are then a set of sequences {(ajk,Njk)}, with j=1, …, m, k=1, …, nj. The log-
likelihood for the whole set of experiments is simply the sum of the log-likelihood for the individual
cracks; writing Lj(bj,qj,σj) for the log-likelihood for the jth crack gives

                                                                      1 ⎛ a0                                                          − (b j − 1)q j ( Nk − N0 ) ⎞
                                                                                                         nj         1− b j      1− b j
                                                                                                                             − ak
L j (b j , q j , σ j ;{(a jk , N jk )}) = −b j ln a jk − n j ln σ j −        ⎜
                                                                      2 k =1 ⎜                           ∑                   (b j − 1)σ j ( N − N0 )1 / 2
                                                                                                                                                                 ⎟ (29)
                                              k =1                           ⎝                                                                                   ⎠


L=    ∑ L (b , q ,σ
      j =1
             j        j   j        j ;{( a jk , N jk )}).                                                                                                          (30)

The global log-likelihood can be used to investigate explicit parametric models for the parameters, or
simply as a way to pool data. Estimation by maximum likelihood proceeds exactly as above; the qj and
σj are obtained as ordinary least-squares estimators from equations like (26) and (27), one for each
crack, and substituted back into the log-likelihood to yield

                                      m           nj                   m                                  m
L(b1 , b2 , ..., bm ) = −          ∑ ∑j =1
                                                  k =1
                                                         ln a jk −     ∑
                                                                       j =1
                                                                              n j ln[σ j (b j )] −
                                                                                                     2   ∑n .
                                                                                                         j =1
                                                                                                                j                                                  (31)

       When the cracks are all assumed to be independent with distinct parameters the estimators from
the joint log-likelihood are precisely those obtained by estimating from each separately as outlined

Session 5. Aviation

      If a common value of b is used and the qj and σj are assumed to absorb most of the experimental
variability, the joint log- likelihood reduces to

            m       nj               m                              m
L(b) = −b   ∑∑
            j =1 k =1
                         ln a jk −   ∑
                                     j =1
                                            n j ln[σ j (b)] −
                                                                2   ∑n .
                                                                    j =1
                                                                           j                       (32)


       These models allow one to describe the uncertainties in the fatigue-crack growth of Paris-
Erdogan’s law via crack growth equation with a stochastic noise dependent, in general, on the crack
size, and parameters modelled as random variables in order to characterize the random properties,
which seem to vary both from specimen to specimen and during crack growth (see Fig. 4). In other
words, the stochastic fatigue-crack growth parameter and propagation lifetime variability model (with
respect to the parameters B and Q, modelled as random variables, and the stochastic noise V
dependent, in general, on the crack size a) may be given, for example, as

               1     ⎛ 1      1 ⎞
N − N0 =             ⎜ B −1 − B −1 ⎟ + V .                                                         (33)
            (B − 1)Q ⎝ 0
                     ⎜a      a ⎟   ⎠

7.1. Crack Exceedance Probability

      In this case, the probability that crack size a(N) will exceed any given (say, maximum
allowable) crack size a• can be derived and expressed as

                     ⎧ ⎡ ⎛                                       δ ⎫
                •    ⎪ ⎢ ⎜ [a( N 0 )]− (b −1) − [a • ]− (b −1) ⎞ ⎤ ⎪
                                                               ⎟ ⎥⎬ .
Pr{a ( N ) > a } = E ⎨exp − ⎜                                  ⎟                                   (34)
                     ⎪ ⎢ ⎝ σ (b − 1)( N − N 0 + V ) ⎠ ⎥ ⎪
                     ⎩ ⎣                                          ⎦⎭


      Let us assume that a fatigue-sensitive component (outboard longeron, Fig. 6)) has been found
cracked on n=10 aircraft within a warranty period. Here a fleet of ten aircraft have all been inspected.

Table 1. Inspection results

   Aircraft                    Flight hours                           Crack size (mm)
                                   (Ni)                                     (ai)
      1                              3000                                        1
      2                              2300                                       0.5
      3                              2200                                        1
      4                              2000                                        2
      5                              1500                                       0.8
      6                              1500                                       1.5
      7                              1300                                        1
      8                              1100                                        1
      9                              1000                                        1
      10                              800                                        1

           The 7th International Conference “RELIABILITY and STATISTICS in TRANSPORTATION and COMMUNICATION - 2007”

                    Figure 6. Location of cracking on the F/A-18 Y497 upper outboard longeron

       It is assumed that cracks start growing from the time the aircraft entered service. For typical
aircraft metallic materials, an initial discontinuity size (a0) found through quantitative fractography is
approximately between 0.02 and 0.05 mm. Choosing a typical value for initial discontinuity state (e.g.,
0.02 mm) is more conservative than choosing an extreme value (e.g., 0.05 mm). This implies that if
the lead cracks can be attributed to unusually large initiating discontinuities then the available life
       We test a goodness of fit of the data of Table 1 with the Weibull’s fatigue-crack growth
parameter variability model (see Case 1), where

              a1−b − ai1−b
Qi (b) =        0
                                , i = 1(1)n,                                                                  (35)
           (b − 1)( N i − N 0 )

with a common value of b, a0=0.02, and N0=0.

8.1. Goodness-of-Fit Testing
       We assess the statistical significance of departures from the Weibull’s model by performing
empirical distribution function goodness-of-fit test. We use the S statistic [7]. For complete datasets,
the S statistic is given by

            n −1
                      ⎛ ln(Qi +1 (b) / Qi (b)) ⎞
           ∑          ⎜
       i =[ n / 2 ]+1 ⎝           Mi
                                               ⎠ = 0.43,
S (b) = n−1                                                                                                   (36)
                   ⎛ ln(Qi +1 (b) / Qi (b)) ⎞
           ∑i =1 ⎝
                   ⎜           Mi

Session 5. Aviation

where [n/2] is a largest integer ≤ n/2, Qi is the ith order statistic, the values of Mi are given in Table 13
[7]. The rejection region for the α level of significance is {S>Sn;1-α}. The percentage points for Sn;1-α
were given in [7]. The value of b is one that minimizes S(b). For this example, b = 0.87 and

S=0.43 < Sn=10; 1-α=0.95=0.69.                                                                          (37)

Thus, there is not evidence to rule out the Weibull’s model. Using the relationship (4), the inspection
results can be extrapolated from the expected initial crack size, a0, to the time of the next inspection
when the maximum allowable crack size is equal to a•=10 mm as presented in Table 2.

Table 2. Predicted next inspection results

     Aircraft         Maximum allowable                Next inspection time
                      crack size a• (mm)                  (flight hours)
       1                      10                               5626
       2                      10                               5503
       3                      10                               4126
       4                      10                               3033
       5                      10                               3030
       6                      10                               2477
       7                      10                               2438
       8                      10                               2063
       9                      10                               1875
       10                     10                               1500


       The authors hope that this work will stimulate further investigation using the approaches on
specific applications to see whether obtained results with it are feasible for realistic applications.


      This research was supported in part by Grant No. 06.1936 and Grant No. 07.2036 from the
Latvian Council of Science and the National Institute of Mathematics and Informatics of Latvia.


1.    Jones, R., Molent, L. and Pitt, S. Studies in Multi-Site Damage of Fuselage Lap Joints, J. Theor.
      Appl. Fract. Mech.,Vol. 32, 1999, pp. 18-100.
2.    Bogdanoff, J.L. and Kozin, F. Probabilistic Models of Cumulative Damage. New York: Wiley,
3.    Yang, J.N. and Manning, S.D. A Simple Second Order Approximation for Stochastic Crack
      Growth Analysis, Engineering Fract. Mech., Vol. 53, 1996, pp. 677-686.
4.    Wu, W.F., Ni, C.C. and Liou, H.Y. Random Outcome and Stochastic Analysis of Some Fatigue
      Crack Growth Data, Chin. J. Mech., Vol. 17, 2001, pp. 61-68.
5.    Paris, R. and Erdogan, F. A Critical Analysis of Crack Propagation Laws, Journal of Basic
      Engineering, Vol. 85, 1963, pp. 528-534.
6.    Nechval, N.A., Nechval, K.N. and Vasermanis, E.K. Statistical Models for Prediction of the
      Fatigue Crack Growth in Aircraft Service. In: Fatigue Damage of Materials 2003 (ed. by
      A. Varvani-Farahani & C. A. Brebbia). Southampton, Boston: WIT Press, 2003, pp. 435-445.
7.    Kapur, K.C. and Lamberson, L.R. Reliability in Engineering Design. New York: Wiley, 1977.


To top