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

and S (the load):

Stochastic model

Theoritical 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

Solving the problem with OpenTURNS' TUI - Release 0.9.2

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

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

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

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

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

Results

Script

Here is the script used for this example: AxialStressedBeam.py

Here is the archive of the three calculation scripts: Download archive