Added a stiffened panel buckling model

python/doc/bibliography.rst

@@ -186,6 +186,8 @@ Bibliography
 .. [knio2010] Le Maître, O., & Knio, O. M. (2010). *Spectral methods for uncertainty
     quantification: with applications to computational fluid dynamics.* Springer
     Science & Business Media.
+.. [ko1994] William L. Ko, Raymond H. Jackson,
+    *Share Buckling Analysis of a Hat-Stiffend Panel*, NASA Technical Memorandum 4644 (November 1994).
 .. [koay2006] Koay C.G., Basser P.J.,
     *Analytically exact correction scheme for signal extraction from noisy magnitude MR signals*,
     Journal of magnetics Resonance 179, 317-322, 2006.

python/doc/examples/reliability_sensitivity/reliability/plot_stiffened_panel.py

@@ -0,0 +1,143 @@
+"""
+Estimate a buckling probability
+===============================
+"""
+# %%
+#
+# In this example, we estimate the probability that the output of a function
+# exceeds a given threshold with the FORM method, the SORM method and an advanced
+# sampling method.
+
+# We consider the :ref:`stiffened panel model <use-case-stiffened-panel>`.
+
+# %%
+# Define the model
+# ----------------
+
+# %%
+from openturns.usecases import stiffened_panel
+import openturns as ot
+import openturns.viewer as viewer
+
+ot.Log.Show(ot.Log.NONE)
+
+# %%
+# We load the stiffened panel model from the usecases module :
+panel = stiffened_panel.StiffenedPanel()
+distribution = panel.distribution
+model = panel.model
+
+# %%
+# See the input distribution
+distribution
+
+# %%
+# See the model
+model.getOutputDescription()
+
+# %%
+# Draw the distribution of a sample of the output.
+sampleSize = 1000
+inputSample = distribution.getSample(sampleSize)
+outputSample = model(inputSample)
+graph = ot.HistogramFactory().build(outputSample).drawPDF()
+_ = viewer.View(graph)
+
+# %%
+# Define the event
+# ----------------
+
+# %%
+# Then we create the event whose probability we want to estimate.
+
+# %%
+vect = ot.RandomVector(distribution)
+G = ot.CompositeRandomVector(model, vect)
+N0 = 165
+event = ot.ThresholdEvent(G, ot.Less(), N0)
+event.setName("buckling")
+
+# %%
+# Estimate the probability with FORM
+# ----------------------------------
+
+# %%
+# Define a solver.
+
+# %%
+optimAlgo = ot.Cobyla()
+optimAlgo.setMaximumEvaluationNumber(1000)
+optimAlgo.setMaximumAbsoluteError(1.0e-10)
+optimAlgo.setMaximumRelativeError(1.0e-10)
+optimAlgo.setMaximumResidualError(1.0e-10)
+optimAlgo.setMaximumConstraintError(1.0e-10)
+
+# %%
+# Run FORM.
+
+# %%
+startingPoint = distribution.getMean()
+algo = ot.FORM(optimAlgo, event, startingPoint)
+n0 = model.getCallsNumber()
+algo.run()
+n1 = model.getCallsNumber()
+result = algo.getResult()
+standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()
+
+# %%
+# Retrieve results.
+
+# %%
+result = algo.getResult()
+probability = result.getEventProbability()
+print("Pf (FORM)=%.3e" % probability, "nb evals=", n1 - n0)
+
+# %%
+# Importance factors.
+
+# %%
+graph = result.drawImportanceFactors()
+view = viewer.View(graph)
+
+# %%
+# Estimate the probability with SORM
+# ----------------------------------
+
+# %%
+# Run SORM.
+
+# %%
+algo = ot.SORM(optimAlgo, event, startingPoint)
+n0 = model.getCallsNumber()
+algo.run()
+n1 = model.getCallsNumber()
+
+# %%
+# Retrieve results.
+
+# %%
+result = algo.getResult()
+probability = result.getEventProbabilityBreitung()
+print("Pf (SORM)=%.3e" % probability, "nb evals=", n1 - n0)
+
+# %%
+# We see that the FORM and SORM approximations give significantly different
+# results. Use a simulation algorithm to get a confidence interval.
+
+# %%
+# Estimate the probability with PostAnalyticalControlledImportanceSampling
+# ------------------------------------------------------------------------
+
+# %%
+algo = ot.PostAnalyticalControlledImportanceSampling(result)
+algo.setBlockSize(100)
+algo.setMaximumOuterSampling(100)
+algo.setMaximumCoefficientOfVariation(0.1)
+n0 = model.getCallsNumber()
+algo.run()
+n1 = model.getCallsNumber()
+result = algo.getResult()
+Pf = result.getProbabilityEstimate()
+print("Pf (sim) = %.3e" % Pf, "nb evals=", n1 - n0)
+width = result.getConfidenceLength(0.95)
+print("C.I (95%)=[" + "%.3e" % (Pf - 0.5 * width), ",%.3e" % (Pf + 0.5 * width), "]")

python/doc/usecases/use_case_stiffened_panel.rst

@@ -0,0 +1,207 @@
+.. _use-case-stiffened-panel:
+
+Stiffened panel buckling
+========================
+
+
+Introduction
+------------
+
+The following figure presents a stiffed panel subject to buckling on a military aircraft.
+
+This use-case implements a simplified model of buckling for a stiffened panel, detailed in [ko1994]_.
+
+.. figure:: ../_static/stiffened_panel_illustration.jpg
+    :align: center
+    :alt: buckling illustration
+    :width: 100%
+
+    **Figure 1.** Buckling of a stiffened panel.
+
+
+.. figure:: ../_static/stiffened_panel_simulation.png
+    :align: center
+    :alt: buckling simulation
+    :width: 100%
+
+    **Figure 2.** 3D simulation of buckling.
+
+
+.. figure:: ../_static/stiffened_panel_description.png
+    :align: center
+    :alt: stiffened panel geometry
+    :width: 100%
+
+    **Figure 3.** Parameterization of the stiffened panel.
+
+
+This test case is composed of nine random variables:
+
+- :math:`E\sim\mathcal{TN}(110.0e9, 55.0e9, 99.0e9, 121.0e9)` : Young modulus (Pa)
+
+- :math:`nu\sim\mathcal{U}(0.3675, 0.3825)` : Poisson coefficient (-)
+
+- :math:`h_c\sim\mathcal{U}(0.0285, 0.0315)` : Distance between the mean surface of the hat and the foot of the Stiffener (m)
+
+- :math:`\ell\sim\mathcal{U}(0.04655, 0.05145)` : Length of the stiffener flank (m)
+
+- :math:`f_1\sim\mathcal{U}(0.0266, 0.0294)` : Width of the stiffener foot (m)
+
+- :math:`f_2\sim\mathcal{U}(0.00627, 0.00693)` : Width of the stiffener hat (m)
+
+- :math:`t\sim\mathcal{U}(8.02e-5, 8.181e-5)` : Thickness of the panel and the stiffener (m)
+
+- :math:`a\sim\mathcal{U}(0.6039, 0.6161)` : Width of the panel (m)
+
+- :math:`b_0\sim\mathcal{U}(0.04455, 0.04545)` : Distance between two stiffeners (m)
+
+- :math:`p\sim\mathcal{U}(0.03762, 0.03838)` : Half-width of the stiffener (m)
+
+The output of interest is:
+
+- :math:`(N_{xy})_{cr}`: the critical shear force (N)
+
+We assume that the input variables are independent except the :math:`f_1` and
+:math:`f_2` for which we measure a Spearman correlation of :math:`\rho^S_{12}=-0.8`,
+modelled using a :class:`~openturns.NormalCopula`.
+
+The critical load :math:`(\tau_{xy})_{cr}` of a stiffened panel subject to shear load is given by
+
+.. math::
+  (\tau_{xy})_{cr}=k_{xy}\frac{\pi^2D}{b_0^2t_s}
+
+where
+
+- :math:`a` is the width of the panel;
+
+- :math:`b_0` is the width between too consecutive stiffener feet;
+
+- :math:`t_s` is the thickness of the panel main surface;
+
+- :math:`E_s` is the Young modulus of the panel main surface;
+
+- :math:`\nu_s` is the Poisson coefficient of the panel main surface;
+
+- :math:`D=\frac{E_st_s^3}{12(1-\nu_s^2)}` is the bending coefficient of the
+  panel main surface;
+
+- :math:`k_{xy}` is the load factor associated to shear buckling. It is given as
+  a function of :math:`\frac{b_0}{a}` through the empirical relation
+  :math:`k_{xy}=5.35 + 4\left(\frac{b_0}{a}\right)^2`
+
+It is more convenient to use the shear force :math:`N_{xy}` instead of the shear
+stress component :math:`\tau_{xy}`. It leads to the relation:
+
+.. math::
+   N_{xy}=q_1+q_c
+
+where :math:`q_1` abd :math:`q_c` are the shear fluxes in the panel main surface
+and its stiffener. They are given by:
+
+.. math::
+   q_1=\tau_{xy}t_s=2G_sh_0t_s\frac{\partial^2w}{\partial x\partial y}
+
+and
+
+.. math::
+  q_c=\frac{G_ct_cp}{\ell}\left[h-2h_0+\frac{h_c}{2p}(f_1-f_2)\right]\frac{\partial^2w}{\partial x\partial y}
+
+where
+
+- :math:`G_s=\frac{E_s}{2(1+\nu_s)}` is the shear modulus of the panel main
+  surface;
+
+- :math:`\frac{\partial^2w}{\partial x\partial y}` is the torsion strain of the
+  panel main surface;
+
+- :math:`G_c=\frac{E_c}{2(1+\nu_c)}` is the shear coefficient of the stiffener;
+
+- :math:`t_c` is the thickness of the stiffener;
+
+- :math:`h_c` is the distance between the mean surfaces of the stiffener hat and
+  foot;
+
+- :math:`h=h_c+\frac{t_c+t_s}{2}` is the distance between the mean surfaces of
+  the stiffener hat and the panel main surface;
+
+- :math:`f_1` is the width of the foot of the stiffener;
+
+- :math:`f_2` is the width of the hat of the stiffener;
+
+- :math:`p` is the half-widht of the stiffener;
+
+- :math:`R` is the radius of the circular part of the stiffener;
+
+- :math:`\theta` is the angle of the circular part of the stiffener;
+
+- :math:`\ell` is the length of the stiffener flank;
+
+- :math:`d=\frac{\ell-f_2}{2}-R\theta- :math:` is the half-lenght of the straight
+  part of the flank of the stiffener;
+
+- :math:`A=\ell t_c` is the area of the section of an half-ondulation;
+
+- :math:`\bar{A}=A+pt_s+\frac{1}{2}(f_1-f_2)t_c` is the area of the section of
+  the full panel (main surface and stiffener) bounded by :math:`p`;
+
+- :math:`h_0=\frac{1}{2\bar{A}}\left[A(h_c+t_c+t_s)+\frac{1}{2}t_c(f_1-f_2)(t_c+t_s)\right]`
+  is the distance between the mean surface of the panel main surface and the
+  global geometric center of the pa,nel;
+
+It leads to:
+
+.. math::
+  N_{xy}=q_1(1+q_c/q1)=\tau_{xy}t_s\left[1+\frac{1}{4}\frac{G_ct_c}{G_st_s}\frac{\left(2p(h-2h_0)-h_c(f_1-f_2)\right)}{h_0\ell}\right]
+
+and finally, :math:`(N_{xy})_{cr}` is given by:
+
+.. math::
+  (N_{xy})_{cr}=\bigg[5.35 + 4\left(\frac{b_0}{a}\right)^2\bigg]\bigg[\frac{\pi^2}{b_0^2}\frac{E_st_s^3}{12(1-\nu_s^2)}\bigg]\bigg[1+\frac{1}{4}\frac{G_ct_c}{G_st_s}\frac{(2p(h-2h_0)-h_c(f_1-f_2))}{h_0\ell}\bigg]
+
+For industrial constraints, the stiffener and the main surface are cut in the
+same metal sheet, so :math:`E_c=E_s=E`, :math:`\nu_c=\nu_s=\nu`, :math:`t_c=t_s=t`.
+The final expression of the critical shear force is then:
+
+.. math::
+  (N_{xy})_{cr}=\bigg[5.35 + 4\left(\frac{b_0}{a}\right)^2\bigg]\bigg[\frac{\pi^2}{b_0^2}\frac{Et^3}{12(1-\nu^2)}\bigg]\bigg[1+\frac{1}{4}\frac{(2p(h-2h_0)-h_c(f_1-f_2))}{h_0\ell}\bigg]
+
+with:
+
+- :math:`A=\ell t`
+
+- :math:`\bar{A}=A+t\left(p+\frac{f_1-f_2}{2}\right)`
+
+- :math:`h_0=\frac{1}{2\bar{A}}\left[A(h_c+2t)+t^2(f_1-f_2)\right]`
+
+- :math:`h=h_c+t`
+
+
+References
+----------
+
+* [ko1994]_
+
+Load the use case
+-----------------
+
+We can load this model from the use cases module as follows :
+
+.. code-block:: python
+
+    >>> from openturns.usecases import stiffened_panel
+    >>> sp = stiffened_panel.StiffenedPanel()
+    >>> # Load the stiffened panel use case
+    >>> model = sp.model()
+
+API documentation
+-----------------
+
+.. currentmodule:: openturns.usecases.stiffened_panel
+
+.. autoclass:: StiffenedPanel
+    :noindex:
+
+Examples based on this use case
+-------------------------------
+
+.. minigallery:: openturns.usecases.stiffened_panel.StiffenedPanel

python/doc/usecases/usecases.rst

@@ -27,3 +27,4 @@ Contents
    use_case_oscillator
    coles
    use_case_linthurst
+   use_case_stiffened_panel

python/doc/user_manual/usecases.rst

@@ -24,3 +24,4 @@ Use cases from the usecases module
     usecases.fireSatellite_function.FireSatelliteModel
     usecases.wingweight_function.WingWeightModel
     usecases.oscillator.Oscillator
+    usecases.stiffened_panel.StiffenedPanel