Notes Reliability Theory 2004
Notes Reliability Theory 2004
Notes Reliability Theory 2004
Preface
The present collection of notes have been developed during the last 10 years as lecture notes for the
course in ‘Structural Reliability Theory’ at the 4th year in the MSc study in Civil and Structural
Engineering at Aalborg University.
Content:
Page
Note 0 Introduction to risk analysis 1
Note 1+2 Structural reliability 27
Note 3 First order reliability methods 49
Note 4 First order reliability analysis with correlated and non-normal stochastic
variables 65
Note 5 SORM and simulation techniques 83
Note 6 Reliability evaluation of series systems 103
Note 7 Reliability evaluation of parallel systems 119
Note 8 Structural reliability: Level 1 methods 133
Note 9 Time-variant reliability 161
Note 10 Load combinations 179
Note 11 Example: Fatigue / Reliability-based inspection planning 193
Note 12 Reliability updating 217
Contents:
1 Introduction....................................................................................................................................... 1
2 Definition of Risk.............................................................................................................................. 3
3 Framework for Risk Analysis ........................................................................................................... 3
4 Implementation of Risk Analysis...................................................................................................... 4
5 QRA methods.................................................................................................................................... 7
5.1 Event Tree Analysis .................................................................................................................. 7
5.2 Fault Tree Analysis ................................................................................................................. 13
5.3 Risk Matrix ............................................................................................................................. 16
5.4 Decision trees.......................................................................................................................... 17
6 Risk acceptance criteria .................................................................................................................. 18
7 References....................................................................................................................................... 24
1 Introduction
Operational Failure
(Eng. Judgement/
Mechanical Failure data/experience) Structural Failure
(Data/Fault trees etc.) (SRA)
Economic Risk
Decision
Sustainable development related to conservation of the environment, the welfare and safety of the peo-
ple have been subject to increasing concern of the society during the last decades. At the same time
optimal allocations of available natural and financial resources are considered very important. There-
fore methods of risk and reliability analysis in civil engineering developed during the last decades are
1
Note 0 – Risk analysis
becoming more and more important as decision support tools in civil engineering applications. The
decision process is illustrated in figure 1.
Civil engineering facilities such as bridges, buildings, power plants, dams and offshore platforms are
all intended to contribute to the benefit and quality of life. Therefore when such facilities are planned it
is important that the benefit of the facility can be identified considering all phases of the life of the fa-
cility, i.e. including design, manufacturing, construction, operation and eventually decommissioning.
Benefit has different meanings for different people in the society, simply because different people have
different preferences. However, benefit for the society can be understood as
• economically efficient for a specific purpose
• fulfil given requirements with regard to safety of people directly or indirectly involved with and
exposed to the facility
• fulfil given requirements to the effects of the facility on the community and environment
Taking into account these requirements it is seen that the task of the engineer is to make decisions or to
provide the decision basis for others such that it may be ensured that engineering facilities are estab-
lished, operated, maintained and decommissioned in such a way that they will optimise or enhance the
possible benefits to society and individuals of society.
For many years it has been assumed in design of structural systems that all loads and strengths are de-
terministic. The strength of an element was determined in such a way that it exceeded the load with a
certain margin. The ratio between the strength and the load was denoted the safety factor. This number
was considered as a measure of the reliability of the structure. In codes of practice for structural sys-
tems values for loads, strengths and safety factors are prescribed.
As described above structural analysis and design have traditionally been based on deterministic meth-
ods. However, uncertainties in the loads, strengths and in the modeling of the systems require that
methods based on probabilistic techniques in a number of situations have to be used. A structure is usu-
ally required to have a satisfactory performance in the expected lifetime, i.e. it is required that it does
2
Note 0 – Risk analysis
not collapse or becomes unsafe and that it fulfills certain functional requirements. Generally structural
systems have a rather small probability that they do not function as intended, see table 1.
This note gives an introduction to the different aspect of risk analysis and risk acceptance criteria and is
partly based on the JCSS working paper, Faber & Stewart [2].
2 Definition of Risk
Risk is here defined as the expected consequences associated with a given activity. Considering an ac-
tivity with only one event with potential consequences risk R is thus the probability that this event will
occur P multiplied with the consequences given the event occurs C i.e.
R = P ⋅C (1)
n
R = ∑ Pi ⋅ C i (2)
i =1
This definition is consistent with the interpretation of risk used e.g. in the insurance industry (expected
losses) and risk may e.g. be given in terms of DKKs, dollars, number of human fatalities, etc.
Inprovement of knowledge
Risk assessment is used in a number of situations with the general intention to indicate that important
aspects of uncertainties, probabilities and/or frequencies and consequences have been considered in
some way or other. Decision theory provides a theoretical framework for such analyses, see Figure 2.
In typical decision problems encountered the information basis is often not very precise. In many situa-
tions it is necessary to use historical and historical data. The available historical information is often
not directly related to the problem considered but to a somewhat similar situation. Furthermore, an im-
portant part of a risk assessment is to evaluate the effect of additional information, risk reducing meas-
ures and/or changes of the considered problem. It is therefore necessary that the framework for the de-
cision analysis can take these types of information into account and allow decisions to be updated
3
Note 0 – Risk analysis
based upon new information. This is possible if the framework of Bayesian decision theory is used, see
e.g. Raiffa and Schlaifer [3] and Benjamin and Cornell [4].
A fundamental principle in decision theory is that optimal decisions must be identified as those result-
ing in the highest expected utility, see e.g. Ditlevsen and Madsen [5]. In typical engineering applica-
tions the utility may be related to consequences in terms of costs, fatalities, environmental impact etc.
In these cases the optimal decisions are those resulting in the lowest expected costs, the lowest ex-
pected number of fatalities and so on.
• Who are the decision maker(s) and the parties with interests in the activity (e.g. society, cli-
ent(s), state and organizations).
• Which matters might have a negative influence on the impact of the risk analysis and its results.
• What might influence the manner in which the risk analysis is performed (e.g. political, legal,
social, financial and cultural).
Furthermore the important step of setting the acceptance criteria must be performed. This includes the
specification of the accepted risks in regard to economic, public or personnel safety and environmental
criteria. In setting the acceptable risks – which might be considered a decision problem itself, due ac-
count should be taken to both international and national regulations in the considered application area.
However, for risk analysis performed for decision making in the private or inter-company sphere with
no potential consequences for third parties the criteria may be established without the consideration of
such regulations. In these cases the issue of risk acceptance is reduced to a pure matter of cost or re-
source optimisation involving the preferences of the decision maker alone. Risk criteria are discussed
in section 6.
System Definition
The system (or the activity) considered has to be described and all assumptions regarding the system
representation and idealizations stated.
A hazard is typically referred to as a failure event for the considered system or activity. Occurrence of a
hazard is therefore also referred to as a system failure event. System failures may thus represent events
such as collapse of a building structure, flooding of a construction site or explosion in a road or rail
4
Note 0 – Risk analysis
tunnel. Identification of hazards is concerned about the identification of all events, which might have
an adverse consequence to
• People
• Environment
• Economy
Different techniques for hazard identification have developed from various engineering application
areas such as the chemical, nuclear power and aeronautical industries. Examples are:
• Preliminary Hazard Analysis (PHA)
• Failure Mode and Effect Analysis (FMEA)
• Failure Mode Effect and Criticality Analysis (FMECA)
• Hazard and Operability Studies (HAZOP)
• Risk Screening (Hazid sessions)
Analysis of Consequences
Typical consequences are economic consequences, loss of life and effects on the environment. The
estimation of consequences given failure of the system of sub-systems requires a good understanding of
the system and its interrelation with its surroundings and is thus best performed in collaboration with
experts who have hands on experience with the considered type of activity.
5
Note 0 – Risk analysis
Analysis of Probability
Evaluation of probabilities of failure for the individual components and sub-systems may be based on,
in principle, two different approaches: failure rates for e.g. electrical and production systems or meth-
ods for structural reliability for structural systems as buildings and bridges.
Risk analyses are typically made on the basis of information, which is subject to uncertainty. These
uncertainties may be divided in
Analyse Sensitivities
The sensitivity analysis is useful for analysis of the identified risk scenarios and normally includes an
identification of the most important factors influencing the risks associated with the different risk sce-
narios. Also the sensitivity analysis may include studies of “what if” situations for the evaluation of the
importance of various system simplifications performed under the definition of the system.
Risk Treatment
Calculated risks are compared with the accepted risks initially stated in the risk acceptance criteria.
Should the risks not be acceptable in accordance with the specified risk acceptance criteria there are
principally four different ways to proceed.
Risk mitigation: Risk mitigation is implemented by modification of the system such that the source of
risk is removed. For example, the risk of fatalities from a ship collision with a bridge may be mitigated
by traffic lights stopping traffic proceeding onto the bridge whenever a ship navigates under the bridge.
Risk reduction: Risk reduction may be implemented by reduction of the consequences and/or the prob-
ability of occurrence – in practice risk reduction is normally performed by a physical modification of
the considered system.
Risk transfer: Risk transfer may be performed by e.g. insurance or other financial arrangements where
a third party takes over the risk.
Risk acceptance: If the risks do not comply with the risk acceptance criteria and other approaches for
risk treatment are not effective then risk acceptance may be an option.
6
Note 0 – Risk analysis
back of information from the system. Whenever new information is provided the risk analysis may be
updated.
5 QRA methods
Quantitative Risk Analysis (QRA) is used in assessment of the risks. Three calculation methods are:
• Event Tree Analysis (ETA)
• Fault Tree Analysis (FTA)
• Risk matrix
These methods are described in the following partly based on [7].
An Event Tree is constructed by defining an initial event and the possible consequences that result from
this, when the emergency systems function or not. The initial event is usually placed on the left and
branches are drawn to the right, each branch representing a different sequence of events and terminating in
an outcome. The main elements of the tree are event definitions and branch points, or logic vertices.
The initial event is usually expressed as a frequency (events/year) and the subsequent splits as probabilities
(events/demand), so that the final outcomes are also expressed as frequencies (event/year). Each branch of
the Event Tree represents a particular scenario. An example of a simple Event Tree is shown in figure 4
and 5. The fire protection is provided by a sprinkler system. A detector will either detect the rise in tempera-
ture or it will not. If the detector succeeds the control box will either work correctly or it will not - and so
on. There is only one branch in the tree that indicates that all the subsystems have succeeded:
7
Note 0 – Risk analysis
SPRINKLER
INITIATING FIRE SPREADS PEOPLE CAN- RESULTANT
FAILS TO
EVENT QUICKLY NOT ESCAPE EVENT
WORK SCENARIO
P= 0.5
MULTIPLE
1
FATALITIES
P= 0.3
NO
P=0.1 LOSS/DAMAGE 2
P= 0.5
YES
FIRE FIRE CON-
NO 3
STARTS TROLLED
NO
FIRE CON- 4
Frequency = 1/yr
TAINED
P=0.9
Figure 5. Event Tree Analysis for Building Protected by a Sprinkler System.
The results of the Event Tree are outcome event frequencies (probabilities) per year. The outcome frequen-
cies may be processed further to obtain the following results:
Risk to Workforce:
• Annual risk
• Individual risk
• Fatal Accident Rate
Risk to Public:
Physical effects:
• Contours on site map
• Transects on site map
Individual risk:
• Contours on site map
• Transects on site map
• Annual risk at fixed location
Societal risk:
• FN table
• FN curve
• Equivalent annual fatalities
Risk to Workforce
Annual Risk
The annual risk may be expressed as Potential Loss of Life (PLL), where the PLL expresses the probability
of fatalities per year for all the operation personnel. As such the PLL is a risk indicator which is valid for
the whole installation, rather than for an individual. The calculation for a given event i is of the form
8
Note 0 – Risk analysis
Individual Risk
The Individual Risk (IR) expresses the probability per year of fatality for one individual. It is also termed as
Individual Risk Per Annum (IRPA). The IR depends on the location of the individual at a given time and its
contents of work. In practice, for the operating personnel of an installation an Average Individual Risk, AIR
may be estimated for groups of persons taking into account the percentage of time of exposure to the hazard
per year. For all the personnel involved in the annual operation of the installation, the AIR may be derived
from the PLL
PLL
AIR = (5)
Nb ⋅ Q
where
PLL potential loss of life for the installation, (1/year)
Nb number of personnel involved per year, (1/year)
Q is average percentage of exposure
9
Note 0 – Risk analysis
Scenario PLL/year %
Hydrocarbon 2.20E-03 64
Ship collision 7.60E-04 22
Helicopter crash 8.90E-05 3
Structural + earthquake 1.20E-04 3
Flying 2.60E-04 8
TOTAL Installation Risk (PLL) 3.43E-03 100
PLL/year
4.00E-03
3.50E-03
3.00E-03
2.50E-03
2.00E-03 PLL/year
1.50E-03
1.00E-03
5.00E-04
0.00E+00
Hydrocarbon Ship collision Helicopter Structural + Flying TOTAL
crash earthquake Installation
Risk (PLL)
% PLL
Hydrocarbon
Ship collision
Helicopter crash
Structural + earthquake
Flying
10
Note 0 – Risk analysis
The Fatal Accident Rate (FAR) is defined as the potential number of fatalities in a group of people exposed
for a specific exposure time to the activity in question. Generally, the FAR is expressed as a probability of
fatality per 100 million exposure hours for a given activity. It is mainly used for comparing the fatality risk
of activities. The 100 million exposure hours is to represent the number of hours at work in 1000 working
lifetimes.
PLLarea ⋅ 10 8
FARarea = (6)
N area ⋅ 8760
where
FARarea area specific Fatal Accident Rate, (1/year)
PLLarea Potential Loss of Life in an area per year, (1/year)
N area manning in an area calculated as an average value over a typical year of operation
8760 number of hours per year
The installation FAR value is calculated as an average value over a year of operation as
PLLtot ⋅ 10 8
FARinstallation = (7)
N tot ⋅ 8760
where
FARinstallation whole installation FAR, (1/year)
PLLtot total PLL for the installation, (1/year)
N tot number of personnel on the installation
11
Note 0 – Risk analysis
FAR/year
10
9
Helicopter - passengers
8 Helicopter - platform
7 Structural failure
6 Ship collision
5 Shale shaker
4 Turbines
Risers and pipelines
3
Process
2 Well events
1
0
Drilling Personnel FAR/year Production personnel FAR/year
% FAR
Well events
Process
Risers and pipelines
Turbines
Shale shaker
Ship collision
Structural failure
Helicopter - platform
Helicopter - passengers
12
Note 0 – Risk analysis
Risk to Public
FN Curves
Number (N) of Potential Frequency of Scenario Frequency of Incidents with Potential (N) or more
Scenario Fatalities per Year Fatalities per Year
1 1 0.1 0.12021
2 20 0.014 0.01141
3 70 0.0075 0.00713
4 150 0.00023 0.00022
5 300 0.00009 0.00011
6 500 0.00001 0.00001
F/N
1.00E+00
Frequency of Potential Event, F(N) with N or
1.00E-01
more Fatalities per Year
1.00E-02
F(N)
1.00E-03
1.00E-04
1.00E-05
1 10 100 1000
Potential Fatalities per Event (N)
FN curves or F/N plots (generally also called the “Cumulative Frequency Graphs”) are probability versus
consequence diagrams where “F” denotes frequency of a potential event and “N” the number of associated
fatalities.
A Cumulative Frequency Graph shows the probability of N or more fatalities occurring. Such graphs
tend to be of interest when the risk acceptance criterion selected, or, as is more often the case, imposed
by the Regulator, includes an aversion to potential incidents that would result in, say, more than ten
fatalities. In simple terms, risk aversion exists if society regards a single accident with 100 fatalities as
in some sense worse than 100 accidents (e.g. road accidents) with a single fatality each. An example of
a F/N graph is shown in figure 8.
13
Note 0 – Risk analysis
The Fault Tree is both a qualitative and a quantitative technique. Qualitatively it is used to identify the indi-
vidual scenarios (so called paths or cut sets) that lead to the top (fault) event, while quantitatively it is used
to estimate the probability (frequency) of that event.
A component of a Fault Tree has one of two binary states, either in the correct state or in a fault state. In
other words, the spectrum of states from total integrity to total failure is reduced to just two states.
The application of a Fault Tree may be illustrated by considering the probability of a crash at a road junction and
constructing a tree with AND and OR logic gates (figure 9). The Tree is constructed by deducing in turn the
preconditions for the top event and then successively for the next levels of events, until the basic causes are iden-
tified.
Qualitative Analysis
By using the property of the Boolean algebra it is possible first to establish the combinations of basic
(components) failures which can lead to the top (undesirable) event when occurring simultaneously.
These combinations are so called "minimal cut sets" (or "prime impliquant" ) and can be derived from
the logical equation represented by the Fault Tree.
Considering the Fault Tree representing figure 9, six scenarios can be extracted:
• Driving too fast AND Car at main road junction;
• Driver too ill AND Car at main road junction;
• Vision obscured AND Car at main road junction;
• Road too Slippery AND Car at main road junction;
• Brake failure AND Car at main road junction;
• Tyres worn AND Car at main road junction.
These 6 minimal cut sets are in first approach equivalent. However, a common cause failure analysis could
show, for example that the "Road too slippery" increase the probability of "Car at main road junction" be-
cause it is too slippery from both side. Therefore the 4th cut set is perhaps more likely than the others.
14
Note 0 – Risk analysis
Semi-Quantitative Analysis
The second step consists of calculating the probability of occurrence of each scenario. By ascribing prob-
abilities to each basic event we obtain the next figures for our example:
• Driving too fast AND Car at main road junction : 0.1 x 0.01 = 10-3
• Driver too ill AND Car at main road junction : 0.01 x 0.01 = 10-4
• Vision obscured AND Car at main road junction : 0.01 x 0.01 = 10-4
• Road too Slippery AND Car at main road junction : 0.01 x 0.01 = 10-4
• Brake failure AND Car at main road junction : 0.001 x 0.01 = 10-5
• Tyres worn AND Car at main road junction : 0.001 x 0.01 = 10-5
Total = 1.32 10-3
Now it is possible to sort the minimal cut sets in a more accurate way i.e. into three classes: One cut set
at 10-3, three at 10-4 and two at 10-5. Of course, it is better to improve the scenarios with the higher
probabilities first if we want to be efficient.
As by-product of this calculation, the global failure probability 1.32 10-3 is obtained by a simple sum of
all the individual probabilities. But this simple calculation is a conservative approximation, which
works well when the probabilities are sufficiently low (in case of safety, for example). It is less accu-
rate when the probabilities increase and it can even exceed 1 when probabilities are very high. This is
due to cross terms that are neglected. Therefore, this approach must be used with care.
These simple calculations only work on the basis of the above hypothesis. For example, as soon as the
Fault Tree contains repeated events (same events in several location in the Tree) the independence hy-
pothesis is lost. Therefore the calculation becomes wrong and even worse it is impossible to know if
the result is optimistic or pessimistic. On the other hand, the estimation of the top event probability is
less and less accurate (more and more conservative) when the probabilities increase (even if the events
are independent).
15
Note 0 – Risk analysis
Figure 10.
In-
creased Unacceptable risk
probabil-
ity
Further evaluation
or attention re-
quired
Acceptable risk
Increasing consequence →
Figure 11. Risk Matrix.
Further evaluations have to be carried out for the region between acceptable and unacceptable risk, to
determine whether further risk reduction is required or more studies should be performed.
16
Note 0 – Risk analysis
The limit of acceptability is set by defining the regions in the matrix which represent unacceptable and ac-
ceptable risk. The Risk Matrix may be used for qualitative as well as quantitative studies. If probability is
classified in broad categories such as “rare” and “frequent” and consequence in “small”, “medium”, “large”
and “catastrophic”, the results from a qualitative study may be shown in the Risk Matrix. The definition of
the categories is particularly important in case of qualitative use. The categories and the boxes in the Risk
Matrix may be replaced by continuous variables, implying a full quantification. An illustration of this is
shown in Figure 12.
Probability
Consequence
Figure 12. Risk Matrix Presentation with Continuous Variables.
The upper tolerability limit (figures 11 and 12) is almost always defined, whereas the lower limit is related
to each individual risk reducing measure, depending on when the cost of implementing each measure be-
comes unreasonably disproportional to the reduction of risk.
• Risk to safety of personnel for different solutions such as integrated versus separate quarters platform;
• Risk of operations such as exploration drilling;
• Risk of the use of a particular system such as mechanical pipe handling;
• Environmental risk.
17
Note 0 – Risk analysis
According to decision theory the alternative with the lowest expected costs should be chosen, i.e. alter-
native B should be chosen here.
The decision tree is constructed from left to right. Each consequence is associated with probabilities
(summing up to 1) after each node. For each branch the expected cost /utility is determined by multi-
plying probabilities and costs/utilities for that branch.
The “as low as reasonably practicable” (ALARP) principle is sometimes used in the industry as the only
acceptance principle and sometimes in addition to other risk acceptance criteria.
The use of the ALARP principle may be interpreted as satisfying a requirement to keep the risk level as low
as reasonably practicable, provided that the ALARP evaluations are extensively documented.
The ALARP principle is shown in Figure 14.
18
Note 0 – Risk analysis
The risk level should be reduced as far as possible in the interval between acceptable and unacceptable risk.
The common way to determine what is possible is to use cost-benefit evaluations as basis for decision on
whether to implement certain risk reducing measures or not.
The upper tolerability limit is almost always defined, whereas the lower tolerability limit is sometimes
defined, or not defined. The lower limit is individual to each individual risk reducing measure, depend-
ing on when the cost of implementing each measure becomes unreasonably disproportional to the risk
reducing effect.
The ALARP principle is normally used for risk to safety of personnel, environment and assets.
The value for the upper tolerability limit derived from accident statistics, for example, indicate that “a risk of
death around 1 in 1,000 per annum is the most that is ordinarily accepted by a substantial group of workers in
any industry in the UK”.
HSE (Health and Safety Excecutive), [8] suggested the upper maximum tolerable risk level as a line with a slope
of –1 through point n = 500 (number of fatalities), F = 2 x 10-4 (frequency) per year. This line corresponds to n
= 1 at F = 10-1 per year, and n = 100 at F = 10-3 per year. However, in the document [9], HSE quotes that risk of
a single accident causing the death of 50 people or more with the frequency of 5 x 10-3 per year is intolerable.
For the “negligible” level, the HSE recommends a line drawn three decades lower than the intolerable line. This
line corresponds to one fatality, n = 1, in one per ten thousand per year, F = 10-4 per year, and similarly, n = 100
corresponds to one in a million per year, F = 10-6 per year.
Railway transport
For the railway area different railway operators in the UK has suggested the risk criteria in table 2.
19
Note 0 – Risk analysis
Road Transport
The current convention adopted by Department of Transport (DTLR) for the Value of Statistical Life (VOSL) is
about £1,000,000 (2001 prices). This value is used in a cost-benefit evaluation of road improvement schemes in
the UK. There are no risk criteria specific to road transport, but only to transport of hazardous material by road
[8].
It is interesting to note that according to [8], the national scrutiny level of societal risk for road transport could be
defined by a line passing through points defined by 10 fatalities corresponding to 0.7 times per year (7 in 10
years), 100 fatalities corresponding to 0.07 times per year (7 in 100 years), and 1,000 fatalities corresponding to
0.007 times per year (7 in 1,000 years).
Recently developed criteria for a tunnel under a river were based on the HSE’s suggestion in [7], which defines
the upper maximum tolerable risk level as a line with a slope of –1 passing through points N = 1 and F = 10-1 per
year, and N = 100 and F = 10-3 per year. For the negligible level a line drawn three decades lower is suggested
by the HSE. This line corresponds to one fatality, N = 1, in 1 per 10,000 per year (F = 10-4 per year), and simi-
larly, N = 100 corresponds to 1 in a million per year (F = 10-6 per year). For the broadly acceptable level a line
two decades lower that the maximum tolerable line is suggested. This line is defined as N = 1 and F = 10-3 per
year, and N = 100 and F = 10-5 per year. The proposed criteria for the societal risk for the tunnel are presented
in figure 15.
The Dutch criteria are compared with the proposed UK criteria in figure 16. It can be seen that up to N = 100,
the both criteria can be assumed to have similar impact, in spite of the fact that the Dutch maximum tolerable
level is more stringent that the proposed one for the UK. The equalising factor is the ALARP requirement in the
UK which drives the risk towards the broadly acceptable level.
20
Note 0 – Risk analysis
1.E+0
Negligible Level
1.E-1
Maximum Tolerable
Annual Frequency of N or more Fatalities
Level
Broadly Acceptable
1.E-2
1.E-3
1.E-4
1.E-5
1.E-6
1.E-7
0.1 1 10 100 1000
Number of Fatalities N
1.E-1
1.E-2
Annual Frequency of N or more Fatalities
Broadly Acceptable
Negligible Level
1.E-3
Dutch Maxium Tolerable Level
1.E-4
1.E-5
1.E-6
1.E-7
1 10 100 1000
Number of Fatalities N
Figure 16. Comparison of proposed and Dutch criteria for road tunnel users.
21
Note 0 – Risk analysis
Risk criteria were developed for a 6,600 m long tunnel under the River Westerschelde completed in 2002. No
risk criteria for road users were available. So the criteria for population in the vicinity of the road used for trans-
portation of hazardous substances have been used as the basis. These criteria define the upper tolerability level
as the individual risk per annum of 10-6, and the societal risk of 10 fatalities once in 10,000 years, and 100 fatali-
ties once in a million years. Since the risk to road users can be considered as voluntary, a factor of 10 has been
used for the criteria for road users [10], as presented in table 4.
MaximumAnnual Frequency
Number of Fatalities
per km
1 10-1
10 10-3
100 10-5
Table 4. Societal risk criterion used for Westerschelde tunnel.
It is interesting to note that the risk criterion for the overall length of the tunnel is also mentioned, but it is not
clear if this was supposed to be used to “average” risks along the tunnel, as it is also mentioned that risks vary
significantly throughout the tunnel (i.e. exits and entrances and their slopes present a considerably greater risk
than horizontal sections).
Building sector
Buildings are usually designed for a life time of 50 years and in practice most buildings have a much
longer life, also because of maintenance and refurbishment activities which also for a part of the build-
ing sector. A number of 80 years is often mentioned as an average. For some buildings, of course, a
smaller design life is chosen on purpose, especially for some kind of industrial buildings.
22
Note 0 – Risk analysis
Usually they do not belong to one firm and this may make communication quite difficult. It also contri-
butes to the difficulty to make successful innovations in this sector. The following hazards are to be
considered:
Design usually takes care of all those issues, either implicitly or explicitly. The degree in which atten-
tion is paid depends very much on the nature and importance of the building.
Most buildings are designed without any special and explicit risk analysis and risk management. Im-
plicitly, of course, the building codes can be considered as a relatively effective risk management tool.
Special attention is given to the aspect of fire safety, both the structural side and the human safety side.
Most people in fires do not die from fire or structural collapse, but from smoke. The most important
issue over here is a good system of detection, warning, early extinguishing and good possibilities of
evacuation and escape. In this field an increasing interest in explicit risk analysis approach can be ob-
served.
In most buildings codes safety classes are distinguished for various types of buildings and various types
of building elements. A very common system is to divide the consequences of failure into small, me-
dium and large. Some codes make a distinction between economic losses and danger to human lives.
Based on these classes there is a safety differentiation in the code, which usually comes down to the
prescription of different values for the partial factors or the introduction of a set of adjustment factors.
It is well known, but not always fully appreciated, that the reliability of a structure as estimated on the
basis of a given set of probabilistic models for loads and resistances may have limited bearing to the
actual reliability of the structure. This is the case when the probabilistic modeling forming the basis of
the reliability analysis is highly influenced by subjectivity and then the estimated reliability should be
interpreted as being a measure for comparison only. In these cases it is thus not immediately possible to
judge whether the estimated reliability is sufficiently high without first establishing a more formalized
reference for comparison.
Such a reference may be established by the definition of an optimal or best practice structure. The idea
behind the "best practice" reference is that if the structure of consideration has been designed according
to the "best practice" then the reliability of the structure is "optimal" according to agreed conventions
for the target reliability. Typical values for the corresponding target annual failure probability are in the
range of 10-6 to 10-7 depending on the type of structure and the characteristics of the considered failure
23
Note 0 – Risk analysis
mode. Using this approach the target reliability is determined as the reliability of the "best practice"
design as assessed with the given probabilistic model.
The determination of the "best practice" design can be performed in different ways. The simplest ap-
proach is to use the existing codes of practice for design as a basis for the identification of "best prac-
tice" design. Alternatively the "best practice design" may be determined by consultation of a panel of
recognized experts.
In case where the probabilistic modeling does not rest on subjective assessments the most rational ap-
proach is to establish the optimal design on the basis of the economic decision theory.
In tables 1 target failure probabilities are given for ultimate limit states based on recommendations of
JCSS [11]. Note that the values given correspond to a year reference period and the stochastic models
recommended in JCSS [11].
7 References
[1] Melchers, R.E.: Structural Reliability, analysis and prediction. John Wiley & Sons, New York,
1987.
[2] Faber, M.H. & M.G. Stewart: Risk analysis for civil engineering facilities: overview and discus-
sion. JCSS (Joint Committee on Structural Safety), 2001.
[3] Raiffa, H. & R. Schlaifer: Applied Statistical Decision Theory. Harward University Press, Cam-
bridge University Press, Cambridge, Mass., 1961.
[4] Benjamin, J.R. & C.A. Cornell: Probability, Statistics and Decision for Civil Engineers. McGraw-
Hill, NY, 1970.
[5] Ditlevsen, O. & H.O. Madsen: Structural Reliability Methods. Wiley, Chichester, 1996.
[6] Stewart, M.G. & R.E. Melchers: Probabilistic Risk Assessment of Engineering Systems. Chap-
man & Hall, 1997.
[7] Saferelnet: Quantitative risk analysis – report on current practice. 2002.
[8] HSC: Major hazard aspects of the transport of dangerous substances. Advisory Committee on
Dangerous Substances, HMSO, London 1991.
[9] HSE: Reducing Risk, Protecting People. HMSO, 2001.
[10] Worm, E.W. & J. Hoeskma: The Westerschelde Tunnel: Development and application of am in-
tegrated safety philosophy. Safety in Road and Rail Tunnels, 3rd International Conference organ-
ised by University of Dundee and ITC Ltd., Nice, France, 9-11 March, 1998.
[11] JCSS: Probabilistic model code. The Joint Committee on Structural Safety, 2001.
http://www.jcss.ethz.ch/
24
Note 0 – Risk analysis
Opgave 1 – Risikoanalyse
I forbindelse med en eksisterende bro er der opstået tvivl om den har tilstrækkelig bæreevne. Derfor
ønskes et svigttræ og et hændelsestræ opstillet. For at kunne vurdere forskellige alternative handlinger
ønskes også et beslutningstræ opstillet.
Hvis derimod svigtet sker øjeblikkeligt som et totalt kollaps indtræder, så er omkostningerne 10 mill.
kr. På basis af trafikobservationer vides at sandsynligheden for at der ingen køretøjer er på broen er
95% i svigtøjeblikket. Med en sandsynlighed på 5% er der 10 køretøjer og 20 personer på broen på et
tilfældigt tidspunkt. Konsekvenserne af tab af et menneskeliv antages at kunne ækvivaleres med 4 mill.
kr.
Opstil et hændelsestræ og vis at de forventede omkostninger ved svigt af broen er 4.55 mill. kr.
Bestem de forventede omkostninger for hver af de 3 alternativer. Hvilket alternativ bør vælges.
De forventede omkostningerne for et svigt af den eksisterende bro tages fra spørgsmål 2.
25
Note 0 – Risk analysis
26
Note 1+2: Structural reliability
1 Introduction
For many years it has been assumed in design of structural systems that all loads and strengths are
deterministic. The strength of an element was determined in such a way that it exceeded the load
with a certain margin. The ratio between the strength and the load was denoted the safety factor.
This number was considered as a measure of the reliability of the structure. In codes of practice for
structural systems values for loads, strengths and safety factors are prescribed.
These values are traditionally determined on the basis of experience and engineering judgement.
However, in new codes partial safety factors are used. Characteristic values of the uncertain loads
and resistances are specified and partial safety factors are applied to the loads and strengths in order
to ensure that the structure is safe enough. The partial safety factors are usually based on experience
or calibrated to existing codes or to measures of the reliability obtained by probabilistic techniques.
As described above structural analysis and design have traditionally been based on deterministic
methods. However, uncertainties in the loads, strengths and in the modeling of the systems require
that methods based on probabilistic techniques in a number of situations have to be used. A struc-
ture is usually required to have a satisfactory performance in the expected lifetime, i.e. it is required
that it does not collapse or becomes unsafe and that it fulfills certain functional requirements. Gen-
erally structural systems have a rather small probability that they do not function as intended, see
table 1.
Reliability of structural systems can be defined as the probability that the structure under considera-
tion has a proper performance throughout its lifetime. Reliability methods are used to estimate the
probability of failure. The information of the models which the reliability analyses are based on are
27
Note 1+2: Structural reliability
generally not complete. Therefore the estimated reliability should be considered as a nominal meas-
ure of the reliability and not as an absolute number. However, if the reliability is estimated for a
number of structures using the same level of information and the same mathematical models, then
useful comparisons can be made on the reliability level of these structures. Further design of new
structures can be performed by probabilistic methods if similar models and information are used as
for existing structures which are known to perform satisfactory. If probabilistic methods are used to
design structures where no similar existing structures are known then the designer has to be very
careful and verify the models used as much as possible.
The reliability estimated as a measure of the safety of a structure can be used in a decision (e.g. de-
sign) process. A lower level of the reliability can be used as a constraint in an optimal design prob-
lem. The lower level of the reliability can be obtained by analyzing similar structures designed after
current design practice or it can be determined as the reliability level giving the largest utility (bene-
fits – costs) when solving a decision problem where all possible costs and benefits in the expected
lifetime of the structure are taken into account.
In order to be able to estimate the reliability using probabilistic concepts it is necessary to introduce
stochastic variables and/or stochastic processes/fields and to introduce failure and non-failure be-
havior of the structure under consideration.
Typical failure modes to be considered in a reliability analysis of a structural system are yielding,
buckling (local and global), fatigue and excessive deformations.
28
Note 1+2: Structural reliability
can be related to e.g. formation of a mechanism in the structure, exceedance of the material strength
or instability (buckling).
The fundamental quantities that characterize the behavior of a structure are called the basic vari-
ables and are denoted X = ( X 1 ,..., X n ) where n is the number of basic stochastic variables. Typical
examples of basic variables are loads, strengths, dimensions and materials. The basic variables can
be dependent or independent, see below where different types of uncertainty are discussed. A sto-
chastic process can be defined as a random function of time such that for any given point in time the
value of the stochastic process is a random variable. Stochastic fields are defined in a similar way
where the time is exchanged with the space.
The uncertainty modeled by stochastic variables can be divided in the following groups:
Physical uncertainty: or inherent uncertainty is related to the natural randomness of a quantity, for
example the uncertainty in the yield stress due to production variability.
The above types of uncertainty are usually treated by the reliability methods which will be de-
scribed in the following chapters. Another type of uncertainty which is not covered by these meth-
ods are gross errors or human errors. These types of errors can be defined as deviation of an event
or process from acceptable engineering practice.
Generally, methods to measure the reliability of a structure can be divided in four groups, see Mad-
sen et al. [2], p.30:
• Level I methods: The uncertain parameters are modeled by one characteristic value, as for ex-
ample in codes based on the partial safety factor concept.
• Level II methods: The uncertain parameters are modeled by the mean values and the standard
deviations, and by the correlation coefficients between the stochastic variables. The stochastic
variables are implicitly assumed to be normally distributed. The reliability index method is an
example of a level II method.
• Level III methods: The uncertain quantities are modeled by their joint distribution functions.
The probability of failure is estimated as a measure of the reliability.
• Level IV methods: In these methods the consequences (cost) of failure are also taken into ac-
count and the risk (consequence multiplied by the probability of failure) is used as a measure of
the reliability. In this way different designs can be compared on an economic basis taking into
account uncertainty, costs and benefits.
Level I methods can e.g. be calibrated using level II methods, level II methods can be calibrated
using level III methods, etc.
29
Note 1+2: Structural reliability
Level II and III reliability methods are considered in these notes. Several techniques can be used to
estimate the reliability for level II and III methods, e.g.
• simulation techniques: Samples of the stochastic variables are generated and the relative num-
ber of samples corresponding to failure is used to estimate the probability of failure. The simu-
lation techniques are different in the way the samples are generated.
• FORM techniques: In First Order Reliability Methods the limit state function (failure func-
tion) is linearized and the reliability is estimated using level II or III methods.
• SORM techniques: In Second Order Reliability Methods a quadratic approximation to the
failure function is determined and the probability of failure for the quadratic failure surface is
estimated.
In level IV methods the consequences of failure can be taken into account. In cost-benefit analyses
(or RISK analyses) the total expected cost-benefits for a structure in its expected lifetime are maxi-
mized
max W ( z ) = B( z ) − C I ( z ) − C IN ( z ) − C REP ( z ) − Pf ( z )C F (1)
z
where z represents design/decision variables, B is the expected capitalized benefits, C I is the ini-
tial (or construction) costs, C IN is the expected capitalized inspection costs, C REP is the expected
capitalized repair costs and C F is the capitalized failure costs. Cost-optimized inspection strategies
are based on cost-benefit analyses where the costs due to inspection, repair and failure are mini-
mized with e.g. inspection locations and inspection times and qualities as decision variables.
For a detailed introduction to structural reliability theory reference is made to the following text-
books: Melchers [1], Madsen, Krenk & Lind [2], Thoft-Christensen & Baker [3] and Ditlevsen &
Madsen [4].
If failure of one element gives failure of the system, then a union (series system) is used to model
the system failure, E :
m
E = E1 ∪ ... ∪ Em = U Ei (2)
i =1
where Ei is the event modeling failure of element i and m is the number of events.
If failure of all elements are needed to obtain failure of the system, then an intersection (parallel
system) is used to model the system failure, E :
m
E = E1 ∩ ... ∩ E m = I Ei (3)
i =1
30
Note 1+2: Structural reliability
P U Ei = ∑ P ( Ei )
m m
(10)
i =1 i =1
Using the multiplication rule in (13) and considering mutually exclusive events E1 , E 2 ,..., E m the
total probability theorem follows:
P ( A) = P( A E1 )P( E1 ) + P(A E 2 )P( E 2 ) + ... + P(A E m )P( E m )
(15)
= P( A ∩ E1 ) + P( A ∩ E 2 ) + ... + P( A ∩ E m )
where A is an event.
31
Note 1+2: Structural reliability
Using also the total probability theorem in (15) the so-called Bayes theorem follows from:
P (A Ei )P( Ei ) P(A Ei )P( Ei )
P (Ei A) = = m (17)
P( A)
(
∑ P A E j P( E j ) )
j =1
32
Note 1+2: Structural reliability
Consider a structural element with load bearing capacity R which is loaded by the load S . R and
S are modeled by independent stochastic variables with density functions f R and f S and distribu-
tion functions FR and FS , see figure 3. The probability of failure becomes
33
Note 1+2: Structural reliability
∞ ∞
PF = P ( failure ) = P (R ≤ S ) = ∫ P (R ≤ x )P ( x ≤ S ≤ x + dx )dx = ∫ FR ( x ) f S ( x )dx
−∞ −∞
It is noted that it is important that the lower part of the distribution for the strength and the upper
part of the distribution for the load are modeled as accurate as possible.
Figure 4 shows the density function for a Normal distributed stochastic variable with expected value
10 and standard deviation 3.
Relative Frequency
Density Plots (1 Graphs) - [w]
0.1329
X1
Normal (Gauss)
0.1197
0.1066
0.0934
0.0802
0.0670
0.0538
0.0407
0.0275
0.0143
0.0011
0.0000 2.0000 4.0000 6.0000 8.0000 10.0000 12.0000 14.0000 16.0000 18.0000 20.0000
Value of X
Figure 4. Normal distributed stochastic variable with expected value 10 and standard deviation 3.
34
Note 1+2: Structural reliability
σ 2
σ Y = ln + 1 (27)
µ
1
µ Y = ln µ − σ Y2 (28)
2
is the standard deviation and expected value for the Normal distributed stochastic variable
Y = ln X .
0.1185
0.1038
0.0890
0.0743
0.0595
0.0447
0.0300
0.0152
0.0005
0.0000 2.0000 4.0000 6.0000 8.0000 10.0000 12.0000 14.0000 16.0000 18.0000 20.0000
Value of X
Figure 5. Lognormal distributed stochastic variable with expected value 10 and standard deviation
3.
σ
If the coefficient of variation V = is small (less than ≈ 0.25) then the standard deviation and
µ
expected value of Y = ln X can approximately be obtained from
σY ≈V
µ Y ≈ ln µ
The 5 % quantile y0.05 defined such that P (Y = ln X ≤ y0.05 ) = 0.05 can be obtained from
y0.05 = µ Y − 1.645σ Y
where –1.645 is obtained from the standard Normal distribution function such that
Φ (−1.645) = 0.05 . Correspondly the 5% qualtile x0.05 of the LogNormal variable can be obtained
from
35
Note 1+2: Structural reliability
Figure 6 shows the density function for a 2-parameter Weibull distributed stochastic variable with
expected value 10 and standard deviation 3.
Relative Frequency
Density Plots (1 Graphs) - [w]
0.1283
X3
Weibull (min)
0.1156
0.1029
0.0902
0.0775
0.0648
0.0522
0.0395
0.0268
0.0141
0.0014
0.0000 2.0000 4.0000 6.0000 8.0000 10.0000 12.0000 14.0000 16.0000 18.0000 20.0000
Value of X
Figure 6. 2-parameter Weibull distributed stochastic variable with expected value 10 and standard
deviation 3.
36
Note 1+2: Structural reliability
Relative Frequency
Density Plots (1 Graphs) - [w]
0.1278
X4
Weibull (min)
0.1151
0.1024
0.0898
0.0771
0.0644
0.0517
0.0390
0.0264
0.0137
0.0010
0.0000 2.0000 4.0000 6.0000 8.0000 10.0000 12.0000 14.0000 16.0000 18.0000 20.0000
Value of X
Figure 7. 3-parameter Weibull distributed stochastic variable with expected value 10, standard de-
viation 3 and lower threshold γ =3.
where α and β are the parameters. The allowable intervals for the parameters are:
37
Note 1+2: Structural reliability
3(1 + 2 β )(3 − β + 2 β 2 )
Kurtosis: β2 =
(1 + 3β )(1 + 4 β )
where α and β are shape and scale parameters. These are related to µ and σ by:
0.5772
µ=β+ (42)
α
π
σ= (43)
α 6
Figure 8 shows the density function for a Gumbel distributed stochastic variable with expected
value 10 and standard deviation 3.
Relative Frequency
Density Plots (1 Graphs) - [w]
0.1572
X5
Gumbel (max)
0.1415
0.1258
0.1102
0.0945
0.0788
0.0631
0.0475
0.0318
0.0161
0.0004
0.0000 2.0000 4.0000 6.0000 8.0000 10.0000 12.0000 14.0000 16.0000 18.0000 20.0000
Value of X
Figure 8. Gumbel distributed stochastic variable with expected value 10 and standard deviation 3.
1 6
x0.98 = u − ln (− ln (0.98)) = µ 1 − V (0.5772 + ln(− ln(0.98)))
α π
The distribution function for the maximum Y = max{X 1 , X 2 ,..., X n } of n independent Gumbel dis-
tributed stochastic variables X 1 , X 2 ,..., X n becomes (see also section 3):
38
Note 1+2: Structural reliability
6 6
µY = µ + σ ln(n) = µ 1 + V ln(n)
π π
σY = σ
It is seen that
Cov[ X 1 , X 1 ] = Var[ X 1 ] = σ 12 (47)
− 1 ≤ ρ X1 , X 2 ≤ 1 (49)
39
Note 1+2: Structural reliability
i =1 j =1
Finally, it can be shown that if X 1 , X 2 ,L, X n are Normal distributed then Y is also Normal distrib-
uted.
The following general comments can be made in relation to choice of distribution functions.
For extreme loads for example the annual maximum / extreme value of the load (wind velocity /
wind pressure, significant wave height, ...) is the important value. If X 1 , X 2 ,L, X n are independent
stochastic variables with identical distribution function FX then the maximum value
Y = max{X 1 , X 2 ,L, X n }
has the distribution function
FY ( y ) = P (Y ≤ y ) = P (max{X 1 , X 2 ,L, X n } ≤ y )
= P({X 1 ≤ y}∩ {X 2 ≤ y}∩ ... ∩ {X n ≤ y}) = (FX ( y ) )
n
40
Note 1+2: Structural reliability
Figure 9 illustrates the density functions f X ( y ) and f Y ( y ) . It is seen that as n increases the den-
sity function for the maximum value becomes more narrow (smaller coefficient of variation) and
the expected value increases. It can be shown that in general f Y ( y ) approaches one of the so-called
extreme distribution.
For fatigue analysis where a good fit is important for the central part of the distribution of the load
variations (stress ranges) the following distribution types will be relevant for wind velocities and
significant wave heights:
• Normal distribution
• LogNormal distribution
• Weibull distribution. This distribution is used e.g. in Windatlas [10]
A number of methods can be used to estimate the statistical parameters in distribution functions, for
example:
• The Maximum Likelihood method
• The Moment method
• The Least Square method
• Bayesian statistics
41
Note 1+2: Structural reliability
In general the Maximum Likelihood method or Bayesian statistics is recommended. The Maximum
Likelihood method gives a consistent estimate of the statistical uncertainties. In Bayesian statistics
it is possible to take consistently into account subjective / prior information.
The statistical parameters are α , β and γ . xi , i = 1, n are n data values. The optimization prob-
lem to obtain the maximum value of the log-Likelihood function, max ln L(α , β , γ ) can be solved
α , β ,γ
using non-linear optimization algorithms, e.g. NLPQL, [11]. The result is the best estimate of the
statistical parameters α , β and γ .
Since the statistical parameters α , β and γ are determined from a limited number of data, the es-
timates will be subjected with statistical uncertainty. If the number of data is larger than 25-30 α ,
β and γ can be assumed to be asymptotically Normal distributed with expected values equal to the
solution of the optimization problem and with the following covariance matrix, see e.g. Lindley, [4]
σ α2 ρ αβ σ α σ β ρ αγ σ α σ γ
[ ]
−1
Cα , β ,γ = − H αβγ = ρ αβ σ α σ β σ β2
ρ βγ σ β σ γ (54)
ρ αγ σ α σ γ ρ βγ σ β σ γ σγ 2
where H αβγ is the Hessian matrix with second derivatives of the Log-Likelihood function. σ α ,
σ β and σ γ are standard deviations of α , β and γ . ρ αβ is the correlation coefficient between α
and β .The Hessian matrix is determined by numerical differentiation.
3.2 Moment method
The unknown parameters in a given distribution function FX ( x θ ) for a stochastic variable X is
denoted θ = (θ1 ,θ 2 ,...,θ m ) . The theoretical statistical moments are with given θ = (θ1 ,θ 2 ,...,θ m )
m j = ∫ x j f X ( x θ )dx (55)
On the basis of data / observations x$ = ( x$1 , x$ 2 ,..., x$ n ) the empirical moments are
1 n
mˆ j = ∑ xˆi (56)
n i =1
42
Note 1+2: Structural reliability
Requiring that the theoretical moments are equal to the empirical moments the statistical parameters
θ = (θ1 ,θ 2 ,...,θ m ) can be determined.
It is noted that the method does not give an estimate of the statistical uncertainties and that it is not
possible to include prior information. However, bootstrapping can in some situations be used to
estimate the statistical uncertainties.
On the basis f data / observations x$ = ( x$1 , x$ 2 ,..., x$ n ) an empirical distribution function is deter-
mined using e.g. the Weibull – plotting formula:
i
Fˆi = , x = xˆ i (57)
n +1
The solution of this optimization problem is a central estimate of the statistical parameters
θ = (θ1 ,θ 2 ,...,θ m ) .
If the distribution has to be good in the tails of the distribution the summation in (58) can be re-
duced to e.g. the smallest 30% of the data.
Consider a stochastic variable X with distribution function FX ( x θ ) which depends on the statisti-
cal parameters θ = (θ 1 ,θ 2 ,...) . For a Normal distribution the statistical parameters are equal to the
expected value and the standard deviation.
It is assumed that one or more of the statistical parameters are uncertain, and that prior knowledge
on this uncertainty can be expressed in a prior density function for parameters: fθ' (θ ) .
If data is available these can be used to update this prior knowledge. The updated – posterior den-
sity function for the statistical parameters can be determined by
f X (xˆ θ ) fθ' (θ )
fθ (θ xˆ ) =
''
(59)
∫ f X (xˆ θ ) f θ (θ )dθ
'
n
where f X (xˆ θ ) = ∏ f X ( xˆ i θ ) is the probability (Likelihood) for the given data / observations
i =1
The predictive (updated) density function for X given data x$ = ( x$1 , x$ 2 ,..., x$ n ) is determined by
43
Note 1+2: Structural reliability
Prior, posterior and predictive distributions can be established in e.g. the following cases:
• Normal distribution with known standard deviation
• Normal distribution with known expected value
• Normal distribution with unknown expected value and standard deviation
• Lognormal distribution
• Gumbel distribution
• Weibull distribution
• Exponential distribution
It is noted that statistical uncertainty automatically is included in this modeling and that engineering
judgments based on experience can be quantified rationally via prior distributions of the statistical
parameters θ . Further it can be mentioned that in Eurocode 0, Basis of Design, annex D, [1] it is
recommended that Bayesian statistics is used in statistical treatment of data and in design based on
tests.
44
Note 1+2: Structural reliability
COV x0.05
Non-parametric 0.26 21.6
Normal 0.26 22.4
Normal – tail 0.25 22.7
LogNormal 0.28 24.1
LogNormal - tail 0.38 22.8
Weibull-2p 0.27 21.3
Weibull-2p - tail 0.23 22.8
Weibull-3p 0.26 23.3
Weibull-3p - tail
Table b. Statistical data (in MPa).
1.00 1.00
0.90 0.90
0.80 0.80
0.70 0.70
0.60 0.60
0.50 0.50
0.40 0.40
0.30 0.30
0.20 0.20
0.10 0.10
0.00 0.00
0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00
LT20: κ =30% truncation LT20: κ =100% truncation
1.00 1.00
0.90 0.90
0.80 0.80
0.70 0.70
0.60 0.60
0.50 0.50
0.40 0.40
0.30 0.30
0.20 0.20
0.10 0.10
0.00 0.00
0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00
45
Note 1+2: Structural reliability
1.00 1.00
0.90 0.90
0.80 0.80
0.70 0.70
0.60 0.60
0.50 0.50
0.40 0.40
0.30 0.30
0.20 0.20
0.10 0.10
0.00 0.00
0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00
1.00 1.00
0.90 0.90
0.80 0.80
0.70 0.70
0.60 0.60
0.50 0.50
0.40 0.40
0.30 0.30
0.20 0.20
0.10 0.10
0.00 0.00
0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00
46
Note 1+2: Structural reliability
The parameters α and β are determined using available data and are thus subject to statistical un-
certainty. If the parameters are estimated by the Maximum Likelihood technique the uncertainty can
be quantified and included in the stochastic model.
Data obtained by continuous simulations of significant wave heights and wave directions for the
central part of the North Sea covering the period 1979 to 1993 are used. All storm events are identi-
fied and the maximum significant wave height within the eight directional sectors: N, NE, E, S,
SW, W, and NW are determined. The simulated wave heights have been calibrated against available
measurements from the same location. The calibrated statistical parameters and other optimal pa-
rameters are shown in table a together with estimates of the characteristic 100 year wave heights,
H S ,100 . In figure b and c empirical and fitted distribution functions are shown for the omnidirec-
tional, SouthWest, West and NorthWest directions.
αj βj λj γj H S ,100
N 3.06 4.25 m 1.20 4.0 m 7.8 m
NE 2.55 2.93 m 1.40 3.0 m 5.9 m
E 3.23 4.36 m 1.60 4.0 m 7.6 m
SE 3.00 3.90 m 1.07 4.0 m 7.0 m
S 3.53 4.75 m 1.53 5.0 m 8.1 m
SW 4.97 6.23 m 1.47 6.5 m 9.5 m
W 6.03 6.90 m 2.20 6.0 m 8.8 m
NW 4.98 6.25 m 1.80 5.25 m 9.1 m
Omni 5.52 6.64 m 3.73 6.0 m 9.2 m
Table a. Estimated statistical parameters and characteristic significant wave height: H S ,100 .
1.00 1.00
0.80 0.80
0.60 0.60
0.40 0.40
0.20 0.20
0.00 0.00
5.50 6.00 6.50 7.00 7.50 8.00 8.50 9.00 6.00 6.50 7.00 7.50 8.00 8.50
Figure b. Southwest (left) and West (right) empirical and fitted distribution functions.
47
Note 1+2: Structural reliability
1.00 1.00
0.80 0.80
0.60 0.60
0.40 0.40
0.20 0.20
0.00 0.00
5.50 6.00 6.50 7.00 7.50 8.00 8.50 6.00 6.50 7.00 7.50 8.00 8.50 9.00
Figure c. Northwest (left) and Omnidirectional empirical and fitted distribution functions.
6 References
[1] Melchers, R.E.: Structural Reliability, analysis and prediction. John Wiley & Sons, New
York, 1987.
[2] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
[3] Thoft-Christensen, P. and M.J. Baker: Structural Reliability Theory and Its Applications.
Springer Verlag, 1982.
[4] Ditlevsen, O. & H.O. Madsen: Structural Reliability Methods. Wiley, 1996.
[5] JCSS (2002). Joint Committee on Structural Safety: Model code. www.jcss.ethz.ch
[6] DS 410: Code of Practice for Loads for the Design of Structures. DS 1998.
[7] EN 1990 EN 1990 (2000). Eurocode, Basis of Structural Design, EN1990. Draft, December
2000.
[8] ISO 2394. General principles on reliability for structures. 1998.
[9] van Gelder (2000). van Gelder, P.H.A.J.M.: Statistical methods for the Risk-Based Design of
Civil Structures. PhD thesis, Delft University of Technology, Holland.
[10] Windatlas (1989). Troen, I., B. Petersen & E. Lundtag: European Wind Atlas, Risø, Roskilde.
[11] NLPQL (1986). Schittkowski, K.: NLPQL: A FORTRAN Subroutine Solving Non-Linear
Programming Problems. Annals of Operations Research.
[12] Sørensen, J.D. & P. Hoffmeyer: Statistical analysis of data for timber strengths. Report, Aal-
borg University, 2001.
[13] Sørensen, J.D., M. Sterndorff & A. Bloch: Reliability Analysis of Offshore Structures Using
Directional Loads. Proc. ASCE Joint Specialty Conf. on ‘Probabilistic Mechanics and Struc-
tural reliability’, Notre Dame, Indianapolis, July 2000.
48
Note 3: First order reliability methods
3.1 Introduction
In this section the problem of estimating the reliability or equivalently the probability of failure is
considered. Generally, methods to measure the reliability of a structure can be divided into four
groups, see Madsen et al. [3.1], p.30:
• Level I methods: The uncertain parameters are modelled by one characteristic value, as for ex-
ample in codes based on the partial coefficients concept.
• Level II methods: The uncertain parameters are modelled by the mean values and the standard
deviations, and by the correlation coefficients between the stochastic variables. The stochastic
variables are implicitly assumed to be normally distributed. The reliability index method is an
example of a level II method.
• Level III methods: The uncertain quantities are modelled by their joint distribution functions.
The probability of failure is estimated as a measure of the reliability.
• Level IV methods: In these methods the consequences (cost) of failure are also taken into ac-
count and the risk (consequence multiplied by the probability of failure) is used as a measure of
the reliability. In this way different designs can be compared on an economic basis taking into
account uncertainty, costs and benefits.
If the reliability methods are used in design they have to be calibrated so that consistent reliability
levels are obtained. This is further discussed in a later note.
Level I methods can e.g. be calibrated using level II methods, level II methods can be calibrated
using level III methods, etc.
In this note level II and III reliability methods are considered. Several techniques can be used to
estimate the reliability for level II and III methods, e.g.
• Simulation techniques: Samples of the stochastic variables are generated and the relative num-
ber of samples corresponding to failure is used to estimate the probability of failure. The simula-
tion techniques are different in the way the samples are generated. Simulation techniques are de-
scribed in note 5.
• FORM techniques: In First Order Reliability Methods the limit state function (failure function,
see below) is linearized and the reliability is estimated using level II or III methods. FORM
techniques for level II methods are described in this note. FORM techniques for level III meth-
ods are described in note 4.
• SORM techniques: In Second Order Reliability Methods a quadratic approximation to the fail-
ure function is determined and the probability of failure for the quadratic failure surface is esti-
mated. SORM techniques are discussed in note 5.
In section 3.2 basic variables and failure functions are defined. Next, a linear failure function is
considered in section 3.3 and the reliability index β is defined. In section 3.4 non-linear failure
functions are considered. The so-called invariance problem is discussed, and the Hasofer & Lind
reliability index β is defined. A numerical algorithm for determination of the reliability index is
49
Note 3: First order reliability methods
shown. Finally it is shown how a sensitivity analysis of the reliability index with respect to a deter-
ministic parameter can be performed.
The joint density function for the stochastic variables X is denoted f X (x) . The elements in the
vector of expected values and the covariance vector are:
μ i = E [ X i ] , i = 1,K, n (3.1)
Cij
ρ ij = , i, j = 1,K, n (3.3)
σ iσ j
Application of FORM, SORM and simulation methods requires as noted above that it is possible for
given realizations x of the basic variables to state whether the structure (or component/failure
mode) is in a safe state or in a failure state. The basic variable space is thus divided into two sets,
the safe set ω S and the failure set ω F . The two sets are separated by the failure surface (limit state
surface). It is assumed that the failure surface can be described by the equation:
g (x) = g ( x1 ,K, xn ) = 0
Usually the failure function is defined such that positive values of g correspond to safe states and
negative values correspond to failure states, see figure 3.1.
50
Note 3: First order reliability methods
⎧ > 0 , x ∈ωs
g ( x) ⎨ (3.4)
⎩≤ 0 , x ∈ ω f
It is important to note that the failure surface does not define a unique failure function, i.e. the fail-
ure surface can be described by a number of equivalent failure functions. However, whenever pos-
sible, differentiable failure functions should be used. In structural reliability the failure function
usually results from a mechanical analysis of the structure.
If, in the failure function x is replaced by the stochastic variables X , the so-called safety margin M
is obtained:
M = g (X) (3.5)
Pf = P( M ≤ 0) = P( g ( X) ≤ 0) = ∫ω f X (x)dx
f
(3.6)
Example 3.1
In the fundamental case only two basic variables are used, namely the load variable P and the
strength variable S. A failure function can then be formulated as:
g ( s, p ) = s − p (3.7)
The failure surface g ( s, p) = 0 is shown in figure 3.2. The safety margin corresponding to (3.7) is:
M =S−P (3.8)
Instead of the failure function (3.7) the following equivalent failure function can be used:
g ( s, p ) = s 3 − p 3 (3.9)
51
Note 3: First order reliability methods
M = a0 + a1 X 1 + L + an X n (3.10)
where a0 , a1 ,K, a n are constants. The expected value μ M and the standard deviation σ M are:
μ M = a0 + a1μ x + L + an μ x = a0 + aT μ X
1 n
(3.11)
σ M = a T Ca (3.12)
As a measure of the reliability of a component with the linear safety margin (3.10) the reliability
index β can be used:
μM
β= (3.14)
σM
If the basic variables are normally distributed and the safety margin is linear then M becomes nor-
mally distributed. The probability of failure is, see figure 3.3:
52
Note 3: First order reliability methods
⎛ μ ⎞
Pf = P( M ≤ 0) = P(μ M + Uσ M ≤ 0 ) = P⎜⎜U ≤ − M ⎟⎟ = Φ (− β ) (3.15)
⎝ σM ⎠
where Φ is the standard normal distribution function and U is a standard normally distributed vari-
able with expected value zero and unit standard deviation ( μU = 0, σ U = 1) .
Figure 3.3. Illustration of reliability index and probability of failure. ϕ is the standard normal den-
sity function.
Example 3.2
Consider the fundamental case with the linear failure function (3.7). If the stochastic variables P
and S are independent then the reliability index becomes:
μM μ − μP
β= = S
σM σ S2 + σ P2
Assume that P and S are normally distributed with expected values μ P = 2, μ S = 3.5 and standard
deviations σ P = 0.3, σ S = 0.25 .
3.5 − 2
β= = 3.84
0.25 2 + 0.32
g (x) = a0 + a1 x1 + a 2 x2 (3.16)
If normalized stochastic variables U 1 and U 2 with zero expected value and unit standard deviation
are introduced by:
X i − μ Xi
Ui = i = 1, 2 (3.17)
σ Xi
53
Note 3: First order reliability methods
g (u) = a0 + a1 ( μ X + σ X u1 ) + a 2 ( μ X + σ X u 2 )
1 1 2 2
= a0 + a1 μ X + a 2 μ X + a1σ X u1 + a 2σ X u 2
1 2 1 2
g (u) = β − α 1u1 − α 2 u 2
where:
a0 + a1μ X + a2 μ X 2
β= 1
a12σ X2 1 + a22σ X2
2
− aiσ X i
αi = i = 1, 2
a12σ X2 1 + a 22σ X2 2
Figure 3.4. Linear failure function in the x-space and in the normalized u-space.
In figure 3.4 the failure function in the x-space and in the u-space is shown. It is seen that β is the
shortest distance from origo to the failure surface in the normalized space and that the coefficients
α1 and α 2 are elements in a unit vector, α , normal to the failure surface.
54
Note 3: First order reliability methods
A first approximation to obtain an estimate of the reliability index in this case could be to linearize
the safety margin with the point corresponding to the expected values as expansion point:
M ≅ g (μ X ) + ∑
n ∂g
i =1 ∂X i
(X i − μX i
) (3.18)
X =μ X
The reliability index can then be estimated from (3.11) - (3.14). However, as noted above, the fail-
ure surface g (x) = 0 can be defined by many different but equivalent failure functions.
This implies that the reliability index based on the linearized safety margin becomes dependent on
the mathematical formulation of the safety margin. This problem is also known as the invariance
problem.
In 1974 Hasofer & Lind [3.3] proposed a definition of the reliability index which is invariant with
respect to the mathematical formulation of the safety margin.
In this section it is assumed that the stochastic variables X i , i = 1,K, n are independent. Further, it
is implicitly assumed that the variables are normally distributed. The first step in calculation of the
Hasofer & Lind reliability index β HL is to define a transformation from X to stochastic variables
U that are normalized. The normalized variables U i , i = 1,K, n with expected values 0 and stan-
dard deviation 1 are defined by:
X − μ Xi
Ui = i = 1, 2,K, n (3.19)
σ Xi
By this transformation the failure surface in the new u-space is given by, see figure 3.5:
g ( μ X + σ X u1 ,K, μ X + σ X u n ) = g u (u) = 0
1 1 n n
(3.20)
It should be noted that the u-space is rotationally symmetric with respect to the standard deviations.
55
Note 3: First order reliability methods
The Hasofer & Lind reliability index β is defined as the smallest distance from the origin O in the
u-space to the failure surface g u (u) = 0 . This is illustrated in figure 3.6. The point A on the failure
surface closest to the origin is denoted the β -point or the design point. The Hasofer & Lind reli-
ability index defined in the u-space is invariant to different equivalent formulations of the failure
function because the definition of the reliability index is related to the failure surface and not di-
rectly to the failure function. The reliability index is thus defined by the optimization problem:
n
β = min ∑ ui
2
(3.21)
g u ( u ) =0 i =1
If the failure surface is linear it is easy to see that the Hasofer & Lind reliability index is the same as
the reliability index defined by (3.14). The Hasofer & Lind reliability index can thus be considered
a generalization of the Cornell reliability index.
The numerical calculation of the reliability index β defined by (3.21) can be performed in a num-
ber of ways. (3.21) is an optimization problem with a quadratic objective function and one non-
linear constraint. A number of algorithms exist for solution of this type of problem, e.g. the NLPQL
algorithm by Schittkowski [3.4]. Here a simple iterative algorithm will be described. For simplicity
the index u will be omitted on the failure function g (u) in the following.
u ∗ = λ∇g (u ∗ ) (3.22)
u ∗ = u 0 + Δu (3.23)
56
Note 3: First order reliability methods
g (u ∗ ) ≅ g (u 0 ) + ∇g (u 0 ) T (u ∗ − u 0 ) = g (u 0 ) + ∇g (u 0 ) T Δu (3.24)
g (u ∗ ) ≅ g (u 0 ) + ∇g (u 0 ) T (u ∗ − u 0 ) ≅ g (u 0 ) + ∇g (u 0 ) T (λ∇g (u 0 ) − u 0 ) (3.25)
∇g (u 0 ) T u 0 − g (u 0 )
λ= (3.26)
∇g (u 0 ) T ∇g (u 0 )
1. Guess (u 0 )
Set i = 0
2. Calculate g (u i )
3. Calculate ∇g (u i )
If a unit normal vector α to the failure surface at the β point u ∗ is defined by:
∇g (u ∗ )
α=− (3.29)
∇g (u ∗ )
u ∗ = βα (3.30)
57
Note 3: First order reliability methods
It is noted that α is directed towards the failure set. The safety margin corresponding to the tangent
hyperplane obtained by linearizing the failure function at the β point can then be written:
M = β − αT U (3.31)
Further, using that α T α = 1 it is seen from (3.30) that the reliability index β can be written:
β = αT u∗ (3.32)
dβ
= αi (3.33)
du i u =u ∗
i.e. the components in the α vector can be considered measures of the relative importance of the
uncertainty in the corresponding stochastic variable on the reliability index. However, it should be
noted that for dependent (correlated) basic variables the components in the α -vector cannot be
linked to a specific basic variable, see the next section.
An important sensitivity measure related to α i is the socalled omission sensitivity factor ς i sug-
gested by Madsen [3.5]. This factor gives the relative importance on the reliability index by assum-
ing that stochastic variable no. i , i.e. it is considered a deterministic quantity. If variable no. i is
applied to the value ui0 , then the safety margin in the normalized space is written:
n
M i′ = β − α i ui0 − ∑α jU j (3.34)
j =1
j ≠i
β − α i ui0
β i′ = (3.35)
1 − α i2
β i′ 1 − α i ui0 / β
ςi = = (3.36)
β 1−α 2 i
58
Note 3: First order reliability methods
1
ςi = (3.37)
1 − α i2
It is seen that if α i < 0.14 , then ς i − 1 < 0.01 , i.e. the error in the reliability index is less than 1% if
a variable with α < 0.14 is fixed. The omission sensitivity factor can be generalized to non-normal
and dependent stochastic variables, see Madsen [3.5].
In this section it is assumed that the stochastic variables are normally distributed. The normalized
variables U defined by the linear transformation (3.19) are thus also normally distributed. If the
failure function in the u-space is not too non-linear, then the probability of failure Pf can be esti-
mated from:
Pf = P ( M ≤ 0) ≅ P ( β − α T U ≤ 0) = Φ (− β ) (3.38)
where Φ is the standard normal distribution function. The accuracy of (3.38) is further discussed in
section 5.
Example 3.4
1 pl 3
u max =
48 ei
where e is the modulus of elasticity and i the moment of inertia. p, l, e and i are assumed to be out-
comes of stochastic variables P, L , E and I with expected values μ and standard deviations σ .
μ[⋅] σ [⋅]
P 2 kN 0.6 kN
L 6m ~0m
E 2 ⋅10 7 kN/m 2 3 ⋅10 6 kN/m 2
I 2 ⋅10 −5 m 4 2 ⋅10 −6 m 4
59
Note 3: First order reliability methods
u max 1
≥
l 100
P−2
U1 = → P = 0.6U 1 + 2
0.6
E − 2 ⋅10 7
U2 = → E = (0.3U 2 + 2)10 7
0.3 ⋅10 7
I − 2 ⋅10 −5
U3 = → I = (0.2U 3 + 2)10 −5
0.2 ⋅10 5
∂g u
a1 = = −2160
∂u1
∂g
a 2 = u = 1440(0.2 u 3 + 2)
∂u 2
∂g
a3 = u = 960(0.3 u 2 + 2)
∂u 3
60
Note 3: First order reliability methods
Iteration 1 2 3 4 5
u1 1.00 1.29 1.90 1.91 1.90
u2 1.00 -1.89 -2.20 -2.23 -2.25
u3 1.00 -1.32 -1.21 -1.13 -1.12
β 1.73 2.64 3.15 3.15 3.15
a1 -2160 -2160 -2160 -2160
a2 3168 2500 2532 2555
a3 2208 1376 1286 1278
∑ ai
2 19.58·106 12.81·106 12.73·106 12.83·106
∑ aiui 3216 -9328 -11230 -11267
gu (u) 14928 1955 3.5 8.1
λ -0.598·10-3 -0.881·10-3 -0.882·10-3 -0.879·10-3
The reliability index is thus β = 3.15 and the corresponding α -vector is:
(
( p ∗ , e ∗ , i ∗ ) = 0.6 ⋅1.90 + 2, (−0.3 ⋅ 2.25 + 2) ⋅10 7 , (−0.2 ⋅1.12 + 2) ⋅10 −5 )
= (3.14, 1.33 ⋅10 7 , 1.78 ⋅10 −5 )
The omission sensitivity factor ς 3 corresponding to a fixed variable u3 = 0 is, see (3.37):
1
ς3 = = 1.07
1 − (−0.36) 2
* * *
Another very important sensitivity measure is the reliability elasticity coefficient defined by:
dβ p
ep = (3.39)
dp β
where p is a parameter in a distribution function (e.g. the expected value or the standard deviation)
or p is a constant in the failure function. From (3.39) it is seen that if the parameter p is changed by
1%, then the reliability index is changed by e p %. dβ dp is determined as follows:
g (u, p) = 0 (3.40)
61
Note 3: First order reliability methods
If the parameter p is given a small increment then β and the β -point change, but (3.40) still has to
be satisfied, i.e.:
n ∂g ∂ui ∂g
∑ + =0 (3.41)
i =1 ∂ui ∂p ∂p
dβ dp is determined from:
dβ d n
= ∑ ui
2
dp dp i =1
(3.42)
1 n ∂u
= ∑ ui i
β i =1 ∂p
dβ 1 n − β ∂g ∂ui
= ∑
dp β i =1 ∇g ∂ui ∂p
(3.43)
1 ∂g
=
∇g ∂p
i.e. dβ dp can be estimated on the basis of a partial differentiation of the failure function with re-
spect to the parameter p. ∇g is already determined in connection with calculation of β .
What is the reliability elasticity coefficient el for the length l in example 3.4? Using (3.43), dβ dl
is:
dβ 1 ∂g
=
dl ∇g ∂l
1
= (−200 p ⋅ l )
∑ ai
2
= −1.05
and thus:
l
el = −1.05 = −2.00
β
i.e. if the length is increased by 1%, then the reliability index decreases approximately by 2%.
62
Note 3: First order reliability methods
3.5 References
[3.1] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
[3.2] Cornell, C.A.: A Probability-Based Structural Code. ACI-Journal, Vol. 66, 1966, pp. 974-985.
[3.3] Hasofer, A.M. & N.C. Lind: An Exact and Invariant First Order Reliability Format. ASCE,
Journ. Eng. Mech. Div, 1974, pp. 111-121.
[3.4] Schittkowski, K.: NLPQL: A FORTRAN Subroutine Solving Constrained Non-Linear Pro-
gramming Problems. Annals of Operations Research, Vol. 5, 1985, pp. 485-500.
[3.5] Madsen, H.O.: Omission Sensitivity Factors. Structural Safety, Vol. 5, No. 1, 1988, pp. 33-45.
[3.6] Thoft-Christensen, P. & M.B. Baker: Structural Reliability Theory and Its Applications.
Springer Verlag, 1982.
[3.7] Ditlevsen, O. & H.O. Madsen: Bærende Konstruktioners sikkerhed. SBI-rapport 211, Statens
Byggeforskningsinstitut, 1990 (in Danish).
63
Note 3: First order reliability methods
64
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
4.1 Introduction
In note 3 it was described how a first order reliability analysis can be performed for uncorrelated
and normally distributed stochastic variables. The reliability method which is also named the "First
Order Reliability Method" (FORM) results in a reliability index β . In this note it is described how
a reliability index β can be determined when the stochastic variables are correlated and non-
normally distributed.
Let the stochastic variables X i , i = 1,K, n be normally distributed with expected values
µ X , K, µ X , standard deviations σ X1 ,K,σ X n and with correlation coefficients ρ ij , i, j = 1,K, n.
1 n
Further, let a failure function g (x) be given. In order to determine a reliability index for this failure
mode a transformation from correlated to uncorrelated stochastic variables is added to the procedure
described in section 3.4. This transformation can be performed in several ways, e.g. by determining
eigenvalues and eigenvectors, see Thoft-Christensen & Baker [4.1]. Here Choleski triangulation is
used. The procedure described in the following requires that the correlation coefficient matrix ρ is
positive definite.
The first step is to determine normalized variables Yi , i = 1,K, n with expected value 0 and standard
deviation 1:
Xi − µX
Yi = i
, i = 1,K, n (4.1)
σX i
It is easy to see that Y will have a covariance matrix (and correlation coefficient matrix) equal to
ρ.
The next step is to define a transformation from Y to uncorrelated and normalized variables U
with expected values 0 and standard deviations 1. The transformation is written:
Y = TU (4.2)
where T is a lower triangular matrix (i.e. Tij = 0 for j > i ). It is seen that the covariance matrix
C Y for Y can be written:
65
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
T11 = 1
Example 4.1
Let the three normalized stochastic variables Y = (Y1 , Y2 , Y3 ) have the correlation coefficient matrix:
1 0 .5 0 .2
ρ = 0.5 1 0.4
0.2 0.4 1
1 0 0
T = 0.5 0.87 0
0.2 0.34 0.92
Y1 = U 1
Y2 = 0.5U 1 + 0.87U 2
Y3 = 0.2U 1 + 0.34U 2 + 0.92U 3
* * *
X = µ X + DTU (4.5)
where D is a diagonal matrix with standard deviations in the diagonal. Using (4.5) the failure func-
tion can be written g (x) = g (µ X + DTu) and a reliability index β can be determined as shown in
section 3.4.
66
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
Example 4.2
A failure mode is modelled by a failure function with three normally distributed variables
X1, X 2 , X 3 :
g (x) = x1 − x2 x32
where µ X = 25.0 , σ X = 0.25 , µ X = 4.0 , σ X = 0.2 , µ X = 2.0 and σ X = 0.1. The variables are
1 1 2 2 3 3
correlated as the variables in example 4.1. The standardized normalized and uncorrelated u-
variables are obtained from example 4.1 and (4.5) as:
X 1 = µ X + σ X U1
1 1
X 2 = µ X + σ X (0.5U 1 + 0.87 U 2 )
2 2
g (u) = 25.0 + 0.25u1 − (4.0 + 0.2(0.5u1 + 0.87u 2 ) )(2.0 + 0.1(0.2u1 + 0.34u 2 + 0.92u 3 ) )
2
The failure function can be used to find β as explained in section 3.4 by the iteration scheme used
in example 3.4.
The solution is β = 3.86 ( Pf = 5.67 ⋅ 10 −5 ) , u ∗ = {1.051, 2.426, 2.812} and α = {0.27, 0.63, 0.73}.
* * *
Generally the stochastic variables are not normally distributed. In order to determine a measure of
the reliability of a component (failure mode) with non-normally distributed variables it is natural, as
for normally distributed variables, to establish a transformation to standardized (uncorrelated and
normalized) normally distributed variables and to determine a Hasofer & Lind reliability index β .
A simple transformation from X i to U i can be defined by the identity:
Φ (U i ) = FX i ( X i ) (4.6)
x1 = FX−1 (Φ (u1 ) )
1
M (4.7)
xn = F −1
Xn (Φ(u n ) )
67
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
(
g ( x1 , K , x n ) = g FX−1 (Φ (u1 ) ), K , FX−1 (Φ (u n ) ) = 0
1 n
) (4.8)
In the algorithm for determination of β (see section 3.4) the gradient of the failure function with
respect to ui is needed. From (4.8):
∂g ∂g ∂xi ∂g ϕ Φ FX ( xi )
−1
( ( ))
= = i
(4.9)
∂u i ∂xi ∂u i ∂xi f X ( xi ) i
ln x − µ L
FX ( x) = Φ (4.10)
σL
where:
σ 2
σ L = ln 2 + 1 and µ L = ln µ − 12 σ L2
µ
x = exp(σ L u + µ L ) (4.11)
* * *
Example 4.4 Gumbel Variable
For a Gumbel distributed variable X with expected value µ and standard deviation σ the distribu-
tion function is:
where:
π 0.5772
a= and b=µ−
6σ a
1
x =b− ln[− ln Φ (u )] (4.13)
a
68
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
* * *
(
u1 = Φ −1 FX ( x1 ) 1
)
M (4.14)
(
u n = Φ −1 FX ( x n ) n
)
When the transformation defined above is applied in connection with the β -algorithm in section
3.4 it is also known under the name of principle of normal tail approximation. In the normal tail
approximation a normal distribution with parameters µ i′ and σ i′ is determined for each non-normal
stochastic variable such that the distribution function values and the density function values are the
same at a point xi′ :
x ′ − µ i′
Φ i = FX i ( xi′ ) (4.15)
σ ′
i
1 xi′ − µ i′
ϕ = f X i ( xi′ ) (4.16)
σ i′ σ i′
σ i′ =
(
ϕ Φ −1 (FX i ( xi′ ) ) ) (4.17)
f X i ( xi′ )
xi − µ i′
ui = (4.19)
σ i′
69
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
∂g ∂g (x) ∂xi
=
∂ui ∂xi ∂ui
∂g (x)
= σ i′ (4.21)
∂xi
=
( (
∂g (x) ϕ Φ FX ( xi′ )
−1
i
))
∂xi f X ( xi′)
i
At the β -point, u ∗ , and the corresponding point x ∗ in the x-space the gradient estimated by (4.9) is
equal to the gradient estimated by (4.21) if xi′ = xi∗ , i = 1,2,K, n . This indicates that if the current
guess of the β -point in the algorithm u i is used as u´ in (4.17) - (4.21) and if the points u1 , u 2 ,K
converge to u ∗ then the transformation defined by (4.7) is equivalent to the transformation defined
by the normal tail approximation, see Ditlevsen [4.2] for further details.
Example 4.5
Consider the safety margin:
M = g ( X) = X 1 − 2 X 22
where:
X 1 : is log-normally distributed with expected value µ1 = 10 and standard deviation σ 1 = 3 (or
LN(10.0, 3.0)). From (4.10) ( µ L , σ L ) = (2.26, 0.294) is obtained.
X 2 : is Gumbel distributed with expected value µ1 = 1 and standard deviation σ 1 = 0.1 (or
EX1(1.0, 0.1)). From (4.12) (a, b) = (12.8, 0.955) is obtained.
The transformation from the physical x-space to the standard normal u-space is found from (4.11)
and (4.13):
2
1
g (u) = exp(σ L u1 + µ L ) − 2 b − ln[− ln Φ (u2 )]
a
By application of the β -iteration scheme explained in section 3.4 β can be found as β = 4.040
and u ∗ = {– 2.587, 3.103}, α = {– 0.640, 0.768}.
* * *
In this section two techniques are described, which can be used to determine a reliability index
when the stochastic variables are dependent and non-normally distributed, namely methods based
on the Rosenblatt transformation, see [4.3] and the Nataf transformation, see [4.4]. It should be
noted that if all the stochastic variables are normally and log-normally distributed then the tech-
nique described in section 4.2 can be used because the log-normal variables can easily be trans-
formed to normal variables, see example 4.6.
70
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
Example 4.6
Consider 3 stochastic variables X i , i = 1, 2, 3 with expected values µ[⋅] , standard deviations σ [⋅]
and coefficients of variation V [⋅] as shown in this table:
1 sym.
ρ = ρ X X
2 1
1
ρ X X ρX 1
3 1 3X2
X 1 is assumed to be normally distributed, but X 2 and X 3 are log-normally distributed. Two new
variables are defined by Yi = ln X i , i = 2, 3. They become normally distributed. The expected
values and standard deviations of the normally distributed variables X 1 , Y2 and Y3 become, see
example 4.3,
µ[⋅] σ [⋅]
X1 µ X1 σ X1
Y2 µY2 = ln µ X 2 − 12 σ Y22 σ Y = ln(V X2 + 1)
2 2
Y3 µY3 = ln µ X 3 − 1σ 2
2 Y3 σ Y3 = ln(V X2 + 1)
3
The new correlation matrix ρ' of correlation coefficients between X 1 , Y2 and Y3 can be obtained
from the definition of the covariance between two stochastic variables:
1 sym.
ρ X X VX
ρ' = 2 1 2
1
σY 2
ρ X X VX ln(1 + ρ X 2X3
VX VX )
3 1 3 2 3
1
σ Y 3
σY σY 2 3
Example 4.7
Consider a normally distributed variable X 1 and two log-normally distributed variables X 2 and
X 3 with the statistic parameters:
71
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
1 (sym.)
ρ = 0 .2 1
0.5 0.3 1
From example 4.6 the following parameters are obtained for X 1 , Y2 = ln X 2 and Y3 = ln X 3
µ[⋅] σ [⋅]
X 1 10.0 2.0
Y2 1.50 0.472
Y3 1.94 0.05
and:
1 (sym.)
ρ' = 0.21 1
0.50 0.37 1
It is seen that the absolute value of the correlation coefficients become higher (which will always be
the case). Furthermore, it is seen from the example and the expressions in the ρ' -matrix that the
difference between ρ ij′ and ρ ij vanishes for small coefficients of variation V, which is also the rea-
son why the difference between ρ ij′ and ρ ij is sometimes neglected.
From this example it is concluded that a failure function of normally and log-normally distributed
stochastic variables can be transformed to a failure function of normally distributed variables. The
failure function in the u-space can then be obtained from ρ' and the transformation explained in
section 5.2. Next the reliability index β can be obtained as usual.
* * *
For dependent stochastic variables X i , i = 1,K, n the Rosenblatt transformation, see [4.3], can be
used to define a transformation to the u-space of uncorrelated and normalized normally distributed
variables U i , i = 1,K, n . The transformation is defined as, see also (4.7):
x1 = FX−1 (Φ (u1 ) )
1
x2 = F −1
X 2 | X1 (Φ(u2 ) | X 1 = x1 )
(4.22)
M
xn = FX−1| X L X
n 1 n −1
(Φ(un ) | X 1 = x1 ,K, X n−1 = xn−1 )
72
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
where FX | X L X
i 1 i −1
(xi | X 1 = x1 ,K, X i −1 = xi −1 ) is the distribution function of Xi given
X 1 = x1 , K, X i −1 = xi −1 :
xi
∫ f X LX i −1 X i
( x1 , L, xi −1 , t )dt
(xi | X 1 = x1 ,K, X i −1 = xi −1 ) =
1
−∞
FX | X L X (4.23)
i 1 i −1
f X L X ( x1 ,K, xi −1 )
1 i −1
f X LX ( x1 ,K, xi ) is the joint density function of X 1 ,K, X i . The transformation starts for given
1 i
u1 ,K, u n by determination of x1. Next x2 is calculated using the value of x1 determined in the
first step. x3 ,K, xn are then calculated in the same stepwise manner.
(
u1 = Φ −1 FX ( x1 ) ) 1
u2 = Φ −1
(F (x X 2 | X1 2 )
| X 1 = x1 )
(4.24)
M
un = Φ −1 FX ( n | X 1L X n−1
(xn | X 1 = x1,K, X n−1 = xn−1 ))
The Rosenblatt transformation is very useful when the stochastic model for a failure mode is given
in terms of conditional distributions. For example, this is often the case when statistic uncertainty is
included. Examples 4.8 and 4.9 show how the Rosenblatt transformation can be used.
1 2πk
4kb4 kγ H S2π 3
4
Sηη (ω ) = 5 exp − b γ a
(a)
ω (k pTZ ) 4
π ωk pTZ
where γ = 3 , k b = 1.4085 , k p = 0.327 exp(−0.315γ ) + 1.17 and k γ = 1 − 0.285 ln(γ ) . H S is the sig-
nificant wave height and TZ is the zero crossing period. The superscript a is:
(k p TZ 2ωπ − 1) 2
a = exp −
2σ 2
a
where:
73
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
2π
0.07 for ω < k T
p Z
σa =
2π
0.09 for ω ≥ k T
p Z
The distribution function of the maximum wave elevation H m within a given time period [0, T] can
be estimated from, see Davenport [4.6]:
1 hm 2
FH (hm ) = exp − ν 0T exp − 2
(b)
m
σ
where :
m2
ν0 = (c)
m0
σ = m0 (d)
1 ∞
mi = ∫
(2π ) i 0
ω i Sηη (ω )dω (e)
H S and TZ are usually modelled as stochastic variables. Here H S is modelled by a Rayleigh dis-
tribution with the parameter s:
1 h 2
FH (h) = 1 − exp − , h ≥ 0 (f)
2 s
S
t γ
(h)
FT |H S (t | H S = h) = 1 − exp − (g)
Z
k ( h)
where:
74
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
The distribution function for H m given H S and TZ is given by (b). The distribution function for
TZ given H S is given by (g), and the distribution function for H S is given by (f). (j) can then be
estimated by FORM using the failure function:
g = hm − H m ( H S ,T Z ) (k)
Φ (U1 ) = FH ( H S )
S
Φ (U 2 ) = FT Z |H S (TZ | H S ) (l)
Φ (U 3 ) = FH m |H S ,T Z ( H m | H S , TZ )
The reliability index β for (k) is determined by the algorithm in section 3.4 and:
P ( H m > hm ) ≅ Φ (− β ) (m)
Example 4.9
Consider a failure function with two stochastic variables X 1 and X 2 : (Madsen et al. [4.7], p. 77)
g (x) = 18 − 3 x1 − 2 x2 (a)
75
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
∂ 2 FX X ( x1 , x2 )
f X X ( x1 , x2 ) = 1 2
= ( x1 + x2 + x1 x2 ) exp[−( x1 + x2 + x1 x2 )] , x1 > 0, x2 > 0 (c)
1 2
∂x1∂x2
Realisations u1 and u 2 of standard normal variables U 1 and U 2 are obtained from the Rosenblatt
transformation as:
(
u1 = Φ −1 FX ( x1 ) 1
)
(d)
u2 = Φ −1
(F X 2 | X1 ( x2 | x1 ) )
where:
∞
f X ( x1 ) = ∫ f X X ( x1 , x 2 )dx 2 = exp(− x1 ) , x1 > 0
1 1 2
(e)
0
x1
FX ( x1 ) = ∫ f X ( x1 )dx1 = 1 − exp(− x1 ) , x1 > 0
1 1
(f)
0
FX 2 |X1
( x 2 | x1 ) = 0
(g)
f X ( x1 )
1
For the transformation from the x-space to the u-space the formulas become:
(h)
x2 = F −1
X 2 | X1 (Φ(u2 ) | X 1 = x1 )
1 − (1 + x 2 ) exp[−( x 2 + x1 x 2 )] = Φ (u 2 ) (i)
76
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
The β -optimization problem includes a local and a global minimum. The β -point (which is also
the global minimum) is u1∗ = {2.78, 0.1} with β1 = 2.78 and Pf ≈ 2.68 ⋅10 −3 . Further, the local mi-
nimum point u ∗2 = {− 1.30, 3.25} is identified with β 2 = 3.50 .
* * *
An alternative way to define the transformation from the u-space to the x-space is to use the Nataf
transformation, see [4.4] and [4.8]. This transformation is in general only an approximate transfor-
mation. The basic idea is to establish the marginal transformations defined in section 4.3 (as if the
stochastic variables were independent) and to use a correlation coefficient matrix ρ e in an y-space,
which is obtained from the correlation coefficient matrix ρ in the x-space by multiplying each cor-
relation coefficient by a factor F, which depends on distribution types and the statistical parameters.
To describe the Nataf transformation it is thus sufficient to consider two stochastic variables X i
and X j .
X i = FX−1 (Φ (Yi ) )
i
(4.25)
Xj =F −1
Xj (Φ(Y )) j
The stochastic variables Yi and Y j have an (equivalent) correlation coefficient ρ ije , which in the
Nataf transformation is determined such that dependence between X i and X j is approximated as
well as possible.
Xk − µX
Zk = k
k = i, j (4.26)
σX k
77
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
FX−1 (Φ ( y k ) ) − µ X
zk = k k
k = i, j (4.27)
σX k
Further, from (4.2) it is seen that uncorrelated variables U i and U j can be introduced by:
yi = ui
(4.28)
y j = ρ u + 1 − (ρ ) u j
e
ij i
e 2
ij
ρ ij can then be related to the (unknown) equivalent correlation coefficient ρ ije by:
∞ ∞
ρij = ∫ ∫ zi z jϕ 2 ( yi , y j , ρije )dyi dy j
−∞ −∞
FX−1 (Φ ( yi ) ) − µ X FX (Φ ( y j ) ) − µ X
−1
∞ ∞
= ∫ ∫ i i j j
ϕ 2 ( yi , y j , ρije )dyi dy j (4.29)
−∞ −∞ σX σX
(( ))
i j
where ϕ 2 ( ) is the two-dimensional normal density function. From (4.29) ρ ije can be determined by
iteration.
Based on ρ ije the following approximate joint density function f Xei X j ( xi , x j ) is obtained:
f X i ( xi ) f X j ( x j )
f Xei X j ( xi , x j ) = ϕ 2 ( yi , y j , ρ ije ) (4.30)
ϕ ( yi )ϕ ( y j )
where y i = Φ −1 FX ( xi ) . ( i
)
(4.29) has been solved for ρ ije by der Kiureghian & Liu [4.8] for a number of distribution functions
and approximations for the factor:
ρ ije
F= (4.31)
ρ ij
has been obtained. With ρ = ρ ij and Vi = σ X i µ X i examples of approximations for F are shown in
table 4.1.
78
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
e
For n = 2 it should be checked that ρ12 ≤ 1 . For n > 2 the corresponding requirement is that ρ e is
positive definite. In der Kiureghian & Liu [4.8] or Ditlevsen & Madsen [4.9], approximations for F
are also shown for Gamma, Frechet, Uniform, Rayleigh and Gumbel distributions.
Xi Xj F
normal log-normal Vj ln(1 + V j2 )
Example 4.10
Consider the same problem as in example 4.9 but use the Nataf transformation instead of the Ro-
senblatt transformation. The correlation coefficient between X 1 and X 2 is:
ρ=
1
σX σX
[E[ X X 1 2 ] − µX µX 1 2
]
1 2
1 ∞ ∞
= 0∫ 0∫ x1 x2 f X ( x1 x2 )dx1dx2 − µ X µ X
σX σX
1 2
1X 2
1 2
= ∫0 ∫0 x1 x2 ( x1 + x2 + x1 x2 ) exp[− ( x1 + x2 + x1 x2 )]dx1dx2 − 1
∞ ∞
= −0.40366
where:
∞ ∞
µ X = µ X = E[ x1 ] = ∫ x1 f X ( x1 )dx1 = ∫ x1 exp(− x1 )dx1 = 1
2 1 1
0 0
∞ ∞
σ X2 = σ X2 = E[ x12 ] − µ X2 = ∫ x12 f X ( x1 )dx1 − µ X2 = ∫ x12 exp(− x1 )dx1 − 1 = 1
2 1 1 1 1
0 0
79
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
ρ e = Fρ = −0.566
The transformation form (u1 , u 2 ) to ( x1 , x2 ) is given by (4.25) and (4.2) (or (4.28) for two stochas-
tic variables)
x1 = − ln[1 − Φ (u1 )]
[ (
x 2 = − ln 1 − Φ ρ e u1 + 1 − ( ρ e ) 2 u 2 )]
Using the failure function in example 4.9 the two β -points are determined as:
* * *
As described in note 3 three important sensitivity measures can be used to characterize the sensitiv-
ity of the reliability index with respect to parameters and the stochastic variables, namely:
α -vector
The elements in the α -vector characterize the importance of the stochastic variables. From the lin-
earized safety margin M = β − α T U it is seen that the variance of M is:
σ M2 = α12 + α 22 + L + α n2 = 1 (4.32)
For independent stochastic variables α i2 thus gives the percentage of the total uncertainty associ-
ated with U i (and X i ). If for example X 2 , X 3 and X 4 are dependent, then α 22 + α 32 + α 42 gives the
percentage of the total uncertainty which can be associated with X 2 , X 3 and X 4 altogether.
1 ∂g p
ep = (4.33)
∇g ∂p β
80
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
For parameters p in the distribution function for X , which is related to standardized variables U by
U = T(X, p), e p is obtained as:
1 ∂T(x ∗ , p ) p
∗ T
e p = (u ) (4.34)
β ∂p β
where u ∗ and x ∗ are the β -points in the u-space and the x-space.
1
ξi = (4.35)
1 − α i2
gives a measure of the change in the reliability index if stochastic variable no. i is fixed. This sto-
chastic variable is assumed to be independent of the other stochastic variables. As described in
Madsen [4.10], the omission sensitivity factor can be generalized to dependent stochastic variables.
4.6 References
[4.1] Thoft-Christensen, P. & M.B. Baker: Structural Reliability Theory and Its Applications.
Springer Verlag, 1982.
[4.2] Ditlevsen, O.: Principle of Normal Tail Approximation. ASCE, Journal of Engineering
Mechanics Division, Vol. 107, 1981, pp. 1191-1208.
[4.3] Rosenblatt, M.: Remarks on a Multivariate Transformation. The Annals of Mathematical
Statistics, Vol. 23, 1952, pp. 470-472.
[4.4] Nataf, A.: Determination des distribution dont les marges sont données. Comptes redus de
l'Academie des Sciences, Vol. 225, 1962, pp. 42-43.
[4.5] Barltrop, N.D.P. & A.J. Adams: Dynamics of fixed marine structures. Butterworth -
Heinemann, London, 1991.
[4.6] Davenport, A.G.: Note on the Distribution of the Largest Value of a Random Function
with Application to Gust Loading. Proc. Institution of Civil Engineers, London, Vol. 28,
1964, pp. 187-196.
[4.7] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
[4.8] Der Kiureghian, A. & P.-L. Liu : Structural Reliability Under Incomplete Probability In-
formation. ASCE, Journal of Engineering Mechanics, Vol. 112, No. 1, 1986, pp. 85-104.
[4.9] Ditlevsen, O. & H.O. Madsen: Bærende Konstruktioners sikkerhed. SBI-rapport 211, Sta-
tens Byggeforskningsinstitut, 1990 (in Danish).
[4.10] Madsen, H.O.: Omission Sensitivity Factors. Structural Safety, Vol. 5, No. 1, 1988, pp. 33-
45.
81
Note 4: First order reliability analysis with correlated and non-normal stochastic variables
82
Note 5: SORM and simulation techniques
First Order Reliability Methods can be expected to give reasonable results when the failure func-
tions are not too non-linear. FORM techniques are described in notes 3 and 4. If the failure func-
tions in the standardized u-space are rather non-linear then Second Order Reliability Methods
(SORM) techniques, where a second order approximation of the failure function is established, can
be used. These techniques are described in section 5.1.
Other techniques, which can be used for such types of problems, are simulation techniques. Simula-
tion methods, which are described in sections 5.2 - 5.7, can also be powerful when the failure func-
tions in the u-space have more than one β -point, i.e. there are several local, probable failure re-
gions.
In simulation methods realisations (outcomes) x̂ of the stochastic variables X are generated for
each sample. When simulation methods are used to estimate Pf the failure function is calculated for
each realisation x̂ and if the realisation is in the failure region, then a contribution to the probability
of failure is obtained. In section 5.2 different techniques to generate realisations of stochastic vari-
ables are described. In the literature a large number of simulation methods are described. Sections
5.3 to 5.7 contain a description of some of the most important methods, namely:
Finally in section 5.8 it is described how sensitivity measures can be obtained by simulation.
Compared with a FORM estimate of the reliability of a component (or failure mode) an improved
estimate can be obtained by using a second order approximation of the failure surface at the β -
point u∗ in the u-space:
where D the Hessian matrix of second order partial derivatives of the failure surface at the β -
point:
83
Note 5: SORM and simulation techniques
∂2 g
Dij = , i, j = 1,2, K , n (5.2)
∂ui ∂u j ∗
u =u
In the following it is described how a second order reliability index can be determined.
The β -point and the gradient vector can be written, see (3.29) and (3.30):
u ∗ = βα ∇g (u ∗ ) = − ∇g (u ∗ ) α (5.3)
y = Ru (5.4)
1 ~
β − yn + ∗
y T RDR T ~
y=0 (5.6)
2 ∇g (u )
where ~
y = ( y1 , y2 ,K, yn −1 , yn − β )T .
The solution of (5.6) with respect to y n using up to second order terms in y1 , y 2 ,K, y n−1 gives the
hyperparabolic surface:
Aij =
1
∗
{RDR } T
ij i, j = 1,2,K, n − 1 (5.8)
2 ∇g (u )
y´= H v (5.9)
where the columns in H are the eigenvectors of A . (5.7) can then be written:
n −1
y n = β + ∑ λi vi2 (5.10)
i =1
84
Note 5: SORM and simulation techniques
where λi , i = 1,2,K, n − 1 are the eigenvectors in A . The eigenvectors and eigenvalues can e.g. be
found by Jacobi-iteration or subspace-iteration for large problems, where only the largest eigenval-
ues are important, see e.g. [5.11].
The probability of failure Pf estimated using the second-order approximation of the failure surface
is:
∞ ∞ ∞
PfSO = ∫ L ∫ ϕ (v1 )Lϕ (vn−1 ) ∫ ϕ ( yn )dyn dv1 L dvn−1 (5.11)
−∞ −∞ β + Σλi vi2
The approximation is illustrated in figure 5.1, which also shows the first-order approximation (see
(3.38)) to the exact probability of failure Pf = P( g (U) ≤ 0).
PfFO = Φ (− β ) (5.12)
Figure 5.1. Illustration of first and second order approximations of the failure surface.
It should be noted that due to the rotational symmetry of the normal density function the points in
the area close to the β -point (which is the point closest to origo) has the largest probability density.
Therefore, the largest contributions to the probability of failure come from this area. Further, it is
noted that the u-dimensional normal density function for uncorrelated variables ϕ n ( ) ∝ exp(− r 2 2)
decreases fast with the distance r from origo. If the failure surface is rather non-linear then a second
order approximation of the failure surface can be expected to give a much better estimate of the
probability of failure than the first-order approximation. Finally it should be noted that for β → ∞
the first (and second) order estimates of the probability converge to the exact result: PfFO → Pf .
n −1 − 12
PfSO ≈ Φ (− β ) ∏ (1 + 2 βλ j ) (5.13)
j =1
Improvements to (5.13) have been suggested by for example Tvedt [5.2] and [5.3].
85
Note 5: SORM and simulation techniques
β SO = −Φ −1 ( PfSO ) (5.14)
The approximation in (5.13) - (5.14) assumes that the matrix I + 2 βA is positive definite.
A necessary tool in simulation techniques for estimation of the probability of failure is to simulate
outcomes of stochastic variables with an arbitrary distribution. For this a method to generate uni-
formly distributed numbers is first described. Next it is shown how the inverse method can be used
to generate outcomes of stochastic variables with a general distribution. Finally methods to generate
outcomes of normally distributed variables are described.
In this subsection a stochastic variable V which is uniformly distributed between 0 and 1 is consid-
ered. The distribution function is:
v if 0 ≤ v ≤1
FV (v) = (5.15)
0 else
The most widely used techniques to simulate (generate) pseudo-random numbers of V is the multi-
plicative congruential generators, see Hammersley & Handscomb [5.4] and the XOR generator, see
Ditlevsen & Madsen [5.5]. In multiplicative congruential generators the pseudo-random numbers
are determined sequentially by:
(
vi = 89069vi −1 + 1 modulo 232 ) (5.17)
86
Note 5: SORM and simulation techniques
The numbers generated by (5.16) are not completely independent. It can be shown that the correla-
tion between successive numbers lies in the interval, see Hammersley & Handscomb [5.4]:
1 6c c a 1 6c c a
a − am (1 − m ) − m , a − am (1 − m ) + m (5.18)
Numerical investigations have shown that if the multiplicative congruential generator is used to
generate outcomes of stochastic vectors then the generated vectors are not uniformly distributed in
the n -dimensional space. An algorithm which generates numbers much more random in the n-
dimensional space is the so-called XOR random number generator, see Ditlevsen & Madsen [5.5].
The method is illustrated in figure 5.2. It is seen that the distribution function for X̂ with outcomes
generated by this procedure is:
( )
FXˆ ( x) = P( Xˆ ≤ x) = P FX−1 (V ) ≤ x = P(V ≤ FX ( x) ) = FX ( x) (5.20)
Example 5.1
Let X be exponential distributed with distribution function:
FX ( x) = 1 − exp(−λx)
87
Note 5: SORM and simulation techniques
1
xˆ = − ln(1 − vˆ)
λ
* * *
U 1 = − 2 ln V1 cos(2πV2 )
(5.21)
U 2 = − 2 ln V1 sin( 2πV2 )
where V1 and V2 are independent stochastic variables uniformly distributed between 0 and 1.
V1 + V2 + L + Vn − a → U for n → ∞ (5.22)
where V1 ,V2 ,K are independent equidistributed random variables uniformly distributed between 0
1
and 1 (expected value µV = 1
2
and variance σ V2 = ∫ ( x − 12 ) 2 dx = 121 ).
0
A reasonable choice is a = n2 and n = 12. Then U becomes approximately normal with expected
value 0 and standard deviation 1.
88
Note 5: SORM and simulation techniques
X = µ X + DTU (5.23)
where the elements in U are uncorrelated with zero mean and unit standard deviation. Using the
techniques described above to generate outcomes of normally distributed variables and (5.23) reali-
sations of X can be generated.
* * *
In the following sections different simulation methods to estimate the probability of failure are de-
scribed:
Pf = P( g (U) ) ≤ 0 (5.24)
1 N
Pˆ f = ∑ I [ g (uˆ j )] (5.25)
N j =1
where N is the number of simulations and û j is sample no. j of a standard normally distributed
stochastic vector U The indicator function I [ g (u)] is defined by:
Pˆ f (1 − Pˆ f )
s= (5.26)
N
Confidence intervals for the estimate of the probability of failure can be determined using that P̂f
becomes normally distributed for N → ∞.
The idea in importance sampling is to concentrate the sampling in the area of the total sample space
which has the largest contribution to the probability of failure. In this way the standard error of the
estimate of Pf can be reduced significantly. Pf is written:
f U (y )
Pf = ∫ L ∫ I [ g (u)] f U (u)du1 L du n = ∫ L ∫ I [ g (y )] f S (y )dy1 L dyn (5.27)
f S (y )
89
Note 5: SORM and simulation techniques
where f S (y) is the sampling density function and f U (y ) = ϕ ( y1 )Lϕ ( y n ) is the standard normal
density function for U .
In theory, if the sampling density f S is chosen to be proportional to f U in the failure region then
the standard error on Pf would be zero. Unfortunately, this choice is not possible because Pf is not
known beforehand. In the following it is shown how f S can be chosen reasonably.
1 N f U (yˆ j )
Pˆ f = ∑ I [ g (yˆ j )] (5.28)
N j =1 f S (yˆ j )
where f S (y) is the sampling density function from which the sample vectors ŷ j are generated.
2 2
1 N f U (yˆ j )
1 N f U (yˆ j )
s= ∑ I [ g (yˆ j )] − ∑ I [ g (yˆ j )] (5.29)
N ( N − 1) j =1 f S (yˆ j ) N j =1 f S (yˆ j )
g ( x1 , x2 ) = x2 − x1
where f X (x) is the joint density function of the stochastic variables modelling the load and the
strength.
In importance sampling the simulations are concentrated in the area which contributes most to the
probability of failuire. Pf is estimated by (5.28):
1 N f X (yˆ j )
Pf = ∑ I [ g (yˆ j )]
N j =1 f Y (yˆ j )
90
Note 5: SORM and simulation techniques
where f Y (y) is the sampling density function from which the sample vector ŷ j is generated. Fig-
ure 5.3 shows the general difference between crude Monte Carlo simulation and importance sam-
pling.
Example 5.3
Consider the example from Madsen et al. [5.6], where the cross-section of a reinforced concrete
beam is analysed. n = 7 stochastic variables are used. The failure function is written
x5 x32 x42
g ( x ) = x 2 x3 x 4 − − x1
x 6 x7
91
Note 5: SORM and simulation techniques
VARIABLE DIST. µ V
x1 bending moment N 0.01 MNm 0.3
x2 eff. depth of reinforcement N 0.3 m 0.05
x3 yield stress of reinforcement N 360 MPa 0.1
x4 area of reinforcement N 226·10 −6 m 2 0.05
x5 factor N 0.5 0.1
x6 width of beam N 0.12 m 0.05
x7 compressive strength of concrete N 40 MPa 0.15
Table 5.1. Statistical data. µ is the expected value and v = σ / µ is the coefficient of variation. N
indicates normal (Gauss) distribution.
The statistical data are shown in table 5.1. The stochastic variables are assumed to be independent.
A transformation to normalized stochastic variables (with expected value 0 and standard deviation
1) U i , i = 1,2,K,7 is established by:
X i = σ iU i + µ i , i = 1,K,7
(σ 5 u5 + µ 5 )(σ 3u3 + µ 3 ) 2 (σ 4 u 4 + µ 4 ) 2
g (u) = (σ 2 u 2 + µ 2 )(σ 3u3 + µ 3 )(σ 4 u 4 + µ 4 ) − − (σ 1u1 + µ1 )
(σ 6 u 6 + µ 6 )(σ 7 u 7 + µ 7 )
In importance sampling Pf is estimated by (5.28) with yˆ = uˆ + u ∗ and f S (yˆ ) = f U (yˆ − u ∗ ), i.e. the
samples are concentrated around the point u ∗ . û is a sample generated from the standard normal
vector U . In this example u ∗ is chosen as (see the next section)
u ∗ = (2.5, − 1, − 2, − 1, 0, 0, 0)
92
Note 5: SORM and simulation techniques
If the β -point has been determined before simulation techniques are used importance sampling can
be very effective with the β -point as the point around which the samplings are concentrated, see
figure 5.4. Such a technique is described in this section. The sampling density function f S in (5.28)
is the normal density of uncorrelated variables with expected values u ∗i , i = 1,2,K, n and common
standard deviations, σ .
Pf is estimated by:
1 N f U (σuˆ j + u ∗ ) n
Pf = ∑ I [ g (σuˆ j + u )]
∗
σ (5.30)
N j =1 f U (uˆ j )
where f U (u) is the standardized normal density function and û j is a sample generated from stan-
dardized normal variables.
Figure 5.5. Different standard deviations of the sampling density, σ 1 < σ 2 < σ 3 .
The standard error is estimated by (5.29). The efficiency of the importance sampling can be ex-
pected to be dependent on the choice of standard deviation of the sampling density, see figure 5.5.
It should be noted that if a failure mode has multiple β -points importance sampling based on only
one β -point is not efficient. In this case more general methods have to be used, see section 5.7.
93
Note 5: SORM and simulation techniques
In this technique the space is separated into two disjoint regions D1 and D2 , see figure 5.6. It is
assumed that D1 is selected such that no failure occurs in this region. Here D1 is chosen as the re-
gion inside a sphere with radius β . The probability of being in D1 is:
p1 = P ∑ U 12 ≤ β 2 = χ 2 (n, β 2 )
n
(5.31)
i =1
1 − p1 N
Pˆ f = ∑ I [ g (uˆ j )] (5.32)
N j =1
where û j is sample no. j from D2 (simulated from a standard normally distributed stochastic vec-
tor U = (U 1 ,K,U n ) , but only those samples outside D1 are used).
The standard error is:
Pˆ f (1 − Pˆ f )
s = (1 − p1 ) (5.33)
N
The standard error is thus reduced by a factor (1 − p1 ) when compared with crude Monte Carlo
simulation. Usually this is a significant reduction. However, it should be taken into account that it is
more difficult to generate the samples to be used. If the samples are generated by taking the samples
from simulation of normal distributed variables with û > β then on average 1−1p samples should
1
94
Note 5: SORM and simulation techniques
be generated before one sample is outside the β -sphere. So only in cases where the failure function
require much more computational work than the generation of the samples û it can be expected that
this technique is efficient.
Example 5.4
Consider an example where the failure surface in standardized coordinates can be written:
Pf = Φ (−3.305) = 0.000228
The simulation results are shown in table 5.3 with standard errors in ( ). It is seen that importance
sampling and Crude Monte Carlo simulation by excluding the β -sphere are much better than crude
Monte Carlo simulation. Further it is seen that in this example σ = 1 is the best choice for impor-
tance sampling.
In this section some other simulation methods are described, namely directional sampling, Latin
hypercube simulation and adaptive simulation techniques.
95
Note 5: SORM and simulation techniques
Directional simulation
Instead of formulating the reliability problem in rectangular coordinates it is possible to formulate it
in polar coordinates. Directional simulation methods are based on such a formulation and were first
suggested by Deak [5.7] in connection with evaluation of the multinormal distribution function.
U = RA (5.34)
where the radial distance R > 0 is a stochastic variable and A is a unit vector of independent sto-
chastic variables, indicating the direction in the u-space.
It is now assumed that the origin u = 0 is in the safe area (g (0) > 0 ) and that the failure region de-
fined by {u : g (u) ≤ 0} is star shaped with respect to the u = 0 , i.e. every half-line starting form
u = 0 only crosses the failure surface once.
96
Note 5: SORM and simulation techniques
(
P (g ( RA) ≤ 0 A = a ) = ∫r (a ) f R ( s A = a)ds = 1 − χ n2 r (a) 2
∞
) (5.36)
where χ n2 ( ) is the χ n2 distribution with n degrees of freedom. r (a) is the distance from the origin
u = 0 to the failure surface, i.e. g (r (a)a ) = 0 in the a direction.
1 N 1 N
(
Pˆ f ≈ E[ Pˆ f ] = ∑ pˆ j = ∑ 1 − χ n2 (r (aˆ j ) 2 )
N j =1 N j =1
) (5.37)
The basic idea in this method is to assure that the entire range of each variable is sampled, in order
to obtain an efficient estimate of the probability of failure. The range of each variable is divided
into m intervals. The probability of an outcome in each interval should be equal.
In the simulation procedure the samples are generated in such a way that an interval of each vari-
able will be matched just one time with an interval from each of the rest of the variables. In figure
5.8 the Latin hypercube method is illustrated by an example with n = 2 stochastic variables and m =
7 intervals.
97
Note 5: SORM and simulation techniques
1. For each variable generate one point from each of the intervals. uˆ ij , j = 1,2,K, m thus repre-
sents the m points for variable i.
2. The first point uˆ 1j in the Latin hypercube sample is generated by sampling one value uˆij1 from
each axis ui . The second point is generated in the same way, except that the values uˆij1 are de-
leted from the sample. In this way m points are generated.
1 m
Pˆ f = ∑ I [ g (uˆ j )]
m j =1
There is no simple form for the standard error of this simulation method but in general the standard
1 times the standard error of crude Monte Carlo simulation.
error is of the magnitude mN
The initial sampling density is suggested to be standard normal with standard deviation 1 but with
the expected value point moved to a point û ( 0) in or close to the failure region. This can be difficult,
but based on the initial knowledge of which variables represents load variables and which variables
represents strength variables such a point can be selected (for strength variables û (0) should be
negative and for load variables û (0) should be positive). The initial density is used until a sample
point is generated in the failure domain.
When multiple points in the failure region are generated the sampling density is modified such that
the regions around the points with the largest probability density are emphasized. The simplest ap-
proach is to locate the expected value point at the point in the failure region with the largest prob-
ability density.
Another approach is to use a so-called multi-modal sampling density which generates samples
around a number of points in the failure region, but emphasizes the region around a point in propor-
tion to the probability density at the point. This allows us to emphasize more than one point and is
98
Note 5: SORM and simulation techniques
closer to the ideal sampling density (which is proportional to the probability density at each point in
{ }
the failure domain). Let uˆ (1) , uˆ ( 2 ) ,K, uˆ ( k ) be the set of points in the failure region, which are used
to construct the multi-modal sampling density. The corresponding multi-modal density is:
k
hUk (u) = ∑ w j f U( j ) (u) (5.38)
j =1
where f U( j ) (u) is the density function of a normally distributed stochastic vector with uncorrelated
variables, standard deviations 1 and expected value point equal to uˆ ( j ) . The weights are determined
by:
f U (uˆ ( j ) )
wj = (i )
(5.39)
∑ik=1 f U (uˆ )
An estimate of the probability of failure can now be obtained on the basis of N simulations where
the importance sampling technique is used :
ˆ 1 N f U (uˆ ( j ) )
Pf = ∑ j ( j ) I g (uˆ ( j ) )
N j =1 hU (uˆ )
[ ] (5.40)
In many cases it is very interesting to know how sensitive an estimated probability of failure is with
respect to a change of a parameter p. p is here assumed to be the expected value or the standard de-
99
Note 5: SORM and simulation techniques
viation of a stochastic variable. The transformation from the basic stochastic variables X to stan-
dardized normal variables is written:
X = T(U, p) (5.41)
1 N
[
Pˆ f = ∑ I g (T(uˆ j , p ))
N j =1
] (5.43)
∂Pf
By direct differentiation the gradient ∂p
of Pf with respect to p can be estimated by introducing a
small change ∆p in p and calculating:
∂Pf ∆Pˆ f
∂p
≈
∆p
=
1 1 N
[( 1 N
)]
∑ I g T(uˆ , p + ∆p) − ∑ I g T(uˆ , p)
∆p N j =1
j
N j =1
j
[( )] (5.44)
The two terms in (5.44) are estimated separately. This estimate of ∆Pf can be expected to be both
inaccurate because it is the difference between two "uncertain" estimates and time consuming be-
cause two sets of samples has to be generated.
∂Pf
Alternatively, ∂p
can be written:
∂Pf ∂ ∂
= ∫ I [g (T(u, p ) )]f U (u)du = ∫ I [g (x)]f X ( p ) (x)dx
∂p ∂p ∂p
∂f ( x)
= ∫ I [g (x)] X ( p ) dx (5.45)
∂p
∂f ( x) 1
= ∫ I [g (x)] X ( p ) f X ( p ) (x)dx
∂p f X ( p ) ( x)
where f X ( p ) (x) is the density function of X with the parameter p. Corresponding to (5.43) and
(5.45) the following estimates can be obtained by simulation:
1 N
[
Pˆ f = ∑ I g (xˆ j )
N j =1
] (5.46)
∂f Xˆ ( p ) (xˆ j )
1 N
Pˆ f = ∑ I g (xˆ )
N j =1
j
[ ∂p
] 1
f X ( p ) (xˆ j )
(5.47)
The samples x̂ j are generated from the density function f X ( p ) (x) using for example the inverse
100
Note 5: SORM and simulation techniques
simulation method. The advantage of this formulation is that the same samples can be used to esti-
∂Pˆ f
mate both P̂f and ∂p
. This increases the accuracy and reduces the computational effort compared
with direct differentiation. Similar formulations can be derived for other simulation types.
5.9 References
[5.1] Breitung, K.: Asymptotic approximations for multinormal integrals. Journal of the Engi-
neering Mechanics Division, ASCE, Vol. 110, 1984, pp. 357-366.
[5.2] Tvedt, L.: Two second order approximations to the failure probability. Veritas report
DIV/20-004-83, Det norske Veritas, Norway, 1983.
[5.3] Tvedt, L.: Second order reliability by an exact integral. Lecture notes in engineering, Vol.
48, Springer Verlag, 1989, pp.377-384.
[5.4] Hammersley, J.M. & D.C. Handscomb: Monte Carlo methods. John Wiley & sons, New
York, 1964.
[5.5] Ditlevsen, O. & H.O. Madsen: Bærende Konstruktioners sikkerhed. SBI-rapport 211, Sta-
tens Byggeforskningsinstitut, 1990 (in Danish).
[5.6] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
[5.7] Deak, I.: Three digit accurate multiple normal probabilities. Numerical Mathematik, Vol.
35, 1980.
[5.8] Melchers, R.: Simulation in time-invariant and time-variant reliability problems. Proceed-
ings, IFIP WG 7.5, Munich, September 1991, pp. 39-82.
[5.9] McKay, M.D., Beckman, R.J. & W.J. Conover: A comparison of three methods for select-
ing values of input variables in the analysis of output from a computer code. Technomet-
rics, Vol. 21, No. 2, 1979.
[5.10] Karamchandani, A.: New methods in systems reliability. Department of Civil Engineering,
Stanford University, Report No. RMS-7, Ph.D. Thesis, May 1990.
[5.11] Bathe, K.-J. Finite element procedures in engineering analysis. Prentice-Hall, 1982.
101
Note 5: SORM and simulation techniques
102
Note 6: Reliability evaluation of series systems
6.1 Introduction
So far, in the previous notes, only reliabilities of individual failure modes or limit states have been
considered. In this note it is described how the individual limit states interact on each other and how
the overall systems reliability can be estimated when the individual failure modes are combined in a
series system of failure elements.
In section 6.2 a series system is defined, followed by section 6.3 where it is explained how the
FORM-approximation of the reliability of a series system is obtained and how the correlation be-
tween failure elements are interpreted. In section 6.4 it is described how the multi-dimensional nor-
mal distribution function needed for the series system reliability estimation can be evaluated using
bounds and approximations. Finally, section 6.5 introduces sensitivity analysis of series systems.
A failure element or component, see figure 6.1, can be interpreted as a model of a specific failure
mode at a specific location in the structure.
The combination of failure elements in a series system can be understood from the statically deter-
minate (non-redundant) truss-structure in figure 6.2 with n structural elements (trusses). Each of the
n structural elements is assigned 2 failure elements. One with a failure function modelling material
yielding failure and one with a failure function modelling buckling failure.
For such a statically determinate structure it is clear that the whole structural system fails as soon as
any structural element fails, i.e. the structure has no load-carrying capacity after failure of one of
the structural elements. This is called a weakest link system and is modelled as a series system. The
series system which then becomes the systems reliability model consists of 2n failure elements
shown in figure 6.3.
103
Note 6: Reliability evaluation of series systems
Figure 6.3. Weakest link system modelled as a series system of failure elements.
It is in this connection important to notice the difference between structural components and failure
elements and the difference between a structural system and a systems reliability model.
If failure of one failure element is defined as systems failure the reliability of the series system can
be interpreted as the reliability of failure. That also includes the case of statically indeterminate
structures where failure of more than one failure element cannot be accepted.
Consider a structural system where the system reliability model is a series system of m failure ele-
ments. Each of the failure elements is modelled with a safety margin:
M i = g i ( X) , 1,2,K, m (6.1)
The transformation between the standard normal stochastic U -variables and the stochastic variables
X can be obtained as explained in note 4 and is symbolically written as X = T(U) . Furthermore, it
is known from notes 3 and 4 that the FORM probability of failure for failure element i can be writ-
ten:
Pf = P ( M i ≤ 0) = P (g i ( X) ≤ 0 ) = P( g i (T(U ) ) ≤ 0 )
(6.2)
i
≈ P ( β i − α Ti U ≤ 0) = Φ (− β i )
The series system fails if just one of the elements fails, i.e. the probability of failure of the series
system is:
Thus, if all the failure functions as in (6.2) are linearized at their respective β -points the FORM
approximation of PfS of a series system can be written:
{
PfS ≈ P U − α Ti U ≤ − β i }
m
(6.4)
i =1
i =1
{
} m
i =1
{ }
PfS ≈ 1 − P I − α Ti U > − β i = 1 − P I α Ti U < − β i = 1 − Φ m (β; ρ)
m
(6.5)
104
Note 6: Reliability evaluation of series systems
where Φ m is the m-dimensional normal distribution function (see the following section 6.4). It has
been used that the correlation coefficient ρ ij between two linearized safety margins
M i = β i − α Ti U and M j = β j − α Tj U is:
ρ ij = α Ti α j (6.6)
From (6.5) a formal or so-called generalized series systems reliability index β S can be introduced
from:
or:
In figure 6.4 the exact failure domain, which is the union of the individual element failure domains
is hatched. Furthermore, the reliability indices β i , i = 1,2,3 and the safety margins linearized at their
corresponding β -points u ∗i , i = 1,2,3 are shown. It is seen that (6.7) or (6.8) is an approximation
when the failure functions are non-linear in the u-space.
* * *
105
Note 6: Reliability evaluation of series systems
cosθ ij = α Ti α j = ρ ij
where θ ij is the angle between the α-vectors α i and α j or simply between the linearized safety
margins. I.e., the correlation coefficients ρ ij can be comprehended as a measure of the angle be-
tween the linearized safety margins and hereby as a measure of the extent of the failure domain.
* * *
106
Note 6: Reliability evaluation of series systems
107
Note 6: Reliability evaluation of series systems
* * *
From the previous section it is obtained that if β i and ρ ij , i, j = 1,2,L, m are known the problem is
to evaluate the m-dimensional normal distribution function Φ m (β; ρ) in (6.8) for the FORM ap-
proximation of β S .
β1 β 2 βm
Φ m (β; ρ) = ∫−∞ ∫−∞ L ∫−∞ ϕ m (x;ρ)dx1 dx 2 K dx m (6.9)
1 1
ϕ m (x;ρ) = 12
exp(− x T ρ −1 x) (6.10)
(2π ) m 2 ρ 2
The multi-dimensional integral in (6.9) can only in special cases be solved analytically and will for
even small dimensions, say five, be too costly to evaluate by numerical integration. Instead, so-
called bounds methods are used for hand calculations and so-called asymptotic approximate meth-
ods are used for computational calculations.
In the following, so-called simple bounds and Ditlevsen bounds will be introduced as bounds for
the reliability of series systems.
Simple Bounds
Simple bounds can be introduced as:
m m
max P( M i ≤ 0) ≤ PfS ≤ ∑ (P( M i ≤ 0) ) (6.11)
i =1 i =1
108
Note 6: Reliability evaluation of series systems
where the lower bound corresponds to the exact value of PfS if all the elements in the series system
are fully correlated.
m m
∑
− Φ −1 Φ (− β i ) ≤ β S ≤ min β i
i =1
(6.12)
i =1
When the failure of one failure element is not dominating in relation to the other failure elements
the simple bounds are generally too wide and therefore often of minor interest for practical use.
Ditlevsen Bounds
Much better bounds are obtained from the second-order bounds called Ditlevsen bounds [6.4]. The
derivation of the Ditlevsen bounds can be seen in [6.1], [6.4], [6.6], [6.7] or [6.8]. The bounds are:
m i −1
PfS ≥ P( M 1 ≤ 0) + ∑ max P( M i ≤ 0) −
∑ P(M i ≤ 0 I M j ≤ 0),0 (6.13a)
i =2 j =1
{P(M i I M j ≤ 0)}
m m
PfS ≤ ∑ P( M i ≤ 0) − ∑ max
j <i
≤0 (6.13b)
i =1 i =2
m i −1
Φ (− β S ) ≥ Φ (− β1 ) + ∑ maxΦ(− β i ) − ∑ Φ 2 (− β i ,− β j ; ρ ij ),0 (6.14a)
i =2 j =1
Φ (− β S ) ≤ ∑ Φ (− β i ) − ∑ max{Φ 2 (− β i ,− β j ; ρ ij )}
m m
(6.14b)
i =1 i =2 j <i
The numbering of the failure elements influences the bounds. However, experience suggests that it
is a good choice to arrange the failure elements according to decreasing probability of failure, i.e.
P( M 1 ≤ 0) ≥ P( M 2 ≤ 0) ≥ L ≥ P( M m ≤ 0) . The Ditlevsen bounds are usually much more precise than
the simple bounds in (6.11) - (6.12), but require the estimation of Φ 2 (− β i ,− β j ; ρ ij ) in (6.14).
∂ 2 Φ 2 ( β i , β j ; ρ ij ) ∂Φ 2 ( β i , β j ; ρ ij )
= (6.15)
∂β i ∂β j ∂ρ ij
Therefore:
109
Note 6: Reliability evaluation of series systems
ρ ij ∂Φ 2 ( β i , β j ; t )
∫0
Φ 2 ( β i , β j ; ρ ij ) = Φ 2 ( β i , β j ;0) +
∂t
t = z dz
(6.16)
ρ
= Φ ( β i )Φ ( β j ) + ∫ ϕ 2 ( β i , β j ; z )dz
ij
From figure 6.8 it is seen that P( M i ≤ 0I M j ≤ 0) equals the probability contents in the hatched
angle BAE. P is greater than the probability content in the angle BAD and in the angle CAE. How-
ever, P is less than the sum of the probability contents in the angles BAD and CAE. This observa-
tion makes it possible to derive simple bounds for P ij = Φ 2 (− β i ,− β j ; ρ ij ) .
The probability contents pi and p j in the angles CAE and BAD, respectively, are:
β i − ρ ij β j β j − ρ ij β i
γi = γj= (6.18)
1 − ρ ij2 1 − ρ ij2
110
Note 6: Reliability evaluation of series systems
max( pi , p j ) ≤ Φ 2 (− β i ,− β j ; ρ ij ) ≤ pi + p j (6.19)
0 ≤ Φ 2 (− β i ,− β j ; ρ ij ) ≤ min( pi , p j ) (6.20)
These bounds are easy to use and Pij can be approximated as the average of the lower and the upper
bounds. If the gap between the lower and the upper bounds is too wide, a more accurate method,
such as numerical integration of (6.16) should be used.
PfS ≥ P( D1 ) + P( D2 ) − P( D2 I D1 ) + P( D3 ) − P( D3 I D1 ) − P( D3 I D2 )
from which it is seen that the hatched domain in figure 6.9 is the difference between the lower Dit-
levsen bound and the exact PfS .
PfS ≤ P( D1 ) + P( D2 ) + P ( D3 ) − P( D2 I D1 ) − P( D3 I D1 )
From which it is seen that the dotted domain in figure 6.9 is the difference between the upper Dit-
levsen bound and the exact PfS .
* * *
111
Note 6: Reliability evaluation of series systems
g1 (u) = exp(u1 ) − u 2 + 3
g 2 (u) = u1 − u 2 + 5
g 3 (u) = exp(u1 + 4) − u 2
g 4 (u) = 0.1u12 − u 2 + 4
The reliability indices β i with the corresponding Pf , α -vectors and β -points are shown in table
i
6.1.
112
Note 6: Reliability evaluation of series systems
From table 6.1 the correlation matrix ρ can be obtained from (6.6):
1.000 sym.
0.878 1.000
ρ=
0.712 0.961 1.000
0.962 0.714 0.492 1.000
Simple Bounds
From (6.12) the simple bounds of β S can be obtained as:
Ditlevsen Bounds
For Ditlevsen bounds it is necessary to evaluate Φ 2 (− β i ,− β j ; ρ ij ) , i, j = 1,…,4, for j < i , which
can be done approximately by (6.17) - (6.20). In the following matrix γ i and γ j from (6.18) are
shown. ( γ i from (6.18) is shown in the lower triangle and γ j is shown in the upper triangle)
From (6.17) -(6.20) it is then possible to obtain the following table with bounds of
Φ 2 (− β i ,− β j ; ρ ij ) .
113
Note 6: Reliability evaluation of series systems
{
Φ(− β S ) ≥ 2.276 ⋅10−4 + max 2.035 ⋅10−4 − 7.95 ⋅10 −5 , 0 }
{
+ max 5.738 ⋅10 −5 − (1.40 + 4.71) ⋅10−5 , 0 }
+ max { 3.174 ⋅10 −5
− (3.09 + 0.776 + 0.0927) ⋅10 −5 , 0 }
= 3.52 ⋅10− 4
corresponding to:
3.36 ≤ β S ≤ 3.39
If instead the average approximations of Φ 2 (− β i ,− β j ; ρ ij ) in the bottom row of table 6.2 are
used only approximations of the bounds are obtained (i.e, there is no guarantee that β S is within
the bounds):
3.36 ≤ β S ≤ 3.37
If Φ 2 (− β i ,− β j ; ρ ij ) is calculated exactly from (6.16) the following exact bounds are obtained:
3.381 ≤ β S ≤ 3.383
It is seen that the Ditlevsen bounds in this case are narrow. This will often be the case.
* * *
114
Note 6: Reliability evaluation of series systems
Consider again example 4.9 where the failure function in the u-space was found as shown in figure
6.11.
Instead of estimating the probability of failure as Pf1 = Φ(− β1 ) = Φ(−2.78) = 2.68 ⋅10 −3 , the prob-
ability of failure is estimated as Pf = P( M 1 ≤ 0 U M 2 ≤ 0) where M 1 and M 2 are safety margins
from linearization at the β -points u1∗ and u ∗2 , respectively, (see figure 6.11). The safety margins
are written M 1 = β1 − α 1T U and M 2 = β 2 − α T2 U .
With β 1 = 2.784 , β 2 = 3.501 ( Pf = 2.31 ⋅10 −4 ) and the α -vectors α 1 = (0.999, 0.036)
2
and
α 2 = ( −0.370, 0.929) , the correlation coefficient is ρ12 = α α 2 = −0.337 . The probability of failure
T
1
is then obtained as Pf = 1 − Φ 2 ( β1 , β 2 ; ρ12 ) , which is:
Φ 2 (− β 1 ,− β 2 ; ρ12 ) is estimated from (6.17) -(6.20). From (6.18) it can be obtained that γ 1 = 4.2101
and γ 2 = 4.715 , which by use of (6.17) results in p1 = 3.25 ⋅ 10 −9 and p 2 = 2.960 ⋅ 10 −9 . An average
estimate from (6.20) is then obtained as Φ 2 (− β 1 ,− β 2 ; ρ12 ) = 1.48 ⋅ 10 −9 . Pf then is Pf = 2.68 ⋅10 −3 +
2.32 ⋅ 10 −4 − 1.48 ⋅ 10 −9 = 2.91 ⋅ 10 −3 , which corresponds to β S = 2.758 . Compared to the exact result
β S = 2.755 obtained by numerical integration with formula (c) in example 4.9 inserted into (3.6)
this is a satisfactory estimate.
* * *
115
Note 6: Reliability evaluation of series systems
m β − ρt
Φ m (β; ρ) = ∫−∞∞ ϕ (t )∏ Φ i dt
(6.21)
i =1
1 − ρ
m β − ρt
PfS = 1 − ∫−∞∞ ϕ (t )∏ Φ i dt
(6.22)
i =1
1 − ρ
when the correlation coefficients are not all equal an approximation of the probability of failure can
be obtained by using an average correlation coefficient ρ as ρ in (6.22). ρ is determined from:
m i −1
r 2
ρ=
m(m − 1)
∑∑ ρ ij (6.23)
i =1 j =1
The approximation based on the average correlation coefficient can be considered as the first term
in a Taylor expansion of PfS at the average correlation coefficient point with respect to the correla-
tion coefficients.
Using (6.22) with ρ = ρ , an approximation of PfS is obtained. The approximation will in many
cases be conservative.
Example 6.7
Consider the series system of example 6.5 again. The average correlation coefficient becomes:
1
ρ = (0.878 + 0.712 + 0.961 + 0.962 + 0.714 + 0.492) = 0.786
6
with β = (3.51, 3.54, 3.85, 4.00) in (6.22) the average correlation coefficient approximation be-
comes PfS = 4.28 ⋅10 −4 corresponding to β S = 3.33 , which from example 6.5 is seen to give a
conservative estimate of the series system reliability.
* * *
116
Note 6: Reliability evaluation of series systems
are used. Two of these methods are the Hohenbichler approximation, see [6.5], and the approxima-
tion by Gollwitzer and Rackwitz [6.3]. These methods are in general very precise and make it pos-
sible to calculate Φ m within reasonable computer time.
From (6.8) it can be shown that the sensitivity of β S with respect to a model parameter p can be
found as:
dβ S 1 m ∂Φ m (β; ρ) dβ i i −1 ∂Φ m (β; ρ) dρ ij dρ ij
= ∑ + 2 ∑ (6.24)
dp ϕ ( β ) i =1 ∂β i
S
dp j =1 ∂ρ ij dp
However, to get an estimate of the sensitivity of a systems reliability index β S it is often sufficient
to use:
dβ S 1 m ∂Φ m (β; ρ) dβ i
≈ ∑ (6.25)
dp ϕ ( β S ) i =1 ∂β i dp
where dβ i dp can be obtained as already described in note 4 and ∂Φ m (β, ρ) ∂β i can be determined
either numerically by finite differences or by the semi-analytical methods described in [6.9] where
also details of sensitivity analysis can be found.
6.6 References
[6.1] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
[6.2] Hohenbichler, M. & R. Rackwitz: First-Order Concepts in Systems Reliability. Structural
Safety, Vol. 1, No. 3, pp. 177-188, 1983.
[6.3] Gollwitzer, S. & R. Rackwitz: Comparison of Numerical Schemes for the Multinormal
Integral. Springer Verlag. In Proceedings of the first IFIP WG 6.5 Working Conference on
Reliability and Optimization of Structural Systems, P. Thoft-Christensen (ed.) pp. 157-174,
1986.
[6.4] Ditlevsen, O.: Narrow Reliability Bounds for Structural Systems. Journal of Structural Me-
chanics, Vol. 7, No. 4, pp. 453-472. 1979.
[6.5] Hohenbichler, M.: An Asymptotic Formula for the Probability of Intersections. Berichte
zur Zuverlassigkeitstheorie der Bauwerke, Heft 69, LKI, Technische Universität München,
pp. 21-48, 1984.
[6.6] Thoft-Christensen, P. & M.J. Baker: Structural Reliability Theory and Its Applications.
Springer Verlag, 1982.
[6.7] Thoft-Christensen, P. & Y. Murotsu: Application of Structural Systems Reliability Theory.
Springer Verlag, 1986.
[6.8] Ditlevsen, O. & H.O. Madsen: Bærende Konstruktioners sikkerhed. SBI-rapport 211, Sta-
tens Byggeforskningsinstitut, 1990 (in Danish).
117
Note 6: Reliability evaluation of series systems
118
Note 7: Reliability evaluation of parallel system
7.1 Introduction
In this note it is described how the reliability of a system can be evaluated when more than one fail-
ure element have to fail before the whole system is defined to be in a state of failure. This is per-
formed by introduction of parallel systems in section 7.2, followed by sections 7.3 and 7.4 where
the FORM approximation of the reliability of a parallel system and reliability evaluation techniques
are introduced, respectively. In section 7.5 it is described how the parallel systems are combined
into a systems reliability model of a series system of parallel systems and finally, in section 7.6, it is
shown, how the corresponding reliability evaluations can then be performed.
The introduction and the necessity of parallel systems for the reliability modelling of some struc-
tural systems can be illustrated by considering the statically indeterminate (redundant) truss-
structure in figure 7.1 with N structural elements (trusses). Two failure elements are assigned to
each of the N structural elements, one with a failure function modelling material yielding failure
and one with a failure function modelling buckling failure.
For such a statically indeterminate (redundant) structure it is clear that the whole structural system
will not always fail as soon as one of structural element fails, because the structure has a load-
carrying capacity after failure of some of the structural elements. This load-carrying capacity is ob-
tained after a redistribution of the load effects in the structure after the element failure. Failure of
the entire redundant structure will then often require failure of more than one structural element. (It
is in this connection very important to define exactly what is understood by failure of the structural
system). Clearly the number of systems failure modes in a redundant structure is generally high.
Each of these system failure modes can be modelled by a parallel system consisting of generally n
elements, where n is the number of failure elements which have to fail in the specific systems fail-
ure mode before the entire structure is defined to be in a state of failure. The parallel system with n
elements is shown in figure 7.2.
119
Note 7: Reliability evaluation of parallel system
Since a redistribution of the load effects has to take place in a redundant structural system after fail-
ure of one or more of the structural elements it becomes very important in parallel systems to de-
scribe the behaviour of the failed structural elements after failure has taken place. If the structural
element has no strength after failure the element is said to be perfectly brittle. If the element after
failure has a load-bearing capacity equal to the load at failure, the element is said to be perfectly
ductile.
In figure 7.3 a perfectly brittle and a perfectly ductile element are shown with an example of the
behaviours and the symbols used for perfectly brittle and perfectly ductile elements, respectively.
Figure 7.3. Perfectly brittle and perfectly ductile elements with symbols.
Clearly all kinds of structural components and material behaviours cannot be described as perfectly
brittle or perfectly ductile. All kinds of combinations in between exist, i.e. some, but not all, of the
failure strength capacity is retained. One of these modellings are the elastic-residual model shown
in figure 7.4.
Before the reliability modelling in a parallel system of failure elements can be performed the struc-
tural behaviour of the considered failure mode must be clarified. More specifically the failure of the
120
Note 7: Reliability evaluation of parallel system
structural elements and consequences with determination of residual load-carrying capacity and
load redistribution in each step in the structural element failure sequence must be described. Then
the failure functions of the failure elements in the parallel system can be formulated. Failure func-
tion no. 1 models failure in parallel system element no. 1 without failure in any other elements.
Failure function no. 2 models failure in parallel system element no. 2 with failure in the structural
element corresponding to failure element no. 1 (i.e. after redistribution of loads). Failure function
no. 3 then models failure of parallel system element no. 3 with failure in the structural elements
corresponding to failure element nos. 2 and 1, etc. etc.
The obtained failure functions can then be used in the reliability evaluations of the parallel system
without further consideration of the structural system and structural behaviour.
Consider a fibre bundle with n perfectly ductile fibres modelled by a parallel system. The strength
Ri , i = 1,2,K, n of the individual fibres is identically normal distributed N( µ , σ ) with a common
correlation coefficient ρ . The fibre bundle is loaded by a deterministic load S = nS e , where S e is
the constant load on each fibre. The reliability indices of the fibre are the same for all
fibres and equal to:
µ − Se
β=
σ
The strength R of the ductile fibre bundle is obtained as the sum of the individual fibre strengths,
i.e. R is normally distributed with:
µ R = nµ and σ R2 = nσ 2 + n(n − 1) ρσ 2
The reliability index of the parallel system (fibre bundle) then is:
µR − S nµ − n( µ − βσ ) n
βP = = =β
σR nσ 2 + n(n − 1)σ 2 ρ 1 + ρ (n − 1)
It is also possible to obtain β P of a ductile fibre bundle when the fibres are not correlated by a
common correlation coefficient ρ . This can e.g. be done by use of the average correlation coeffi-
cient defined in (6.23) and used in the above expression, see [7.4].
Another case of a fibre bundle is the Daniels system [7.7] of n perfectly brittle fibres. The strengths
of the n fibres are r1 , r2 ,K, rn , where r1 ≤ r2 ≤ L ≤ rn . The strength of the fibre bundle then is:
rs = max{nr1 , (n − 1)r2 ,K,2rn −1 , rn }
Now, let ri , i = 1, 2,K, n be realizations of independent random variables Ri with identical distri-
bution functions. Similarly, rs is the realization of Rs . Daniels showed that Rs is normally distrib-
uted N( µ R , σ R ) for n → ∞ , where:
s s
121
Note 7: Reliability evaluation of parallel system
µ R = nr0 [1 − FR (r0 )]
s
and σ R2 = nr02 FR (r0 )[1 − FR (r0 )]
s
where r0 is the maximum point of the function r[1 − FR (r )] . The result is valid under the condition
that r0 is unique and r[1 − FR (r )] = 0 for r → ∞ .
For a closer description also for small values of n, see [7.8 p. 249].
* * *
After the failure functions of the failure elements in a parallel system have been formulated it is
possible to estimate the reliability by FORM from the following description.
Consider a parallel system of n failure elements each modelled with a failure function and a safety
margin:
M i = g i ( X), i = 1, 2,K, n (7.1)
The transformation between the standard normal U -variables and the stochastic variables X can be
obtained as explained in note 4 and is symbolically written as X = T(U).
The parallel system fails if all of the elements fail, i.e. the probability of failure of the parallel sys-
tem is defined as the intersection of the individual failure events:
Then a so-called joint β -point is introduced as the point in the failure domain (defined from (7.2))
closest to the origin, see figure 7.5. The n A out of the n failure functions which equal zero at u∗ are
then linearized at u∗ :
M i = β iJ − α Ti U, i = 1, 2,K, n A (7.3)
where:
− ∇ u g i (T(u ∗ ))
αi = and β iJ = α Ti u∗ (7.4)
∇ u g i (T(u∗ ))
122
Note 7: Reliability evaluation of parallel system
use of the joint β -point and not the individual β -points as in calculation of an element reliability
index β .
n
{ } = P I {− α }
A nA
PfP ≈ P I β iJ − α Ti U ≤ 0 T
i U ≤ − β iJ = Φ n (−β J ; ρ) (7.5)
i =1 i =1 A
where Φ nA is the n A -dimensional normal distribution function and the correlation coefficient ρ ij
between two linearized safety margins M i = β iJ − α Ti U and M j = β jJ − α Tj U is:
ρ ij = α Ti α j (7.6)
From (7.5) a formal generalized parallel systems reliability β P can be introduced by:
PfP = Φ n A (−β J ; ρ) = Φ (− β P ) (7.7)
as:
β P = −Φ −1 ( PfP ) = −Φ −1 (Φ n (−β J ; ρ) )
A
(7.8)
The joint β -point is from its definition determined as the solution of the following optimization
problem:
min γ = 12 u T u
u (7.9)
s.t. g i (u) ≤ 0 , i = 1, 2, K , n
The solution of the joint β -point problem can be obtained by a general non-linear optimization
algorithm as NLPQL [7.1] or the problem specific algorithm JOINT3 described in [7.2].
123
Note 7: Reliability evaluation of parallel system
In figure 7.6 the exact failure domain as the intersection of the individual element failure domains is
hatched. Furthermore, the n A = 2 active safety margins linearized at the joint β -point u∗ are
shown.
It is seen that (7.7) or (7.8) is an approximation when the failure functions are non-linear in the u-
space or if so-called secondary joint β -points exist (a secondary β -point is shown in figure 7.6 as
u 2 ). For high reliability levels the approximation in (7.8) including the n A active constraints of
(7.9) is often sufficiently accurate.
* * *
The formulation in (7.9) requires that at least one of the failure functions is greater than zero in the
origin. If this is not the case the problem can be converted to a series system problem by writing the
safe domain as a union. For further explanation and inclusion of the secondary joint β -points for a
more precise estimation, see [7.3].
In some references a cruder and older formulation of the FORM parallel system reliability is util-
ized. The failure domain is estimated as the intersection of the linearized failure functions at the
individual β -points, i.e. only the individual β -point optimization problems are solved and not the
joint β -point problem in (7.9).
124
Note 7: Reliability evaluation of parallel system
125
Note 7: Reliability evaluation of parallel system
From figure 7.8 it is seen that 3.0 ≤ β P ≤ ∞ corresponding to the fully positive correlated and the
fully negative correlated cases, respectively.
* * *
The result from the previous section is that if β iJ and ρ ij , 1, 2,K , n A are known the problem is to
evaluate the n A -dimensional normal distribution function Φ n (−β J ; ρ) in (7.8) for the FORM ap-
A
P
proximation of β . As described in note 6, this can generally not be performed by numerical inte-
gration within a reasonable computing time for higher dimensions. Instead bounds or approximate
methods are used.
In the following, simple bounds and a second order bound will be introduced as bounds for the reli-
ability of parallel systems.
Simple Bounds
If only the active constraints of (7.9) are assumed to influence the reliability of the parallel system
the simple bounds can be introduced as:
( )
nA
0 ≤ PfP ≤ min P( M iJ ≤ 0) (7.10)
i =1
where M iJ , i = 1, K , n A are the linearized safety margins at the joint β -point. The upper bound
corresponds to the exact value of PfP if all the n A elements are fully correlated with ρ ij = 1 .
If all correlation coefficients ρ ij between the n A elements are higher than zero, the following sim-
ple bounds are obtained:
126
Note 7: Reliability evaluation of parallel system
nA nA
∏ P ( M i ≤ 0) ≤ Pf ≤ min P ( M i ≤ 0)
J P J
(7.12)
i =1 i =1
max β iJ ≤ β P ≤ −Φ −1 ∏ Φ (− β iJ )
n A n A
(7.13)
i =1 I =1
The simple bounds will in most cases be so wide that they are of little practical use.
In (7.15) it is seen that the probability of failure of a parallel system of two elements
Φ 2 (− β iJ ,− β jJ , ρ ij ) is necessary. These probabilities are the same as the probabilities used in the
Ditlevsen bounds for series systems, see note 6. In note 6 both a method by numerical integration
(6.16) and a bounds method (6.17) - (6.20) are described. Hereby the tools for evaluation of the
bounds are described.
More refined and complicated bounds can also be developed, see [7.4], but will not be shown here.
g1 (u) = exp u1 − u 2 + 1
g 2 (u) = u1 − u 2 + 1
g 4 (u) = 0.1u12 − u 2 + 2
127
Note 7: Reliability evaluation of parallel system
Figure 7.9. Four failure functions for a parallel system and joint β -point, u∗ .
It is seen directly from figure 7.9 that n A = 2 and the joint β -point is the intersection between g 3
and g 4 . The joint β -point can be found to be u ∗ = (−1.23; 2.16) . The α -vectors are found from
(7.4) as α 3 = (−0.908; 0.420) and α 4 = (0.233; 0.971) , i.e the correlation coefficient from (7.6) is
ρ 34 = 0.18 . From (7.3) β J = (2.02; 1.81) .
The second order lower bound will in this two-dimensional case be exact if Φ 2 (− β 3J , − β 4J ; ρ 34 ) s
evaluated to be exact. The result is:
β P = 2.92
If instead the bounds technique from note 6, (6.17)-(6.20) is used, the bounds are obtained as 2.84
≤ β P ≤ 3.04 or by taking the average of the bounds in (6.19) β P = 2.92 .
* * *
128
Note 7: Reliability evaluation of parallel system
It is clear that a real redundant structural system generally has many failure modes, i.e. different
sequences of element failure. Each sequence can then be modelled by a parallel system. If one of
these parallel systems fails then the whole system fails, i.e. the overall systems reliability model is a
series system of the failure modes or parallel systems. This is schematically shown in figure 7.10.
It is also possible to formulate the systems reliability model as a parallel system of series systems,
see [7.5].
It is seen in figure 7.11 that the truss structure becomes statically determinate if any of the elements
1,2,3,4,5 or 6 is removed (fails). It is furthermore seen that the structure fails if any pair of the ele-
ments 1,2,3,4,5 and 6 fails. The structure also fails if one of the elements 7,8,9 or 10 fails. The sys-
tems reliability model is then a series system with 19 elements where 15 of the elements are parallel
systems each with two failure elements. The elements in the series system are: {1,2}, {1,3}, {1,4},
{1,5}, {1,6}, {2,3}, {2,4}, {2,5}, {2,6}, {3,4}, {3,5}, {3,6}, {4,5}, {4,6}, {5,8}, {7}, {8}, {9} and
{10}.
* * *
129
Note 7: Reliability evaluation of parallel system
The FORM estimate of the generalized systems reliability index β S is written as in note 6, see (6.1)
- (6.8):
β S = −Φ −1 (1 − Φ n (β P ; ρ P ) )P
(7.17)
where β P is an n P -vector of generalized reliability indices for the individual parallel systems calcu-
lated as in (7.8) and ρ P is a matrix of the corresponding approximate correlation coefficients
between the parallel systems.
For approximation of the coefficients in the correlation matrix ρ P each of the parallel systems is
approximated by a failure element with a linear safety margin, see [7.6]:
M P = β iP − α iP
i
( ) T
U , i = 1, 2, K , n (7.18)
where the vectors α iP , i = 1, 2,K, n are determined such that the sensitivity of β P with respect to
changes in the joint β -point: ∇ u ∗ β P are equivalent when obtained from (7.18) (formulated as
(β ) u ) and when obtained from (7.8). Furthermore, a normalization is performed for calculation
i
P T ∗
of correlations:
a iP
α iP = , i = 1, 2, K, n (7.19)
a iP
n ∂Φ n (−β J , ρ i )
T i
i dα k i
i
1 Ai ∗
aij = ∑ α kj + ∗ u
P
(7.20)
Ai
ϕ (− β iP ) k =1 i
du j ∂β kJ
In (7.20) the influence on β P in (7.18) of the correlations ρ i are neglected. n Ai is the number of
active constraints in the ith parallel system. dα k du ∗j is obtained from differentiation of (7.4):
dα k − I ∂∇ g
T
∇ u g ik ∇ u g ik
= + u ik
(7.21)
du ∗j ∇ u g ik ∇ u g ik
3 ∂u ∗j
The elements in the matrix of correlation coefficients between the parallel systems are then calcu-
lated from:
= (α mP ) α nP
T
ρ mn
P
(7.22)
130
Note 7: Reliability evaluation of parallel system
Now β S can be estimated from (7.17). For further explanations and details of reliability estimation
of series systems of parallel systems, see [7.6].
The sensitivities for evaluation of the obtained systems reliability indices in (7.17) or (7.9) can in
principle be obtained as explained in section 7.5. The sensitivity evaluation of a generalized reli-
ability index of series system of parallel systems or of a parallel system, however, requires much
more numerical effort and several perturbation analyses of optimality conditions of the included
optimization problems, see [7.6].
7.8 References
[7.2] Enevoldsen, I. & J. D. Sørensen: Optimization Algorithms for Calculation of the Joint De-
sign Point in Parallel Systems. Structural Optimization, Vol. 4, pp. 121-127, 1992.
[7.3] Hohenbichler, M., S. Gollwitzer, W. Kruse & R. Rackwitz: New Light on First- and Sec-
ond Order Reliability Methods. Structural Safety, Vol. 4, No. 4, pp. 267-284. 1987
[7.5] Ditlevsen, O. & H.O. Madsen: Bærende konstruktioners sikkerhed. SBI-rapport 211, Sta-
tens Byggeforskningsinstitut, 1990 (in Danish).
[7.7] Daniels, H. E.: The statistical Theory of the Strength of Bundles of Threads. Proceedings
of the Royal Society, Vol. A183, pp. 405-435, 1945.
[7.8] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
131
Note 7: Reliability evaluation of parallel system
132
Note 8: Structural reliability: Level 1 approaches
1 Introduction
During the last two decades calibration of partial safety factors in level 1 codes for structural
systems and civil engineering structures has been performed on a probabilistic basis in a number of
codes of practice, see e.g. OHBDC (Ontario Highway Bridge Design Code) [1], NBCC (National
Building Code of Canada) [2], Ravindra & Galambos [3], Ellingwood et al. [4] and Rosenblueth &
Esteva [5].
The calibration is generally performed for a given class of structures, materials and/or loads in such
a way that the reliability measured by the first order reliability index β estimated on the basis of
structures designed using the new calibrated partial safety factors are as close as possible to the
reliability indices estimated using existing design methods. Procedures to perform this type of
calibration of partial safety factors are described in for example Ravindra & Lind [6], Thoft-
Christensen & Baker [7].
A code calibration procedure usually includes the following basic steps, see e.g. Nowak [8]:
A first guess of the partial safety factors is obtained by solving an optimization problem where the
objective is to minimize the difference between the reliability for the different structures in the
class considered and a target reliability level. In order to ensure that all the structures in the class
considered have a satisfactory reliability, constraints are imposed on the reliability for the whole
range of structures. In this note it is shown how this optimization problem can be formulated and
solved. Next, the partial safety factors determined in this way are adjusted taking into account
current engineering judgement and tradition.
In section 2 the partial safety factor method is briefly described. In section 3 it is shown how partial
safety factors can be determined for a single failure mode using the results from a first order
reliability method (FORM). In section 4 a general procedure for estimating partial safety factors is
described. This procedure can be used to calibrate partial safety factors for a class of structures. In
section 5 the ‘design value method’ is presented and illustrated by an example. Finally section 6,
describes the calibration of partial safety factors in the Danish structural codes from 1999, [24].
133
Note 8: Structural reliability: Level 1 approaches
In the partial safety factor method single structural elements are usually considered and it has to be
verified that the design resistance Rd is larger than the design load effect S d for the structural
element considered:
Rd > S d (1)
This requirement has to be verified for a number of different load combinations, see below in
equation (5) and in section 6. The design value of the load effect is determined on the basis of
permanent actions, variable actions and accidental loads.
Gd = γ G Gc (2)
where
γG partial safety factor
Gc characteristic value for permanent actions, typically the 50 % quantile
Design values for variable actions are determined by, see figure 1:
Qd = γ Q Qc (3)
where
γQ partial safety factor
Qc characteristic value for variable actions, typically the 98 % quantile
134
Note 8: Structural reliability: Level 1 approaches
The design value of the load effect S d is obtained considering the following load combination
Design value for resistances Rd are usually determined using design values for material and
geometrical parameters.
135
Note 8: Structural reliability: Level 1 approaches
In code calibration based on first order reliability methods (FORM) it is assumed that the limit state
function can be written
g (x, p, z ) = 0 (7)
Using FORM (First Order Reliability Methods) the reliability index β is determined. The
corresponding estimate of the probability of failure is
Pf = Φ (− β ) (8)
If the partial safety factors and if the number of design variables is N = 1 then the design (modeled
by z ) can be determined from the design equation
G (x c , p, z, γ ) ≥ 0 (9)
x c = ( xc1 ,..., xcn ) are characteristic values corresponding to the stochastic variables X .
γ = (γ 1 ,..., γ m ) are m partial safety factors. The partial safety factors γ are usually defined such
that γ i ≥ 1, i = 1,..., m . In the most simple case, m = n .
The design equation is closely connected to the limit state function (7). In most cases the only
difference is that the state variables x are exchanged by design values x d obtained from the
characteristic values x c and the partial safety factors γ .
The characteristic values are for load variables usually the 90 %, 95 % or 98 % quantiles of the
distribution function of the stochastic variables, e.g.
where FX ( xi ) is the distribution function for X i . The design values for load variables are then
i
obtained from
136
Note 8: Structural reliability: Level 1 approaches
The characteristic values are for strength variables usually the 10 %, 5 % or 2 % quantiles of the
distribution function of the stochastic variables. The design values for strength variables are then
obtained from
xci
xdi = (11)
γi
For geometrical variables usually the median (50 % quantile) is used and the design values are
If m = n = 2 , x1 is a load variable and x 2 is a strength variable, the design equation can be written:
x
( ) ( )
G ( xc1 , xc 2 ), p, z, (γ 1 , γ 2 ) = G ( xd 1 , xd 2 ), p, z = G γ 1 xc1 , c 2 , p, z
γ2
(13)
A reliability analysis by FORM with the limit state function (7) gives the reliability index β and
the β -point x * . Partial safety factors can then be obtained from
xci
γi = for strength variables
xi*
xi*
γi = for load variables
xci
If more than one variable load type are important then e.g. the Turkstra rule can be used to model
the combined effect, see e.g. Thoft-Christensen & Baker [7]. Let X 1 ,..., X v model v different
variable load variables. The variables modeling permanent loads are denoted X v +1 ,..., X v + p and the
remaining stochastic variables are denoted X v + p +1 ,..., X n . The design equation is written
x
( )
x
G x c , p, z, γ = G γ 1ψ 1 xc1 ,..., γ vψ v xcv , γ v+1 xc ,v +1 ,..., γ v + p xc ,v + p , c ,v + p +1 ,..., cn , p, z = 0
γ γ
(14)
v + p +1 n
where ψ i ≤ 1 is a load combination factor for the variable load X i . Usually v different load
combinations are investigated where in combination j , ψ j = 1 and ψ i < 1 for j ≠ i .
g = x2 − x1
where x1 is a load variable and x 2 is a strength variable. The design equation becomes:
xc 2
G = xd 2 − xd 1 = − γ 1 xc1
γ2
* * * * * * * *
137
Note 8: Structural reliability: Level 1 approaches
1. Definition of the scope of the code, i.e. the class of structures to be considered is defined.
2. Definition of the code objective. The code objective may be defined at any higher level than the
level of the reliability method used in the code. In a level 1 reliability method (which uses a
single characteristic value of each uncertain quantity and partial safety factors) the objective
may be to obtain on average the same reliability (measured by the target reliability index β as
obtained by a reliability method on a higher level.
5. Definition of a measure of closeness between code realizations and the code objective.
6. Determination of the "best" code format, i.e. calculation of the 'optimal' partial safety factors
which gives the closest fit to the objective measured by the closeness criteria.
In general, the target reliability index can be determined by calibration to the reliability level of
existing similar structures. Alternatively or supplementary the target reliability indices can be
selected on the basis of e.g. the recommended minimum reliability indices specified in ISO [19] or
138
Note 8: Structural reliability: Level 1 approaches
NKB [10]. The maximum probability of failure (or equivalently the minimum reliability) are
assumed to be related to the consequences of failure specified by safety classes and failure types.
The following safety classes are considered, see NKB [10] and DS409 [24]:
Less serious: 1- and 2-storey buildings, which only occasionally hold persons, for instance
stock buildings, sheds, and some agricultural buildings, small pylons, roofs and
internal walls.
Serious: Buildings of more than two stories and hall structures which only occasionally
hold people, small 1- and 2-storey buildings often used for people, for example
houses, offices or productions buildings, tall pylons, scaffolds and moulds,
external walls, staircases and rails.
Very serious: Buildings of more than two stories, hall structures, and stages which will often
hold many persons and e.g. be used for offices, sports or production. 1- and 2-
storey buildings with large spans often used by many persons, stands,
pedestrian bridges, road bridges, railroad bridges.
The following failure types are considered (see NKB [10] and DS409 [24]):
Failure type I: Ductile failures where it is required that there is an extra carrying capacity
beyond the defined resistance, i.e. in the form of strain hardening.
Failure type III: Failures such as brittle failure and instability failure.
For ultimate limit states NKB [10] recommends the maximum probabilities of failure shown in
table 1 based on a reference period of 1 year. The corresponding minimum reliability indices are
shown in table 2.
As explained above calibration of partial safety factors is generally performed for a given class of
structures, materials or loads in such a way that the reliability measured by the first order reliability
index β estimated on the basis of structures designed using the new calibrated partial safety factors
is as close as possible to the target reliability index or to the reliability indices estimated using
existing design methods, see Thoft-Christensen & Baker [7], Ditlevsen & Madsen [12], Östlund
[13], Shinozuka et al. [14], Vrouwenvelder [15] and Hauge et al. [16]. Procedures to perform this
type of calibration of partial safety factors are described in e.g. Thoft-Christensen & Baker [7].
139
Note 8: Structural reliability: Level 1 approaches
In the following this procedure is described and extended in some directions. For each failure mode
the limit state function is written, see (7)
g (x, p, z ) = 0 (15)
Using FORM (First Order Reliability Methods) the reliability index β can be determined.
If the number of design variables is N = 1 then the design can be determined from the design
equation, see (9)
G (x c , p, z, γ ) ≥ 0 (16)
If the number of design variables is N > 1 then a design optimization problem can be formulated:
min C ( z)
s.t. ci ( z ) = 0 , i = 1,..., me
(17)
ci ( z ) ≥ 0 , i = me + 1,..., m
z il ≤ z i ≤ z iu , i = 1,..., N
C is the objective function and ci , i = 1,..., m are the constraints. The objective function C is often
chosen as the weight of the structure. The me equality constraints in (17) can be used to model
design requirements (e.g. constraints on the geometrical quantities) and to relate the load on the
structure to the response (e.g. finite element equations). Often equality constraints can be avoided
because the structural analysis is incorporated directly in the formulation of the inequality
constraints. The inequality constraints in (17) ensure that response characteristics such as
displacements and stresses do not exceed codified critical values as expressed by the design
equation (16). The inequality constraints may also include general design requirements for the
design variables. The lower and upper bounds, z il and z iu , to z i in (17) are so-called simple
bounds. Generally, the optimization problem (17) is non-linear and non-convex.
The application area for the code is described by the set I of L different vectors p i , i = 1,..., L . The
set I may e.g. contain different geometrical forms of the structure, different parameters for the
stochastic variables and different statistical models for the stochastic variables.
The partial safety factors γ are calibrated such that the reliability indices corresponding to the L
vectors p are as close as possible to a target probability of failure Pft or equivalently a target
( )
reliability index β t = −Φ −1 Pft . This is formulated by the following optimization problem
min W (γ ) = ∑ w j (β j (γ ) − β t )
L 2
(18)
γ j =1
where w j , j = 1,..., L are weighting factors ( ∑ Lj=1 w j = 1 ) indicating the relative frequency of
appearance / importance of the different design situations. Instead of using the reliability indices in
(18) to measure the deviation from the target, for example the probabilities of failure can be used.
Also, a nonlinear objective function giving relatively more weight to reliability indices smaller than
the target compared to those larger than the target can be used. β j (γ ) is the reliability index for
140
Note 8: Structural reliability: Level 1 approaches
combination j obtained as described below. In (18) the deviation from the target reliability index is
measured by the squared distance.
The reliability index β j (γ ) for combination j is obtained as follows. First, for given partial safety
factors γ the optimal design is determined by solving the design equation (16) if N = 1 or by
solving the design optimization problem (17) if N > 1 . Next, the reliability index β j (γ ) is
estimated by FORM on the basis of (15) using the design z from (16) or (17).
It should be noted that, following the procedure described above for estimating the partial safety
factors two (or more) partial safety factors are not always uniquely determined. They can be
functionally dependent, in the simplest case as a product, which has to be equal to a constant.
In the above procedure there is no lower limit on the reliability. An improved procedure which has a
constraint on the reliability and which takes the non-uniqueness problem into account can be
formulated by the optimization problem
(
min W (γ ) = ∑ w j (β j (γ ) − β t ) + δ ∑ γ i − γ ji
* 2
)
L 2 m
γ j =1 i =1
s.t. β j (γ ) ≥ β tmin , j = 1,..., L (19)
γ il ≤ γ i ≤ γ iu , i = 1,..., m
where w j , j = 1,..., L are weighting factors ( ∑ Lj=1 w j = 1 ). δ is a factor specifying the relative
importance of the two terms in the objective function. β j (γ ) is the reliability index for combination
j obtained as described above. γ *ji is an estimate of the partial safety factor obtained by
considering combination j in isolation. The second term in the objective function (19) is added due
to the non-uniqueness-problem and has the effect that the partial safety factors are forced in the
x*
direction of the "simple" definition of partial safety factors. For load variables: γ = . If only one
xc
x *ji
combination is considered then γ *ji = where x *ji is the design point. Experience with this
xc , ji
formulation has shown that the factor δ should be chosen to be of magnitude one and that the
calibrated partial safety factors are not very sensitive to the exact value of δ . The constraints have
the effect that no combination has a reliability index smaller than β tmin .
This type of code calibration has been used in Burcharth [17] for code calibration of rubble mound
breakwater designs. These structures are known to have reliabilities which vary considerably. The
reason is that the structures are used under widely different conditions.
As discussed above a first guess of the partial safety factors is obtained by solving these
optimization problems. Next, the final partial safety factors are determined taking into account
current engineering judgement and tradition.
Example 1
In this example partial safety factors are determined for one failure mode in one application ( L = 1) .
Consider the limit state function:
141
Note 8: Structural reliability: Level 1 approaches
g = zR − G − Q
where z is a design variable, R a resistance, G a permanent load and Q a variable load. The
stochastic model in table 3 is used.
If the target reliability index is chosen to β t =3.8 and V =0.4 then z =15.6 m 2 . The corresponding
β -point in basic variable space is (r * , g * , q * ) = (0.76, 2.04, 9.83) .
In table 4 results are shown for other coefficients of variation, V and target reliability indices. It is
seen that the partial safety factors for the variable load become rather large compared with the two
other partial safety factors, especially for large values of V .
βt V z qc / µ Q r* g* q* γR γG γQ
3.8 0.2 12.8 1.52 0.70 2.08 5.82 1.10 1.04 1.28
4.3 0.2 14.5 1.52 0.68 2.08 6.53 1.13 1.04 1.43
4.8 0.2 11.3 1.52 0.65 2.08 7.30 1.18 1.04 1.60
3.8 0.3 13.4 1.78 0.74 2.05 7.86 1.04 1.02 1.47
4.3 0.3 15.5 1.78 0.71 2.05 8.96 1.08 1.02 1.68
4.8 0.3 17.9 1.78 0.68 2.05 10.2 1.12 1.02 1.91
3.8 0.4 15.6 2.04 0.76 2.04 9.83 1.01 1.02 1.61
4.3 0.4 18.3 2.04 0.73 2.04 11.3 1.05 1.02 1.84
4.8 0.4 21.4 2.04 0.70 2.04 13.0 1.10 1.02 2.12
Table 4. Partial safety factors obtained by direct reliability-based calibration.
* * * * * * * *
142
Note 8: Structural reliability: Level 1 approaches
In the Eurocodes [18] and ISO [19] the so-called design value format is proposed to estimate partial
safety factors. According to that format the design value xd of an uncertain variable X is estimated
from
FX ( xd ) = Φ (−αβ t ) (20)
where FX is the distribution function for X and β t is the target reliability index, e.g. β t =3.8. α
is the α -coefficient associated with the type and importance of the stochastic variable considered.
The following values are recommended:
When the design value have been estimated the partial safety factor is determined from:
xc
γ =θ for strength variables
xd
x
γ =θ d for load variables
xc
Example 2
Example 1 is considered again but now the partial safety factors are determined using the design
value method. The result is shown in table 5. When compared to the results in table 4 it is seen that
the partial safety factors for resistance and permanent loads are larger than those obtained in
example 1 whereas the partial safety factors for the variable load are smaller. If the design value and
corresponding reliability index are determined using the partial safety factors in table 5, it is seen
that for V =0.4 the reliability indices are almost the same as the target reliability indices. However,
for smaller values of V the reliability indices become larger than the target reliability indices when
using the design value method.
143
Note 8: Structural reliability: Level 1 approaches
βt V qc / µ Q rd gd qd γR γG γQ z β (γ R , γ G , γ Q )
3.8 0.2 1.52 0.63 2.53 5.31 1.22 1.27 1.16 12.4 4.19
4.3 0.2 1.52 0.60 2.60 5.85 1.28 1.30 1.28 14.0 4.67
4.8 0.2 1.52 0.56 2.67 6.39 1.38 1.34 1.40 16.2 5.25
3.8 0.3 1.78 0.63 2.53 6.48 1.22 1.27 1.21 14.3 4.03
4.3 0.3 1.78 0.60 2.60 7.26 1.28 1.30 1.36 16.4 4.50
4.8 0.3 1.78 0.56 2.67 8.10 1.38 1.34 1.52 19.4 5.08
3.8 0.4 2.04 0.63 2.53 7.65 1.22 1.27 1.25 16.1 3.90
4.3 0.4 2.04 0.60 2.60 8.67 1.28 1.30 1.42 18.8 4.39
4.8 0.4 2.04 0.56 2.67 9.81 1.38 1.34 1.60 22.4 4.94
Table 5. Partial safety factors obtained using the design value method. z is the design value
obtained using the values of γ R , γ G and γ Q in columns 7-9. β (γ R , γ G , γ Q ) is the corresponding
reliability index.
* * * * * * * *
144
Note 8: Structural reliability: Level 1 approaches
This section describes the main steps in the probabilistic code calibration performed for the Danish
Structural Codes (1999). It is based on [20]. First, the reliability level is evaluated in a number of
typical, simple structures designed according to the current (old) Danish structural codes (1982) and
with a reasonable stochastic model for the uncertain quantities. The reliability analyses show a non-
uniform reliability level for different materials and actions. Next, new partial safety factors in a
slightly modified code format are calibrated such that the safety level is the same in the new code as
in the current codes, i.e. it is assumed that the reliability level in the old structural codes is
satisfactory. Using the optimized partial safety factors a more uniform reliability level is obtained
for different types of materials / structures and for different types of loads. The calibrations are
performed with the assumption that characteristic values for actions and strengths are the same in
the old and the new codes, except for changes in some of the quantile percentages used.
Three main load combinations have to be checked by design. Load combination 2 consists of four
combinations and load combination 3 consists of three combinations:
γ m = γ 0γ 1γ 2γ 3γ 4γ 5 (21)
where
γ0 takes into account the consequences of failure see table 6
γ1 takes into account the type of failure, see table 7
γ2 takes into account the possibility of unfavorable differences from the characteristic value
of the material parameter, see table 8. The values in the table are indicated by ? since they
are determined on the basis of a calibration, see the following subsections
γ3 takes into account the uncertainty in the computational model, see table 10
145
Note 8: Structural reliability: Level 1 approaches
γ4 takes into account the uncertainty in connection with determination of the material
parameter in the structure on the basis of the controlled material parameter, see table 11
γ5 takes into account the amount of control at the working place (in excess of the statistical
quality control), see table 12
For other quantile values than the 5 % quantile, the γ 2 values should be multiplied with
exp((1,65 − kγ )δ )
The stochastic model used for calibration of partial safety factors is shown in table 13 and is partly
based on the following references: SAKO [21], DGI [22] and Foschi et al. [23].
146
Note 8: Structural reliability: Level 1 approaches
For timber the Lognormal distribution is used. Also the Weibull distribution has been considered,
but since the partial safety factors in the old codes implicitely were based on Lognormal
distributions for the strengths this distribution is also used in the following calibrations. The
Lognormal distribution is considered to be reasonable for laminated timber. For single lumber
members theoretical considerations and statistical analysis of available data indicate that a Weibull
distribution should be considered. A Weibull distribution usually results in significantly smaller
reliability indices than the Lognormal distribution. Note, that in the new Danish codes all quantile
values for material strengths are defined as 5 % values.
For the variable action, two action types are used, namely imposed action and environmental action
(e.g. wind and snow).
6.4 Probabilistic calibration of partial safety factors for load combination 2.1 and 2.3
147
Note 8: Structural reliability: Level 1 approaches
The following six load cases with different ratios between permanent and variable actions are
considered for beams and columns of concrete, steel and timber and for the footings on sand and
clay:
where Gc , Qc,i and Qc, e are characteristic values for permanent actions, variable imposed action
and environmental action, respectively. For the gravity wall only load case 2.1 is considered.
Therefore, in total 98 different structures are used in the calibration. The limit state functions for the
considered failure modes are described in [20].
Each structure is first designed according to the old structural codes (1982). As described in section
6.2, the material partial safety factors are determined as a product of a number of factors, where γ 2
models the physical uncertainty related to a given type of material strength parameter. In table 14 is
shown the partial safety factors γ 2 for materials and γ f for actions corresponding to the old code.
148
Note 8: Structural reliability: Level 1 approaches
Next, reliability indices are determined for each structure using the partial safety factors in table 14.
In table 15 the average and the standard deviation of the reliability indices are shown for each of the
9 groups of structures considered and in figure 3 the distribution of reliability indices is shown.
Table 16 shows the annual probability of failure corresponding to some typical reliability indices. It
is seen that:
• The reliability for concrete and steel structures is larger than for glued laminated timber
structures.
• For the geotechnical problems the reliability levels are slightly lower than for concrete, steel and
timber structures.
• The reliability indices for the concrete beam (where the reinforcement strength is important) are
significantly larger than for the concrete column (where the concrete strength is important).
From table 15 it is also seen that the average reliability index for all example structures are 4.79. In
the next section new partial safety factors are calibrated such that the average reliability index will
remain equal to 4.79, i.e the target reliability index is β t =4.79. In table 17 the average reliability
index is shown for the six different load cases described in section 6.4.1. For each load case is
shown the parameter α (=0 if all load is permanent and = 1 if all load is variable) defined by:
Qc
α= (22)
Qc + Gc
It is seen that
• the largest reliability indices are obtained for α =0, i.e. when all the load is permanent load
• the smallest reliability indices are obtained for α =1, i.e. when all the load is variable load
• relatively high reliability indices are obtained for α =0.2 – 0.5
In summary the reliability analysis using the old partial safety factors shows a rather non-uniform
reliability level and therefore some changes in the partial safety factors can be expected if a
homogenous reliability level is to be obtained. Especially, it can be expected that the partial safety
factors for timber and variable actions are increased and the partial safety factors for concrete and
steel are decreased.
149
Note 8: Structural reliability: Level 1 approaches
Relative Frequency
1.20
1.08
0.96
0.84
0.72
0.60
0.48
0.36
0.24
0.12
0.00
3.50 3.75 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25
beta
6.4.3 Calibration of new partial safety factors for load combination 2.1 and 2.3
The code format used is basically the same as in the current codes, except for two changes:
• Separate partial safety factors for imposed variable actions and for environmental variable
actions (wind and snow) are introduced.
• In load combination 2.3 (permanent action dominating) also partial safety factors for the
variable actions are introduced.
150
Note 8: Structural reliability: Level 1 approaches
The partial safety factor γ G 21 for permanent action in load combination 2.1 is chosen to 1.0 (as in
the current code). The partial safety factors corresponding to load combination 2.1 (variable action
dominating) and 2.3 (permanent action dominating) are:
Structures of concrete, steel and timber are checked by both load combination 2.1 and 2.3.
Geotechnical problems are checked by load combination 2.1 only.
The total set of partial safety factors to be calibrated is: γ Q 21,i , γ Q 21,e , γ G 23 , γ Q 23,i , γ Q 23,e , γ a , γ c ,
γ s , γ t , γ ϕ and γ cu . The optimal partial safety factors are obtained by minimizing the deviation
between the target reliability index β t and the reliability indices obtained from the 98 example
structures. Each reliability index is determined in the following way. First, the considered structure
is designed using the deterministic design equation with characteristic values for the variables and
the actual guess on the partial safety factors. Next with the obtained design a reliability analysis is
performed, now treating the uncertain quantities as stochastic variables.
where β i (γ ) is the reliability index for example structure no i designed with partial safety factors
γ . L = 98 is the number of example structures. ω i =1 is assumed for all structures.
Results
Using a target reliability index β t equal to the average of the reliability indices of the example
structures and a standard nonlinear optimization program, the results shown in table 18 are
obtained. The table shows the optimized partial safety factors directly and modified with the factors
from table 14. The modification consists of multiplying the optimized partial safety factors γ 2 with
the relevant values of the γ 1 , γ 3 , γ 4 and γ 5 factors in the old codes, see table 14. Further, also
optimized and modified partial safety factors obtained by rounding and fixing the action partial
safety factors are shown.
151
Note 8: Structural reliability: Level 1 approaches
It is seen that
• the partial safety factor for imposed actions in load combination 2.1 is almost unchanged while
the partial safety factor for environmental actions should be increased from 1.3 to 1.5.
• In load combination 2.3 the partial safety factor for permanent action is unchanged 1.15.
• the partial safety factors for reinforcement, concrete and steel can be decreased.
• the partial safety factor for glued laminated timber should be increased
• the partial safety factors for geotechnical parameters are almost unchanged.
152
Note 8: Structural reliability: Level 1 approaches
Generally, the partial safety factors for actions should be increased while the partial safety factors
for the material strengths can be decreased. Further it is, as expected, seen that the difference in
reliability levels in the example structures are much smaller using the optimized partial safety
factors than using the partial safety factors in the old code (1982).
Relative Frequency
1.20
1.08
0.96
0.84
0.72
0.60
0.48
0.36
0.24
0.12
0.00
3.50 3.75 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25
beta
Figure 4. Reliability indices for the 98 example structures using calibrated partial safety factors.
The material partial safety factors in the new structural codes are based on the above calibration
results. The values in the following table have been chosen.
In order to evaluate the reliability level in load combination 2.1 and 2.3 the following simple, but
representative limit state function is considered:
153
Note 8: Structural reliability: Level 1 approaches
where
R strength (modeled by a stochastic variable)
XR model uncertainty (modeled by a stochastic variable)
z design variable, e.g. a cross-sectional area
G permanent action (modeled by a stochastic variable)
Q variable action (modeled by a stochastic variable)
α factor between 0 and 1, giving the relative importance of the variable action.
As examples only steel and concrete structures are considered. As variable actions both
environmental and imposed actions are used. The following stochastic model, based on table 13, is
used:
The design variable z is determined by considering load combination 2.1 and 2.3. The design
equations can then be written:
( )
z1 Rc / γ m − (1 − α )γ G1 Gc + αγ Q1 Qc = 0 (25)
(
z 3 Rc / γ m − (1 − α )γ G Gc + αγ Q Qc = 0
3 3
) (26)
-imposed γ Q = 1.3
1
γ Q = 1.0
3
The reliability index β is determined as function of α . The result is shown in figure 5. It is seen
that an almost uniform distribution of the reliability index is obtained as a function of α . Only steel
structures have a larger reliability when permanent actions are dominating (small α ). However, in
practice variable actions will usually be dominating for steel structures.
154
Note 8: Structural reliability: Level 1 approaches
8.00
6.00
nyttelast - stål
nyttelast - beton
naturlast - stål
naturlast - beton
β 4.00
2.00
0.00
Figure 5. Reliability index β for different combinations of variable actions and material type.
r = γ m (γ G (1 − q) + γ Q q ) (27)
where
zRc
r= (28)
(1 − α )Gc + αQc
αQc
q= (29)
(1 − α )Gc + αQc
It is seen that r can be considered as a ‘total safety factor’, namely as a product of the material
partial safety factor and the weighted action partial safety factor. q is seen to be a measure of the
characteristic variable action compared to the total characteristic action.
155
Note 8: Structural reliability: Level 1 approaches
2.40
naturlast - beton
2.20
1.80
r nyttelast - stål
1.60
1.40
1.20
1.00
In figure 6 the ‘total safety factor’ r is shown as a function of q for different combinations of
variable actions and material type. It is noted that the smallest value of the ‘total safety factor’ r is
obtained for values of q between 0.2 and 0.3. A comparison with figure 5 shows that it is not in this
interval that the smallest reliability index is obtained; in fact the largest reliability indices are
obtained in this interval. The reasons for this are among others that
• the variable actions have a larger coefficient of variation than the permanent actions
• the ‘total safety factor’ r is based on characteristic values which have some ‘safety’ included
since they are obtained as quantile values in the distribution functions for actions and strengths.
156
Note 8: Structural reliability: Level 1 approaches
The final partial safety factors for actions in DS409 [24] are given in table 24.
Load combination
Serv. Ultimate Accidental
Load type 1 2.1 2.21) 2.3 2.4 3.1 3.2 3.3
Permanent action
Self weight
Gc 1.0 1.0 0.8 0.9 1.0 1.0 1.0 1.0
0.25 Gc free action - - - 1.0 - - - -
Weight of soil and 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
ground water
Variable action
Variable action
imposed action - 1.3 1.3 1.0 1.0-1.3 ψ ψ ψ
environmental action - 1.5 1.5 1.0 1.0-1.3 - - ψ
Other variable actions - ψ ψ ψ 1.0-1.3 ψ ψ ψ
Horizontal mass load - 1.0 1.0 1.0 - - - 0.25
Accidental action – - - - - - 1.0 - -
impact, etc.
Accidental load – fire - - - - - - - 1.0
Table 24. Load combinations, partial safety factors and load combination factors ψ . The partial
1)
safety factors for load combination 2.2 are for normal safety class. For low and high safety classes
the partial safety factors for variable action have to be multiplied by γ 0 , see table 6.
The final partial safety factors for materials can be determined using (21) in section 6.2 with γ 2
selected according to the coefficient of variation δ for the material considered. In table 25 γ 2
values are shown for different values of δ . The γ 2 values are determined such that the average
reliability level is obtained for the relevant value of δ . The values for δ =0.05 and 0.15 shown in
the table are derived in section 6.4.3.
157
Note 8: Structural reliability: Level 1 approaches
7 References
[1] OHBDC (Ontario Highway Bridge Design Code), Ontario Ministry of Transportation and
Communication, Ontario, 1983.
[2] NBCC (National Building Code of Canada), National Research Council of Canada, 1980.
[3] Ravindra, M.K. & T.V. Galambos: Load and Resistance Factor Design for Steel. ASCE,
Journal of the Structural Division, Vol. 104, N0. ST9, pp. 1337-1353, 1978.
[4] Ellingwood, B., J.G. MacGregor, T.V. Galambos & C.A. Cornell: Probability Based Load
Criteria: Load Factors and Load Combinations. ASCE, Journal of the Structural Division,
Vol. 108, N0. ST5, pp. 978-997, 1982.
[5] Rosenblueth, E. & L. Esteva : Reliability Basis for Some Mexican Codes. ACI Publication
SP-31, pp. 1-41, 1972.
[6] Ravindra, M.K. & N.C. Lind : Theory of Structural Code Calibration. ASCE, Journal of the
Structural Division, Vol. 99, pp. 541-553, 1973.
[7] Thoft-Christensen, P. & M.B. Baker: Structural Reliability Theory and Its Applications.
Springer Verlag, 1982.
[8] Nowak, A.S.: Probabilistic Basis for Bridge Design Codes. Proc. ICOSSAR'89, pp. 2019-
2026, 1989.
[9] Melchers, R.E.: Structural Reliability, Analysis and Prediction. John Wiley & Sons, 1987.
[10] Recommendation for Loading- and Safety Regulations for Structural Design. NKB-report No.
36, 1978.
[11] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
[12] Ditlevsen, O. & H.O. Madsen: Bærende konstruktioners sikkerhed. SBI-rapport 211, Statens
Byggeforskningsinstitut, 1990 (in Danish).
[13] Östlund, L.: General European Principles of Codes Concerning Reliability. Proc.
ICOSSAR'89, pp. 1943-1948, 1989.
[14] Shinozuka, M. & H. Furuta & S. Emi: Reliability-Based LRFD for Bridges : Theoretical
Basis. Proc. ICOSSAR'89, pp. 1981-1986, 1989.
[15] Vrouwenvelder, A.C.W.M. & A.J.M. Siemes: Probabilistic Calibration Procedure for the
Derivation of Partial Safety Factors for the Netherlands Building Codes. HERON, Vol. 32,
No. 4, pp. 9-29, 1987.
[16] Hauge, L.H., R. Loseth & R. Skjong: Optimal Code Calibration and Probabilistic Design.
Proc. OMAE'92, Vol. II, pp. 191-199, 1992.
[17] Burcharth, H.F.: Development of a Partial Safety Factors System for the Design of Rubble
Mound Breakwaters. PIANC Working Group 12, December 1991.
158
Note 8: Structural reliability: Level 1 approaches
[18] Eurocode 1: Basis of design and actions on structures - Part 1: Basis of design. ENV 1991-1,
1994.
[20] Sørensen, J.D., S.O. Hansen & T. Arnbjerg Nielsen. Calibration of Partial Safety Factors for
Danish Structural Codes. 2000.
[21] SAKO: Probabilistic Calibration of Partial Safety Factors in the Eurocodes. 1999.
[22] Danish Geotechnical Institute: Partial factors of safety in geotechnical engineering. Report
No. 1, 1993.
[23] Foschi, R.O. & B.R. Folz & F.Z. Yao: Reliability-based design of wood structures. Structural
Research Series, Report no. 34, Department of Civil Engineering, University of British
Colombia, Vancouver, Canada, 1989.
[24] DS409:1998 – Code of practice for the safety of structures. Danish Standard, 1999.
159
Note 8: Structural reliability: Level 1 approaches
160
Note 9: Time-variant reliability
1 Introduction
In the previous lectures it has been assumed that all variables could be considered either to be time-
invariant stochastic variables or deterministic parameters. However, loads such as wave-loads,
snow-loads and wind-loads are usually modeled as time-varying stochastic processes. In this case
we are usually interested in determining the probability that the load within a given period of time
exceeds a given threshold, the so-called barrier crossing problem. Further, it is of interest to deter-
mine the distribution of the maximum and minimum values of the process. This note is partly based
on an earlier lecture note by S. Engelund.
2 Stochastic processes
A stochastic process is an indexed set of random variables { X (t ), t ∈ [0, T ]} defined in the sample
space Ω . The index variable t is here assumed to be time, defined on the time interval [0, T ] . Fig-
ure 1 shows a realization of a stochastic process.
At time t the stochastic variable X (t ) with realization x(t ) is described by the distribution function
At times t1 and t 2 the joint distribution function for the two random variables X (t1 ) and X (t 2 ) is,
see figure 1
161
Note 9: Time-variant reliability
FX ( x1 , x2 ,..., xn ; t1 , t 2 ,..., t n ) =
(3)
P({X 1 = X (t1 ) ≤ x1} ∩ {X 2 = X (t 2 ) ≤ x2 } ∩ ... ∩ {X n = X (t n ) ≤ xn })
∂ n FX ( x1 , x2 ,..., xn ; t1 , t 2 ,..., t n )
f X ( x1 , x2 ,..., xn ; t1 , t 2 ,..., t n ) = (4)
∂x1∂x2 ...∂xn
∞
µ X (t ) = E[ X (t )] = ∫ xf X ( x, t )dx (5)
−∞
∞ ∞
R XX (t1 , t 2 ) = E[ X (t1 ) X (t 2 )] = ∫ ∫ x1 x2 f X ( x1 , x2 ; t1 , t 2 )dx1dx2 (6)
− ∞ −∞
σ X2 (t ) = C XX (t , t ) = R XX (t , t ) − µ X2 (t ) (8)
C XX (t1 , t 2 )
ρ XX (t1 , t 2 ) = (9)
σ X (t1 )σ X (t 2 )
If all finite dimensional distribution functions FX (x) , FX ( x1 , x2 ) ,…. are invariant to a linear trans-
lation of the time origin then the process is called strictly stationary. If this invariance assumption
only holds for FX (x) and FX ( x1 , x2 ) then the process is called weakly stationary. For a stationary
process FX ( x; t ) becomes independent on time and FX ( x1 , x2 ; t1 , t 2 ) only becomes dependent on the
time difference τ = t1 − t 2 . Similarly, R XX (t1 , t 2 ) , C XX (t1 , t 2 ) and ρ XX (t1 , t 2 ) will only be dependent
on τ = t1 − t 2 .
For a stationary stochastic process, the spectral density is related to the covariance function by the
Wiener-Khintchine equations:
1 ∞
S X (ω ) = ∫ C XX (τ ) exp(−iωτ )dτ (10)
2π −∞
162
Note 9: Time-variant reliability
∞
C XX (τ ) = ∫ S X (ω ) exp(iωτ )dω (11)
−∞
where ω is the circular frequency in radians per second. From (8) then follows that the variance of
the stationary process is:
∞
C X2 = C XX (0) = ∫ S X (ω )dω (12)
−∞
If measurements of a stationary stochastic process are made, then usually only one realization be-
comes available. In that case the expected value is estimated by:
1T
µ= ∫ x(τ )dτ (13)
T0
If this time average approaches µ X for T → ∞ the process is ergodic in mean value. Similarly if
1 T −τ
R (τ ) = ∫ x(t + τ ) x(t )dτ (14)
T −τ 0
approaches R XX (τ ) for T → ∞ the process is ergodic in correlation. If this property holds for all
moments then the process is called ergodic.
f X ( x1 , x2 ,..., xn ; t1 , t 2 ,..., t n ) =
1 1 n
[ ] (x
exp − ∑ ( xi − µ X (ti ) ) C −1 − µ X (t j ) )
(15)
(2π )
n/2 ij j
C 2 i , j =1
and [C −1 ] ij denotes element i, j in the inverse covariance matrix C −1 . A Gaussian process is thus
completely determined by µ X (t ) and C XX (t1 , t 2 ) . Therefore a stationary Gaussian process is strictly
stationary. The derivative process is also Gaussian:
d
X& (t ) = X (t ) (17)
dt
E[ X& ] = 0 (18)
163
Note 9: Time-variant reliability
d 2 R XX (t )
E[ X& 2 ] = − 2
= − R XX
''
(0) (19)
dt t =0
&
E[ X X ] = 0 (20)
Consider a stationary Gaussian process with mean value µ X and standard deviation σ X . Since
X (t ) is a stationary process the mean value of X& is µ X& = 0 , see (18). The standard deviation of X&
is denoted σ & . The joint density function of X and X& is then
X
x&
1 x − µX
2 2
1
f XX& ( x, x& ) = exp − +
(21)
2πσ X σ X& 2 σ X σ X&
3 Barrier Crossing
In the following a number of different methods by which estimates of (22) can be obtained are pre-
sented.
3.1 Simulation
Monte Carlo simulation of stochastic processes has attracted much attention in the recent years.
Partly because the development of more efficient computers the method has become more attractive
and partly because it often is the only available method to determine the reliability of complicated
nonlinear structural systems. The most commonly used method for simulating Gaussian processes is
the so-called spectral representation method proposed by Borgman [1].
M −1
X M (t j ) = ∑ 2 S X (ω k )∆ω cos(ω k t j + Θ k ) (23)
k =0
164
Note 9: Time-variant reliability
where S X (ω ) is the one-sided spectrum of the stochastic process and ω k = k∆ω . The phases, Θ k ,
are stochastic variables, independent and uniformly distributed in the interval [0;2π ] . The process
X M (t ) is asymptotically Gaussian as M becomes large due to the central limit theorem. Further, it
2π
is important to notice that the process X M (t ) is periodic with the period . It is evident that for
∆ω
longer time histories and finer spectral resolution the computation time becomes excessive. Fortu-
nately, this problem can be overcome by performing the summation in (23) by Fast Fourier Trans-
formation (FFT). The failure probability now can be determined by simulating a large number of
realizations of X (t ) and determining the relative number of times X (t ) exceeds the threshold
value, ξ .
NC
Pf = (24)
N
where N C denotes the number of realizations which exceeds the threshold value and N denotes the
number of realizations of X (t ) .
The simulation method is not restricted to Gaussian processes. It is, however, more complicated to
simulate Non-Gaussian processes. The major disadvantage of the method is the fact that it requires
a very large number of simulations in order to determine an out-crossing probability if the out-
crossings events are rare. In that case the method is very inefficient even if the Fast Fourier Trans-
formation is applied to perform the summation.
Let pk denote the probability of exactly k out-crossings in the interval [0; T ] . It is then evident that
the probability of no out-crossings or the complementary first passage probability is
p0 = 1 − Pf
∞
= 1 − ∑ pk
k =1
∞ ∞ ik
= 1 + ∑ pk ∑ (− 1)
k =1 i =1 i
∞ ( −1) i ∞ k
= 1+ ∑ ∑ i! pk (25)
i =1 i! k =1 i
∞(−1) i ∞
= 1+ ∑ ∑ k (k − 1)...(k − i + 1) pk
i =1 i! k =1
∞ ( −1) i
= 1+ ∑ mi
i =1 i!
1 1
= m0 − m1 + m2 − m3 + ...
2 6
165
Note 9: Time-variant reliability
m0 = 1
∞ (26)
mi = ∑ k (k − 1)...(k − i + 1) pk for i ≥ 1
k =1
k
= 0 for i > k (27)
i
(25) is the so-called Rice's '' in- and exclusion '' series (see Rice [4] which provides an exact solu-
tion to the barrier crossing problem. Of course, the moments mi (i = 1,2,...) must exist and the se-
ries in (25) must converge in order to make (25) a valid representation. The series provides lower
and upper bounds for the survival probability upon truncation after an odd or even term, respec-
tively. The computational effort involved in evaluating Pf (t ) according to this method, however, is
extensive. Further an increasing number of terms has to be taken into account as m1 increases.
Normally the series is truncated after the first term. This provides an upper limit for the failure
probability
Pf ≤ m1 (28)
where m1 is nothing but the mean number of out-crossings. It is evident that Pf (t ) can only be ap-
proximated by m1 if the out-crossing probability is very small, i.e. Pf <<1.
Let the process N + (t , ξ ) be a process that increases by one each time the process X (t ) exceeds the
threshold ξ and let N + (0, ξ ) = 0 . Obviously N + (t , ξ ) is a counting process which counts the num-
ber of exits of X (t ) across ξ .
If it is now assumed that the probability of having two or more out-crossings in ] t , t + ∆t ] is negli-
gible compared to the probability of having exactly one out-crossing, if ∆t is sufficiently small, and
further that the out-crossings in ] t , t + ∆t ] are independent of the previous out-crossings in ]0, t ] ,
then N + (t ) is a Poisson process. The probability that the number of out-crossings N + (t , ξ ) is equal
to n can be determined as
(
P N + (t , ξ ) = n =) 1
n!
(λ (t , ξ ) )n exp(−λ (t )) (29)
λ (t , ξ ) = E[ N + (t , ξ )] = m1 (30)
( )
Pf (t ) = 1 − P N + (t , ξ ) = 0 = 1 − exp(− m1 ) (31)
166
Note 9: Time-variant reliability
For broad-banded processes the correlation length is of the magnitude equal to the zero up-crossing
period. In this case the maxima between succeeding zero-upcrossings are virtually uncorrelated.
Hence, the out-crossings from the safe domain related to these maxima will also be independent and
(31) is valid.
For narrow-banded processes, the out-crossings in case of low to medium barrier levels tend to oc-
cur in clumps, see figure 3. In this case the crossing events are highly correlated, and (31) is no
longer appropriate. However, at higher barrier levels only the highest peak in a clump is likely to
imply an out-crossing. This suggests that the out-crossings tend to become independent as ξ → ∞ .
Actually, this hypothesis can be formally proved for Gaussian processes, see Cramer and Leadbetter
[3].
By (25) and (31) one determines the probability that X (t ) at some time crosses the threshold, ξ . It
has not been taken into account that the process might start in the failure region, i.e. X (0) > ξ . By
taking the initial condition into account the failure probability can be defined as
dPf (T )
= f1 (t ) P( X (0) < ξ ) (33)
dT
where f1 (t ) is the probability density function of the time to the first barrier crossing conditional on
X (0) < ξ . No exact solutions for f1 (t ) are available even for very simple problems. Hence, it is
necessary to determine some approximation by which the failure probability can be determined.
167
Note 9: Time-variant reliability
In order to determine the mean number of exits of X (t ) across the level ξ it is convenient to con-
sider the stochastic process Y (t ) given by
Y (t ) = H ( X (t ) − ξ ) (34)
where H (.) is Heavisides step function. By differentiation of Y (t ) the derivative process Y& (t ) can
be determined by
where δ (.) denotes the dirac delta function. In (35) it has been assumed that X (t ) is a differenti-
able process. For a realization of X (t ) the corresponding realizations of Y (t ) and Y& (t ) are shown
in figure 4. It is seen that Y& (t ) consists of a series of unit pulses which occurs each time an out-
crossing of X (t ) occurs. The number of out-crossings, N (T , ξ ) , within the time interval ]0, T ] can
be determined by integrating the absolute value of Y& (t )
T T
N (T , ξ ) = ∫ Y& (τ ) dτ = ∫ X& (τ ) δ ( X (τ ) − ξ )dτ (36)
0 0
T
E [N (T , ξ )] = ∫ E[ X& (τ ) δ ( X (τ ) − ξ )]dτ
0
T ∞ ∞
= ∫ ∫ ∫ x& (t ) δ ( x(t ) − ξ ) f XX& ( x, x& , τ )dx dx& dτ (37)
0 −∞ −∞
T ∞
= ∫ ∫ x& (t ) f XX& (ξ , x& , τ ) dx&dτ
0 −∞
168
Note 9: Time-variant reliability
where f XX& is the joint density function of X and X& . It should be noted that by deriving (37) both
the up-crossings and down-crossings have been taken into account. However, for a stationary proc-
ess it is reasonable to assume that any positive crossing is followed by a negative crossing:
[ ] [
E N + (T , ξ ) = E N − (T , ξ ) = ] 1
2
E [N (T , ξ )] (38)
where N − (T , ξ ) counts the number of down-crossings of X and X& of the level ξ . This implies
that
[ ]
T∞
E N + (T , ξ ) = ∫ ∫ x&f XX& (ξ , x& , τ ) dx&dτ
00 (39)
= m1
It is often convenient to consider the rate of out-crossings pr unit time, ν + (t , ξ ) which is defined by
∞
ν + (t , ξ ) = ∫ x& f XX& (ξ , x& , t ) dx& (40)
0
which is the so-called Rice's formula, see [4]. For stationary processes the out-crossing intensity
does not depend on t i.e. ν + (t , ξ ) = ν + (ξ ) .
From (28) follows that an upper bound of the probability of failure in the time interval ]0, T ] is
T
Pf (T ) ≤ m1 = ∫ ν + (t , ξ )dt (41)
0
If ν + (t , ξ ) = ν + (ξ ) then
Pf (T ) ≤ ν + (ξ )T (42)
Higher order factorial moments and factorial moments of the number of out-crossing of a given safe
domain by a vector process can be determined on the basis on the so-called Belyaev's formula, see
[1]. This formula, however, can only be solved analytically in a few special cases and a numerical
solution is generally a non-trivial task.
169
Note 9: Time-variant reliability
We have now determined the mean number of out-crossings of X (t ) without taking into account
the initial conditions. The mean number of N + (t ) given X (0) < ξ is often approximated by the
unconditional mean value, m1 . By using (11) one then obtains
(
Pf (T ) = 1 − (1 − Pf (0)) exp - E[ N + (T, ξ )] ) (43)
It has, however, been shown that a better approximation for the mean number of out-crossings
given X (0) < ξ is given by
[
E N + (T, ξ ) X(0) < ξ ≈ ] E[ N + (T, ξ )]
1 − Pf (0)
(44)
whereby
E[ N + (T, ξ )]
Pf (T ) = 1 − (1 − Pf (0)) exp - (45)
1 − P (0)
f
This expression has been shown to yield very accurate results even for relatively low threshold lev-
els, where the out-crossings are not independent.
Let X (t ) be a stationary Gaussian process with density function given by (21). For a given thresh-
old ξ the out-crossing intensity now can be determined on the basis of Rice's formula (40):
∞
ν + (ξ ) = ∫ x&f XX& (ξ , x& ) dx&
0
x&
1 ξ − µ X
2 2
∞ 1
= ∫x exp − +
σ & dx
& &
0 2πσ X σ X& 2 σ X X
(46)
1 ξ − µ 2
∞ 2
=
1
exp − X
∫ x& exp − 1 x&
dx&
2πσ X σ X& 2 σX 0 2 σ &
X
σ X& 1 ξ − µ
2
= exp − X
2πσ X 2 σX
1 σ X&
ν + (µ X ) = (47)
2π σ X
170
Note 9: Time-variant reliability
Example 1
Consider a stationary Gaussian process with
µX = 1 σ X = 0 .3 σ X& = 0.2
ξ =2
0.2 1 2 − 1 2
ν + (2) = exp − =0.00041
2π 0.3 2 0.3
0.2
ν + (1) = =0.11
2π 0.3
* * * * * * * *
Figure 5. Realization of narrow-banded stochastic process and density function of local extremes.
First, consider the simple case of a stationary narrowband Gaussian process, X (t ) . A realization
of a narrow-band process is shown in figure 5. For an ideally narrow-band process the rate of zero-
crossings is equal to the rate of occurrence of maxima. Further the rate of crossings of the level ξ is
equal to the rate of occurrence of maxima above ξ . Therefore, the ratio ν + (ξ ) / ν + (0) may be inter-
preted as the complementary distribution function of the local maxima, Ξ
ν + (ξ ) 1 ξ − µ
2
FΞ (ξ ) = 1 − P(Ξ > ξ ) = 1 − + = 1 − exp − X
, ξ ≥ µX (48)
ν ( 0) 2 σX
171
Note 9: Time-variant reliability
ξ − µX 1 ξ − µ
2
f Ξ (ξ ) = exp− X
, ξ ≥ µX (49)
σX2
2 σX
Example 2
Using the same data as in example 1 and assuming that the process is narrow-banded, the density
function of local maxima (peaks) becomes
ξ −1 1 ξ − 1 2
f Ξ (ξ ) = exp − ξ ≥1
0.32 2 0.3
and the expected number of peaks between 1.5 and 1.6 becomes
1 1.5 − 1 2 1 1.6 − 1 2
P (1.5 ≤ Ξ ≤ 1.6) = exp − − exp − =0.249-0.135=0.114
2 0.3 2 0.3
* * * * * * * *
Example 3
Consider a stationary stochastic process where the joint density function of X and X& is:
c( x − 1) (1 − x& )
2
for ( x, x& ) ∈ [− 1;1] × [− 1;1]
f XX& ( x, x& ) =
0 else
1 1
From the condition that ∫ ∫ f XX& ( x, x& )dxdx& = 1 the constant c is determined:
−1−1
172
Note 9: Time-variant reliability
∫ ∫ c( x − 1) (1 − x& )dxdx& = 1
1 1 2
⇒
−1−1
11
4c ∫ ∫ ( x − 1) (1 − x& )dxdx& = 1
2
⇒
00
c = 1 .5
Thus
The density function is illustrated in figure 6. The marginal density functions are
−1
∫ ( x − 0) 1.5( x − 1) dx =
1
µ X = ∫ x1.5( x − 1) dx = 0
1 1 2
2
σX = 2
−1 −1 10
−1 −1 6
The out-crossing intensity of the level ξ =0.8 is determined using Rice’s formula, see (40):
∞ 1
ν + (ξ = 0.8) = ∫ x&f XX& (ξ = 0.8, x& ) dx& = ∫ x& 1.5(0.8 − 1)2 (1 - x& ) dx& = 0.0100
0 0
If the process is approximated by a Gaussian process with the same expected values and standard
deviations then out-crossing intensity becomes, see (46):
2
1 0.8 − 0
1
6 exp −
ν + (0.8) = = 0.0089
2π 1
2 1
10 10
and the expected number of peaks above 0.8 becomes, see (48):
2
1 0.8 − 0
P(0.8 ≤ Ξ) = exp − = 0.041
2 1
10
* * * * * * * *
173
Note 9: Time-variant reliability
In the following is considered a normalized process X m with expected value equal to zero and unit
standard deviation. A realization xm is obtained from:
x − µX
xm = (50)
σX
For non-narrowband Gaussian processes an expression for the distribution of local maxima can be
derived on the basis of Rice's formula, (40). Using the fact that the occurrence of a maxima of X (t )
implies a down-crossing of X& (t ) of the level ξ = µ X , and by introducing the so-called irregularity
factor
Rice [4] has derived the following expression for the density function of the local maxima:
xm x 2 αxm
fX ( xm ) = 1 − α ϕ
2
+ αxm exp − m Φ
2
(52)
1−α 1−α 2
2
m
where Φ (.) denotes the standard Normal distribution and ϕ (.) denotes the standard Normal density
function. The irregularity factor α takes on values in the interval between zero and one. It can be
shown that when α =1 (an ideally narrow-band process) (52) gives the Rayleigh distribution, eq.
(49). When α is approximately equal to zero, the density function of the local extremes, eq. (52),
tends to the Gaussian density function with zero mean and standard deviation σ X . This shows that
maxima occur randomly and with equal probability of being above and below zero.
6 Global Extremes
FT ( xm ) = FX ( xm ) N = (1 − (1 − FX ( xm ))) N
m
m
m
m
(53)
xm x 2 αxm
1 − FX ( xm ) = 1 − Φ + α exp − m Φ
2
(54)
1−α 1−α 2
2
m
174
Note 9: Time-variant reliability
x2
1 − FX ( xm ) ≈ α exp − m (55)
2
m
( )
Φ ( z ) ≈ 1 − ϕ ( z ) z −1 − z −3 + ... (56)
x2
y = N m (1 − FX ( xm )) = N exp − m (57)
2
m
and using the fact that the largest of N m observed maxima is located around the 1 / N m fractile,
which implies that the variable y is of order unity for increasing N m , we obtain
Nm
y
FT ( xm ) = 1 −
N m
y
= exp N m log1 −
N m
(58)
≈ exp(− y )
x2
= exp − N exp − m
2
The mean value and the standard deviation of the maximum value in the interval [0, T ] now can be
determined on the basis of (58). It is found that
0.577
µ max = 2 log N + (59)
2 log N
π 1
σ max = (60)
6 2 log N
In the Danish codes of practice for wind loads the characteristic wind load is determined using (59)
and (60), see [5].
175
Note 9: Time-variant reliability
A stochastic process {X (t )} is considered where the length of the single events is very small and
can be considered as spikes. The size of the spikes are modeled as independent stochastic variables
X with distribution function FX (x) . The times t1 , t 2 ,... of the events are assumed to be modeled by
a Poisson process with intensity λ (t ) , see section 3.3. In figure 7 is shown the case with a time
varying threshold ξ (t ) . The events where the ‘spikes’ exceed the threshold can be considered as a
new Poisson process with intensity
t
Pf (t ) = 1 − exp − ∫ λξ (τ )dτ (62)
0
t
Pf (t ) = 1 − exp − ∫ λ (τ )[1 − FX (ξ )]dτ (63)
0
When the threshold is constant the probability of no failure in the time interval [0, T ] is equal to the
distribution function FT (ξ ) of the maximum value of the stochastic process. Therefore
T
FT (ξ ) = P max {X (t )} ≤ ξ = 1 − Pf (T ) = exp − ∫ λ (τ )[1 − FX (ξ )]dτ (64)
t∈[ 0,T ] 0
176
Note 9: Time-variant reliability
In figure 8 is shown a so-called square-wave Poisson process. The process {X (t )} is constant be-
tween the changes in size at time t1 , t 2 ,... . The times where the process changes its value are as-
sumed to be modeled by a Poisson process with intensity λ (t ) , see section 3.3. The threshold is
assumed to be constant equal to ξ .
It can be shown that the distribution function FT (ξ ) of the maximum value of the stochastic process
in the time interval [0, T ] is, see e.g. Madsen et al. [6]
T
FT (ξ ) = FX (ξ ) exp − [1 − FX (ξ )]∫ λ (τ )dτ (66)
0
7 References
[1] Belyaev, Y. K.: On the Number of Exits Across the Boundary of a Region by a Vector Sto-
chastic Process, Theor. Probab. Appl., 1968, 13, pp. 320-324.
[2] Borgman, L. E.: Ocean Wave Simulation for Engineering Design, J. Wtrwy. and Harb. Div.,
ASCE, 95, 1969, pp. 557-583.
[3] Cramer, H. & M.R. Leadbetter: Stationary and Related Stochastic Processes, Wiley, New
York, 1967.
[4] Rice, S. O.: Mathematical Analysis of Random Noise, in: Selected Papers on Noise and Sto-
chastic Processes, Ed.: N. Wax, Dover Publications, New York, 1954.
[5] DS410:1998 – Code of practice for loads for the design of structures. Danish Standard, 1999.
[6] Madsen, H.O., S. Krenk & N.C. Lind: Methods of Structural Safety. Prentice-Hall, 1986.
177
Note 9: Time-variant reliability
178
Note 10: Load combinations
1 Introduction
In this note the load combination problem is considered. The situation is considered where two or
more variable loads act on a structure. How is for example the annual maximum combined load
described and how are characteristic and design values determined? These are questions considered
in this note, which is partly based on [1] and [2].
2 Exact model
X (t ) = X 1 (t ) + X 2 (t ) (1)
with derivative X& (t ) = X& 1 (t ) + X& 2 (t ) . Figure 1 shows a realization of the stochastic processes. Note
that in the time interval considered the maximum value for X (t ) does not occur at the same time as
the maximum values for X 1 (t ) or X 2 (t ) . The maximum value of X (t ) in the time interval [0, T ] is
179
Note 10: Load combinations
denoted X max,T . T could for example be 1 year. The distribution function for X max,T is equal to 1
minus the probability that the maximum value exceeds a threshold ξ in [0, T ] :
FX max,T
(ξ ) = 1 − P(max X (t ) > ξ , t ∈ [0, T ])
= 1 − P( X (0) > ξ ) − P(one or more out - crossings of the barrier ξ X (0) ≤ ξ )
∞
≅ 1 - P( X (0) > ξ ) − ∑ P(n out - crossings of the barrier ξ )
n =1
∞
≥ 1 - P( X (0) > ξ ) − ∑ nP(n out - crossings of the barrier ξ ) (2)
n =1
= 1 - P( X (0) > ξ ) − ν X (ξ )T
≈ 1 -ν X (ξ )T
∞
= 1 − T ∫ x&f XX& (ξ , x& )dx&
0
∞
where ν X (ξ ) is the out-crossing intensity determined by Rice’s formula (ν X (ξ ) = ∫ x&f XX& (ξ , x& )dx& ).
0
To calculate ν X (ξ ) we need the joint density junction f XX& ( x, x& ) . First, it is seen that the distribu-
tion and density functions for X can be obtained from the convolution integrals:
FX ( x ) = P ( X 1 + X 2 ≤ x )
= ∫ f X1 ( x1 ) f X 2 ( x 2 )dx1dx2
x1 + x2 ≤ x
∞ x − x1
(3)
= ∫ ∫ f X 2 ( x2 )dx2 f X1 ( x1 )dx1
− ∞ −∞
∞
= ∫ FX 2 ( x − x1 ) f X1 ( x1 )dx1
−∞
∂FX ( x)
f X ( x) =
∂x
∞ ∂F ( x − x )
X2 1
= ∫ f X1 ( x1 )dx1 (4)
−∞ ∂x
∞
= ∫ f X 2 ( x − x1 ) f X1 ( x1 )dx1
−∞
Similarly:
∞ ∞
f XX& ( x, x& ) = ∫ ∫ f X1 X& 1 ( x1 , x&1 ) f X 2 X& 2 ( x − x1 , x& − x&1 )dx1dx&1 (5)
−∞ −∞
180
Note 10: Load combinations
∞ ∞ ∞
ν X (ξ ) = ∫ x& ∫ ∫ f X X& ( x1 , x&1 ) f X
1 1 2X2
& (ξ − x1 , x& − x&1 )dx1dx&1dx&
0 − ∞ −∞
∞ ∞ ∞
= ∫ ∫ ∫ ( x&1 + x& 2 ) f X X& ( x1 , x&1 ) f X
1 1 2X2
& (ξ − x1 , x& 2 )dx& 2 dx1dx&1
0 −∞ − x&1
(6)
∞ ∞
= ∫ ∫ x&1 f X X& ( x1 , x&1 ) f X
1 1
&
2X2
(ξ − x1 , x& 2 )dωdx1 +
−∞ ω
∞ ∞
∫ ∫ x& 2 f X X& ( x1 , x&1 ) f X
1 1
&
2X2
(ξ − x1 , x& 2 )dωdx1
−∞ ω
Figure 2. Domain ω .
∞
ν X (ξ ) ≤ ∫ ∫ x&1 f X X& ( x1 , x&1 ) f X
1 1
&
2X2
(ξ − x1 , x& 2 )dω1dx1 +
−∞ ω1
∞
∫ ∫ x& 2 f X X& ( x1 , x&1 ) f X
1 1 2X2
& (ξ − x1 , x& 2 )dω 2 dx1
−∞ ω 2
∞ ∞ ∞
= ∫ ∫ ∫ x&1 f X X& ( x, x&1 ) f X
1 1 2X2
& (ξ − x, x& 2 )dx& 2 dx&1dx + (7)
−∞ 0 − ∞
∞ ∞ ∞
∫ ∫ ∫ x& 2 f X X& ( x, x&1 ) f X
1 1 2X2
& (ξ − x, x& 2 )dx& 2 dx&1dx
−∞ −∞ 0
∞ ∞
= ∫ ν X ( x) f X (ξ − x)dx + ∫ ν X (ξ − x) f X ( x)dx
1 2 2 1
−∞ −∞
181
Note 10: Load combinations
∞
ν X (ξ ) ≤ ∫ ν X ( x1 ) f X
1 2 +X3
(ξ − x1 )dx1 +
−∞
∞
∫ ν X 2 ( x2 ) f X1 + X 3 (ξ − x 2 )dx 2 + (8)
−∞
∞
∫ ν X 3 ( x3 ) f X1 + X 2 (ξ − x3 )dx3 +
−∞
Example 1
Consider two independent stochastic processes with joint density functions:
1
f X X& ( x2 , x& 2 ) = 24
(
9 − x22 (1 − x& 2 )
2
) for (x 2 , x& 2 ) ∈ [−3;3] × [−1;1]
2 2
0 else
1
2 0.5 cos( x1 )(1 − x&1 )dx&1 for ( x1 , x&1 ) ∈ [−0,5π ;0.5π ] × [−1;1]
= 0∫
0 else
0.5 cos( x1 ) for x1 ∈ [−0,5π ;0.5π ]
=
0 else
∞
f X ( x2 ) = ∫ f X
2
&
2X2
( x2 , x& 2 )dx& 2
−∞
1
= 36
9 − x22 ( ) for x2 ∈ [−3;3]
0 else
∞
ν X (ξ ) = ∫ x&1 f X X& (ξ , x&1 )dx&1
1 1 1
0
1
cos(ξ ) for ξ ∈ [−0.5π ;0.5π ]
= 12
0 else
182
Note 10: Load combinations
∞
ν X (ξ ) = ∫ x& 2 f X
2 2X2
& (ξ , x& 2 )dx& 2
0
1
= 288
9 −ξ 2 ( ) for ξ ∈ [−3;3]
0 else
π π
An upper bound for the out-crossing intensity for X (t ) = X 1 (t ) + X 2 (t ) for 3 − ≤ξ ≤ 3+ can
2 2
then be determined from, see integration limits in figure 4:
∞ ∞
ν X (ξ ) ≤ ∫ ν X ( x) f X (ξ − x)dx + ∫ ν X (ξ − x) f X ( x)dx
1 2 2 1
−∞ −∞
( ) ( )
π /2 1 1 π /2 1 1
= ∫ cos( x ) 9 − (ξ − x) 2 dx + ∫ 9 − (ξ − x) 2 cos( x )dx
ξ −3 12 36 ξ −3 288 2
=−
7 cos(ξ − 3) 7 sin (ξ − 3) 7 π 2 − 4ξπ + 4ξ 2
− −
( − 44)
, 3−
π
≤ξ ≤ 3+
π
288 864 6912 2 2
ν X (2) ≤ −
7 cos(2 − 3) 7 sin (2 − 3) 7 π 2 − 8π − 28
− − = 0.0375
( )
288 864 6912
* * * * * * * *
183
Note 10: Load combinations
Example 2
Consider two independent, stationary Gaussian stochastic processes X 1 (t ) and X 2 (t ) with statisti-
cal parameters:
µX =1 1
σX =3 σ X& = 2
1 1
µX = 0
2
σX =2
2
σ X& = 1
2
X (t ) = 2 X 1 (t ) + X 2 (t )
µ X = 2 ⋅1 + 0 = 2
σ X = 2 2 ⋅ 32 + 2 2 = 6.325
σ X& = 2 2 ⋅ 2 2 + 12 = 4.123
4.123 1 12 − 2 2
ν (12) =
+
exp − =0.0297
2π 6.325 2 6.325
Note, that in this case with Gaussian processes and a linear combination of the processes, the out-
crossing intensity can be determined without evaluation of the integrals in (7).
* * * * * * * *
184
Note 10: Load combinations
In the Ferry Borges-Castanheta load model, [3] it is assumed that the variable load processes in the
load combination problem can be approximated by ‘square-wave’ processes, see figure 5 for the
case of two load processes X 1 (t ) and X 2 (t ) . The following description is based on ISO [2].
185
Note 10: Load combinations
For each process (load) three different stochastic variables are defined:
1. An arbitrary point in time stochastic variable for load no j : X *j with distribution function
FX j ( x j ) .
2. The maximum value for load no j : X j ,max,T during the reference time T with the distribution
rj
function: FX j , max,T
( x j ) = [ FX ( x j )] .
j
• for X 1 the combination load is denoted X 1C and is equal to the arbitrary point in time vari-
able X 1* . The distribution function thus becomes: FX1C ( x1 ) = FX1 ( x1 ) .
These stochastic variables and quantiles of them can be used in reliability analyses and in defining
design values if the partial safety factor method is used, see next two sections. The distribution
functions are illustrated in figure 6.
Note that the combination load X C is not the same as the characteristic value of a stochastic vari-
able, xc .
186
Note 10: Load combinations
As can be seen from section 2 it is generally very difficult to obtain an exact / analytic expression
for the out-crossing intensity (and the probability of failure). Therefore a number of approximations
have been suggested. In practice (and in codes) the Turkstra rule, [4] is often used in load combina-
tion problems to estimate the probability of failure and to establish design values to be checked in a
level 1 safety approach.
Z1 = max{X 1 (t )} + X 2 (t * ) + ... + X r (t * )
T
Z 2 = X 1 (t * ) + max{X 2 (t )} + ... + X r (t * )
T (10)
M
Z r = X 1 (t * ) + X 2 (t * ) + ... + max{X r (t )}
T
where t * is an arbitrary point in time and max{ X j (t )} is the maximum value of X j (t ) in the time
T
interval [0, T ] . X max,T is then approximated by X max,T = max{Z1 , Z 2 ,..., Z r }. This stochastic variable
can be used in evaluation of the reliability of the structure considered.
If the partial safety factor method is used the Turkstra rule is often applied together with the Ferry
Borges-Castanheta load model.
It is assumed that the load effect (e.g. a cross sectional force) can be written as a function of two (or
more loads):
S = S ( X1, X 2 ) (11)
In a level 1 code, two combinations for the design load effect corresponding to these two combina-
tions are generally considered:
1. S d ,1 = S ( xd 1,1 , x d 1, 2 ) = S (γ 1 xc1 , γ 2ψ 2 xc 2 )
2. S d , 2 = S ( xd 2,1 , x d 2, 2 ) = S (γ 1ψ 1 xc1 , γ 2 xc 2 )
where γ 1 , γ 2 are partial safety factors, ψ 1 , ψ 2 are load combination factors and xc1 , xc 2 are charac-
teristic values (usually 98 % quantiles in the distribution function for the annual maximum load
which is FX ,max,T =1year ( x j ), j = 1,2 .
j
187
Note 10: Load combinations
Note that this definition of load combination factors are different from that used in the Danish codes
[5] where design values for non-dominating variable loads are calculated as ψ j xc j i.e. without the
partial safety factor γ j .
S d = max{S d ,1 , S d , 2 } (12)
In the following it is shown how the above model can be used in a reliability analysis and how it
can be used to determine load combination factors.
the load effect modeled by (11). The result is a reliability index β1 and design-point values: x1*,1 and
x1*, 2 .
Load combination 2: Similarly, X 1 and X 2 are modeled as stochastic variables with distribution
functions FX1C ( x1 ) and FX 2 ,max,T ( x2 ) . A reliability analysis is performed using a given limit state
function and the load effect modeled by (11). The result is a reliability index β 2 and design-point
values: x2*,1 and x2*, 2 .
The two cases can be considered as two failure modes and a series system reliability index can be
estimated.
Note that if there is a non-zero probability that the load X j is equal to zero during some of the time
intervals τ j , then the distribution function for X j has to include this probability, i.e.
It is assumed that a target reliability index β t and characteristic values are given. The partial safety
factors and load combination factors can now be estimated.
Load combination factors calculated from reliability analyses: reliability analyses are performed
such that for both load combinations the reliability index becomes equal to β t . From the corre-
188
Note 10: Load combinations
sponding design-point values the partial safety factors and load combination factors can be esti-
mated:
Load combination factors calculated from design-value format: In Eurocodes, Basis of Design [5]
and ISO [2] it is recommended that design values of dominating variable loads are determined
from:
FX j , max,T
( xdj , j ) = Φ (−α S β t ) (15)
where α S =-0.7 is recommended. If the Turkstra rule is applied design values for non-dominating
loads can be determined from:
where the factor 0.4 is chosen as a reasonable value. Partial safety factors and load combination
factors can then be obtained from (note the similarity with (14):
Example 3
It is assumed that X 1 and X 2 are modeled as Gumbel distributed stochastic variables with statisti-
cal parameters:
Characteristic values are assumed to be 98% quantiles in the distribution functions for annual
maximum loads. The partial safety factors are determined from:
xdj , j FX
−1
Φ 0.7 β t 1−V j
π
6
( ( ))
0.577 + ln − ln Φ 0.7 β t { ( [ ( )])}
γj = = =j , max,T
(19)
xcj FX−1 (0.98) 6
1−Vj {0.577 + ln(− ln[0.98])}
j , max,T
π
The load combination factors are determined from:
x
ψ 1 = 1C = −1
( (
FX−1 Φ 0.28 β t )) = F −1
X 1, max,T (Φ(0.28β ) ) = 1 − V
t r1
1
π
6
{0.577 + ln(− ln[Φ(0.28β )]) + ln r }
t
1
(Φ(0.7 β )) F (Φ(0.7 β ))
1C
−1
{0.577 + ln(− ln[Φ(0.7 β )])}
t t
x d 1,1 FX 6
1, max,T X 1, max,T
1 − V1 t
π
(20)
189
Note 10: Load combinations
ψ2 =
x2C
=
FX−1 Φ 0.28 β t
=
( (
FX−1 )) (Φ(0.28β ) ) = 1 − V
t r1
2
6
π
{0.577 + ln(− ln[Φ(0.28β )]) + ln r }
t
1
( ( )) (Φ(0.7 β ))
2C 2 , max,T
π
(21)
* * * * * * * *
Example 4. Load combination factors for imposed loads and wind load
The reference period is assumed to be T = 1 year and the corresponding target reliability index
β t =4.3.
If the design-value format is used (18)-(21) give the partial safety factors and load combination fac-
tors in table 2.
In load combination 1 (imposed load dominating) Q1 and Q2 have the distribution functions
FQ 1, max
(
(q1 ) and FQ2 C (q 2 ) = FQ2 , max (q 2 ) ) 1 / r1
. In load combination 2 (wind load dominating) Q1 and Q2
have the distribution functions FQ (q1 ) = FQ 1C
( 1, max
(q1 ) )
1 / r1
and FQ 2 , max
(q 2 ) . FQ1, max and FQ2 , max refer to
the statistical data in table 1. The design values for the two load combinations corresponding to
β t =4.3 and the associated load combination factors obtained by (14) are shown in table 3. The load
combination factors are rather large, whereas compared to the results in [6], example 1 and 2 the
partial safety factors are small. Therefore, it is of interest to calculate modified load combination
factors where the partial safety factors in [6], example 1 are used, see table 4 below. The load com-
bination factors in table 4 are larger but comparable to those obtained by the design-value method,
see table 2.
190
Note 10: Load combinations
* * * * * * * *
The reference period is assumed to be T = 1 year and the corresponding target reliability index
β t =4.3. The statistical data are shown in table 5. Note that snow load only occurs in the period
November – March. The same limit state function as in example 1 is used.
In load combination 1 (snow load dominating) Q1 and Q2 have the distribution functions FQ 1, max
(q1 )
and FQ (q2 ) = FQ
2C
( 2 , max
(q2 ) )
1 / 24
.
In load combination 2 (wind load dominating) Q1 and Q2 have the distribution functions
FQ (q1 ) = FQ
1C
( 1, max
(q1 ) ) 1 / 10
and FQ ( 2 , max
(q2 ) )
150 / 360
.
The design values for the two load combinations corresponding to β t =4.3 and the associated load
combination factors obtained by (14) are shown in table 6. Note that the load combination factor for
snow is much larger than the load combination factor for wind.
191
Note 10: Load combinations
* * * * * * * *
6 References
[1] Thoft-Christensen, P. and M.J. Baker: Structural Reliability Theory and Its Applications.
Springer Verlag, 1982.
[3] Ferry Borges, J. & M. Castanheta: Structural Safety. 2nd edition. Laboratorio Nacional de En-
genharia Civil, Lisbon, 1972.
[4] Turkstra, C.J. & H.O. Madsen: Load combinations in codified structural design. J. Struct.
Div., ASCE, Vol. 106, No. St. 12, 1980.
[5] DS409:1998 – Code of practice for the safety of structures. Dansk Standard, 1999.
[6] Sørensen, J.D.: Structural reliability: Level 1 approaches. Note, Institute of Building Technol-
ogy and Structural Engineering, Aalborg University, April 2000.
192
Note 11: Example: Fatigue / Reliability-based inspection planning
1 Introduction
This note gives an introduction to the main steps in a probabilistic fatigue analysis and inspection
planning for welded joints. As example tubular joints in fixed offshore platforms of the steel jacket
type are considered, but the probabilistic modeling is general and can be used for other types of
structures. Initially the fatigue loading is described, here as an example wave loading. Next stress
analysis is considered. Based on a spectral analysis the stress spectra for critical points (hot spots) in
the joint can then be determined using an influence matrix approach. From the stress spectra stress
ranges can be determined and the number of stress cycles can be estimated, e.g. by the Rainflow
counting method.
Two models for the fatigue strength are described, namely the classical SN approach and a fracture
mechanics approach where the size of the crack is compared with a critical crack length, e.g. the
thickness of the tubular member. The basic steps in a reliability analysis with respect to fatigue and
in reliability-based inspection planning is described and illustrated. Part of this note is based on EFP
[1].
2 Fatigue loading
The most important load for fatigue failure of welded offshore structures is wave loading. Current is
insignificant because the time variation is very slow compared with wave loading. The fatigue load
due to wind excitation can contribute by 10-15 % of the total fatigue load but usually it is of minor
importance. In this section we therefore concentrate on wave loading.
The statistical properties of sea waves are most often modeled using so-called short-term sea states.
The duration of a sea state is normally taken as 3 hours. Within each sea state the wave elevation is
assumed modeled by a stationary, Gaussian stochastic process {η (t )} . The wave elevation η (t ) is
assumed Normal distributed with expected value µη = 0 and standard deviation σ η . The auto-
spectrum of {η (t )} can be modeled by a number of different spectra, e.g.
• Pierson-Moskowitz
• JONSWAP
4π 3 H S2 1
4
Sηη (ω ) = exp − 16π
3
(1)
TZ ω 5 T ω
Z
193
Note 11: Example: Fatigue / Reliability-based inspection planning
where ω is the cyclical frequency, H S is the significant wave height and TZ is the zero up-crossing
period. The parameters H S and TZ are constant within each sea state. In figure 1 a typical wave
spectrum is shown.
TZ [sec] 0-1 1-2 2-3 3-4 4-5 5-6 6-7 7-8 8-9 9-10 10-11 11-12
H S [m]
10.5-11.0 +
10.0-10.5 + +
9.5-10.0 1 +
9.0-9.5 1 +
8.5-9.0 + 1 +
8.0-8.5 1 1 1 +
7.5-8.0 1 2 1 +
7.0-7.5 2 2 1 +
6.5-7.0 + 2 3 1 +
6.0-6.5 1 5 4 2 1
5.5-6.0 3 7 5 1 1
5.0-5.5 1 9 11 5 2 1
4.5-5.0 3 18 13 7 2 1
4.0-4.5 1 11 25 10 7 2 1
3.5-4.0 3 22 30 8 4 1
3.0-3.5 1 20 35 25 5 3 2
2.5-3.0 3 51 42 18 3 1 1
2.0-2.5 15 70 30 15 3 1
1.5-2.0 5 71 58 20 10 2 1
1.0-1.5 23 91 38 10 3 1
0.5-1.0 7 32 16 6 3 1
0.0-0.5 1 1 2 1
Table 1. Representative scatter diagram for central North Sea. Numbers are probabilities in parts per
thousand. +: 0.0<probability<0.0005.
Long-term observations of the sea are usually performed by observing the sea surface for 20 min-
utes every third hour. For each observation H S and TZ are estimated. The relative number of pairs
of H S and TZ can be represented in so-called scatter diagrams, see table 1. Based on the observa-
tions it is also possible to fit the long-term distribution functions for H S , e.g. by a Weibull distribu-
tion
h − H γ
FH (h) = 1 − exp − 0
, h ≥ H0 (2)
S
H C − H 0
194
Note 11: Example: Fatigue / Reliability-based inspection planning
From table 1 it is seen that H S and TZ are dependent. Based on the observations a long-term distri-
bution function for TZ given H S can be fitted, for example by a two-parameter Weibull distribu-
tion
t
k 2 ( hS )
FT |H
(t Z | hS ) = 1 − exp − Z
(3)
Z S
k1 (hS )
where k1 (hS ) and k 2 ( hS ) are functions of hS . In [2] the following models are obtained based on
data from the Northern North Sea ( hS in meters).
Generally the distribution functions for H S and TZ are dependent on the wave direction Θ . If eight
directions (N, NE, E, SE, S, SW, W, NW) with probabilities of occurrence PΘ , i = 1,2,...,8 are i
γi
h−H
FH (h, Θ i ) = 1 − exp − 0 ,i
h ≥ H 0,i , i = 1,2,...,8 (6)
S H −H
C ,i 0 ,i
The parameters in (3)-(5) can be considered independent of the direction. Together with the pa-
rameters in (6) for the 8 directions the probabilities PΘ , i = 1,2,...,8 for waves in the eight direc-
i
Measurements of the directional characteristics of the wave elevation show a variation of both the
mean direction and a spread with frequency. The spreading of the waves can result in a significant
reduction in the wave loading. The directional spectra are assumed modeled by
Γ(s + 1)
Ψ ( Θ) =
1
[cos(0.5(Θ − Θ ))]2 s (8)
2 π Γ(s + 0.5)
195
Note 11: Example: Fatigue / Reliability-based inspection planning
3 Stress Analysis
Above it is described how the wave load can be described by the spectral density Sηη (ω ) and the
distribution functions FH (h) and FT |H (t Z | hS ) . Next, it is of interest to calculate the spectral den-
S Z S
sity Sσσ (ω ) for the stresses in a critical hot spot. One way to calculate Sσσ (ω ) is to perform a sto-
chastic response analysis to find the cross-spectral density functions S S S (ω ) for the cross-sectional
k l
forces in a given structural element and then calculate Sσσ (ω ) as described below. Details of such
an analysis can be found in e.g. Langen & Sigbjørnson [3].
S S S (ω ) = H Fη (ω ) H F* η (ω ) Sηη (ω )
k l i j
(9)
where H Fη (ω ) is the transfer function from wave elevation to cross-sectional force no i and * de-
i
In order to illustrate the procedure the K-joint in figure 2 is considered. The cross-sectional forces
on the joint can be determined using a beam model of the structure. These forces will be in equilib-
rium. A local stress analysis of the joint can therefore be performed by fixing one of the cross-
sections (see figure 2) and applying the cross-sectional forces from the beam model as external
loads on the joint. The cross-sections where the forces are determined should be located in some
distance from the joint in order to be able to apply the cross-sectional loads as distributed line loads
on a shell element model of the joint, i.e. the stress distribution is unaffected by the joint.
The local fatigue inducing hot spot stress σ in a critical point, namely the principal stress perpen-
dicular to the crack, see figure 3 is estimated by
196
Note 11: Example: Fatigue / Reliability-based inspection planning
N
σ = ∑α k Sk (10)
k =1
where N is the number of cross-sectional forces applied as loads to the joint (=18 in figure 4 where
each cross-section has 6 degrees of freedom). α k is the coefficient of influence giving the stress in
the critical point for a unit load S k .
Based on the cross spectral densities for the cross-sectional forces the auto spectral density of the
fatigue hot spot stress σ can be determined from, see also figure 4
N N N N
Sσσ (ω ) = ∑ ∑ α kα l S S S (ω ) = ∑ ∑ α kα l H F η (ω ) H F* η (ω ) Sηη (ω )
k l k l
(11)
k =1l =1 k =1l =1
197
Note 11: Example: Fatigue / Reliability-based inspection planning
For computational reasons it is more convenient to calculate the cross spectral densities S S S (ω ) of
k l
the load effects first. Next, when the auto-spectral density of a stress is required this can be calcu-
lated using (11). If the result of the spectral analysis had been the auto-spectral density of the fa-
tigue stress, a new spectral analysis would be required whenever the fatigue stress in a new location
is needed. This would be rather unfortunate, as a full spectral analysis is very time consuming.
The location of the most critical hot spots is usually not known in advance. Therefore 8 (or 12)
points located as shown in figure 2 are investigated. The auto-spectral density functions are deter-
mined for each location and a fatigue analysis is performed as described in the following sections.
This is done for the 8 points in the brace and for the corresponding 8 points on the intersection in
the chord.
198
Note 11: Example: Fatigue / Reliability-based inspection planning
4 Fatigue strength
4.1 SN approach
Figure 5. Experimental SN results for circular K-joint. From EC3 background document [7].
Assuming that the fatigue damage is accumulated linearly in an interaction free manner the damage
accumulation law attributed to Palmgren [5] and Miner [6] can be applied. Failure occurs when the
accumulated damage exceeds 1, i.e. the failure criteria is
ni
∑ ≥1 (12)
i N (∆σ i )
where ni is the number of stress cycles at a particular stress range level ∆σ i and N (∆σ i ) is the
number of constant amplitude stress cycles at that stress range level which leads to failure. The
summation in (12) is over the number of different stress range levels. N (∆σ i ) is usually deter-
mined on the basis of experiments and therefore has a random character, see figure 5.
is assumed and the material parameters m and K are fitted to experimentally obtained data.
It is seen from (12) and (13) that a limit state function can be written as
1
g =∆− ∑ ni ∆σ i
m
(14)
K i
199
Note 11: Example: Fatigue / Reliability-based inspection planning
or
1
g =∆− ∑ ∆σ i
m
(15)
K i
if the summation is over all individual stress range cycles. ∆ is a stochastic variable modeling the
uncertainty related to application of the Miner rule for linear accumulation of damage from each
stress cycle. Usually, ∆ has an expected value equal to 1.
The most simple and generally applicable crack growth equation is due to Paris & Erdogan [8]:
da
= C (∆K ) m , ∆K > 0 (16)
dn
where a is the crack size (depth), n is the number of stress cycles, ∆K is the stress intensity factor
range in a stress cycle. C and m are material constants.
da
According to (16) a plot of log versus log(∆K ) should be linear but a typical plot obtained ex-
dn
perimentally would be more like the one shown in figure 6.
The agreement between (16) and experiments is seen to be reasonable in region II (almost linear)
whereas (16) overestimates the crack growth rate in region I and underestimates the crack growth
rate in region III. ∆K th is a threshold stress intensity range below which the crack will not grow.
K IC is the value of the stress intensity factor at which the crack becomes unstable and brittle frac-
ture takes place. The stress intensity factor (SIF) can be shown to have the form:
∆K = Y (a)∆σ πa (17)
where
200
Note 11: Example: Fatigue / Reliability-based inspection planning
∆K is a factor accounting for a redistribution of the hot spot fatigue stresses. The reason for this
redistribution is the influence of the crack itself and other local geometry boundary conditions.
da
dn
( )
= CY (a) m (∆σ ) m πa
m
(18)
( 2−m ) / 2 2 − m m
2 /( 2 − m )
a + Cπ ∆σ n
m/2
for m ≠ 2
a(n) = 0 2 (19)
a exp(Cπ∆σ 2 n ) for m = 2
0
For offshore joints it is generally not sufficient to model cracks as being one-dimensional. This is
because both the crack depth a and the crack length c influence the geometry function Y = Y (a, c) .
201
Note 11: Example: Fatigue / Reliability-based inspection planning
Consider a flat plate with a semi-elliptical surface crack under tension or bending fatigue loads, see
figure 7. The depth of the crack is a and its length is 2c , while the thickness of the plate is t .
Shang-Xian [9] assumed that the growth rates at the deepest point A and the end point B of the
crack follow independently the Paris & Erdogan equations:
da
= C a (∆K a ) m with a(0) = a0 (20)
dn
dc
= Cc (∆K c ) m with c(0) = c0 (21)
dn
The variation in the three-dimensional stress field is accounted for by the constants Ca and Cc ,
while ∆K a and ∆K c denote respectively the ranges of the stress intensity factor at the deepest point
A and the summit B, see figure 7.
From the two coupled equations, the differential equation of the shape change is derived as
m
dc Cc ∆K c
= with c(a0 ) = c0 (22)
da Ca ∆K a
together with
dn 1
= with n(a0 ) = 0 (23)
da C a (∆K a ) m
The statistics of the amplitude or stress-ranges and the corresponding number of stress-ranges in a
given time internal must be obtained in order to assess the fatigue damage.
If the fracture mechanics approach (see section 4.2) is used, crack growth is governed by Paris' law.
In order to illustrate how fatigue cracks can be counted a one-dimensional crack model is used in
the following. Integration of (18) gives for constant stress-range amplitudes ∆σ
ac
da
∫ = C ∆σ m n (24)
a0 (Y (a) πa )
m
where a0 and ac are the initial and the final (critical) crack size, respectively. Y (a) is the geometry
function, ∆σ is the constant amplitude stress-range and n is the number of stress cycles during the
considered time interval [0, T ] . A generalization to variable stress-range amplitudes can be obtained
by using instead of ∆σ m the equivalent stress range to power m , E[∆σ m ]
1 n
E[∆σ m ] = ∑ ∆σ i
m
(25)
n i =1
202
Note 11: Example: Fatigue / Reliability-based inspection planning
neglecting any sequence effects. ∆σ m is treated as a stochastic variable and E[.] denotes the expec-
tation operation.
If SN curves (see section 4.1) are used to model the fatigue strength it is seen from (15) that also in
this case the damage accumulation is governed by (25).
For offshore structures the expectation (25) must be performed for a given sea state because the
state term statistics of the stresses are conditional on the sea states. Therefore an expectation over
all sea states must be performed:
∑ ∆σ i = ∑ P( H S ,i , TZ ,i , Θ i ) (∆σ ( H S ,i , TZ ,i , Θ i ) )
1 n nS
E[∆σ m ] = m m
(26)
n i =1 i =1
where n S is the number of non-zero boxes in the scatter diagram and P ( H S ,i , TZ ,i , Θ i ) is the prob-
ability having ( H S ,i , TZ ,i , Θ i ) corresponding to the i th box. ∆σ ( H S ,i , TZ ,i , Θ i ) is the corresponding
stress range.
Figure 8 shows three different sample curves of stress histories. The first case corresponds to con-
stant amplitude loading, where the stress-ranges are the same for all stress cycles. The second case
corresponds to a stationary ideal narrow band Gaussian process. Again the stress cycle is easily de-
fined in terms of the stress process between two constitutive up-crossings of the mean value. The
third case, which is the more general case with broad banded stress variation, is not quite as obvi-
ous. In this case one has to use counting methods.
In section 4.3.1 narrow band stress spectra are considered. Next broad band spectra are considered.
In section 4.3.2 and 4.3.3 it is shown how the range counting and the Rainflow counting methods
can be used to estimate E[∆σ m ] and the expected number of stress cycles n .
For a narrow-banded Gaussian process, the stress-ranges are Rayleigh distributed. The mean value
in (25) is then
( )(
E[∆σ m ] = 2 2
m
m0 ) Γ(1 + m / 2)
m
(27)
203
Note 11: Example: Fatigue / Reliability-based inspection planning
where m0 is the zero'th spectral moment of the stress spectrum, Sσσ (ω ) ( m0 is the standard de-
viation of the stress process). Generally, the i th moment is defined by (if Sσσ (ω ) is a two-sided
spectrum):
∞
mi = 2 ∫ ω i Sσσ (ω )dω (28)
0
The number of stress cycles n in the time interval [0, T ] is estimated from
1 m2
n = ν 0T = T (29)
2π m0
where ν 0 is the mean zero-crossing rate and m2 is given by (28). Note that m2 is equal to the
standard deviation for the derivative process {σ& (t )} .
In the range counting method a half stress cycle is defined of the difference between successive
local extremes, see figure 9. The range counting method uses only local information. Therefore in-
formation on larger stress cycles can be lost if small stress reversals are superimposed on the larger
stress cycles. The method gives a lower bound on the fatigue damage.
The mean number of stress cycles in a time interval [0, T ] is equal to the mean number of local
maxima in the time interval
1 m4
n = ν mT = T (30)
2π m2
where m4 is given by (28).
204
Note 11: Example: Fatigue / Reliability-based inspection planning
Using a double envelope process to model the stress process it can be shown that, see [10]
( )(
E[∆σ m ] = α m 2 2
m
m0 ) Γ(1 + m / 2)
m
(31)
m2 ν
α= = 0 (32)
m0 m4 ν m
(31) deviates from (27) by the factor α m . In the limit where the process is narrow banded ( α =1)
(31) and (27) are identical.
Rainflow counting is considered to give the most accurate predictions of the fatigue life when com-
pared to actual fatigue life results. Rainflow counting is widely used. Material hysterises loops are
sometimes used to justify its use. Rainflow counting is illustrated in figure 10 where the largest cy-
cles are extracted first and the smaller cycles are considered to be superimposed on the larger cy-
cles, see [11] and [12].
The Rainflow counting method counts the number of stress cycles by converting a realization of the
stress process {σ (t )} to a point process of peaks and troughs as shown in figure 11. The peaks are
identified by even numbers and the troughs by odd numbers. The following rules are imposed on
"rain dropping on the roofs", so that cycles and half cycles are defined, see Wirshing & Sheheta
[13].
205
Note 11: Example: Fatigue / Reliability-based inspection planning
more positive than that at the start of the rain path under consideration (e.g. in figure 11, path
[2-3], path [4-5] and path [6-7]).
3. If the rain flowing down a roof intercepts a flow from the previous path, the present path is
stopped, (e.g. in figure 4.6, path [3-3a], path [5-5a], etc.)
4. A new path is not started until the path under consideration is stopped.
Half-cycles of trough-originated range magnitudes hi are projected distances on the x-axis (e.g. in
figure 11, [1-8], [3-3a], [5-5a] etc.). If the realization of {σ (t )} is sufficiently long, any trough-
originated half-cycle will be followed by another peak originated half-cycle of the same range.
Figure 11. Illustration of Rainflow cycle counting applied to sample of {σ (t )} (from [13]).
Due to the complexity of the Rainflow algorithm it is very difficult to derive a density function f ∆σ
for the stress ranges and to estimate the number of stress cycles.
However, based on extensive computer simulations, Dirlik [14] has derived empirical expressions
for f ∆σ :
D s 2
+ D2 s exp − s D3 s
2
1 s
f ∆σ ( s ) = 1 exp − +
8m R 2 2 m exp − (33)
2 m0 Q 2Q m 2 m R2 8m
0 0 0 0 0
where
D1 = 2( β − α 2 ) /(1 + α 2 ) (34)
D3 = 1 − D1 − D2 (37)
Q = 1.25(α − D3 − D2 R) / D1 (38)
206
Note 11: Example: Fatigue / Reliability-based inspection planning
m1 m2
β= (39)
m0 m4
n = ν mT (40)
From (19) and (24) it is seen that if Y (a ) =1 and n = νT then failure can be modeled by the limit
state function
2−m
g = a0( 2−m ) / 2 + Cπ m / 2 ∆σ mνT − ac( 2−m ) / 2
2
T [years] β α1 α2 α3 α4
2.5 5.50 0.49 -0.03 0.77 0.41
5.0 4.38 0.52 -0.02 0.74 0.43
7.5 3.78 0.53 -0.02 0.72 0.45
10.0 3.32 0.54 -0.02 0.70 0.46
12.5 2.99 0.55 -0.03 0.69 0.46
15.0 2.72 0.56 -0.02 0.68 0.47
17.5 2.50 0.57 -0.02 0.67 0.48
20.0 2.31 0.58 -0.02 0.66 0.48
22.5 2.15 0.58 -0.01 0.66 0.48
25.0 2.00 0.59 -0.01 0.66 0.49
Table 3. Reliability index and sensitivity coefficients α1 ,…, α 4 at different times.
The results of a reliability analysis is shown in table 3. It is seen that the reliability index β de-
creases from 5.50 to 2.00 when T increases from 2.5 year to 25 years. Further it is also seen, that
a0 , ln C and ∆σ are the most important stochastic variables in this example.
207
Note 11: Example: Fatigue / Reliability-based inspection planning
On the basis of a deterministic analysis the joints and elements given in table 4 have been identified
as being critical with respect to fatigue. The numbers given in table 4 refers to the numbering given
in the report "Structure Identification, Platforms" [15].
Joint Element
154050 111404
152050 111424
253050 52501
253050 52006
Table 4. Critical joints and elements.
Both the SN-approach and the fracture mechanics approach are based on the evaluation of an
equivalent stress range. The equivalent stress range is determined on the basis of the assumption
that the stress history in a given sea-state can be modeled as a narrow band Gaussian process. This
implies that for a given sea state the stress range, ∆σ , follows a Rayleigh distribution
∆σ 2
F∆σ (∆σ H S , TZ , Θ ) = 1 − exp − (42)
8m0 (H S , TZ , Θ )
where H S , TZ and Θ denotes the significant wave height, the wave period and the direction, m0 is
the 0th order spectral moment of the stress history.
(
F∆σ (∆σ ) = ∑ P(H S ,i , TZ ,i , Θ i )F∆σ ∆σ H S ,i , TZ ,i , Θ i )
nS
(43)
i =1
where nS denotes the number of considered sea states and P(H S ,i , TZ ,i , Θ i ) is the probability of hav-
ing significant wave height H S ,i , wave period TZ ,i and direction Θ i in sea state no. i .
∆σ λ
F∆σ (∆σ ) = 1 − exp −
k
(44)
The Weibull distribution may be fitted so that it corresponds to the long-term stress distribution in
the 95 and 99 % quantiles. Using this approximation it can be found that
208
Note 11: Example: Fatigue / Reliability-based inspection planning
[ ] m
E ∆σ m = k m Γ1 + (45)
λ
m 2 ,i
N = T ∑ P(H S ,i , TZ ,i , Θ i )
nS
(46)
i =1 m 0 ,i
where m0,i = m0,i (H S ,i , TZ ,i , Θ i ) and m2,i = m2,i (H S ,i , TZ ,i , Θ i ) are spectral moments for sea state no.
i.
The spectral moments are determined by a spectral analysis, see [1]. The spectral moments deter-
mined are deterministic constants. However, due to the uncertainty related to the transfer function
and the uncertainty related to the stress concentration factors the spectral moments must be modeled
as stochastic variables. The uncertainty related to the spectral moments is introduced by
m0 = U H U SCF
2
m0∗ (47)
where U H and U SCF are stochastic variables describing the random variation of the transfer func-
tion and the stress concentration factors and where m0∗ is the spectral moment determined by the
spectral analysis where the mean values of the transfer function and stress concentration factors
have been used.
The uncertainty related to the spectral moments is taken into account by determining the joint dis-
tribution of the parameters k and λ in the Weibull model.
The joint distribution of these variables is determined by simulation. For a given outcome of the
variables U H and U SCF the parameters k and λ are determined. On the basis of a large number of
simulations (here 1000 simulations), the joint distribution of the variables is determined.
209
Note 11: Example: Fatigue / Reliability-based inspection planning
For the joints and elements considered here the distributions of the parameters k and λ are given
in table 5.
5.2 SN-Approach
The limit state function for reliability analysis using the SN-approach is written, see section 4.1:
1 n
g =∆−D= ∑ ∆σ i
m
(48)
K i =1
where ∆ is the damage limit corresponding to failure and D models the accumulated damage dur-
ing the considered time interval, see (15).
In table 6 and 7 the input parameters to the limit state function are given.
In figure 12 the reliability index as a function of time is given for the four different joints and ele-
ments for the hot spot 270°. In figure 13 the squared alpha values of the three stochastic variables
are shown.
6
Joint 152050, element 111424, hotspot 270
0
0 5 10 15 20 25 30
Time [years]
210
Note 11: Example: Fatigue / Reliability-based inspection planning
∆
ln(K)
k
The sensitivity measures are only shown for one of the joints, but they are very similar to the three
other joints. The values do not change within time. It is seen that approximately 60% of the uncer-
tainty originate from the mean value of the Weibull distribution modeling the long term stress, k,
30% from the logarithm of the material parameter, ln(K), and 10% from the damage limit corre-
sponding to failure, ∆.
A reliability level corresponding to a safety index of 3.7 will assure the required safety. From figure
12 it is seen that this level is exceeded after approximately 5 years. Therefore it will be necessary to
perform an inspection after 5 years in order to check the condition of the joint.
a
Y (a ) = 1.08 − 0.7 (49)
t
where t is the thickness of the plate.
This expression is valid for a crack in a plate. For cracks in a tubular joint the stress intensity factors
can be determined by multiplying by a correction factor determined by Smith and Hurworth [17] as
a a
M k (a ) = 1.0 + 1.24 exp − 22.1 + 3.17 exp − 357 (50)
t t
211
Note 11: Example: Fatigue / Reliability-based inspection planning
The probabilistic models for m and C A have been chosen such that the difference between the
distribution functions for the time to fatigue failure using the fracture mechanical approach and the
SN-approach is minimized.
In figure 14 the results of the reliability analyses are shown. It is seen that the results of the analyses
using the SN-approach are virtually identical to the results obtained on the basis of the fracture me-
6.0
Computation based on SN-curve
5.0 Computation based on fracture mechanics
4.0
Reliability index
3.0
2.0
1.0
0.0
0.0 5.0 10.0 15.0 20.0 25.0 30.0
Time [years]
chanics model. The analysis has been carried out for joint 253050, element 52006, hot spot 270o,
i.e. the element with the largest mean stress range.
Figure 14. Results of reliability analysis.
212
Note 11: Example: Fatigue / Reliability-based inspection planning
6 Inspection Planning
The inspection planning is performed on the basis of the assumption that no cracks are detected by
the inspections. Further, an inspection is carried out when the reliability with respect to the consid-
ered event decreases below a given acceptable limit. For a more detailed description of the method
see [18].
Each of these inspection methods are associated with a given probability of detection depending on
the size of the crack and depending on whether the inspection is carried out above or below water.
As mentioned it is assumed that no crack is detected by the inspection, i.e. the crack is smaller than
the detectable crack size. This event can be modeled by the following limit state function
h = a − ad (51)
The distribution of the size of the detectable crack or the so-called POD-curve (Probability of De-
tection) depends on the inspection method. The distribution of the detectable crack can be modelled
by an upward bounded Exponential distribution
Fa (a d ) = P0 (1 − exp(− a / λ ))
d
(52)
The inspection planning is carried out using an inspection method with the following parameters.
λ = 2,67 mm -1
P0 = 1.0
In figure 15 the results of the inspection planning are shown. On the basis of the probabilistic mod-
eling of the fatigue crack growth described an inspection plan has been determined in accordance
with the methodology described in [19]. It is found that 4 inspections performed at year 4, 8, 15 and
23 will ensure as sufficient safety for the considered joint throughout its anticipated service life of
30 years.
213
Note 11: Example: Fatigue / Reliability-based inspection planning
6.0
5.0
4.0
Reliability index
3.0
2.0
1.0
1. Inspection 2. Inspection 3. Inspection 4. Inspection
0.0
0.0 5.0 10.0 15.0 20.0 25.0 30.0
Time [years]
7 References
[1] Engelund, S. & L. Mohr: Reliability-Based Inspection Planning. EFP’97 report, COWI, 1999.
[4] Madsen, H.O.: Probability-based fatigue inspection planning. MTD Publication 92/100, The
Marine Technology Directorate Limited, UK, 1992.
[2] Thoft-Christensen, P. & M.J. Baker: Structural reliability theory and its applications. Springer
verlag 1984.
[9] Shang-Xiang, W.: Shape change of surface crack during fatigue growth. Engineering fracture
mechanics, Vol. 22, No. 5, 1985, pp. 897-913.
[5] Palmgren, A.: Die Lebensdauer von Kugellagern. Zeitschrift der Vereines Deutches Ingenieu-
re, Vol. 68, No. 14, 1924, pp. 339-341.
[6] Miner, M.A.: Cumulative damage in fatigue. J. Applied Mech., ASME, Vol. 12, 1945, pp.
A159-A164.
[7] EC3 Background document 9.03: Background information on fatigue design rules for hollow
sections. Part 1: Classification method – statistical evaluation. 1991.
[8] Paris, P.C. & F. Erdogan: A critical analysis of cracks. in Fracture toughness testing and its
applications, ASTM STP 381, 1965, pp. 30-82.
[11] Barltrop, N.D.P. & A.J. Adams: Dynamics of fixed marine structures. Butterworth - Heine-
mann, London, 1991.
214
Note 11: Example: Fatigue / Reliability-based inspection planning
[12] Longuet-Higgins, M.S.:On the joint distribution of periods and amplitudes of sea waves.
Journal of Geophysical Research, Vol. 80, 1975,pp. 2688-2694.
[10] Madsen, H. O., Krenk, S. & Lind, N. C.: Methods of Structural Safety. Prentice Hall, Engle-
wood Cliffs, 1986.
[13] Wirshing, P.H. & A.M. Shebeta: Fatigue under wide band random stresses using the rain-flow
method. J. Eng. Materials and Technology, ASME, July 1977, pp. 205-211.
[14] Dirlik, T.: Application of computers in fatigue analysis. University of Warwick, PhD Thesis,
1985.
[16] Kirkemo, F. and H.O. Madsen: Probability Based Fatigue Inspection Planning and Repair.
A.S. Veritas Research Report No.: 87-2004, 1988.
[17] Smith, I.J. and S.J. Hurworth: Probabilistic fracture mechanic evaluation of fatigue from weld
defects in butt welded joints. Fitness for purpose validation of welded constructions, London
17-19, 1981.
[18] Friis Hansen, P. and H.O. Madsen: Stochastic Fatigue Crack Growth Analysis Program. Dan-
ish Engineering Academy, 1990.
[19] Sørensen, J.D., S. Engelund, A. Bloch & M.H. Faber: Formulation of Decision Problems –
Offshore structures. EFP’97, Aalborg University and COWI, 1999.
215
Note 11: Example: Fatigue / Reliability-based inspection planning
216
Note 12: Reliability updating
1 Introduction
When new information it can be used to update the stochastic models and the estimates of the reli-
ability (probability of failure). In this note it is described how this updating can be performed when
the new information consists of
1. Observation of events described by one or more stochastic variables. The observation is mod-
eled by an event margin and the failure event by a safety margin. Updated / conditional prob-
abilities of failure can then be obtained, see section 2.
2. Samples / measurements of a stochastic variable X . Updating can in this case be performed
using Bayesian statistics, see section 3.
H = h(X) (1)
is introduced. The event function h corresponds to the limit state function. The actual observations
are considered as realizations (samples) of the stochastic variable H. This type of information can
e.g. be
• Inspection events such as measurements of the chloride content in concrete structures or meas-
urements of crack sizes in steel structures exposed to fatigue loading. The event margin can in-
clude the uncertainty related to the measurement.
• Proof loading where a well defined load is applied to a structure and the level of damage (usu-
ally no damage is observed) is observed.
• Repair events where a certain type of repair or maintenance has been performed.
• No-failure events where the ‘simple’ observation that the structure / component considered is
well-functioning after some time in use.
P( g (X) ≤ 0 ∩ h(X) ≤ 0)
PFU = P( g (X) ≤ 0 h(X) ≤ 0) = (2)
P(h(X) ≤ 0)
217
Note 12: Reliability updating
where M = g (X) is the safety margin related to the limit state function g (x) and X = ( X 1 ,..., X n )
are stochastic variables. In (2) it is used that the probability of an event A given an event B (denoted
P ( A ∩ B)
P( A B) ) is equal to . It is seen that P( g (X) ≤ 0 ∩ h(X) ≤ 0) is the probability of a parallel
P ( B)
system with two elements. (2) can be evaluated by simulation or FORM/SORM methods, see Mad-
sen et al. [1].
Other observations can be modeled by equality events {H = 0} , i.e. it is observed that the observed
quantity is equal to some limit. In this case the updated probability of failure can be estimated by,
see Madsen et al. [1], Madsen [2] and Schall and Rackwitz [3]
∂
P(g(X) ≤ 0 ∩ h(X) ≤ z)
P(g(X) ≤ 0 ∩ h(X) = 0) ∂z z=0
PF = P(g(X) ≤ 0 h(X) = 0) =
U
= (3)
P(h(X) = 0) ∂
P(h(X) ≤ z)
∂z z=0
Equation (3) can also be evaluated by FORM/SORM methods and can easily be generalized if more
than one event is observed. In most software packages for reliability analysis efficient algorithms
are available for solving this problem.
Example 2.1
Consider the following limit state function modeling crack growth of a crack with initial length a0
to the critical length ac . m and C are the parameters in Paris’ law. ν is the number of stress cycles
per year with stress range ∆σ . T is the time in years.
2−m
g = a0( 2−m ) / 2 + Cπ m / 2 ∆σ mνT − ac( 2−m ) / 2
2
β (T ) = −Φ −1 (P( g ( X, T ) ≤ 0) − P( g ( X, T − 1) ≤ 0) )
is shown in figure 1. It is seen that the annual reliability index becomes smaller than β = 4 at year
14.
218
Note 12: Reliability updating
It is now assumed that the fatigue critical detail is inspected at time T1 =14 years with the result that
no cracks were detected. The reliability of the inspection is described by the following POD-curve
where a d models the smallest detectable crack size:
a
POD(a) = Fa (a) = 1 − exp −
d
b
b is a parameter modeling the expected value of a d . It is assumed that b =1 mm. The inspection
event (no cracks detected) is modeled by the event margin
2−m
h = a d( 2−m ) / 2 − a0( 2−m ) / 2 − Cπ m / 2 ∆σ mνT ≤ 0
2
β U (T ) = −Φ −1 (P( g ( X, T ) ≤ 0) h( X) ≤ 0) − P( g ( X, T − 1) ≤ 0) h( X) ≤ 0) ) , T > T1
6.0
5.5
5.0
4.5
β 4.0
3.5
3.0
2.5
2.0
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
time
Figure 1. Reliability indices as function of time in years. Full line: β (T ) and broken line: β U (T ) .
219
Note 12: Reliability updating
If possible Bayesian techniques should be used for parameter estimation because Bayesian estima-
tion gives an estimate of the statistical uncertainty related to the estimated parameters and because
updating of the model when new information becomes available is easy.
If observations of one (or more) of the stochastic variables X are available, the probabilistic model
can be updated and thereby also the probability of failure. Consider a stochastic variable X with
density function f X ( x) . If q denotes a vector of parameters defining the distribution for X , the
density function of the stochastic variable X can be written
f X ( x, q ) (4)
If X is normally distributed then q could contain the mean and the standard deviation of X .
If the parameters q are uncertain then f X ( x, q) can be considered as a conditional density function
: f X ( x Q) and q denotes a realization of Q. The initial density function for the parameters Q is
denoted f Q′ ( q) and is denoted the prior density function.
It is assumed that n observations (realizations) of the stochastic variable X are available making up
a sample xˆ = ( xˆ1 , xˆ 2 ,..., xˆ n ) . The realizations are assumed to be independent. The updated density
function f Q′′ (q xˆ ) of the uncertain parameters Q given the realizations is denoted the posterior den-
sity function and is given by, see textbook on Bayesian statistics, e.g. Box & Tiao [4] and Lindley
[5]
f N (xˆ q) f Q′ (q )
f Q′′ (q xˆ ) = (5)
∫ f N (xˆ q) f Q′ (q )dq
N
where f N (xˆ q) = ∏ f X ( xˆ i q) is the probability density at the given observations assuming that the
i =1
distribution parameters are q . The integration in (5) is over all possible values of q .
The updated density function of the stochastic variable X given the realization x̂ is denoted the
predictive density function and is defined by,
Given the distribution function for the stochastic variable X , the prior distribution is often chosen
such that the posterior distribution will be of the same type as the prior distribution (a so-called con-
jugated prior). In the literature a number of prior, posterior and predictive distribution functions can
be found, see e.g. Raiffa & Schlaifer [6], Aitchison & Dunsmore [7] and Rackwitz & Schrupp [8].
Some of these are shown in appendix A.
By use of the Bayesian method presented here, both the physical uncertainty related to the consid-
ered variable as well as the statistical uncertainty related to the model parameters can be quantified.
However, as mentioned the probabilistic model must also be formulated such that the measurement
uncertainty and the model uncertainty are taken into account.
220
Note 12: Reliability updating
Due to the ability of Bayesian statistics to incorporate engineering judgement and experience into
the statistical modeling, different reassessment engineers may reach different results due to the use
of different statistical models. This may obviously be a serious obstacle for the use of such methods.
Also in order to avoid this obstacle it is necessary to define and agree on a common basis for such
analyses thus ensuring that reliability based reassessment analyses are performed on a consistent
and comparable basis
Example 3.1
Consider a Normal distributed stochastic variable X with expected value µ and known standard
deviation σ =3. Prior knowledge on µ is modeled by a Normal distributed prior distribution with
expected value µ ' =10 and standard deviation σ ' =4.
Further it is assumed that n =5 observations are available: x̂ =(11.3, 12.4, 9.5, 10.7, 11.1) imply-
ing that x =10.9. The posterior distribution the becomes Normal with expected value µ ' ' and stan-
dard deviation σ ' ' , see (A-4)-(A-5):
The predictive distribution also becomes Normal distributed with expected value µ ' ' and standard
deviation σ ' ' ' , see (A-8)
nx σ ' 2 + µ ' σ 2
µ '' = 2 2
=10.9 and σ ''' = σ '' 2 +σ 2 =3.26
nσ ' +σ
In figure 2 are shown the prior and posterior distributions for the expected value µ .
0.25
0.20
0.15
0.10
0.05
0.00
0.00 4.00 8.00 12.00 16.00 20.00
µ
Figure 2. Prior (full line) and posterior (broken line) distributions for the expected value µ .
221
Note 12: Reliability updating
0.10
0.08
0.05
0.03
0.00
0.00 4.00 8.00 12.00 16.00 20.00
x
Figure 3. Predictive distribution for the stochastic variable X .
4 References
[1] Madsen, H.O. , Krenk, S. & Lind, N.C., Methods of Structural Safety, Prentice Hall, Engle-
wood Cliffs, N.J. 1986.
[2] Madsen, H.O. : Model Updating in Reliability Theory. Proc. ICASP5, 1987, pp. 564-577.
[3] Schall, G. & R. Rackwitz : On the integration of multinormal densities over equalities. Be-
richte zur Zuverlassigkeitstheorie der Bauwerke. Technical University Munich, Heft 83, 1988.
[4] Box, G. E. P., Tiao, G. C., Bayesian Inference in Statistical Analysis, Wiley, New York,
1992.
[5] Lindley, D.V. Introduction to Probability and Statistics from a Bayesian Viewpoint, Vol. 1+2,
Cambridge Univ. Press, Cambridge 1976.
[6] Raiffa, H. & Schlaifer, R., Applied Statistical Decision Theory, MIT Press, Cambridge, 1968.
[7] Aitchison, J. & I.R. Dunsmore : Statistical Prediction Analysis. Cambridge University Press,
Cambridge 1975.
[8] Rackwitz, R. & K. Schrupp : Quality Control, Proof Testing and Structural Reliability. Struc-
tural Safety, Vol. 2, 1985, pp. 239-244.
222
Note 12: Reliability updating
1 1 x − µ 2
f X ( x µ ,σ ) = f N ( x µ ,σ ) = exp − (A-1)
σ 2π 2 σ
1 1 µ − µ' 2
f µ' ( µ ) = f N ( µ µ ' , σ ') = exp − (A-2)
σ ' 2π 2 σ '
where
nx σ ' 2 + µ ' σ 2
µ '' = (A-4)
nσ '2 +σ 2
σ '2 σ 2
σ '' 2 = (A-5)
nσ ' 2 +σ 2
1 n
x= ∑ x$i (A-6)
n i =1
where
223
Note 12: Reliability updating
1 1 x − µ 2
f X ( x µ ,σ ) = f N ( x µ ,σ ) = exp − (A-9)
σ 2π 2 σ
( )
1 2 (ν ' +2 )/ 2
ν ' ω ' /σ
2 2 1 ν ' ω '
f σ' (σ ) = f iγ 2 (σ ω ' , ν ' ) = exp − (A-10)
Γ ((ν '+2) / 2) 1 2 σ2
(ν ' ω ')
1/ 2
2
where
1 ν
2
ω= ∑ ( x$ i − µ ) (A-14)
ν i =1
where
∞
Γ( a) = ∫ exp( − t ) t a −1 dt (A-16)
0
224
Note 12: Reliability updating
∞
Γ(a , b) = ∫ exp( − t )t b −1dt (A-17)
0
1 1 x − µ 2
f X ( x µ ,σ ) = f N ( x µ ,σ ) = exp − (A-19)
σ 2π 2 σ
σ
f µ' ,σ ( µ , σ ) = f N µ ' x ' , f iγ 2 (σ υ ' , ν ') (A-20)
n'
σ
f µ'',σ ( µ , σ x$ ) = f N µ ' x ' ' , f iγ 2 (σ υ ' ' , ν ' ') (A-21)
n' '
where
( )
υ ′ ′ = ν ' υ '+ n' x ' 2 +νυ + nx 2 − n' ' x ' ' 2 / ν ' ' (A-24)
225
Note 12: Reliability updating
1 n
x= ∑ x$i (A-26)
n i =1
1 n 2
υ= ∑ ( x$ i − x ' ' ) (A-27)
n − 1 i =1
µ − x' '
Fµ'' σ ( µ σ , x$ ) = Φ (A-29)
σ / n' '
n' '+1
f X ( x xˆ ) = f S ( x x ' ' ,υ ' ' ,ν ' ') (A-30)
n' '
1 ln x − µ 2
1
f X ( x µ ,σ ) = exp − (A-31)
xσ 2π 2 σ
The Lognormal stochastic variable X can be treated in the same way as a normal distributed vari-
able because Y = ln X is normal distributed with standard deviation
σ 2
σ Y = ln + 1 (A-32)
µ2
1
µ Y = ln µ − σ Y2 (A-33)
2
226
Note 12: Reliability updating
f u' (u) = exp(n' αu − t ' exp(αu) )αt ' ( n') / Γ (n' ) (A-36)
f u'' (u x$ ) = exp(n' ' αu − t ' ' exp(αu))αt ' ' ( n'') / Γ ( n' ' ) (A-37)
where
n
t = ∑ exp(−αxi ) (A-40)
i =1
u
Fu'' ( u x$ ) = ∫ exp(n' ' αz − t ' ' exp(αz) )αt ' ' ( n '') / Γ (n' ' ) dz (A-41)
−∞
− ( n '' +1)
1 α
f X ( x x$ ) = n' ' 1 + exp( −αx ) exp( −αx ) (A-42)
t' ' t' '
− n ''
1
FX ( x x$ ) = 1 + exp( −αx ) (A-43)
t' '
227
Note 12: Reliability updating
x −ε k
F X ( x ε , u, k ) = 1 − exp − (A-44)
u
k −1 x −ε k
k x −ε
f X ( x ε , u, k ) = exp − (A-45)
u u u
f u' (u) =
1
Γ (n'−1 / k )
( )
u − ( n ' k ) exp − u − k t ' kt ' ( n '−1/ k ) (A-46)
f u'' (u x$ ) =
1
Γ ( n ' ' −1 / k )
( )
u − ( n '' k ) exp − u − k t ' ' kt ' ' ( n ''−1/ k ) (A-47)
where
n
t = ∑ ( x$i − ε ) k (A-50)
i =1
Fu'' (u x$ ) = 1−
(
Γ n ' ' −1 / k , t ' ' u − k ) (A-51)
Γ ( n ' ' −1 / k )
228
Note 12: Reliability updating
− ( n '' −1/ k )
(x − ε)k
F X ( x x$ ) = 1. 1 + (A-53)
t' '
f X ( x λ ) = λ exp( − λ ) (A-54)
where
n
ν ' ' = ν '+ ∑ x$i (A-58)
i =1
229
Note 12: Reliability updating
where
B ( a , b) = Γ ( a ) Γ ( b) / Γ ( a + b ) (A-62)
x
B( x , a , b) = ∫ t ( a −1) (1 − t ) (b −1) dt (A-63)
0
230