Damage Tolerance Monte Carlo Life Estimation Suite (DAMOCLES)

Monte Carlo fatigue crack growth for safety-critical metallic structure, with the NASGRO equation, Newman-Raju stress intensity solutions, inspection planning against probability-of-detection curves, and the FAA AC 33.14-1 hard alpha rotor assessment. Every physics and statistics module is verified against an exact or published reference, and the package reproduces the FAA's own calibration test case inside the acceptance bands the Advisory Circular publishes for qualifying probabilistic rotor codes.

What it does

You describe a component as distributions instead of point values: initial flaw size, stress per cycle (constant amplitude, a rainflow counted block, or an explicitly ordered mission), fracture toughness, growth rate coefficient. For each study DAMOCLES produces:

• probability of failure over the service life, with an exact Clopper-Pearson confidence interval and the equivalent reliability index beta
• the same number with a scheduled inspection program applied: residual risk, risk reduction, and expected detections per part, given an NDT capability stated as a POD curve (lognormal a50/a90 or the digitized AC 33.14-1 curves)
• the probability-of-failure history over cycles, for setting intervals
• total and first-order Sobol indices on log life, so the input scatter that actually drives the answer is identified before money is spent characterising one that does not
• for titanium rotors: the zone-based AC 33.14-1 hard alpha assessment with the published anomaly exceedance curves and ultrasonic POD

The sampling engine offers plain Monte Carlo, Latin hypercube and scrambled Sobol sequences, plus importance sampling with likelihood reweighting for failure probabilities of 1e-5 and below, where plain sampling wastes nearly every draw.

Why the calibration test case matters

AC 33.14-1 Section 3 requires each engine manufacturer to calibrate its probabilistic tools against the Appendix 1 test case, a titanium ring disk with fully specified geometry, loading, material and anomaly distribution: "Test case results in the ranges from 1.27E-09 to 1.93E-09 (for the 'no inspection' case) and from 8.36E-10 to 1.53E-09 (for the 'with in-service inspection' case) are considered acceptable." That is the qualification gate the FAA wrote for codes like DARWIN. DAMOCLES runs the test case in closed form (the test case physics is deterministic given anomaly size, so the risk reduces to a one-dimensional integral over the exceedance curve) and again by Monte Carlo through its own sampling engine, and lands inside both bands. The run is a unit test, so the claim is re-verified on every commit.

One detail discovered on the way: the 2017 Change 1 tabulation of the anomaly distribution (Table A3-1) differs from the original 2001 figure A3-7 by factors of several in the 1,000 to 5,000 square mil region that dominates the test case risk. The acceptance band predates Change 1, so the calibration must run on the 2001 curve. DAMOCLES bundles both vintages, digitized from the official PDFs, and documents which one each assessment uses.

The data

Everything numeric ships with provenance:

• NASGRO 4.0 equation fits for 2024-T3 sheet, 7075-T651 plate, 2014-T6 sheet and 7475-T7351 plate, read from DOT/FAA/AR-05/15 (the FAA fatigue crack growth database report, Forman et al., 2005), with Walker fits where published. Stored in the printed US customary units and converted to SI on load.
• The generic Ti-6Al-4V average-air Paris fit specified by the AC 33.14-1 test case (MCIC-HB-01R via the AC).
• The complete AC 33.14-1 Change 1 Table A3-1 hard alpha exceedance table (184 rows, four billet/forging inspection levels), verbatim, plus the 2001 figure A3-7 curve digitized from the official PDF's vector art.
• The five AC 33.14-1 Appendix 5 POD curves (UT at three calibration levels, FPI, eddy current), digitized the same way; the AC publishes them only as graphs.

None of this is design data. The sources are public; check them, and substitute your own basis values for anything that matters.

How the numbers hold up

Every claim is enforced by the test suite on each commit (CI on Linux and Windows, Python 3.10 to 3.14). Measured results:

Check	Reference	Result
AC 33.14-1 calibration, no inspection	acceptance 1.27e-9 to 1.93e-9 events/cycle	1.44e-9, PASS
AC 33.14-1 calibration, with inspection at 10,000 cycles	acceptance 8.36e-10 to 1.53e-9	1.28e-9, PASS
Quadrature vs the Monte Carlo engine on the inspected test case	internal cross-check	agree within 0.5%
Crack growth life integration	closed-form Paris solution, two exponents	within 0.5%
Failure probability estimator	exact normal R-S problem, pof = 2.035e-4	CI brackets truth
Importance sampling	same reference, equal sample count	estimator CoV 5x below plain MC
Sobol sensitivity indices	Ishigami closed-form values	within 0.03 absolute
Tolerance factors kB, kA	published MMPDS/CMH-17 values (n=10: 2.355, 3.981)	match to 0.005
B-basis confidence	simulation across 2,000 datasets	95% coverage confirmed
NASGRO rates, 2024-T3	the FCGD spline fits of the same data	within a factor of 1.8
Newman-Raju limit cases	published F = 1.04 and 1.11 limits, NASA TM-85793	match within 1%
Rainflow counting	ASTM E1049-85 worked example	exact cycle table
Hoop stress of the AC test disk	AC quoted bore stress 572.4 MPa	574.2 MPa analytic

Quickstart

Needs Python 3.10+.

pip install damocles

The worked examples and the test suite live in the repository:

git clone https://github.com/Krataios14/DAMOCLES
cd DAMOCLES
pip install -e .[dev]
damocles examples/ti64_disk_bore.yaml --sensitivity --plot out/
damocles examples/ordered_overload_mission.yaml
python examples/ac3314_test_case.py
python -m pytest -q

The first command runs a corner crack study at a compressor disk bore on the AC Ti-6-4 material and prints:

  P(failure), no inspection : 9.295e-03
    95% CI                  : [8.879e-03, 9.725e-03]
    reliability index beta  : 2.35

  inspections at            : 4,000, 8,000, 12,000, 16,000
  P(failure), inspected     : 8.221e-05
    risk reduction          : 99.1%
    expected detections/part: 2.712e-01

  target P(failure)         : 1.0e-03  -> MEETS target

  variance drivers (total Sobol index on log-life):
    initial_flaw   total=0.554
    paris_c        total=0.371
    stress_range   total=0.077
    toughness      total=0.000

Read the bottom block before trusting the top one: here toughness scatter is irrelevant and flaw size dominates, so the productive next step is flaw characterisation, not more toughness coupons. A second worked example covers a 2024-T3 fuselage skin panel on the FCGD NASGRO fit with eddy current inspections.

Using it as a library

from damocles import (Lognormal, Normal, DamageToleranceStudy,
                      InspectionPlan, PODCurve, NewmanRajuCornerCrack,
                      material_growth_law)

study = DamageToleranceStudy(
    "disk bore",
    variables={
        "initial_flaw": Lognormal(mean=1.2e-4, cov=0.6),   # metres
        "stress_range": Normal(mean=480.0, std=25.0),      # MPa
        "toughness":    Normal(mean=64.5, std=4.0),        # MPa sqrt(m)
    },
    geometry=NewmanRajuCornerCrack(thickness=0.025, width=0.08),
    growth_law=material_growth_law("ti-6al-4v-ac33.14"),
    service_cycles=20_000,
    inspection_plan=InspectionPlan.at_interval(
        4_000, 20_000, PODCurve.from_a50_a90(0.4e-3, 1.0e-3)),
    n_samples=200_000, method="lhs", seed=42)
print(study.run(sensitivity=True).summary())

The pieces work standalone: estimate_pof takes any limit state over named random variables (it is a general structural reliability engine), grow and grow_spectrum integrate crack growth for arbitrary sample arrays, rainflow counts a stress history per ASTM E1049, b_basis and a_basis turn coupon data into allowables, and ac3314.run_test_case is the calibration run. CustomGeometry accepts any Y(a) callable for FE-derived solutions.

Ordered overload retardation

Load-interaction calculations need the original reversal order rather than merged rainflow classes. Build a SpectrumSequence from that ordered history and use the opt-in Willenborg integrator:

from damocles import (ParisLaw, SpectrumSequence, ThroughCrack,
                      grow_spectrum_retarded)

history = [300.0, 0.0, 180.0, 0.0, 180.0, 0.0]  # MPa
sequence = SpectrumSequence.from_history_ordered(history)

life = grow_spectrum_retarded(
    a0=[1.0e-3],
    sequence=sequence,
    geometry=ThroughCrack(),
    law=ParisLaw(c=1.0e-9, m=3.0),
    k_ic=50.0,       # MPa sqrt(m)
    s_yield=350.0,   # MPa, plane-stress plastic-zone estimate
)
print(life.cycles_to_failure[0])

grow_spectrum_retarded repeats the supplied sequence until fracture or a run-out limit is reached and returns life in cycles. grow_spectrum remains the faster default when load interaction is intentionally excluded.

Version 0.3 also routes the same ordered calculation through DamageToleranceStudy, YAML, and the CLI. Ordered mission records are immutable, bounded, and retained exactly as supplied: they are never sorted, merged, or binned. Each record represents one full cycle, or a residual half cycle from an ordered reversal history. For example:

variables:
  initial_flaw: {dist: deterministic, value: 1.0e-3}
  stress_scale: {dist: deterministic, value: 1.0}
  toughness: {dist: deterministic, value: 100.0}
spectrum:
  type: ordered
  cycles:
    - {delta_sigma: 1500.0, stress_ratio: 0.0}
    - {delta_sigma: 1000.0, stress_ratio: 0.0}
retardation:
  model: willenborg
  yield_strength: 1000.0
  max_cycles: 100.0
  max_blocks: 50
  max_work: 1000000
service_cycles: 100.0
analysis: {samples: 1000, method: lhs, seed: 7}

Set retardation: {model: none} to run the same ordered mission with no load-interaction correction. Omitting retardation has exactly that meaning. Omitted max_cycles and max_blocks are derived from service_cycles and the sequence length; explicit values must cover the whole service interval. The hard ceilings are 100,000 cycle records per block, 100,000,000 elapsed cycles, 1,000,000 blocks, and 100,000,000 prospective sample-cycle evaluations. A smaller max_work may be supplied. All limits are checked before Monte Carlo sampling or geometry evaluation begins.

Inspection times in an ordered mission are strictly increasing elapsed-cycle checkpoints before the end of service. Every service horizon, explicit cycle limit, and inspection time must coincide exactly with a cumulative ordered record endpoint. A full-cycle record cannot be stopped at its half-cycle; a residual half-cycle record may end there. Checkpoints see the post-record crack, and a sample still alive at that endpoint retains its actual crack size for POD evaluation—even if that crack would be critical under a higher load later in the mission. If the opening load of the next record causes failure at the same elapsed endpoint, that endpoint is overwritten with the new record's critical size: the sample has failure life equal to the checkpoint and is not available for detection there. Once failure occurs, that and later checkpoints carry the cycle-specific critical size at which it was detected. DAMOCLES rejects non-boundaries before sampling instead of partly executing or interpolating a record. Endpoint validation accepts only a factory-prepared integer-tick index registered to the exact sequence identity and its canonical structural fingerprint. At each configuration, inspection, or growth trust boundary it recomputes a private endpoint snapshot and rejects changed scalar fields, fingerprints, or tables before arithmetic. A checkpoint batch therefore costs O(N + M) for N records and M checkpoints; a study crossing k such boundaries costs O(kN + M), with constant-time lookup after each validation.

Interpreting the study output

Field	Meaning
pof, ci_low, ci_high	failure probability at end of service, exact binomial CI
reliability_index()	equivalent beta, for comparing designs
inspection.pof_inspected	residual risk with the inspection plan applied
inspection.risk_reduction	1 - inspected/unmitigated
inspection.mean_detections	expected detections per part over the life
sensitivity[name]["total"]	total Sobol index of that input on log life
pof_curve_cycles, pof_curve	the failure probability history
loading, life_unit, retardation	explicit loading and result-unit provenance
prospective_work	preflight upper bound on sample-cycle evaluations for an ordered mission

Run-outs (cracks below threshold that never grow) carry infinite life and are reported as such, not silently dropped.

Design allowables

damocles.allowables computes one-sided lower tolerance bounds from measured samples: B-basis (90% coverage at 95% confidence) and A-basis (99/95), using exact normal-theory factors from the noncentral t distribution, on raw or log-transformed data. Use it to turn your own coupon results into study inputs. The MMPDS/CMH-17 definitions apply; the factors match the published tables and the confidence level is verified by simulation in the test suite.

If you need the toughness input itself predicted from composition with finite-sample guarantees, that is a different problem and a different tool: FTQS (Fracture Toughness Qualification Suite). DAMOCLES consumes toughness as a distribution; FTQS is one defensible way to get one for alloys you have not tested.

Limitations, stated plainly

• LEFM with no load interaction by default: constant amplitude or repeating-block spectra, no retardation. Overload retardation (Willenborg model) is available as an opt-in via grow_spectrum_retarded; the existing grow_spectrum is unchanged.
• The Willenborg option uses the original plane-stress plastic-zone model and the historical zero-effective-R convention. It is empirical and should be supported by representative spectrum-test data before use in a safety-critical life assessment.
• Ordered mission studies are deliberately cycle-by-cycle and capped before execution. Sobol sensitivity is not offered on this expensive path; use explicitly bounded parameter sweeps. The no-retardation ordered baseline is a numerical baseline, not validation of the Willenborg physical model.
• One dominant crack per part. No multi-site damage, no continuing damage after repair (detected parts leave the fleet).
• The AC 33.14 module covers the Appendix 1 ring disk class of problem: axisymmetric stress, the AC's prescribed crack shapes, no surface transition of subsurface cracks. It is a calibration-grade implementation, not a general rotor risk code.
• The AC POD curves and the 2001 anomaly curve are digitized from the official PDFs' vector art (the FAA never published them as numbers). Estimated digitization accuracy is half a percent in POD and one to two percent in size.
• NASGRO 4's threshold reparametrisation is not fully public; DK1 is read as the long-crack R = 0 threshold per the NASA primary source, which is the slightly conservative interpretation, and is documented in nasgro.py.
• Inputs are independent random variables. Correlation requires a custom limit state.
• This is not a certified tool and nothing here is design data. It is built to certification-grade verification standards (every number traced, every claim tested), which is a different statement.

Repository layout

src/damocles/random_vars.py   input distributions
src/damocles/sampling.py      MC / LHS / Sobol sample generation
src/damocles/reliability.py   pof estimation, exact CIs, importance sampling
src/damocles/fracture.py      growth laws, geometry factors, life integration
src/damocles/newman_raju.py   NASA TM-85793 surface and corner crack solutions
src/damocles/nasgro.py        Forman-Mettu equation with Newman closure
src/damocles/spectrum.py      ASTM E1049 rainflow, spectrum blocks, ordered sequences
src/damocles/retardation.py   Willenborg overload retardation model
src/damocles/inspection.py    POD curves, inspection plans, risk arithmetic
src/damocles/ac3314.py        AC 33.14-1 hard alpha assessment + calibration
src/damocles/materials.py     cited material database (data/materials.json)
src/damocles/sensitivity.py   Sobol indices, Saltelli/Jansen estimators
src/damocles/allowables.py    A- and B-basis tolerance bounds
src/damocles/study.py         YAML-driven studies
src/damocles/cli.py           command line entry point
docs/theory.md                the equations and assumptions
docs/verification.md          claim -> reference -> test matrix
examples/                     disk bore, skin panel, AC test case, coin
tests/                        external-reference and numerical regression tests

The original repository was a high school Monte Carlo toy that dropped hexagonal coins through a slotted grating. That problem is still here, in examples/coin_grating.py, now solved on the reliability engine and checked against quadrature; everything else was torn down and rebuilt.

Tests

python -m pytest -q

The verification matrix in docs/verification.md maps each capability to its reference and its tests.

References

• FAA AC 33.14-1, Damage Tolerance for High Energy Turbine Engine Rotors, 2001, and Change 1, 2017. Anomaly distributions, POD curves, the calibration test case and its acceptance bands.
• Forman, Shivakumar, Cardinal, Williams, McKeighan. Fatigue Crack Growth Database for Damage Tolerance Analysis. DOT/FAA/AR-05/15, 2005. NASGRO and Walker fits with the underlying data.
• Newman, Raju. Stress-Intensity Factor Equations for Cracks in Three-Dimensional Finite Bodies Subjected to Tension and Bending Loads. NASA TM-85793, 1984. Also TP-1578 (1979) and TM-83200 (1981), and Eng. Fract. Mech. 15 (1981) 185-192.
• Mettu, Shivakumar, Beek, Yeh, Williams, Forman, McMahon, Newman. NASGRO 3.0: A software for analyzing aging aircraft. NASA NTRS 19990028759. The NASGRO equation and threshold form.
• Newman. A crack opening stress equation for fatigue crack growth. Int. J. Fracture 24 (1984). The closure function.
• Leverant, McClung, Millwater, Enright et al. Turbine Rotor Material Design. DOT/FAA/AR-00/64, 2000. DARWIN's analysis of the AC test case.
• Willenborg, Engle, Wood. A Crack Growth Retardation Model Using an Effective Stress Concept. AFFDL-TM-71-1-FBR, 1971. https://doi.org/10.21236/ADA956517
• AFGROW Damage Tolerance Design Handbook, Section 5.2.1.2. Equations 5.2.3 through 5.2.6 and the cycle-by-cycle zone update procedure.
• Broek, D. Elementary Engineering Fracture Mechanics, 4th ed. Springer, 1986. Chapter 12 reviews retardation models.
• ASTM E1049-85, Standard Practices for Cycle Counting in Fatigue Analysis. The rainflow rules and the worked example used as an acceptance test.
• MIL-HDBK-1823A, Nondestructive Evaluation System Reliability Assessment, 2009. The lognormal POD model.
• MMPDS / CMH-17 for the A- and B-basis definitions.

License

MIT.