Merge pull request #12 from ladisk/TB3

Tb3

FLife/freq_domain/tovo_benasciutti.py

@@ -4,7 +4,7 @@
 
 class TovoBenasciutti(Narrowband):
     """Class for fatigue life estimation using frequency domain 
-    method by Tovo and Benasciutti[1, 2].
+    method by Tovo and Benasciutti[1, 2, 3].
       
     References
     ----------
@@ -14,10 +14,13 @@ class TovoBenasciutti(Narrowband):
     [2] Denis Benasciutti and Roberto Tovo. Spectral methods for lifetime
         prediction under wide-band stationary random processes. International
         Journal of Fatigue, 27(8):867{877, 2005
-    [3] Janko Slavič, Matjaž Mršnik, Martin Česnik, Jaka Javh, Miha Boltežar. 
+    [3] Denis Benasciutti and Roberto Tovo. Comparison of spectral methods for fatigue 
+        analysis of broad-band Gaussian random processes. Probabilistic Engineering Mechanics,
+        21(4), 287-299, 2006
+    [4] Janko Slavič, Matjaž Mršnik, Martin Česnik, Jaka Javh, Miha Boltežar. 
         Vibration Fatigue by Spectral Methods, From Structural Dynamics to Fatigue Damage
         – Theory and Experiments, ISBN: 9780128221907, Elsevier, 1st September 2020
-        
+
     Example
     -------
     Import modules, define time- and frequency-domain data
@@ -72,18 +75,22 @@ def __init__(self, spectral_data):
 
     def _calculate_coefficient(self, method='method 2'):
         """Calculate weigthing parameter b for the Tovo-Benasciutti method. Parameter b is 
-            defined by Tovo and Benasciutti [1,2].
+            defined by Tovo and Benasciutti [1,2,3].
         
         :param method:  string
             - 'method 1': `b` weighting parameter `b` is defined by Tovo[1].
             - 'method 2': `b` weighting parameter `b` is defined by Tovo and Benasciutti [2].
-                          (This is the improved method)
+                          (This is the 2005 improved method [2])
+            - 'method 3': `b` weighting parameter `b` is defined by Tovo and Benasciutti [3].
+                          (This is the 2006 improved method [3])
         :return b: float
         """
         if method == 'method 1': 
             b = self._calculate_coefficient_method_1()
         elif method == 'method 2': 
             b = self._calculate_coefficient_method_2()
+        elif method == 'method 3': 
+            b = self._calculate_coefficient_method_3()
         else: 
             raise Exception('Unrecognized Input Error')
         return b
@@ -103,7 +110,7 @@ def _calculate_coefficient_method_1(self):
         return b
 
     def _calculate_coefficient_method_2(self):
-        """Calculate weigthing parameter b for improved Tovo-Benasciutti method. Parameter b is 
+        """Calculate weigthing parameter b for the 2005 improved Tovo-Benasciutti method. Parameter b is 
             defined by Tovo and Benasciutti [2].
         
         :return b: float
@@ -115,6 +122,19 @@ def _calculate_coefficient_method_2(self):
 
         return b
 
+    def _calculate_coefficient_method_3(self):
+        """Calculate weigthing parameter b for the 2006 improved Tovo-Benasciutti method. Parameter b is 
+            defined by Tovo and Benasciutti [3].
+        
+        :return b: float
+        """
+        alpha075 = self.spectral_data.alpha075
+        alpha2 = self.spectral_data.alpha2
+
+        b = (alpha075**2-alpha2**2) / (1-alpha2**2)
+
+        return b
+
     def get_PDF(self, s, method='method 2'):
         """Returns cycle PDF(Probability Density Function) as a function of stress s.
 
@@ -125,7 +145,8 @@ def get_PDF(self, s, method='method 2'):
             - 'method 1': weighting parameter `b` is defined by Tovo[1].
               
             - 'method 2': weighting parameter `b` is defined by Tovo and Benasciutti [2].
-              (This is the improved method)
+
+            - 'method 3': weighting parameter `b` is defined by Tovo and Benasciutti [3].
 
         :return: function pdf(s)
         """
@@ -141,7 +162,7 @@ def pdf(s):
         return pdf(s)
 
     def get_life(self, C, k, method='method 2', integrate_pdf=False):
-        """Calculate fatigue life with parameters C, k, as defined in [3].
+        """Calculate fatigue life with parameters C, k, as defined in [4].
 
         :param C: [int,float]
             S-N curve intercept [MPa**k].
@@ -151,7 +172,7 @@ def get_life(self, C, k, method='method 2', integrate_pdf=False):
         
             - 'method 1': weighting parameter `b` is defined by Tovo[1].
             - 'method 2': weighting parameter `b` is defined by Tovo and Benasciutti [2].
-              (This is the improved method)
+            - 'method 3': weighting parameter `b` is defined by Tovo and Benasciutti [3].
 
         :param integrate_pdf:  boolean
             If true the the fatigue life is estimated by integrating the PDF, 

FLife/tools.py

@@ -62,7 +62,7 @@ def pdf(x):
     return pdf
 
 
-def random_gaussian(freq, PSD, T, fs, rg=None, **kwargs):
+def random_gaussian(freq, PSD, T, fs, rg=None, random_amplitude=True, **kwargs):
     """
     Stationary Gaussian realization of random process, characterized by PSD.
     
@@ -78,6 +78,9 @@ def random_gaussian(freq, PSD, T, fs, rg=None, **kwargs):
         Sampling frequency [Hz].
     :param rg: numpy.random._generator.Generator
         Initialized Generator object
+    :param random_amplitude: Boolean
+        If true, Rayleigh distributed amplitude is used in addition 
+        to uniformly distributed phase. Defaults to True
     :return: t, signal
         Time and stationary Gaussian realization of random process
     
@@ -92,6 +95,7 @@ def random_gaussian(freq, PSD, T, fs, rg=None, **kwargs):
     Dover Publications, 2005
     """ 
     # time and frequency data
+    var = np.trapz(PSD,freq)
     N = int(T * fs)
     M = N//2 + 1
     t = np.arange(0,N) / fs # time vector
@@ -100,20 +104,27 @@ def random_gaussian(freq, PSD, T, fs, rg=None, **kwargs):
     freq_new = np.arange(0, M, 1) / N  * fs # frequency vector
     PSD_new = np.interp(freq_new, freq, PSD, left=0, right=0)
 
-    ampl_spectra = np.sqrt(PSD_new * N * fs / 2)  # amplitude spectra modulus
-
     if rg == None:
         rg = np.random.default_rng()
+
     if isinstance(rg, np.random._generator.Generator):
+
+        phase = rg.uniform(0, 1, len(PSD_new))
+
+        ampl_spectra = np.sqrt(PSD_new * N * fs / 2)  # amplitude spectra modulus
+        if random_amplitude == True:
+            ampl_spectra = np.array([rg.rayleigh(size=1, scale=x)[0] for x in ampl_spectra]) # Rayleigh distributed amplitudes
+
         ampl_spectra_random = ampl_spectra * np.exp(
-            1j * rg.uniform(0, 1, len(PSD_new)) * 2 * np.pi
-        )  # amplitude spectra,  random phase
+            1j * phase * 2 * np.pi)  # amplitude spectra,  random phase
+
     else:
         raise ValueError(
             '`rg` must be initialized Generator object (numpy.random._generator.Generator)!'
         )
 
     signal = np.fft.irfft(ampl_spectra_random, n=N)  # time signal
+    signal = signal/np.std(signal) * var**.5
     return t, signal
 
 

tests/test_basic.py

@@ -17,6 +17,7 @@ def test_data():
         'Dirlik': 1595.031603,
         'Tovo Benasciutti 1': 1099.438159,
         'Tovo Benasciutti 2': 1657.218508,
+        'Tovo Benasciutti 3': 1842.281870,
         'Zhao Baker 1': 1473.240454,
         'Zhao Baker 2': 1564.902582,
         'Narrowband': 1064.200649,
@@ -81,6 +82,7 @@ def test_data():
         'Rice': quad(sd.get_peak_PDF, a=-np.Inf, b=np.Inf)[0],
         'Tovo Benasciutti 1': quad(tb.get_PDF, a=0, b=np.Inf, args=('method 1',))[0],
         'Tovo Benasciutti 2': quad(tb.get_PDF, a=0, b=np.Inf, args=('method 2',))[0],
+        'Tovo Benasciutti 3': quad(tb.get_PDF, a=0, b=np.Inf, args=('method 3',))[0],
         'Dirlik':  quad(dk.get_PDF, a=0, b=np.Inf)[0],
         'Zhao Baker 1': quad(zb.get_PDF, a=0, b=np.Inf, args=('method 1',))[0],
         'Zhao Baker 2': quad(zb.get_PDF, a=0, b=np.Inf, args=('method 2',))[0],
@@ -99,6 +101,7 @@ def test_data():
         'Alpha 0.75': a075.get_life(C=C, k=k),
         'Tovo Benasciutti 1': tb.get_life(C=C, k=k, method='method 1'),
         'Tovo Benasciutti 2': tb.get_life(C=C, k=k),
+        'Tovo Benasciutti 3': tb.get_life(C=C, k=k, method='method 3'),
         'Dirlik': dk.get_life(C=C, k=k),
         'Zhao Baker 1': zb.get_life(C=C, k=k),
         'Zhao Baker 2': zb.get_life(C=C, k=k, method='method 2'),
@@ -129,6 +132,7 @@ def test_data():
         'Narrowband': nb.get_life(C = C, k=k, integrate_pdf=True),
         'Tovo Benasciutti 1': tb.get_life(C = C, k=k, method='method 1', integrate_pdf=True),
         'Tovo Benasciutti 2': tb.get_life(C = C, k=k, integrate_pdf=True),
+        'Tovo Benasciutti 3': tb.get_life(C = C, k=k, method='method 3', integrate_pdf=True),
         'Dirlik': dk.get_life(C = C, k=k, integrate_pdf=True),
         'Zhao Baker 1': zb.get_life(C = C, k=k, integrate_pdf=True),
         'Zhao Baker 2': zb.get_life(C = C, k=k, method='method 2', integrate_pdf=True),