# Axial stressed beam

## Description

This example is a simple beam, restrained at one side and stressed by a traction load F at the other side :

The geometry is supposed to be deterministic: The diameter D is fixed to D=20 mm.

It is considered that failure occurs when the beam plastifies, i.e. when the axial stress gets bigger than the yield stress :

Therefore, the state limit G used here is :

Two independent random variables R and S are considered:

R (the strength):

### Theoretical results

This two dimensional stochastic problem can be solved by calculating directly the failure probability:

If R and S are independant, then:

The numerical application gives :

Pf = 0.0292

## Model definition

• Run Python
```# Main OpenTURNS module
from openturns import *
# optional matplotlib viewer
from openturns.viewer import View
```
• Constant definition:
```# dimension
dim = 2
```
• Random generator initialization:
```RandomGenerator.SetSeed(0)
```
• Analytical model definition:
```limitState = NumericalMathFunction(['R', 'F'], ['G'], ['R-F/(_pi*100.0)'])
```
• Test of the limit state function:
```x = [300., 75000.]
print 'x=', x
print 'G(x)=', limitState(x)
```

Output:

```x= [300.0, 75000.0]
G(x)= [61.2676]
```
• Stochastic model definition:
```# Create a first marginal : LogNormal distribution 1D, parameterized by
# its mean and standard deviation
R_dist = LogNormal(300., 30., 0., LogNormal.MUSIGMA)
R_dist.setName('Yield strength')
R_dist.setDescription('R')
# Graphical output of the PDF
R_pdf = R_dist.drawPDF()
View(R_pdf).save('pdf_R.png')

# Create a second marginal : Normal distribution 1D
F_dist = Normal(75000., 5000.)
F_dist.setDescription('F')
# Graphical output of the PDF
F_pdf = F_dist.drawPDF()
View(F_pdf).save('pdf_F.png')

# Create a copula : IndependentCopula (no correlation)
aCopula = IndependentCopula(dim)
aCopula.setName('Independent copula')

# Instanciate one distribution object
myDistribution = ComposedDistribution([R_dist, F_dist], aCopula)
myDistribution.setName('myDist')

# We create a 'usual' RandomVector from the Distribution
vect = RandomVector(myDistribution)

# We create a composite random vector
G = RandomVector(limitState, vect)

# We create an Event from this RandomVector
myEvent = Event(G, Less(), 0.0)

```

Output:

## Using Monte Carlo simulations

• Resolution options:
```cv = 0.05
NbSim = int(1e5)

algoMC = MonteCarlo(myEvent)
algoMC.setMaximumOuterSampling(NbSim)
algoMC.setBlockSize(1)
algoMC.setMaximumCoefficientOfVariation(cv)
# For statistics about the algorithm
initialNumberOfCall = limitState.getEvaluationCallsNumber()

```
• Perform the analysis:
```algoMC.run()
```
• Results:
```result = algoMC.getResult()
probability = result.getProbabilityEstimate()
print 'MonteCarlo result=', result
print 'Number of executed iterations =', result.getOuterSampling()
print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall
print 'Pf = ', probability
print 'CV =', result.getCoefficientOfVariation()
MCcvgraph = algoMC.drawProbabilityConvergence()
View(MCcvgraph).save('montecarlo_convergence.png')
```

Output:

```MonteCarlo result= probabilityEstimate=3.029830e-02 varianceEstimate=2.294261e-06 standard deviation=1.51e-03 coefficient of variation=5.00e-02 confidenceLength(0.95)=5.94e-03 outerSampling=12806 blockSize=1
Number of executed iterations = 12806
Number of calls to the limit state = 12806
Pf =  0.030298297673
CV = 0.0499923113245
```

## Using FORM analysis

• Resolution options:
```myCobyla = Cobyla()
# Resolution options:
eps = 1e-3
myCobyla.setMaximumIterationsNumber(100)
myCobyla.setMaximumAbsoluteError(eps)
myCobyla.setMaximumRelativeError(eps)
myCobyla.setMaximumResidualError(eps)
myCobyla.setMaximumConstraintError(eps)

# For statistics about the algorithm
initialNumberOfCall = limitState.getEvaluationCallsNumber()

# We create a FORM algorithm
# The first parameter is a NearestPointAlgorithm
# The second parameter is an event
# The third parameter is a starting point for the design point research

algoFORM = FORM(myCobyla, myEvent, myDistribution.getMean())
```
• Perform the analysis:
```algoFORM.run()
```
• Results:
```result = algoFORM.getResult()
print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall
print 'Pf =', result.getEventProbability()

# Graphical result output
importanceFactorsGraph = result.drawImportanceFactors()
View(importanceFactorsGraph).save('importance_factors.png')
```

Output:

```Number of calls to the limit state = 98
Pf = 0.0299850318913
```

## Using Directional sampling

• Resolution options:
```cv = 0.05
NbSim = 100000

algoDS = DirectionalSampling(myEvent)
algoDS.setMaximumOuterSampling(NbSim)
algoDS.setBlockSize(1)
algoDS.setMaximumCoefficientOfVariation(cv)
# For statistics about the algorithm
initialNumberOfCall = limitState.getEvaluationCallsNumber()
```
• Perform the analysis:
```algoDS.run()
```
• Results:
```result = algoDS.getResult()
probability = result.getProbabilityEstimate()
print 'Number of executed iterations =', result.getOuterSampling()
print 'Number of calls to the limit state =', limitState.getEvaluationCallsNumber() - initialNumberOfCall
print 'Pf = ', probability
print 'CV =', result.getCoefficientOfVariation()
DScvgraph = algoDS.drawProbabilityConvergence()
View(DScvgraph).save('directionalsampling_convergence.png')
```

Output:

```Number of executed iterations = 481
Number of calls to the limit state = 8830
Pf =  0.0297305755073
CV = 0.0499171711488
```

## Script

Here is the script used for this example: wiki:ExampleAxialStressedBeam:AxialStressedBeam-v1.4.py