Monte Carlo simulation of failure probability of a hydrogen tank¶

This notebook is an element of the free risk-engineering.org courseware. It can be distributed under the terms of the Creative Commons Attribution-ShareAlike licence.

Author: Eric Marsden eric.marsden@risk-engineering.org

---

This notebook contains an introduction to the use of Python and the NumPy library for Monte Carlo simulation applied to a mechanical design problem, used for estimating failure probability. See the associated course materials for an introduction to the use of stochastic simulation methods and to download this content as a Jupyter/Python notebook.

Acknowledgement. The problem description in this notebook is taken from UQWorld, who themselves extracted the problem from the article by Bichon et al in 2011 cited below.

Problem statement¶

The 5-dimensional liquid hydrogen tank problem is a reliability analysis benchmark problem (Bichon et al., 2011). The problem consists in quantifying the failure probability of a liquid hydrogen fuel tank on a space launch vehicle. The structure of the tank is subjected to stresses caused by ullage (empty space) pressure, head pressure (caused by the weight of the fuel), axial forces due to acceleration, and bending and shear stresses caused by the weight of the fuel.

Three different failure modes exist for the tank:

• Von Mises stress failure ($PVM$)
• Isotropic strength failure ($PIS$)
• Honeycomb buckling failure ($PHB$)

The limit state function of the overall system is then defined as:

$$g(x) = \min(PVM, PIS, PHB)$$

where:

$$PVM = \frac{8400 \cdot t_{plate}}{\sqrt{N_x^2 + N_y^2 - N_xN_y + 3N_{xy}^2}} - 1$$$$PIS = \frac{8400 \cdot t_{plate}}{|Ny|} - 1$$$$PHB = 0.847 + 0.96 \cdot x_1 + 0.986 \cdot x_2 - 0.216 \cdot x_3 + 0.077 \cdot x_1^2 + 0.11 \cdot x_2^2 + 0.007 \cdot x_3^2 + 0.378 \cdot x_1x_2 - 0.106 \cdot x_1x_3 - 0.11 \cdot x_2x_3$$

and:

$x_1 = 4 * (t_{plate} - 0.075)$

$x_2 = 20 * (t_h - 0.1)$

$x_3 = -6000 * (\frac{1}{N_{xy}} + 0.003)$

The failure event is defined as $g(x) ≤ 0$ (this means that one or more of PVM, PIS and PHB have failed).

Source article: B. J. Bichon, J. M. McFarland, and S. Mahadevan, “Efficient surrogate models for reliability analysis of systems with multiple failure modes”, Reliability Engineering & System Safety, 96(10):1386-1395, 2011. DOI: 10.1016/j.ress.2011.05.008.

Input variables¶

We will make the following assumptions concerning the distribution of our input variables, which we assume are independent:

Name	Distribution	Description
$t_{plate}$	Normal(0.07433, 0.005)	Plate thickness
$t_h$	Normal(0.1, 0.01)	Honeycomb thickness
$N_x$	Normal(13, 60)	Load on tank, x-component
$N_y$	Normal(4751, 48)	Load on tank, y-component
$N_{xy}$	Normal(-648, 11)	Load on tank, xy-component

We want to use a Monte Carlo simulation to estimate the probability of failure, i.e. the probability that $g(x) ≤ 0$.

Method¶

Start by defining a function that generates a random input value (note that the input space is 5-dimensional, here represented using a Python dictionary), then functions that calculate the value of PVM, PIS and PHB for some input.

To get a feeling for the phenomena we are simulating, we can plot histograms of the three failure types PVM, PIS and PHB.

Now we can use a stochastic simulation method (here a naïve Monte Carlo method) to estimate the failure probability. We run a simulation with random inputs a large number of times, and count the proportion of simulations where a failure occurred.

Vectorization¶

The Jupyter magic code %%time used in the cell above leads to the computation time for the cell being printed. As you will have noticed if you’re executing the notebook cell by cell, the element-by-element approach used above (one calculation per loop iteration) is very very slow in Python!

We can make things run faster using NumPy’s vectorization capabilities, its ability to operate on (add, subtract, take the square root, etc.) vectors of values, and not just on single values. Below, we write the functions that calculate PVM, PIS and PHB in a form that allows vectorization (each argument such as tplate can be either a single value or a vector of values).

Now we run the naïve Monte Carlo simulation in vectorized form. We start by generating vectors of values for the 5 input variables, then call the failure calculation functions on these vectors. We can them compare the failure vectors itemwise with a vector full of zeros and count the number of simulations where at least one of the results was negative.

Note that we are running 100 times more simulations here, but the calculation is nevertheless more than 60 times faster. So the vectorization has led to a speedup of a factor of 600.

For some computations of this form, the specialized vector processing hardware present in graphics cards (the GPU) is much faster than the general-purpose hardware present in the CPU. Various software libraries are available to use the GPU instead of the CPU for vector processing, and you can access these in Python using the CuPy array library, which provides a drop-in replacement for the NumPy and SciPy libraries.

Calculations of this form are nonetheless not very fast in Python. Another notebook on the risk-engineering.org website shows how to run this calculation in the Rust programming language, with much higher performance.

A Latin Hypercube Sampling (LHS) approach¶

The LHS method consists of dividing the input space into a number of equiprobable regions, then taking random samples from each region. We can use it conveniently in Python thanks to the pyDOE library, which you will probably need to install on your computer, using a command such as

pip install pyDOE

or if you’re using a Google CoLaboratory notebook, execute a code cell containing

!pip install pyDOE

The lhs function in this library returns an “experimental design” consisting of points in the $[0, 1]^d$ hypercube, where $d$ is the dimension of your problem (it’s 2 in this simple example). You need to scale these points to your input domain.

Using Halton’s low discrepancy sequences¶

A low discrepancy (or quasi-random) sequence is a deterministic mathematical sequence of numbers that has the property of low discrepancy. This means (informally) that there are no clusters of points and that the sequence fills space roughly uniformly. The Halton sequence is a low discrepancy sequence that has useful properties for pseudo-stochastic sampling methods (also called “quasi-Monte Carlo” methods).

Attempting an analytical solution¶

It’s quite easy to specify this problem in analytical terms, using the SymPy library for symbolic calculation in Python. However, this problem is too complicated for the current version (in 2020) of SymPy, and leads to an error. (It also consumes a large amount of memory before failing.)