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

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