Monte Carlo in Python¶

This tutorial builds a Monte Carlo simulation of the bolt fatigue damage with the random samplers of pyflange.stats. It does the same as the sample command, but in a script, so you can change any distribution, add random variables, or post-process the results as you like.

Download the full script: montecarlo.py

Samplers¶

A sampler is a Python generator: an object that returns a new random value every time you call next() on it.

from pyflange.stats import norm_sampler

preload = norm_sampler(2876e3, 0.03)    # normal, mean 2876 kN, CoV 3%
next(preload)                           # a random preload, e.g. 2.91e6
next(preload)                           # another one, e.g. 2.84e6

PyFlange has generic samplers (norm_sampler, lognorm_sampler, weibull_sampler, beta_sampler, sampler for any scipy distribution) and samplers of PyFlange objects, which combine other samplers:

Sampler	Yields
`standard_gap_sampler`	Random `Gap` (AMD1 gap height distribution)
`standard_gap_sampler_ED2`	Random `Gap` (ED2 distributions)
`standard_PolynomialLFlangeSegment_sampler`	Random L-flange segments
`standard_markov_matrix_sampler`	Markov matrices with random load ranges
`standard_bolt_fatigue_curve_sampler`	Random bolt S-N curves
`fatigue_case_sampler`	Random `BoltFatigueAnalysis` objects

You can write your own samplers too: any generator that yields objects of the right type works.

1. Reproducibility¶

The samplers use the global numpy random generator. Seed it to get the same sample at every run:

np.random.seed(1)       # makes the simulation reproducible

2. Flange segment sampler¶

import pyflange.stats as stats
from pyflange.bolts import StandardMetricBolt, RoundNut, standard_hole_diameter

D = 8000*mm
n_bolts = 156
bolt = StandardMetricBolt("M80", "10.9", shank_length=160*mm,
                          shank_diameter_ratio="pitch", stud=True)

fseg_sampler = stats.standard_PolynomialLFlangeSegment_sampler(
    a = 150*mm, b = 122*mm, s = 54*mm, t = 172*mm, R = D/2,
    central_angle = 2*pi / n_bolts,
    Zg = -18044*kN / n_bolts,
    bolt = bolt,
    Do = standard_hole_diameter("M80"),
    washer = None,
    nut = RoundNut("M80"),

    # Random preload before losses: mean 2876 kN, CoV 3%
    preload_sampler = stats.norm_sampler(2876*kN, 0.03),

    # Random preload losses (ref. [5], section 6.4.1)
    ppl_sampler = stats.norm_sampler(0.03, 1/3),        # plastic loss: mean 3%, std 1%
    dft_sampler = stats.norm_sampler(250*um, 0.12),     # coating: mean 250 um, std 30 um

    # Random gap: AMD1 gap height distribution with random angle and shape
    gap_sampler = stats.standard_gap_sampler(
        flange_diameter = D,
        flange_flatness_tolerance = 0.0014,
        gap_angle_sampler = stats.lognorm_sampler(100*deg, 1.0),
        gap_shape_factor_sampler = stats.norm_sampler(1.0, 0.15)),

    # Random inside flange inclination: LN(0.1 deg, CoV 0.5)
    tilt_sampler = stats.lognorm_sampler(0.1*deg, 0.5),

    settlement_factor = 0.5)

For each realization, the sampler draws a gap and a flange inclination, then draws the preload and its losses for half of the bolts within the gap and averages them (ref. [5], 4.4). The yielded segment has the averaged preload after losses as Fv, and PPL = DFT = 0.

The distributions are those of the stochastic model of ref. [5], Table 8, except the gap angle, which the reference does not give.

3. Loads, S-N curve and Miner's rule¶

from pyflange.fatigue import MarkovMatrix

design_mkvm = MarkovMatrix(
    range  = np.array([ 25,  50, 100, 150, 200, 250,  300, 375, 500, 625]) * 1e6,
    mean   = np.array([100, 100, 100, 100, 100, 100,  100, 100, 100, 100]) * 1e6,
    cycles = np.array([3e8, 1e8, 3e7, 8e6, 2e6, 5e5, 1.5e5, 4e4, 8e3, 1e3]),
    duration = 30)

# Load ranges scaled by X_FLS ~ LN(1.0, 0.12)
mkvm_sampler = stats.standard_markov_matrix_sampler(design_mkvm, 1.0, 0.12)

# Bolt S-N curve with DS_ref = X_FAT ~ N(62 MPa, 0.10)
sn_sampler = stats.standard_bolt_fatigue_curve_sampler(bolt.nominal_diameter, 1.0, 0.10, 62*MPa)

# Miner's rule uncertainty: allowable damage ~ LN(1.0, 0.30)
allowable_damage_sampler = stats.lognorm_sampler(1.0, 0.30)

4. Fatigue case sampler¶

The four samplers are combined into a sampler of complete fatigue analyses:

fcase_sampler = stats.fatigue_case_sampler(fseg_sampler,
                                           mkvm_sampler,
                                           sn_sampler,
                                           allowable_damage_sampler)

fcase = next(fcase_sampler)         # one random BoltFatigueAnalysis
print(f"One realization: gap angle {fcase.fseg.gap.angle/deg:.0f} deg, "
      f"gap height {fcase.fseg.gap.height/mm:.2f} mm, damage {fcase.damage:.3g}")

One realization: gap angle 46 deg, gap height 0.19 mm, damage 0.0397

Each realization is a BoltFatigueAnalysis: fcase.fseg is the random segment, fcase.flange_mkvm the random Markov matrix, fcase.fatigue_curve the random S-N curve, fcase.allowable_damage the random allowable damage, and fcase.damage and fcase.fatigue_life the results.

5. Generate the sample¶

Collect the values of interest in a pandas table:

import pandas as pd

rows = []
for i in range(SAMPLE_SIZE):
    fcase = next(fcase_sampler)
    rows.append({
        "GapAngle":        fcase.fseg.gap.angle,
        "GapHeight":       fcase.fseg.gap.height,
        "BoltPreload":     fcase.fseg.bolt_preload,
        "Damage":          fcase.damage,
        "AllowableDamage": fcase.allowable_damage,
        "FatigueLife":     fcase.fatigue_life,
    })
    if (i+1) % 500 == 0:
        print(f"... {i+1} realizations")

sample = pd.DataFrame(rows)

6. Results¶

utilization = sample.Damage / sample.AllowableDamage
pf = np.mean(utilization > 1)       # fraction of failed realizations

print(sample.describe(percentiles=[0.05, 0.5, 0.95]).transpose().to_string())
print(f"Probability of failure over {design_mkvm.duration} years: {pf:.2e}")

sample.to_csv("montecarlo-sample.csv")      # for further processing, e.g. in Excel

                  count          mean           std           min            5%           50%           95%           max
GapAngle         2000.0  1.236226e+00      0.723843  6.096265e-02  2.983545e-01  1.087546e+00  2.675680e+00  3.140577e+00
GapHeight        2000.0  1.143269e-03      0.001238  9.070261e-08  4.384932e-05  7.442454e-04  3.596016e-03  1.128285e-02
BoltPreload      2000.0  2.732521e+06  27989.062233  2.600771e+06  2.688598e+06  2.732288e+06  2.777531e+06  2.916772e+06
Damage           2000.0  1.432734e-01      0.231241  1.589216e-03  1.232016e-02  7.530805e-02  4.771642e-01  4.232725e+00
AllowableDamage  2000.0  9.917610e-01      0.294008  2.884288e-01  5.983693e-01  9.423795e-01  1.534146e+00  2.654757e+00
FatigueLife      2000.0  7.505789e+02   1285.993299  8.399722e+00  5.808260e+01  3.838157e+02  2.471174e+03  2.216393e+04
Probability of failure over 30 years: 1.55e-02

Sample size

With 2000 realizations the estimate of a small probability is very uncertain: 31 failures out of 2000 give a coefficient of variation of about 18% on the estimate. For a target probability of failure of 10^-4, use at least 10⁶ realizations, or a variance reduction technique.

This simulation does not apply a model uncertainty factor (stress multiplication factor) to the stress ranges. Pass a constant or a sampler as the SMF argument of fatigue_case_sampler to include it.