Updated Weibull_3p files

srlife/damage.py

@@ -31,7 +31,7 @@ def __init__(self, pset, *args, cares_cutoff=True, page=False):
         page (bool):    if true, page or store results to disk
         """
         self.cares_cutoff = cares_cutoff
-        self.shear_sensitive = pset.get_default("shear_sensitive", True)
+        self.shear_sensitive = pset.get_default("shear_sensitive", False)
         self.page = page
         self.page_prefix = ""
 
@@ -956,11 +956,13 @@ def calculate_volume_flaw_flattened_eq_stress(
         self.material = material
         self.mandel_stress = mandel_stress
 
+        self.suvals = material.threshold_vol(temperatures) 
         self.mvals = material.modulus_vol(temperatures)
         N = material.Nv(temperatures)
         B = material.Bv(temperatures)
 
         # Temperature average values
+        suavg = np.mean(self.suvals,axis=0)
         mavg = np.mean(self.mvals, axis=0)
         Navg = np.mean(N, axis=0)
         Bavg = np.mean(B, axis=0)
@@ -1047,6 +1049,11 @@ def calculate_volume_flaw_flattened_eq_stress(
                 + (sigma_e_max ** (Navg[..., None, None] - 2))
             ) ** (1 / (Navg[..., None, None] - 2))
 
+            # Subtracting threshold stress 
+            sigma_e_0 -= suavg[...,None, None]
+            sigma_e_0[sigma_e_0 < 0] = 0
+
+
         # Defining area integral element wise
         integral = (
             (sigma_e_0 ** mavg[..., None, None])
@@ -1098,6 +1105,7 @@ def calculate_surface_flaw_flattened_eq_stress(
         self.material = material
         self.mandel_stress = mandel_stress
 
+        suvals = material.threshold_surf(temperatures) 
         mvals = material.modulus_surf(temperatures)
         N = material.Ns(temperatures)
         B = material.Bs(temperatures)
@@ -1106,6 +1114,7 @@ def calculate_surface_flaw_flattened_eq_stress(
         count_surface_elements = np.count_nonzero(surface_elements)
 
         # Temperature average values
+        suavg = np.mean(suvals,axis=0)[:count_surface_elements]
         mavg = np.mean(mvals, axis=0)[:count_surface_elements]
         Navg = np.mean(N, axis=0)[:count_surface_elements]
         Bavg = np.mean(B, axis=0)[:count_surface_elements]
@@ -1198,6 +1207,10 @@ def calculate_surface_flaw_flattened_eq_stress(
                 + (sigma_e_max ** (Navg[..., None, None, None] - 2))
             ) ** (1 / (Navg[..., None, None, None] - 2))
 
+            # Subtracting threshold stress 
+            sigma_e_0 -= suavg[...,None, None, None]
+            sigma_e_0[sigma_e_0 < 0] = 0
+
         # Defining area integral element wise
         integral = (sigma_e_0 ** mavg[..., None, None, None]) * self.ddelta
 
@@ -1262,8 +1275,6 @@ def calculate_volume_flaw_element_log_reliability(
         mavg = np.mean(mvals, axis=0)
         kavg = np.mean(kvals, axis=0)
 
-        # shear_sensitive = True
-
         if self.shear_sensitive is True:
             try:
                 kbar = self.calculate_volume_flaw_kbar(
@@ -1428,19 +1439,21 @@ def calculate_volume_flaw_element_log_reliability(
         pstress = self.calculate_volume_principal_stress(mandel_stress)
 
         # Material parameters
+        suvals = material.threshold_vol(temperatures) 
         svals = material.strength_vol(temperatures)
         mvals = material.modulus_vol(temperatures)
         kvals = svals ** (-mvals)
         N = material.Nv(temperatures)
         B = material.Bv(temperatures)
 
         # Temperature average values
+        suavg = np.mean(suvals,axis=0)
         mavg = np.mean(mvals, axis=0)
         kavg = np.mean(kvals, axis=0)
         Navg = np.mean(N, axis=0)
         Bavg = np.mean(B, axis=0)
 
-        # Only tension
+        # Only tension 
         pstress[pstress < 0] = 0
 
         # Max principal stresss over all time steps for each element
@@ -1463,10 +1476,15 @@ def calculate_volume_flaw_element_log_reliability(
             + (pstress_max ** (Navg[..., None] - 2))
         ) ** (1 / (Navg[..., None] - 2))
 
-        # Log reliability in each element
-        log_reliability = (
-            -kavg * np.sum(pstress_0 ** mavg[..., None], axis=-1) * volumes
-        )
+        # Subtracting threshold stress 
+        pstress_0 -= suavg[...,None]
+        pstress_0[pstress_0 < 0] = 0
+
+        ## Without threshold
+        log_reliability =  (
+                -kavg * np.sum(pstress_0 ** mavg[..., None], axis=-1) * volumes
+                )
+
 
         return log_reliability
 
@@ -1498,6 +1516,7 @@ def calculate_surface_flaw_element_log_reliability(
         """
 
         # Material parameters
+        suvals = material.threshold_surf(temperatures) 
         svals = material.strength_surf(temperatures)
         mvals = material.modulus_surf(temperatures)
         kvals = svals ** (-mvals)
@@ -1508,6 +1527,7 @@ def calculate_surface_flaw_element_log_reliability(
         count_surface_elements = np.count_nonzero(surface)
 
         # Temperature average values
+        suavg = np.mean(suvals,axis=0)[:count_surface_elements]
         mavg = np.mean(mvals, axis=0)[:count_surface_elements]
         kavg = np.mean(kvals, axis=0)[:count_surface_elements]
         Navg = np.mean(N, axis=0)[:count_surface_elements]
@@ -1545,6 +1565,10 @@ def calculate_surface_flaw_element_log_reliability(
             + (surf_pstress_max ** (Navg[..., None, None] - 2))
         ) ** (1 / (Navg[..., None, None] - 2))
 
+        # Subtracting threshold stress 
+        surf_pstress_0 -= suavg[...,None, None]
+        surf_pstress_0[surf_pstress_0 < 0] = 0
+
         # Summing up over last two axes i.e. over the elements of surface stress and
         # surfaces for which normals are there
 
@@ -1586,11 +1610,13 @@ def calculate_volume_flaw_avg_normal_stress(
         self.material = material
         self.mandel_stress = mandel_stress
 
+        suvals = material.threshold_vol(temperatures) 
         mvals = material.modulus_vol(temperatures)
         N = material.Nv(temperatures)
         B = material.Bv(temperatures)
 
         # Temperature average values
+        suavg = np.mean(suvals,axis=0)
         mavg = np.mean(mvals, axis=0)
         Navg = np.mean(N, axis=0)
         Bavg = np.mean(B, axis=0)
@@ -1674,6 +1700,10 @@ def calculate_volume_flaw_avg_normal_stress(
                 + (sigma_n_max ** (Navg[..., None, None] - 2))
             ) ** (1 / (Navg[..., None, None] - 2))
 
+            # Subtracting threshold stress 
+            sigma_n_0 -= suavg[...,None, None]
+            sigma_n_0[sigma_n_0 < 0] = 0
+
         integral = (
             (sigma_n_0 ** mavg[..., None, None])
             * np.sin(self.A)
@@ -1716,6 +1746,7 @@ def calculate_surface_flaw_avg_normal_stress(
         self.material = material
         self.mandel_stress = mandel_stress
 
+        suvals = material.threshold_surf(temperatures) 
         mvals = material.modulus_surf(temperatures)
         N = material.Ns(temperatures)
         B = material.Bs(temperatures)
@@ -1728,6 +1759,7 @@ def calculate_surface_flaw_avg_normal_stress(
         count_surface_elements = np.count_nonzero(surface)
 
         # Temperature average values
+        suavg = np.mean(suvals,axis=0)[:count_surface_elements]
         mavg = np.mean(mvals, axis=0)[:count_surface_elements]
         Navg = np.mean(N, axis=0)[:count_surface_elements]
         Bavg = np.mean(B, axis=0)[:count_surface_elements]
@@ -1822,6 +1854,11 @@ def calculate_surface_flaw_avg_normal_stress(
                 + (sigma_n_max ** (Navg[..., None, None, None] - 2))
             ) ** (1 / (Navg[..., None, None, None] - 2))
 
+            # Subtracting threshold stress 
+            sigma_n_0 -= suavg[...,None, None, None]
+            sigma_n_0[sigma_n_0 < 0] = 0
+
+
         integral = ((sigma_n_0 ** mavg[..., None, None, None]) * self.ddelta) / (
             2 * np.pi
         )

srlife/data/damage/SiC.xml

@@ -1,6 +1,14 @@
 <!-- Weibull strength should be compatible with units of K, MPa, and mm -->
 <models type="ceramic">
   <base type="StandardModel">
+    <threshold_vol>
+      <temperatures>298.15 1073.15 1273.15 1473.15 1673.15 1773.15</temperatures>
+      <values>0.0 0.0 0.0 0.0 0.0 0.0</values>
+    </threshold_vol>
+    <threshold_surf>
+      <temperatures>298.15 1073.15 1273.15 1473.15 1673.15 1773.15</temperatures>
+      <values>0.0 0.0 0.0 0.0 0.0 0.0</values>
+    </threshold_surf>
     <strength_vol>
       <temperatures>298.15 1073.15 1273.15 1473.15 1673.15 1773.15</temperatures>
       <values>507.0 467.0 528.0 570.0 476.0 471.0</values>
@@ -37,6 +45,14 @@
     </fatigue_Bs>
   </base>
   <cares type="StandardModel">
+    <threshold_vol>
+      <temperatures>298.15 1073.15 1273.15 1473.15 1673.15 1773.15</temperatures>
+      <values>0.0 0.0 0.0 0.0 0.0 0.0</values>
+    </threshold_vol>
+    <threshold_surf>
+      <temperatures>298.15 1073.15 1273.15 1473.15 1673.15 1773.15</temperatures>
+      <values>0.0 0.0 0.0 0.0 0.0 0.0</values>
+    </threshold_surf>
     <strength_vol>
       <temperatures>298.15 811.15 1089.15 1366.15 1566.15</temperatures>
       <values>552.94 552.94 552.94 552.94 552.94</values>

srlife/materials.py

@@ -708,6 +708,9 @@ class StandardCeramicMaterial:
 
     def __init__(
         self,
+        su_temperatures,
+        threshold_v,
+        threshold_s,
         s_temperatures,
         strengths_v,
         strengths_s,
@@ -727,6 +730,11 @@ def __init__(
     ):
         super().__init__(*args, **kwargs)
 
+        self.su_temperatures = su_temperatures
+        self.threhold_v = threshold_v
+        self.threhold_s = threshold_s
+        self.su_v = inter.interp1d(s_temperatures, threshold_v)
+        self.su_s = inter.interp1d(s_temperatures, threshold_s)
         self.s_temperatures = s_temperatures
         self.strengths_v = strengths_v
         self.strengths_s = strengths_s
@@ -750,6 +758,18 @@ def __init__(
         self.Bv = inter.interp1d(Bv_temperatures, Bvvals)
         self.Bs = inter.interp1d(Bv_temperatures, Bsvals)
 
+    def threshold_vol(self, T):
+        """
+        Weibull threshold parameter for volume flaws as a function of temperature
+        """
+        return self.su_v(T)
+
+    def threshold_surf(self, T):
+        """
+        Weibull threshold parameter for surface flaws as a function of temperature
+        """
+        return self.su_s(T)
+
     def strength_vol(self, T):
         """
         Weibull strength for volume flaws as a function of temperature
@@ -824,6 +844,13 @@ def load(cls, node):
         Parameters:
           node:    node with model
         """
+        threshold_v = node.find("threshold_vol")
+        su_temps = threshold_v.find("temperatures")
+        su_vals_v = threshold_v.find("values")
+
+        threshold_s = node.find("threshold_surf")
+        su_vals_s = threshold_s.find("values")
+
         strength_v = node.find("strength_vol")
         s_temps = strength_v.find("temperatures")
         svals_v = strength_v.find("values")
@@ -856,6 +883,9 @@ def load(cls, node):
         Bsvals = Bs.find("values")
 
         return StandardCeramicMaterial(
+            np.array(list(map(float, su_temps.text.strip().split()))),
+            np.array(list(map(float, su_vals_v.text.strip().split()))),
+            np.array(list(map(float, su_vals_s.text.strip().split()))),
             np.array(list(map(float, s_temps.text.strip().split()))),
             np.array(list(map(float, svals_v.text.strip().split()))),
             np.array(list(map(float, svals_s.text.strip().split()))),