= 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 :

[[Image(cylindre.png)]]


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 :

[[Image(equation_1.png)]]

Therefore, the state limit ''G'' used here is :

[[Image(equation_2.png)]]


Two independent random variables ''R'' and ''S'' are considered:

''R'' (the strength):

[[Image(equation_3.png)]]

and ''S'' (the load):

[[Image(equation_4.png)]]


=== Stochastic model ===

[[Image(stochastic_model.png)]]

 

=== Theoritical results ===

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

[[Image(equation_5.png)]]


If ''R'' and ''S'' are independant, then:

[[Image(equation_6.png)]]


The numerical application gives :

'''''Pf'' = 0.0292'''



== Solving the problem with OpenTURNS' TUI - Release 0.9.2 to current ==

== Model definition ==

- Run Python

- Load OpenTURNS:
{{{
# Main OpenTURNS module
from openturns import *
# OpenTURNS viewer capability (since release 0.9.2)
# from openturns_viewer import ViewImage
# New module for the viewer (since release 0.13.0)
from openturns.viewer import ViewImage
}}}

- Constant definition:
{{{
stochasticDimension = 2
}}}

- Random generator initialization:
{{{
RandomGenerator().SetSeed(0)
}}}

- Analytical model definition:
{{{
# Analytical construction : Input
inputFunction = Description(stochasticDimension)
inputFunction[0] = "R"
inputFunction[1] = "F"

# Analytical construction : Output
outputFunction = Description(1)
outputFunction[0] = "G"

formulas = Description(outputFunction.getSize())
# Here, _pi is the built-in constant _pi=3.14159265....
formulas[0] = "R-F/(_pi*100.0)"

LimitState = NumericalMathFunction(inputFunction, outputFunction, formulas)
}}}

- Test of the limit state function:
{{{
x = NumericalPoint(stochasticDimension, 0)
x[0]=300.0
x[1]=75000.0
print "x=" , x
print "G(x)=" , LimitState(x)
}}}
Output:
{{{
x= class=NumericalPoint name=Unnamed dimension=2 implementation=class=NumericalPointImplementation name=Unnamed dimension=2 values=[300,75000]
G(x)= class=NumericalPoint name=Unnamed dimension=1 implementation=class=NumericalPointImplementation name=Unnamed dimension=1 values=[61.2676]
}}}

- Stochastic model definition:
{{{
# Mean
mean = NumericalPoint(stochasticDimension, 0.0)
mean[0] = 300.0
mean[1] = 75000.0

# Standard deviation
sigma = NumericalPoint(stochasticDimension, 0.0)
sigma[0] = 30.0
sigma[1] = 5000.0

# Additional parameters for the lognormal distribution:
component = Description(1)
BorneInf = 0.0

# Initialization of the distribution collection:
aCollection = DistributionCollection()

# Create a first marginal : LogNormal distribution 1D, parameterized by its mean and standard deviation
marginal = LogNormal(mean[0], sigma[0], BorneInf, LogNormal.MUSIGMA)
marginal.setName("Yield strength")
component[0] = "R"
marginal.setDescription(component)
# Graphical output of the PDF
pdfgraph1=marginal.drawPDF()
# If the directory name is omitted, the graph will be produced in the current directory.
# If the dimensions are omitted, the default is 400x300 up to release 0.10.0, and 640x480 for the later releases
pdfgraph1.draw("/tmp", "pdf_R", 640, 480)
# Visualize the graph
ViewImage(pdfgraph1.getBitmap())
# Fill the first marginal of aCollection
aCollection.add( Distribution(marginal, "Yield strength") )

# Create a second marginal : Normal distribution 1D
marginal = Normal(mean[1], sigma[1])
marginal.setName("Traction_load")
component[0] = "F"
marginal.setDescription(component)
# Graphical output of the PDF
pdfgraph2=marginal.drawPDF()
pdfgraph2.draw("/tmp", "pdf_F", 640, 480)
ViewImage(pdfgraph2.getBitmap())
# Fill the second marginal of aCollection
aCollection.add( Distribution(marginal, "Traction_load") )

# Create a copula : IndependentCopula (no correlation)
aCopula = IndependentCopula(aCollection.getSize())
aCopula.setName("Independent copula")

# Instanciate one distribution object
myDistribution = ComposedDistribution(aCollection, Copula(aCopula))
myDistribution.setName("myDist")

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

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

# We create an Event from this RandomVector
myEvent = Event(G, ComparisonOperator(Less()), 0.0, "Event 1")
}}}
Output:

[[Image(pdf_R.png)]]

[[Image(pdf_F.png)]]


== Using Monte Carlo simulations ==

- Resolution options:
{{{
cv = 0.05
NbSim = 100000
myAlgo = MonteCarlo(myEvent)
myAlgo.setMaximumOuterSampling(NbSim)
myAlgo.setBlockSize(1)
myAlgo.setMaximumCoefficientOfVariation(cv)
# For statistics about the algorithm
initialNumberOfCall = LimitState.getEvaluationCallsNumber()
}}}

- Perform the analysis:
{{{
myAlgo.run()
}}}

- Results:
{{{
result = myAlgo.getResult()
probability = result.getProbabilityEstimate()
print "MonteCarlo result=" , result.str()
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()
}}}
Output:
{{{
MonteCarlo result= probabilityEstimate=3.029830e-02 varianceEstimate=2.294261e-06 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
# Full digits for validation only
Pf =  0.030298297673
# Full digits for validation only
CV = 0.0499923113245
}}}
[[Image(montecarlo_convergence.png)]]

== Using FORM analysis ==

- Resolution options:
{{{
eps=1E-3

# We create a NearestPoint algorithm
myCobyla = Cobyla()
myCobyla.setSpecificParameters(CobylaSpecificParameters())
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

myAlgo = FORM(NearestPointAlgorithm(myCobyla), myEvent, mean)
}}}

- Perform the analysis:
{{{
myAlgo.run()
}}}

- Results:
{{{
result = myAlgo.getResult()
print "Number of calls to the limit state =" , LimitState.getEvaluationCallsNumber() - initialNumberOfCall
print "Pf =" , result.getEventProbability()

# Graphical result output
importanceFactorsGraph = result.drawImportanceFactors()
importanceFactorsGraph.draw("/tmp", "importance_factors", 640, 480)
ViewImage(importanceFactorsGraph.getBitmap())
}}}
Output:
{{{
Number of calls to the limit state = 98
# Full digits for validation only
Pf = 0.0299850318913
}}}
[[Image(importance_factors.png)]]


== Using Directional sampling ==

- Resolution options:
{{{
cv = 0.05
NbSim = 100000

myAlgo = DirectionalSampling(myEvent)
myAlgo.setMaximumOuterSampling(NbSim)
myAlgo.setBlockSize(1)
myAlgo.setMaximumCoefficientOfVariation(cv)

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

- Perform the analysis:
{{{
myAlgo.run()
}}}

- Results:
{{{
result = myAlgo.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()
DScvgraph.draw("/tmp", "directionalsampling_convergence", 640, 480)
ViewImage(DScvgraph.getBitmap())
}}}
Output:
{{{
Number of executed iterations = 481
Number of calls to the limit state = 11984
Pf =  0.0297306063238
CV = 0.0499171701744
}}}
[[Image(directionalsampling_convergence.png)]]

== Results ==

[[Image(results.png)]]

== Script ==

Here is the script used for this example:
[attachment:wiki:ExampleAxialStressedBeam:AxialStressedBeam.py AxialStressedBeam.py]

Here is the archive of the three calculation scripts:
[attachment:wiki:ExampleAxialStressedBeam:AxialStressedBeam.tar.gz Download archive]