wiki:ExampleAxialStressedBeam

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

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

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

Model definition

  • Run Python
  • Load OpenTURNS:
    # 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.setName('Traction_load')
    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

Results

Script

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

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

Last modified 3 years ago Last modified on 10/21/14 13:19:27

Attachments (19)

Download all attachments as: .zip